Upper plate controls on deep subduction, trench migrations and deformations at convergent margins

Tectonophysics 483 (2010) 80–92

Contents lists available at ScienceDirect

Tectonophysics 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 / t e c t o

Upper plate controls on deep subduction, trench migrations and deformations at convergent margins F.A. Capitanio a,b,⁎, D.R. Stegman c, L.N. Moresi a,b, W. Sharples a a b c

School of Mathematical Sciences, Monash University, Clayton, Victoria, Australia School of Geosciences, Monash University, Clayton, Victoria, Australia School of Earth Sciences, University of Melbourne, Victoria, Australia

a r t i c l e

i n f o

Article history: Received 3 March 2009 Received in revised form 28 July 2009 Accepted 19 August 2009 Available online 6 September 2009 Keywords: Subduction Upper plate deformations Back-arc basins Plate motions

a b s t r a c t Thus far, relatively simplistic models of free subduction, in which the trench and plate motions are emergent features completely driven by the negative buoyancy of the slab, have investigated the dynamics of a single, isolated subducting plate. Here we extend such models to incorporate an overriding plate and present the results of how such an overriding plate feedbacks into the dynamics of free subduction. In this study, we address three fundamental aspects of these dynamics: 1) How does the presence of an overriding plate change the force balance at the convergent margins? 2) How are the forces from deep subduction propagated to the surface? And 3) what controls the stress regime in a system of coupled upper and subducting plates and how is it expressed in the deformations and plate motions? In general, we find that the evolution of subduction zones is strongly controlled by both the interactions between the slab and the upper-lower mantle discontinuity as well as the strength of the upper plate. When either the subducting or upper plates are unable to move, subduction motions are steady state and partitioned entirely into either slab rollback or plate advance, respectively. When conditions favour a quasi-stationary trench, the subducted lithosphere can form into a pile with multiple recumbent folds of slab material atop the lower mantle. Alternating between forwards- and backwards-draping slab, the corresponding horizontal trench motions at the surface are frontward and rearward, respectively, resulting in either a compressive or extensional regime in the back-arc. Timedependent forcing arising from the slab piling behaviour can have a feedback with upper plate and produce strongly non-steady state, intermittent phases of upper plate deformation as those commonly observed on Earth. Two types of discontinuous back-arc strain evolution are identified: (1) periodic, when recurrent phases of strain over finite durations are accommodated by (viscous) stretching/thickening of the plate, and (2) episodic, when upper plate deformation localizes (plastic strain) and allows for punctuated episodes. These phases can include extension, quiescence, and compression, giving rise to a large variety of possible tectonic evolutions. The models presented here provide insight into the dynamics behind the non-steady-state evolution of subduction, which can help unravel seemingly erratic motions of major convergent margins and back-arc deformations around the Pacific and Indian Oceans during the Cenozoic. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Forces associated with subduction of cold and dense oceanic plates control the motions and deformations around convergent margins. However, how the forces propagate to the upper plates from subducting slabs and evolve through time, is still poorly known. Upper plate thickening and/or stretching is often the result of the interactions between lithosphere on both sides of the subduction zone and the mantle underneath (Uyeda and Kanamori, 1979). Variable stretching characterizes back-arc basins, which are found floored by thinned continental lithosphere, as the Japan Sea (Jolivet et al., 1994), ⁎ Corresponding author. School of Mathematical Sciences, Monash University, Clayton, Victoria, Australia. E-mail address: [email protected] (F.A. Capitanio). 0040-1951/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2009.08.020

or by oceanic crust, i.e. the Mariana Trough. Alternatively, thickened continental lithosphere, such as the Andes (Allmendinger et al., 1997), is associated with convergent margins. Numerous models of subduction have shown that the foundering of cold lithosphere backwards into the mantle and the associated trench retreat drive back-arc basin spreading (Faccenna et al., 2001; Garfunkel et al., 1986; Hall et al., 2003; Kincaid and Olson, 1987; Shemenda, 1993). The constraints imposed by the subducting plates, sliding below upper plates are critical for basal tractions (Sleep and Toksöz, 1971; Wdowinski et al., 1989), thickening and density redistribution (Husson and Ricard, 2004) as well as for the frictional interface (Buiter et al., 2001; Hampel and Pfiffner, 2006). Yet, the dynamic interactions of integrated subducting slabs and upper plates are not well understood. Geological evidence reveals strongly episodic upper plate deformations as well as trench motions are strongly episodic (Faccenna et al.,

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

2001; Jolivet et al., 1994). Different styles of episodicity are identified in models of subduction dynamics including pseudo-, quasi- and hyper-episodicity (Clarke et al., 2008) corresponding to phases of extension, quiescence and compression during subduction evolution (Sdrolias and Müller, 2006). It has been shown that such complex evolution of subducting margins on Earth often responds to deep subduction dynamics (Goes et al., 2008). Here, the interactions of slabs with the transition zone and the lower mantle might contribute to episodic rearrangement of plate motions and deformations at the surface (Capitanio et al., 2009b; Christensen and Yuen, 1984; van der Hilst and Seno, 1993; Zhong and Gurnis, 1995). Recent kinematic models (Boutelier and Cruden, 2008; Guillame et al., 2009) suggested that the interaction of slabs with the bottom of the mantle might be the cause of repeated deformation at surface, however we have no understanding of the dynamic controls on recurring strain episodes, often inverted, during the evolution of convergent margins. So far, only a few physical models of the subduction system have included the subducting and upper plate, however these mainly focused on thermal processes (Arcay et al., 2005; Gerya et al., 2002), global plate motions (Becker and O'Connell, 2001; Conrad and Lithgow-Bertelloni, 2002), deep penetration dynamics (Zhong and Gurnis, 1995), collision (De Franco et al., 2007; Kaus et al., 2008) and buoyant plateau subduction (Espurt et al., 2008; Mason et al., 2010this issue; van Hunen et al., 2002). Therefore they do not provide insight in the dynamics of convergent margins. Here we propose a simplified model of free subduction system with a plate subducting in the upper mantle underthrusting an overriding plate. Our focus is on the self-consistent subduction dynamics where driving and resisting forces are balanced with plate motions and deformations. This allows understanding of (1) the force propagation from deep subduction to upper plate at surface and the controls on the stress regime, (2) the feedback between deformations in the upper plate and the subducting plate, and (3) the impact of the dynamics on the tectonic style of plate motions. The outcomes of the modelling allow a comparison with the temporal evolution of the main subduction zones in the Pacific and Indian Oceans during the Cenozoic (Sdrolias and Müller, 2006), where behaviour similar to our models is shown. 2. Numerical method and model setup We model subduction as the viscous flow of an infinite Prandtl number fluid at a very low Reynolds number in a two-dimensional Cartesian geometry. We solve the governing mass and momentum conservation equations within the geodynamics modelling framework, Underworld, which tracks distinct geologic materials with Lagrangian particles embedded within the Eulerian finite element grid. The details of the numerical method, of the software implementation and relevant numerical benchmarks are described in Moresi et al., 2003 and Stegman et al., 2006. The subduction system includes the downgoing plate and the upper plate confined to the upper mantle. The model space is 650 km thick and 4500 km wide for a total of 64⋅ 288 elements yielding a spatial resolution of approximately 10 ⋅ 15 km with about 40 particles in each element. Boundary conditions are free-slip on top, no-slip on the bottom, and the sidewalls use periodic boundary conditions which allow for the unconstrained lateral migration of slabs (Enns et al., 2005). The viscosity of the mantle is Newtonian with ηUM = 1019 Pa s and the other parameters common to all numerical experiments are given in Table 1. The 100 km thick downgoing plate of thickness has a layered rheology structure (Fig. 1). The upper 30 km of the plate has a viscoplastic rheology to represent the brittle crust, with an initial (unyielded) viscosity of ηCRUST / ηUM = 30 and a von Mises plasticity criterion τy = τ0, where |0 = 35 MPa, in agreement with several published works (Gurnis et al., 2004; Krien and Fleitout, 2008; Toth and Gurnis, 1998). Underlying the top layer is a 30 km thick strong

81

Table 1 Parameters common to all models. Upper mantle Density Height Width Viscosity

ρUM hUM wUM ηUM

3300 kg/m3 650 km 4500 km 1019 Pa s

Subducting plate Initial length Thickness Density contrast. Cohesion

lSP hSP Δρ τ0

2200 km 100 km 60 kg/m3 35 MPa

Upper Plate Length Thickness Density contrast. Viscosity Gravitational acceleration

lUP hUP ΔρUP ηUP g

1200 km 50 km 0 kg/m3 3000·ηUM 9.81 m/s2

core with a Newtonian viscosity of ηCORE / ηUM = 3000. The lowermost portion of the plate is also a Newtonian viscous layer with the same viscosity as the top layer. This layered plate setup is based on the strong temperature-dependent rheology of lithospheric forming minerals, leading to a concentration of strength in the stiff, lowtemperature core (Kameyama et al., 1999; Karato and Wu, 1993). This model has been tested within the subduction system (Schellart et al., 2007) and is essential for realistic force propagation to plates at the surface (Capitanio et al., 2009a). The trailing edge of the subducting plate is tapered from 100 km to 0 in order to avoid the occurrence of an unnaturally large stress jump. The properties of the subducting plate are not varied in this work, other than the density contrast (Table 1). In this work we mostly focus on the role of the upper plate in subduction dynamics. We first investigate a suite of models subjected to different far-field boundary conditions to better understand to what extent these boundary conditions alter the dynamics of subduction. This suite includes models without an upper plate, and models with an upper plate and varying boundary conditions, i.e. the plate's trailing edge which can be either free or fixed to the end of the box (i.e. vx = 0). In this setup, the tails of both trailing plates can be free (UP models), or the tail of either the subducting plate (SUBfix) or upper plate (UPfix) is fixed while the other plate remains free. In the models with a free upper plate, the trailing edge is tapered, as in the free subducting plate. These reference models have an upper plate viscosity of 3000 ηUM and the buoyancy is held constant at Δρ = 60 kg m− 3. The large viscosity of the upper plate will limit the deformation allowing assessment of the force propagation from slab to upper plate in the absence of time-dependent strain. In other sets of models we focus on the impact the deformation of the upper plate has on the entire system, by varying the rheology in the upper plate. The rheological properties and thickness are varied to cover the wide array of structures in overriding plates on Earth. In these setups we (1) reduce the upper plate initial (unyielded) viscosity to 300 ηUM and (2) introduce viscoplasticity with a low yield stress (58 MPa). The plastic properties of the subducting plate's crust are not varied here. This simple setup allows for non-steady-state subduction, where the

Fig. 1. Model set up of the free upper plate model (UP). Details in Table 1. Initial conditions with slab perturbation in the computational domain are shown.

82

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

Table 2 Parameters of additional models different from reference values. Free plate models UP_l600 UP_rho30 UP_rho90 UP_h30 UP_h30rho30 UP_h30rho90 UP_h80 UP_h80rho30 UP_h80rho90 UP_rhoNeg UP_h20 UP_h30 UP_h100 UP_h150

lUP = 600 km Δρ = 30 kg/m3 Δρ = 90 kg/m3 hUP = 30 km hUP = 30 km, Δρ = 30 kg/m3 hUP = 30 km, Δρ = 90 kg/m3 hUP = 80 km hUP = 80 km, Δρ = 30 kg/m3 hUP = 80 km, Δρ = 90 kg/m3 hUP =80 km, ΔρUP =−60 kg/m3 hUP = 20 km hUP = 30 km hUP = 100 km hUP = 150 km

Fixed upper plate models (vx =0 in x<500 km) UPfix UPfix_rho30 UPfix_rho90 UPfix_h30 UPfix_h30rho30 UPfix_h30rho90 UPfix_h80 UPfix_h80rho30 UPfix_h80rho90 UPfix_h100 UPfix_h150 UPfix_h30W UPfix_h50W UPfix_h80W UPfix_h30P UPfix_h50P UPfix_h80P

lUP = 1700 km Δρ = 30 kg/m3 Δρ = 90 kg/m3 hUP = 30 km hUP = 30 km, Δρ = 30 kg/m3 hUP = 30 km, Δρ = 90 kg/m3 hUP = 80 km hUP = 80 km, Δρ = 30 kg/m3 hUP = 80 km, Δρ = 90 kg/m3 hUP = 100 km hUP = 150 km hUP = 30 km, ηUP = 300·ηUM hUP = 50 km, ηUP = 300·ηUM hUP = 100 km, ηUP = 300·ηUM hUP = 30 km, τ0 = 58 MPa hUP = 50 km, τ0 = 58 MPa hUP = 80 km, τ0 = 58 MPa

Fixed subducting plate models (vx =0 in x>4000 km) SUBfix SUBfix_rho30 SUBfix_rho90 SUBfix_h30 SUBfix_h30rho30 SUBfix_h30rho90 SUBfix_h80 SUBfix_h80rho30 SUBfix_h80rho90 SUBfix_h100 SUBfix_h150

lSP = 3300 km Δρ = 30 kg/m3 Δρ = 90 kg/m3 hUP = 30 km hUP = 30 km, Δρ = 30 kg/m3 hUP = 30 km, Δρ = 90 kg/m3 hUP = 80 km hUP = 80 km, Δρ = 30 kg/m3 hUP = 80 km, Δρ = 90 kg/m3 hUP = 100 km hUP = 150 km

episodic nature of overriding plate deformation is characterized. This is done by measuring the velocity on the top layer of the model or in the center of plates. Similarly, the force propagated to the surface is measured along the plates as a horizontal force integrated over plate thickness, i.e. Fxx = − ∫hσxxdy. Subduction initiates freely under the pull of an initial instability represented by the tip of the slab that extends to a 150 km depth (250 km for thicker upper plate models). The subsequent evolution is self-sustaining owing to the negative buoyancy of the developing slab. Stresses generated at the convergent margin are largely accommodated by the viscoplastic upper layer of the subducting plate, localizing the deformation and allowing the downgoing plate to slide past the overriding plate. The force from the slab is propagated back to the surface through the stiff core of the subducting plate. The overriding and subducting plates are coupled through viscous flow in the mantle wedge and the “suction” it generates (see Forsyth and Uyeda, 1975). Using both of these upper plate rheologies, we investigate the role of the subducting plate's buoyancy (reduced to the density contrast Δρ), which is varied between 30, 60 and 90 kg m− 3, which results in different slab pull driving forces. We also investigate the effect of thicker continents and thinner oceanic plates through variations in upper plate thickness, ranging between 30, 50, and 80 km. A full description of all the models is provided in Table 2. 3. Results 3.1. Effect of the upper plate Here we compare the reference model with a free upper plate (UP) to a model similar in all ways except no overriding plate is included, and address briefly the differences between these two cases of subduction (Fig. 2). The first observation is that the sinking velocities are essentially the same. This indicates that the dissipation due to bending or upper plate deformation is minor and that most of the potential energy of the system is dissipated in the mantle (Capitanio et al., 2007). This further suggests that the presence of a free upper plate does not alter how the subduction velocity is partitioned at the surface between slab rollback and downgoing plate advance (Fig. 2A). Additionally, strain rates achieved in the wedge are comparable in magnitude (Fig. 3). However, the concentration of the corner flow in a smaller wedge, when the upper plate is present, increases the strain rates below the overriding plate (Fig. 3B–C). This increases the shear stress at the base of the upper plate that is larger than that measured at the top of the mantle when the upper plate is not included (Fig. 3A, red lines). The stress propagated to

Fig. 2. Comparison of subduction models with an upper plate (UP model) and with a single subducting plate (SP). (A) Horizontal velocity partitioning between plate advance towards the trench (negative) and upper plate-mantle motions (positive). (B) Viscosity of subducting plate model (SP). (C) Viscosity of free upper plate model (UP). Viscosity colorscale as in Fig. 1.

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

83

3.2. Effects of far-field boundary conditions

Fig. 3. Stress and strain rates in the subducting plate (SP) and free upper plate (UP) models. (A) Shear stress (red lines) and horizontal stress (black lines) of the model with the subducting plate only (SP, dotted lines) and the free upper plate (UP, solid lines). The bending area is shaded as stresses are reoriented and not parallel to the horizontal. Direction is positive towards the trench. (B) Detail of the strain rates in the trench zone (note the different horizontal scale) of the subducting plate model (SP). Strain localizes in the crustal layer and fades quickly in the wedge. (C) Detail of the strain rates in the trench zone of the free upper plate model (UP). The larger strain propagates below the upper plate, whereas in the plastic layer and downdip it is very similar to the SP model.

the subducting plate is the same in both setups, and almost constant through the stiff core along the plate (Fig. 3A, black lines). However, the stress above subduction zones is sensitive to whether the upper plate is included. It is interesting that while the mantle does not resist the retreating motion of the trench, the upper plate is observed to be in compression around the trench. Furthermore, it appears that stresses occurring in the upper plate are larger than those in the subducting plate. This suggests that the two effects are coupled as stress propagated across the interface between the two plates increases the basal tractions on the upper plate which generates larger compressive stresses. These are consistent features of all the models tested. As a concluding remark, we note that the stresses propagated through the plates are larger than the shear drag at their base. This shows that in our models, the subduction system is driven by slab pull propagation to plates and not by mantle viscous drag.

The subduction motions evolve in different ways according to their far-field boundary conditions (Fig. 4). When the upper plate is fixed (Fig. 4A, UPfix) subduction occurs entirely by plate advance towards a fixed trench, while slab accumulates underneath the trench at depth. In the opposite case, when the subducting plate is fixed (Fig. 4C, SUBfix) it is consumed by trench retreat as the upper plate overrides. Although mantle flow in the three models is significantly different, similar strain rates are achieved (Fig. 4D–F). In the models with a fixed upper plate, the advancing plate drives flow predominantly, hence mantle circulation is limited under the subducting lithosphere, while the mantle under the upper plate is nearly stagnant (Fig. 4D). When the upper plate is free as in UP and SUBfix models, subduction is largely accommodated by slab rollback (60% and 100%, respectively). The flow in the mantle develops throughout the entire mantle (Fig. 4E–F) which influences the sinking of the models. In fact, all the models driven by the slab pull sink at the same rate, varying with the density contrast (Fig. 5A), although slight differences can be noted. The fixed upper plate model has a timedependent sinking velocity depending on the amount of slab piled atop the lower mantle. Hence, although the buoyancy per unit length is the same, the slab pull integrated over its length varies by up to ±25%. The differences between the sinking rates of the free upper plate model (UP and SUBfix) are here attributed to the larger resistance found in the viscous mantle when the subducting plate is fixed and the slab has larger (forced) lateral migration resulting in slower slab motions (~20%). Despite the difference found in the configurations, the sinking velocities scale with the buoyancy (Fig. 5A), suggesting that the Stokes sinker approximation found in free subduction of a single plate (i.e. without upper plate) (Capitanio et al., 2007) is appropriate. The sinking velocities for a given buoyancy are comparable in the different models (Fig. 5A) and the total subduction velocity scales with the sinking velocity which then determines the magnitude of the slab pull. However, the horizontal partitioning of subduction motions varies substantially from completely rollback-dominated to entirely plate advance towards a stationary trench (Fig. 5B). In the models with both plates free (UP models), the amount of convergence achieved by rollback increases from 40%, to 50% and 60% with increasing buoyancy. This is consistent with the observation that most of the work done by slab pull is being dissipated in the mantle, implying that although farfield boundary conditions constrain the mobility of plates controlling plate and trench motions, the underlying dynamics is the same for the different configurations. We have also investigated varying the thickness in models with both plates free (20, 30, 100 and 150 km in models UPh20, UPh30, UPh100 and UPh150, respectively), the length of upper plate (model UPl600) and the density of the upper plate (positively buoyant, model

Fig. 4. Viscosity, strain rates and velocity distribution of basic model setup at ~ 22 m.y. Models shown are fixed upper plate (UPfix, A,D), free upper plate (UP, B,E) and fixed subducting plate (SUBfix, C,F). The configuration of plates at surface and depth varies largely in the setups. The white area in A and C is where the imposed velocity is zero.

84

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

Fig. 5. Motions of the basic three setups versus buoyancy. (A) Sinking velocity of slabs, measured at mid-mantle depth at the timestep of Fig. 4. Velocities are very comparable and scale with buoyancy (Δρ). The range bar in the UPfix model is due to the variable velocity measured as the slab folds at depth (Fig. 4A). (B) Surface horizontal velocities of subducting (red symbols) and upper plates (i.e. trench motions, blue symbols). Plate motions scale with buoyancy although the partitioning is completely different: 100% subducting plate motions in UPfix models and 100% upper plate motions in SUBfix models, as imposed by boundary conditions. Free upper plate model (UP) subduction motions are achieved by 40%, 50% and 60% trench retreats. The variations to sinking velocity in UPfix models are reflected into surface motions, expressed by a range bar.

UPrhoNeg). However, none of these properties influence how the subduction motions are partitioned.

3.3. Coupling at convergent margins 3.3.1. Stress regimes Large stresses are propagated from slabs to plates at the surface which are then transmitted across to upper plates. The stress is measured in both plates at the stage shown in Fig. 4A–C, after some portion of the slab rests atop the lower mantle. Due to its high viscosity, deformation in the upper plate is negligible and the stress within does not generally change significantly through time, except in the case when the upper plate is fixed due to the deep piling of slab, as shown in the previous section. The most important observation is that the stresses transmitted to subducting and upper plates are very similar in all the three upper plate models (Fig. 6A). This further corroborates the view expressed here that although plate motions at the surface can vary considerably, the underlying dynamics are the same.

Within the upper plate, compressive and extensional stresses occur simultaneously depending on proximity to the trench. Near the trench, henceforth the arc region, stresses are compressive with magnitudes of ~70 MPa, comparable with those in the subducting plate. The arc region is in compression in all models, regardless of the stress distribution in the rest of the upper plate. Further from the trench, typically around ~150–400 km from the trench depending upon the far-field boundary conditions, a sharp transition occurs. At distances greater than 300– 400 km from the trench, henceforth the back-arc region, the upper plate is in an extensional stress regime. Although trench migration rates vary considerably between the different reference models, from 0 (UPfix) to 5 cm/year (SUBfix), the stress distribution in the arc region is insensitive to such variations and remains at approximately the same magnitude of compression. This essentially results from the pattern of the mantle flow beneath the upper plates, which diverges at a distance of ~350 km, where tractions change from trench-directed in the arc to trenchopposed in the back-arc. Increasing the buoyancy of the subducting plate (through Δρ) provides a greater slab pull force, and such variations in the model with a fixed upper plate (UPfix) result in corresponding increases of compressional stress within the arc region of the upper plate (Fig. 6B). However, forces propagated in the far-field of the upper plate are almost insensitive to such variations in buoyancy. Of particular interest is the vertical stress above subduction zones, since this has a net effect on topography. Since our upper plates' properties (i.e. density and viscosity) are uniform, the vertical stress is of similar magnitude as the horizontal stress, after subtracting the buoyancy term of the plates (Fig. 7). The subducting plate buoyancy, for the values considered here, mainly influence the peak stress, which from 70 rises to 90 and 140 MPa for increasing buoyancy, whereas the width of the elevated area is almost the same. We have also investigated the role of upper plate thickness using a fixed upper plate thickness varying between 30, 50, and 100 km. In fact, larger upper plate thickness results in larger plate strength as well as a larger contact area between the plates. Vertical stress above the upper plate decreases with increasing thickness, whereas the areal extent with a non-negligible stress increases (Fig. 7B), i.e. plateau-like. The widening of the elevated area results from a larger total force generated from both increasing the buoyancy of the subducting plate and the thickness of the upper plate (Fig. 8). The trends of increasing force propagation versus either increasing buoyancy or thickness are both approximately linear. This indicates that the buoyancy of the downgoing plate and thickness of the upper plate are controlling factors for the vertical stress around trench zones, and modulate both the peak topography and width of the plateau formed above convergent margins. We conclude that (1) the stress above the trench zone does not depend on the motions of the margin, whether it is migrating or fixed. (2) Larger buoyancies of subducting plates increase the basal traction below upper plates and result in larger peak stresses within ~ 350– 400 km from the trench. (3) Thicker upper plates allow for the transmission of larger forces through an increased plate interface, increasing the areal extent over which large vertical stresses are distributed but decreasing the magnitude of the maximum stress. 3.3.2. Time evolution of coupling The time evolution of the upper plate models reflects the different constraints on plates at the surface and the evolution of the slab at depth. We have measured the integrated horizontal force in the subducting plates and in the upper plates, far from the trench. In all models, the slab pull force increases continuously reaching a maximum of ~9·1012 N/m when the slab tip extends to the bottom of the upper mantle (Fig. 9A, stage I). Slab pull is very sensitive to the folding at depth (Fig. 9A, stage II), which occurs in models of fixed upper plate (UPfix, Fig. 9A red line) and, to a lesser extent, models with both plates free (UP models, Fig. 9A black dashed line). In these models, the slab pull evolution is very similar, despite the fact that trench motions are completely different.

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

85

Fig. 6. Horizontal stress in the whole domain versus distance from trench. Stress is measured at the mid-depth of plates at the time of ~ 22 m.y. (Fig. 4). Positive is extension, negative is compression. In the bending area stresses are reoriented and not parallel to the horizontal. (A) Comparison of stress in the three basic setups, Δρ = 60 kg m− 3, and (B) stress for different buoyancy, 30, 60 and 90 kg/m3, in the free upper plate models (UP). Stress domains in the upper plate are subdivided into arc region (<375 km) and back-arc region.

The UPfix model experiences a longer period of folding and associated oscillation of slab pull, whereas the slab tip in the UP model simply readjusts when it reaches the bottom. The increasing piling and support of the slab in the former lead to a slight decay of the magnitude of slab pull force however it remains sufficiently large (~5·1012 N/m) to drive subduction during the time in which ~80% of the subducting plate is consumed. In contrast, when the subducting plate's tail is fixed (SUBfix), slab pull reaches a relatively constant maximum when the subducted slab extends to the length of the upper mantle (~1013 N/m). It then increases up to almost double due to the internal stretching of the subducting plate resulting from the slab pull applied on one end and a fixed trailing edge boundary condition on the other. Thus, in the SUBfix models, the advancing plate motion is resisted by horizontal stresses at the sidewall (equal but opposite direction) which causes larger inner deformations. A different evolution characterizes the force through the upper plate. The stress regime is always in compression around the trench, however, far from the trench, the regime switches to compressive during the initial stage of subduction (Fig. 9B) to extensional in the later stage (Fig. 9B blue and black dashed lines). The magnitude of the forces in the upper plate is around ~ 4.5·1012 N/m, i.e. half the slab pull

propagated to the subducting plate. The initial stage is similar in all the subduction types, and is mainly a consequence of a steeplydipping slab leaning forward during subduction (Fig. 9B, inset) and generating a forward horizontal velocity. Later on, during the rollback phase, the slab is leaning backward (Fig. 9B, inset), i.e. the subduction direction has a rearward horizontal velocity, and upper plate forces are around 3 · 1012 N/m. In UPfix models, there is no trench retreat, however, the piling of slab mass at the bottom of the mantle causes fluctuations of the slab dip, alternating between forward- and rearward-leaning slabs and their corresponding oscillations in the strain regime of the upper plate. As seen before, the force in the upper plate of this model is around twice that of the other models due to the fixed boundary conditions at the tail of the plate. The two plates are viscously coupled through corner flow in the mantle wedge and variations in slab pull forces arising from deep subduction dynamics are propagated to the upper plate and distributed throughout its entire length. Because the viscosity of the upper plate is sufficiently large to resist deformation, it is essentially passive. The upper plate simply responds to forces from the subducting plate by

Fig. 7. Vertical stress in the arc area versus distance from trench. (A) Stress for different slab buoyancy, 30, 60 and 90 kg/m3, in the free upper plate models (UP). (B) Stress for different upper plate thickness, 30 (blue), 50 (black) and 80 km (red).

Fig. 8. Vertical force integrated in the arc area for different densities (Δρ) and upper plate thickness (hUP) tested. The variable density models are those of a 50 km thick free upper plate (UP). The variable thickness models have density contrast of 60 kg/m3.

86

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

Fig. 9. Force evolution during subduction. (A) In-plane slab pull force propagated to the subducting plate, measured 500 km from the trench. Evolution is divided in two phases: pull increases until slab reaches the bottom of the mantle (I) and a phase in which morphology of the slab results in a different evolution (II). (B) In-plane force propagated to the back-arc area in the upper plate. While the pull increases until the slab tip reaches the bottom of the mantle (phase I), the upper plate back-arc is in compression. As the slab further subduct the free upper plate (UP and SUBfix models) is in extension, the fixed upper plate models (UPfix) experience periodic compression according to the folding at depth. Compression (C) and extension (E) in the back-arc are due to the horizontal component of the subduction vector (insets).

accommodating large tectonic stresses, and this allows estimating the magnitude of forces transmitted across from the subducting plate. In the next section we investigate the effect of a more realistic rheology which has a strength comparable to the stresses generated in these models. In those models, the deformation (i.e. strain) in the upper plate affects the evolution, resulting in more non-linear behaviour than seen here, but the plates are subjected to the same magnitude of forces measured here. The forces propagated to the subducting and upper plates depend on different aspects of subduction. The subducting plate is subjected to an in-plane force that depends only on the slab pull. However, the stress regime in the upper plate's back-arc region is determined by the horizontal component of the subduction vector, rather than the subduction velocity. According to the slab folding, this is either frontward or rearward, corresponding to inducing a compressive or extensional regime, respectively.

which can influence the evolution the system. The force evolution in the upper plate is similar for different thicknesses tested (30, 50 and 80 km) however the deformation here critically depends on the upper plate strength. Although the viscosity of the plate is constant, varying its thickness changes the way it resists stretching and compression, as determined by its resistance to stretching Sr∝ηh. Therefore an 80 km thick plate has ~2.6 times the resistance of a 30 km thick plate without considering any thickening or thinning that may occur. In Fig. 10 we show the upper plate velocity measured in the proximity of the trench. These are considered diagnostic of strain rates in the back-arc region, whereas in the near-trench (<~300 km) compression and thickening are constant in all the models at all stages (Fig. 11). The weaker, 30 km thick plate allows for the greatest trench motions (Fig. 10A, dash–dotted line). These are positive in the initial stage, meaning that upper plate stretches allowing initial trench rollback (Fig. 11A–B). However, as slabs accumulate at depth, the slab pull force decreases and upper plate stretching diminishes (Fig. 10A). During the folding at depth, (Fig. 11B–C) the forward-leaning slab induces compression in the upper plate (Figs. 10A, 11C). However, the mass of the slab pile is supported at depth gradually reducing the slab pull and the amplitudes of corresponding advancing motions and upper plate thickening decrease accordingly. As the overriding plate continues to thicken, the subducting plate transitions to rollback and the trench migration begins to drive the upper plate back into extension. A similar evolution occurs in the model with a weaker, 80 km thick upper plate (Figs. 10A, 11E–H). However the larger resistance to stretching of thicker lithosphere impedes trench motions. The periodic folding at depth is very similar to the previous case, but the nearly stationary trench allows for repeated upper plate compression and trench advance, without extension, more similar to the very stiff upper plate of the basic setup (UPfix model). The dynamic interaction between deep subduction dynamics and surface deformation is, thus, completely due to the ability of the upper plate to deform. The strong periodic character emerges from all the

3.4. Upper plate strain regimes In the following sections we describe how the time evolution of convergent margins is controlled by upper plate rheological properties. We first use a more realistic value of upper plate strength and then introduce viscoplasticity. This allows back-arc regions of upper plates to yield during periods of extension and compression, as well as formation of back-arc basins, which are often observed around convergent margins. All the models here have a fixed upper plate. This is representative of plates attached to continents or very large plates, and this is a more typical tectonic setting for subduction over the past 50 Ma on Earth. 3.4.1. Viscous upper plate The models presented here use the same setups of previous sections, however, the upper plates are much weaker as the viscosity is reduced to 300 ηUM (Table 2). In these models, again the coupling with the upper plate increases as slab pull develops and then modulates during slab folding at depth (Fig. 9). However the ability of the weaker plate to stretch allows for additional, albeit minor, trench motions to occur

Fig. 10. Evolution of upper plate velocities and back-arc deformation for different upper plate thicknesses and rheologies. Plate velocity is measured in the arc zone of the fixed upper plate. Positive motions indicate trench rollback/back-arc extension, and negative motions indicate trench advance/back-arc compression. Dots on curves mark the time steps of Figs. 11 and 12. (A) Velocities for purely viscous UPfix models with upper plate thicknesses of 30 km (dash–dotted line), 50 km (solid line), and 80 km (dashed line) and (B) velocities for viscoplastic UPfix models with upper plate thicknesses of 30 km (dash–dotted line), 50 km (solid line), and 80 km (dashed line).

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

87

Fig. 11. Viscosity and velocity distribution of viscous fixed upper plate models. Time steps correspond to the dots on curves in Fig. 10A. The viscosity color scale as in Fig. 1. Note the correspondence between the subduction vector horizontal component and the strain phase in Fig. 10. (A–D) 30 km thick upper plate and (E–H) 100 km thick upper plate. Zero velocity is denoted by white areas.

models, with the magnitude and direction of the motions proportional to the strength of the upper plates (Fig. 10A). While thinner upper plates allow larger retreat motions and subsequent extension, the thicker upper plates only allow for small trench advance velocities and corresponding upper plate compression, intervened by periods of quiescence. Furthermore, the strength of upper plates is important in the early stages of subduction. Thicker or thinner plates will thus result in upper plate compression or extension, respectively. 3.4.2. Viscoplastic upper plate These models are similar to the reference models with an initial (preyielded) viscosity of 3000 ηUM, however viscoplasticity is introduced. The amount of deformation to be taken up by plastic strain depends mostly upon the yield stress in the upper plate (Table 2). This leads to transient upper plate deformations that are, due to the nonlinearity of plasticity, strongly punctuated and episodic, leading eventually to dramatic switches in subduction modality. Fig. 10B shows the evolution of the velocity measured on the upper plate above the trench zone. The motions at the trench are used to indicate the magnitude of strain achieved in the back-arc region of the plates. As in the previous section, variations of upper plate thickness (30, 50 and 80 km) are investigated. Also similarly to the previous models, the integrated strength of the plate varies strongly with the thickness, as thicker plates readily accommodate larger stresses. The thickest upper plate model (80 km) has enough strength to withstand

the forces transmitted through the plate margins, and upper plate deformation is negligible with upper plate velocity close to zero (Fig. 10B, dashed line). In the case for a slightly less thick (50 km) upper plate, the integrated strength is about the same order of magnitude as the tectonic stresses generated by subduction. This results in episodic stretching, with periods of increased extensional deformation occurring between 8–10 m.y. and 15–19 m.y. (Fig. 10B, solid line). As the slab pull continuously increases while sinking towards the bottom of the upper mantle (dots on solid line in Fig. 10B and corresponding frames in 12E–F), the traction propagated to the upper plate exceeds the upper plate strength and plastic weakening occurs (Fig. 12F) resulting in a maximum upper plate velocity of ~ 2 cm/year at 9 m.y. (Fig. 10B). However, as soon as the slab is partly supported by the lower mantle (bottom boundary), slab pull decreases and the deformation in the upper plate subsides, corresponding as well to a period of forward-leaning of the slab (Fig. 12G). A renewed phase of stretching begins at 15 m.y. as the slab leans backwards. However, at this point the slab is partly supported by the slab pile below, and much smaller upper plate motions and deformations occur. The smaller resistance of the thinner (30 km) upper plate can be easily overcome (well exceeding the yield stress) early during the initial stage of slab subduction. This retards the trench motions until the upper plate breaks apart and a small plate can follow the trench, freely retreating (Fig. 12A–D). The subduction modality following this

Fig. 12. Viscosity and velocity distribution of viscoplastic fixed upper plate models. Time steps correspond to the dots on curves in Fig. 10B. Viscosity color scale as in Fig. 1. (A–D) 30 km thick upper plate and (E–H) 50 km thick upper plate. Zero velocity is denoted by white areas.

88

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

event switches to the free upper plate behaviour as seen before, with rollback velocity comparable to the UP models (Fig. 10B). In contrast to the uniform stretching observed for weaker, purely viscous plates, the plastic strain (yielding) is localized at a particular distance from the trench (~300 km), eventually isolating small upper plates of that length. As observed in previous models, this distance represents the transition of stress regimes along the upper plate from always compressional nearby to the trench to either compression or extension (depending on model specifics) in the back-arc region (Fig. 6). This transition is generated from the change in direction from applied basal tractions (trenchwards above the mantle wedge) which localizes deformation and allows spontaneous formation of back-arc extension/ spreading centres at a distinct location with respect to the trench (Fig. 12). 3.4.3. Plate motions As slab pull develops and evolves, it is important to understand how the motions of both plates, which sum up to the convergence rate, respond. During the initial stage of subduction in models with weaker, purely viscous upper plates, subducting plate velocities (Fig. 13A, red lines) are much larger than that of the overriding plates (Fig. 13A, black lines). However, as the slab piles at depth, subduction velocities decrease by about 50%. In general, the velocities of the upper plates are much less than the subduction plates. The deformation rate of the thicker upper plate is negligible with respect to the subduction velocity. Thicker plates respond to the forcing by thickening, whereas thinner plates exhibit mainly extension in the back-arc region. Subducting plate velocities evolve in a similar way in all models, having the same magnitude and only minor differences in the periodicity of the fluctuation. However, similar fluctuations in the upper plate velocities evolve over a different timescale. This is essentially because the subducting plate advance is controlled only by slab pull (and its variations), whereas the upper plate motions are also sensitive to the forwards or backwards leaning direction of the slab (Fig. 9). The deformation rate of the thicker upper plate is negligible

with respect to the subduction velocity. Thicker plates respond to the forcing by thickening, whereas thinner plates exhibit mainly extension in the back-arc region. In the models with a viscoplastic overriding plate, the evolution of subducting plate velocities is similar to the previous models (Fig. 13B). The ability of these plates to fail (yield plastically) only causes significant differences in the 30 km thick upper plates. In this case, upper plate motions are very large and convergence is achieved by large retreat (40–50%). In models with thicker upper plates, trench motions are negligible although periodic thickening and stretching does occur. The 50 km thick upper plate model has an intermediate behaviour, and episodic stretching occurs at rates of ~ 10 times smaller than subduction velocity. It is important to note that the velocity of subducting plates does not indicate directly the state of strain in the upper plate, whether compressive or extensive, although the strain amplitude does. 3.5. Summary of results Our models show that the evolution of convergent margins responds to the feedback between deep subduction dynamics and rheological controls on upper plate deformations. We obtain consistent deformations in the upper plates that include compression in the arc zone, i.e. at distance <~350 km from the trench, and coexist with the extensional or compressive stress regimes in the upper plate's back-arc zone. During its evolution the preferred mode of slab subduction is by trench retreat, however the presence of a fixed upper plate can significantly influence the trench migrations, depending on its mechanical properties. More specifically, a mechanically strong, fixed upper plate prohibits trench retreat, and results in a stationary trench with the slab piling at the base of the upper mantle. This slab folding influences the subducting velocity and slab pull, varying with the piling, but does not directly affect the state of stress of the upper plate. The latter rather depends on the leaning directions of the slab while piling at depth, with a forwardleaning slab inducing an upper plate compressive stress regime, and a vertical/rearward leaning slab generating extensional stresses in the upper plate. The interaction of the force propagated through the margin with the upper plate produces strain regimes that are strongly controlled by plate rheology, whether viscous or viscoplastic, and its integrated strength (which depends on plate thickness). Our models reproduce non-steady-state subduction with upper plate deformation phases that are either (1) periodic when rheology is linear (viscous) or (2) episodic, when upper plate rheology is non-linear (plastic). The strength of the plate, defined by its rheology and thickness, limits the magnitude of extension/compression regimes and retreating/advancing trench motions. 4. Discussion 4.1. Comparison with previous models

Fig. 13. Evolution of upper and subducting plate velocities for different upper plate thickness and rheology. Upper plate velocities as in Fig. 10. Plate velocity is measured in the arc zone of the fixed upper plate (vUP, black lines), and subducting velocity (vSP, red lines) is measured 500 km form the trench. Positive motions indicate trench rollback/ back-arc extension, and negative motions indicate trench advance/back-arc compression. Subducting plate velocities are positive trenchward. (A) Viscous upper plate velocities and (B) viscoplastic upper plate.

The forces transmitted through the convergent margins have been a subject of lively debate and numerous studies. Thus far, no integrated model of subduction dynamics with an upper plate has been proposed. Major differences reside in the treatment of the bending area and the implementation of the margin contact area. The recent improvements in understanding the role of the bending at the trench (Capitanio et al., 2007; Capitanio et al., 2009a) allow models to reproduce self-consistent, fully dynamic subduction (Stegman et al., 2010-this issue) where plate and trench motions, dips and stresses in the plates more closely resemble those observed on Earth (Capitanio et al., 2009a; Goes et al., 2008). The ability of Underworld to track material properties allows a layer of plastic crust to be included above the higher viscosity lithospheric core and allows strain to localize in the layer being wrenched at the convergent margin. We do not investigate the role of friction at the

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

interface between plates in this work, as we are more interested in the first order force balance at margins. However, very low friction interfaces and low yield stress (~30–50 MPa), as those used here, are appropriate over the long timescales (~50 Ma) and facilitate more realistic subduction in dynamic models, as constrained by geoid and topography (Krien and Fleitout, 2008; Moresi and Gurnis, 1996; Toth and Gurnis, 1998). Another important limitation of our model is the choice of the 2D model space and lack of a free surface. Stress propagation through the contact is not affected by this approximation, although the tractions in the mantle wedge might be. The results presented are more applicable to a wide plate (see Dvorkin et al., 1993), where flow in the center, far from edges, is essentially poloidal. Furthermore, toroidal flow around slab edges acts to reduce the tractions to the upper plate arising from corner flow circulation (which is the only flow present in 2D models). However, if the mantle wedge or uppermost mantle has a nonNewtonian rheology, the toroidal flow will be more concentrated locally near slab edges and a greater lateral extent of the plate will have an approximately 2D circulation. The free upper surface might imply numerical difficulties due to entrainment at trenches (Schmeling et al., 2008); however, in problems such as the one treated here, the entrainment of crust, which is negatively buoyant, is a feature of the model. Furthermore, Kaus et al. (2008), showed that in a setup similar to ours the deeper lithosphere dynamics is not constrained at all when upper layers are viscoplastic. The force evolution of our models, in the subducting and the upper plates, is similar to that of Gurnis et al. (2004), as well as the horizontal stress distribution at the surface across the trench zone that finds a good agreement with the results of De Franco et al. (2007), corroborating the outcomes of our work. Laboratory models (Faccenna et al., 2001) already showed the rearrangement of pull force as the slab interacts with the 660 km discontinuity, as found here, although no upper plate is included in these models. Furthermore, analog models including an upper plate have reproduced intermittent back-arc compression behaviour due to the folding (Boutelier and Cruden, 2008; Guillame et al., 2009) as well as arc thickening (Espurt et al., 2008). However these models are strongly controlled by the velocity imposed at plate boundaries, so only periodic back-arc compression has been reproduced. The setup we present here is similar to that of Arcay et al. (2008). However a comparison is difficult because in the case of Arcay et al. (2008), forces and dissipation are strong functions of the preimposed and constant velocities. In our models the forces are approximately constant, as these are imposed through buoyancy, as well as the dissipation partitioning. 4.2. Upper plate stress and strain regimes Several published works have addressed some aspects of the models presented here usinga analytical, numerical and laboratory modelling. The mantle wedge tractions are well-known to determine the stress state in upper plates (Sleep and Toksöz, 1971), and generate thickening of the upper plate (Wdowinski et al., 1989). The models presented here achieve a strong (viscous) coupling between the downgoing slab and overriding plate through the corner flow circulation in the mantle wedge. As previously stated, the strength of this coupling is possibly overestimated due to the choice of adopting a uniform, linear viscosity throughout the upper mantle. However, the key aspect of this investigation is to understand how the forces of the subducting plate are expressed in upper plate deformation. In particular, the model setup adopted here has several advantages over previous efforts (Clarke et al., 2008) in which a weak arc zone between the upper plate and subducting plate accommodated most of the deformation. The presence of a weak arc zone in Clarke et al., 2008, has the undesirable effect of decoupling the plates. As a consequence, it substantially reduces stress transfer across to the

89

upper plate. It also allowed the subducting plate to rapidly form into a steeply-dipping (nearly vertical) slab, greatly enlarging the mantle wedge and causing the viscous coupling between the slab and upper plate arising from corner flow to become negligible. The models presented here overcome these challenges and provide mantle wedge geometry (with its associated and important corner flow) throughout the entire time evolution, which adjusts self-consistently according to the subduction dynamics. An interesting outcome of our model is that the compression above trench zone is the same in the three different reference models irrespective of whether the trenches are stationary or retreating. The viscous coupling between the downgoing plate and the upper plate varies entirely with the sinking velocity of the slab, which is controlled by the subducting plate buoyancy (i.e. the thermal age). This is in agreement with the analysis of Peterson and Seno, 1984, which demonstrates that plate age provides the strongest correlation with seismic moment release rates at convergent margins amongst a range of subduction parameters. Secondly, the seismic moment release rates of the plate-bounding fault are also found to correlate with upper plate velocity (Peterson and Seno, 1984; Scholz and Campos, 1995), although to a lesser extent. The correlation with upper plate motions is likely smaller because of the dominant effect of far-field boundary conditions on upper plate motions. In this case, the subducting plate age gives the best indication of the strength of coupling between the slab and upper plate, in agreement with the major conclusions of Peterson and Seno, 1984. This is because it is through the sinking velocity of the slab that mantle wedge circulation and the resultant tractions along the base of the upper plate are generated. However, this is only a first order approximation as the mantle wedge flow also depends strongly on the extent of wedge hydration (which depends on slab dip angle), as well as the non-Newtonian rheology. Although in the models without an upper plate the trench motions are always slab-driven, our models show that only rarely during the evolution of a subduction zone are all motions (the upper plate and trench) dynamically correlated to subducting plate velocity. This is reflected in the poor statistical correlation of upper plate velocities versus subduction (R2 < 0.37) (Lallemand et al., 2008), suggesting that only a few subduction zones might meet the conditions of the lab. An example might possibly be the Tonga subduction zone, where a back-arc spreading center allows for the fastest retreating trench motions on Earth, comparable with laboratory models. Another case is the Western Mediterranean subduction zone, where there are two stages of back-arc opening, one with associated oceanic crust and one with extremely thinned continental crust. Apparently, the upper plate did not resist the retreating motions of the free falling slab (Faccenna et al., 2001), in agreement with our model results. The compressive stress found in the arc region of the upper plate is independent of the regime present in the back-arc region. The first depends mainly on the coupling at the trench and the traction in the mantle wedge, whereas the second is more sensitive to the upper plate rheology and the direction of slab subduction in the mantle underneath. These two mechanisms are independent, and can therefore occur in any combination including simultaneous compression and extension. For instance, in the Mediterranean, thrusting in the Apennines and extension in the back-arc coexisted during the whole evolution of subduction (Malinverno and Ryan, 1986). Similarly, this model applies to the Pannonian basin-Carpathian arc and other settings ((Malinverno and Ryan, 1986) and ref. therein). The control on the upper plate regime exerted by trench motions has been long recognized (Uyeda and Kanamori, 1979); however, we showed here that these conditions change throughout the evolution of subduction. The models presented here clearly demonstrate that extension, quiescence or compression may arise in the back-arc region while the trench is continuously retreating, depending upon the deep dynamics. Furthermore, if the upper plate strongly resists deformation, the trench remains nearly stationary regardless of whether the back-arc region experiences compression or extension,

90

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

and in such cases the stress regime cannot be understood from trench motion alone. We suggest that the subduction vector, possibly defined by surface motions, and its angle with the dip of the slab (i.e. Boutelier and Cruden, 2008; Guillame et al., 2009), potentially offers a more diagnostic tool to understand the non-steady-state subduction. The upper plate thickness also has an impact on the total force propagated to upper plates and influences the dynamics in two ways. First, a thick overriding lithosphere reduces the wedge volume, resulting in increased corner flow circulation and larger tractions under the upper plate. Secondly, thicker upper plates are also stronger, and therefore distribute a given magnitude of compressive stress over a larger area, producing a wider, but less elevated area of thickened plate, i.e. plateaulike. This is in agreement with the results of Husson and Ricard, 2004, where the difference along the Andean topography, cordillera versus plateau, can be fitted by a combination of basal traction, in-plane compression and internal stress due to heterogeneous density distribution. However, these authors adopt a constant plate thickness, and hence strength, leading to a somewhat counterintuitive result in which lower stresses are required to build the Altiplano–Puna plateau, compared with the cordilleras in Peru and South Chile Andes. Alternatively, one can adjust the magnitude of the resistive force at the plate boundary to balance the forces arising from gravitational spreading of topographically elevated regions and the presumed tractions based on global mantle circulation models (Iaffaldano et al., 2006). Similarly, we find that a larger force is required to sustain topographically high plateaus, but we suggest the origin of this force, and inherent deformations, can arise selfconsistently due to the natural variation in the thickness of the South American continental lithosphere, in agreement with Allmendinger and Gubbels, 1996. The main difference is that the surface tractions in global mantle circulation models (which are then applied to elastic shells that represent plates) result from double-sided downwellings (zones of convergence) in contrast to the single-sided subduction modeled here in which the motions of the trench, upper plate and subducting plate are all emergent properties of the model. Thus, if the models presented here are more appropriate for convergent margin dynamics, then the upper plate

plays the most significant role in the evolution of the Andes (Sobolev and Babeyko, 2005). 4.3. Intermittent plate motions and deformation at Earth's convergent margins Recent compilations of Cenozoic motions and oceanic back-arc basin evolution showed a strong intermittent behaviour in subducting and upper plate migrations as well as basin formation (Clarke et al., 2008; Sdrolias and Müller, 2006). With the understanding provided in the previous section we now analyse these features trying to recognize the evolution characteristic and how these can be reconciled with the dynamics showed here. Cenozoic motions of major Pacific–Indian Ocean convergent margins, suggest that subduction on Earth is often far from the steady-state behaviour indicated by the laboratory models. Both subducting plate and upper plate motions have strong fluctuations, like those found here. In Fig. 14 we show the absolute velocity of subducting (vSP) and upper plates (vUP) from the database of Sdrolias and Müller, 2006. The motions of overriding plates account for back-arc basin spreading and thus represent the trench motions. These appear to be rather independent of downgoing plate motions. In the case of the Andes, Cascadia, Central America and Aleutians, upper plate motions are almost constant, whereas the downgoing plate rates vary considerably, as in the case of Farallon slab detachment (Grand et al., 1997) conferring a dramatic rearrangement of Pacific subducting plate velocities below Cascadia and Central America (Fig. 14F, G), without altering the rate of the upper plate. On the other hand, subduction zone motions such as Japan–Kuril, Izu–Bonin–Mariana, Tonga–Kermadec and Sumatra–Java, show the opposite behaviour. In these zones, downgoing plate rates are rather smooth, with no large variations, whereas the upper plate motions are affected by episodic and/or periodic intermittent evolution. In all these zones, the intermittent back-arc deformation phase elapses 15 to 20 m.y. with variable trench motions ~2 to 5 times smaller than subducting plate rates, in agreement with our models. The motions of

Fig. 14. Subducting plate (vSP, red lines) and upper plate velocities (vUP, black lines) of major subduction zones in the Pacific and Indian Oceans (Sdrolias and Müller, 2006). In grey, major deformation events and names of associated basins, from Folguera et al., 2002; Grand et al., 1997; Jolivet et al., 1994; Schellart et al., 2006; and Sdrolias and Müller, 2006. Solid grey is for extension and hatched for compression.

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

the upper plate are comparable with subducting plate motions only when back-arc basins are spreading (or the upper plate has been strongly attenuated), for example in the Tonga–Kermadec back-arc system, Japan Sea, or Wharton basin where trench retreat rates are larger than ~5 cm/year. In other cases, when the upper plate is thicker (i.e. the Andes) or fixed by far-field boundary conditions such as being attached to a large continent, the partitioning between plate motions increases, in agreement with the models presented here. The motion of the subducting plate alone cannot predict the evolution of the margin, except when the upper plate and trench are both stationary in which case the subducting plate velocity and convergence rate are one and the same. Otherwise, subducting plate motions can be constant but the convergence rate can be essentially episodic with discrete periods of back-arc basin opening (i.e. Tonga– Kermadec, Japan–Kuril). The most prominent feature of the Japan–Kuril margin's Cenozoic evolution is the opening of the Japan Sea from ~28 to 12 Ma (Jolivet et al., 1994). This follows a phase of relatively stationary trench, during which the slab piled at depth and eventually sank into the lower mantle (Goes et al., 2008). Large strike–slip motions characterized the seafloor spreading of the Japan Sea (Jolivet et al., 1994), which likely helped avoid significant stretching of the continental margins. We thus suggest that an enhanced slab pull due to the sinking slab pile might have been responsible for upper plate yielding, followed by trench retreat and a flat-lying slab at depth (Fukao et al., 2001). Subsequently, the deformation at Japan margin turned into compressive during the last ~7 Ma (Jolivet et al., 1994). Depending on the choice of reference frame, trench motions are either slowly advancing or retreating (Funiciello et al., 2008), while subducting plate speed increases again. A dynamic correlation with the downgoing plate is difficult to see, as periodically during the Cenozoic the subducting plate motions decreased while upper plate motions increased, and vice-versa, summing up to a roughly constant convergence (Sdrolias and Müller, 2006). It is therefore difficult to establish a correlation between the motions and the dynamics of convergence, unless assuming that the underlying driving mechanisms are independent, as shown in our models. The Tonga–Kermadec subduction zone recorded compression in the Norfolk basin (38–34 Ma) subsequent to its initiation (~ 48 Ma) (DiCaprio et al., 2009; Sdrolias and Müller, 2006). This event is rapidly followed by the opening of the South Fiji basin (34–25 Ma) and after a period of quiescence, by the Lau basin behind the Tonga arc (7–0 Ma) ((Sdrolias and Müller, 2006) and ref. therein). Strong change in strain regimes during the early stage of subduction and following intermittent deformation characterizes this subduction zone. A detailed plate reconstruction of this area, with tomographic interpretation of subducted slab and estimates of the amount of consumed slab (Schellart et al., 2006), helps unravel the regional tectonics. The initial stage of compression occurred ~10–14 m.y. after subduction initiated (DiCaprio et al., 2009; Gurnis et al., 2004) and subsequent transition into an extensional strain regime are consistent with subduction models which have a strong upper plate. The pattern of advancing motions preceding back-arc spreading is commonly found in all the subduction zones studied by Sdrolias and Müller, 2006. This is consistent with the upper plate strain transitions during slab piling as seen in the models presented here. This stage of subudction is accompanied with very slow trench motions for which piling atop the upper-lower mantle interface can be inferred (Schellart et al., 2006). The penetration of the lower mantle by folded slab material may provide the extra slab pull needed to induce upper plate failure and opening of the large South Fiji basin. Upper plate failure resulted in back-arc spreading and allowed the Tonga–Kermadec trench to rollback, shallowing the slab's dip angle, similar to that suggested for the Japan Sea case. The resistance offered by a strong upper plate reduces rollback motions to rates below what slab sinking would otherwise produce. This leads to a second episode of this scenario, with another slab pile flushing into the lower mantle, and finishing with the Lau basin opening (Schellart et al., 2006). This

91

evolution is very similar to our models of a thinner viscous upper plate constraining the trench motions during this time with strong feedback with deep interactions. In this case the second phase of piling during the Lau basin formation probably happened during the rollback phase between ~20 and 10 Ma, as the lower trench motion recorded in this time is more compatible with those of the folding of slab. The upper plate deformation evolution of this subduction zone can be classified as periodic, essentially explained by a thin (weak) upper plate. There has been a substantial amount of previous effort directed towards understanding the development of the Andes, and it is simply beyond the scope of this work to comprehensively cover this subject. However, we would like to briefly address the nature of recurrent upper plate deformations in the Andes. Here compression and extension have been recorded throughout the Meso–Cenozoic evolution of the chain (Glodny et al., 2006; Oncken et al., 2006) and episodes of slab shallowing and deepening have been associated with the extension and compression, respectively (Folguera et al., 2002). However a clear correlation with the convergence motions has not been found, and some have associated the Oligocene–Miocene extension with an increase of convergent motions (Jordan et al., 2001), whereas others relate Oligocene extensional phases with a decrease in convergence (Charrier et al., 2002). The correlation between the deep slab folding, the dip and episodicity of strain in the Andes has already been investigated by several laboratory models (Espurt et al., 2008; Guillame et al., 2009; Heuret et al., 2007), the results of which are essentially similar to what is found here, although motions are completely prescribed. From the point of view of the dynamics we study here, a very thick, stiff upper plate hampers almost completely the trench motions, but still undergoes phases of compression and extension responding to the interaction of the slab and lower mantle. In fact, the low ratio of upper versus subducting plate velocities suggests that this interaction and continuous deeper penetration occurred throughout the Cenozoic time (Goes et al., 2008), compatible with the alternating phases of shallow slab/extension and steeper slab/compression found here. Indeed, a correlation between the building of Andes with the variation of shear stress at the subducting margin can be established (Lamb and Davis, 2003), as this relates to the dynamics of subduction, however, we favour the explanation given here, where the upper plate deformations are the result of the balance of subduction forces, not only those acting near the surface. This will be an interesting area of further work. 5. Conclusions Numerical models of subduction systems with an overriding plate allow us to gain new insights into the dynamics of convergent margins and their control on the evolution of plate tectonics on Earth. The subduction of a single downgoing plate is strongly dominated by the slab pull, with a typical steady-state evolution. We found the same dynamics when the upper plate interacts with subducting plates, with negligible differences in the stress patterns. However, this simple dynamic setting allows for very different plate motions and deformation styles at the surface, often showing non-stationary evolutions. The periodic slab folding atop the deep mantle transition zone confers similar non-stationary evolution of stress regimes in the upper plate, whereas the deformations achieved here respond to the mechanical strength of the upper plates, i.e. their rheology and thickness. Both compression and extension occur as a result of thicker and thinner upper plates, respectively. These might evolve into periodic back-arc deformations, when upper plates accommodate strain viscously, or into strongly punctuated, episodic strain, where oceanic basin might form, if plastic strain localizes. These same characteristics are found in the Cenozoic motions and deformations of major subduction zones, suggesting that the interactions of deep dynamics and the upper plate constraints have a strong influence on the tectonics. We conclude that the upper plates are an integral part of

92

F.A. Capitanio et al. / Tectonophysics 483 (2010) 80–92

the subduction dynamic system, and their impact upon the diversification of tectonic styles at plate margins cannot be neglected. Acknowledgments We would like to thank Maria Seton for kindly providing the dataset used in this paper. This research was supported under the Australian Research Council's Discovery Projects funding scheme (project numbers DP0663258, DP0878501, DP0666472 and DP0987374). References Allmendinger, R.W., Gubbels, T., 1996. Pure and simple shear plateau uplift, Altiplano– Puna, Argentina and Bolivia. Tectonophys. 259, 1–13. Allmendinger, R.W., Jordan, T.E., Kay, S.M., Isacks, B.L., 1997. The evolution of the Altiplano–Puna plateau of the central Andes. Ann. Rev. Earth Plan. Sc. 25, 139–174. Arcay, D., Tric, E., Doin, M.-P., 2005. Numerical simulations of subduction zones — effect of slab dehydration on the mantle wedge dynamics. Phys. Earth Plan. Int. 149 (1–2), 133–153. Arcay, D., Lallemand, S., Doin, M.P., 2008. Back-arc strain in subduction zones: statistical observations versus numerical modeling. Geoc. Geophy. Geosys. 9 (5). doi:10.1029/ 2007GC001875. Becker, T.W., O'Connell, R.J., 2001. Predicting plate velocities with mantle circulation models. Geoc. Geophy. Geosys. 2. Boutelier, D., Cruden, A.R., 2008. Impact of regional mantle flow on subducting plate geometry and interplate stress: insights from physical modelling. Geophys. J. Int. 174, 719–732. 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.A., Morra, G., Goes, S., 2007. Dynamic models of downgoing plate-buoyancy driven subduction: subduction motions and energy dissipation. Earth Plan. Sci. Let. 262, 284–297. Capitanio, F.A., Morra, G., Goes, S., 2009a. Dynamics of plate bending at the trench and slab–plate coupling. Geochem. Geophys. Geosyst. 10 (Q04002). doi:10.1029/ 2008GC002348. Capitanio, F.A., Faccenna, C., Funiciello, R., 2009b. Opening of Sirte Basin: result of slab avalanche? Earth Plan. Sci. Let. doi:10.1016/j.epsl.2009.06.019. Charrier, R., et al., 2002. Evidence for Cenozoic extensional basin development and tectonic inversion south of the flat-slab segment, southern Central Andes, Chile (33˚–36˚S.L.). J. South Am. Earth Sci. 15, 117–139. Christensen, U.R., Yuen, D.A., 1984. The interaction of a subducting lithosphere slab with a chemical or phase boundary. J. Geophys. Res. 89 (B6), 4389–4402. Clarke, S.R., Stegman, D.R., Müller, R.D., 2008. Episodicity in back-arc tectonic regimes. Phys. Earth Plan. Int. 171, 265–279. Conrad, C.P., Lithgow-Bertelloni, C., 2002. How mantle slabs drive plate tectonics. Science 298 (5591), 207–209. De Franco, R., Govers, R., Wortel, M.J.R., 2007. Numerical comparison of different convergent plate contacts: subduction channel and subduction fault. Geophy. J. Int. 171 (1), 435–450. DiCaprio, L., Müller, R.D., Gurnis, M., Goncharov, A., 2009. Linking active margin dynamics to overriding plate deformation: synthesizing geophysical images with geological data from the Norfolk Basin. Geoc. Geophy. Geosys. 10 (1). doi:10.1029/2008GC002222. Dvorkin, J., Nur, A., Mavko, G., Ben-Avraham, Z., 1993. Narrow subducting slabs and the origin of back arc basins. Tectonophys. 227, 63–79. Enns, A., Becker, T.W., Schmeling, H., 2005. The dynamics of subduction and trench migration for viscosity stratification. Geophy. J. Int. 160 (2), 761–769. Espurt, N., et al., 2008. Flat subduction dynamics and deformation of the South American plate: insights from analog modeling. Tectonics 27. doi:10.1029/2007TC002175. Faccenna, C., Funiciello, F., Giardini, D., Lucente, P., 2001. Episodic back-arc extension during restricted mantle convection in the Central Mediterranean. Earth Plan. Sci. Let. 187, 105–116. Folguera, A., Ramos, V.A., Melnick, D., 2002. Strain partitioning along the arc in the Neuquen Andes (36–39 degrees S) for the last 30 million years. Re. Geol. de Chile 29 (2), 227–240. Forsyth, D., Uyeda, S., 1975. On the relative importance of driving forces of plate motion. Geophys. J. R. Astr. Soc. 43, 163–200. Fukao, Y., Widiyantoro, S., Obayashi, M., 2001. Stagnant slabs in the upper and lower mantle transition region. Rev. Geophys. 39 (3), 291–323. Funiciello, F., et al., 2008. Trench migration, net rotation and slab–mantle coupling. Earth Plan. Sci. Let. doi:10.1016/j.epsl.2008.04.006. Garfunkel, Z., Anderson, C.A., Schubert, G., 1986. Mantle circulation and the lateral migration of subducted slab. J. Geophys. Res. 91 (B7), 7205–7223. Gerya, T., Stockhert, B., Perchuk, A.L., 2002. Exhumation of high-pressure metamorphic rocks in a subduction channel: a numerical simulation. Tectonics 21 (6). doi:10.1029/2002TC001406. Glodny, J., et al., 2006. Long-Term Geological Evolution and Mass-Flow Balance of the South-Central Andes. Frontiers in Earth Sciences, III. Springer, Heidelberg. Goes, S., Capitanio, F.A., Morra, G., 2008. Evidence of lower mantle slab penetration phases in plate motions. Nature 451. doi:10.1038/nature06691. Grand, S.P., van der Hilst, R.D., Widiyantoro, S., 1997. High resolution global tomography : a snapshot of convection in the Earth. Geol. Soc. Am. Today 7 (4), 1–7. Guillame, B., Martinod, J., Espurt, N., 2009. Variations of slab dip and overriding plate tectonics during subduction: insights from analogue modelling. Tectonophys. 463, 167–174. Gurnis, M., Hall, C.E., Lavier, L.L., 2004. Evolving force balance during incipient subduction. Geoc. Geophy. Geosys. 5 (7). doi:10.1029/2003GC000681.

Hall, C.E., Gurnis, M., Sdrolias, M., Lavier, L.L., Müller, R.D., 2003. Catastrophic initiation of subduction following forced convergence across fracture zones. Earth Plan. Sci. Let. 212 (1–2), 15–30. Hampel, A., Pfiffner, A., 2006. Relative importance of trenchward upper plate motion and friction along the plate interface for the topographic evolution of subductionrelated mountain belts. Geol. Soc. London Spec. Publ. 253, 105–115. 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 Plan. Sci. Let. 256, 473–483. Husson, L., Ricard, Y., 2004. Stress balance above subduction: application to the Andes. Earth Plan. Sci. Let. 222, 1037–1050. Iaffaldano, G., Bunge, H.P., Dixon, T.H., 2006. Feedback between mountain belt growth and plate convergence. Geology 34 (10), 893–896. Jolivet, L., Tamaki, K., Fournier, M., 1994. Japan Sea, opening history and mechanism: a synthesis. J. Geophys. Res. 99 (B11), 22237–22259. Jordan, T.E., et al., 2001. Extension and basin formation in the southern Andes caused by increased convergence rate: a mid-Cenozoic trigger for the Andes. Tectonics 20 (3), 308–324. Kameyama, M., Yuen, D.A., Karato, S.I., 1999. Thermal–mechanical effects of lowtemperature plasticity (the Peierls mechanism) on the deformation of a viscoelastic shear zone. Earth Plan. Sci. Let. 168, 159–172. Karato, S.I., Wu, P., 1993. Rheology of the upper mantle — a synthesis. Science 260 (5109), 771–778. Kaus, B.J.P., Steedman, C., Becker, T.W., 2008. From passive continental margin to mountain belt: insights from analytical and numerical models and application to Taiwan. Phys. Earth Plan. Int. 171, 235–251. Kincaid, C., Olson, P., 1987. An experimental study of subduction and slab migration. J. Geophys. Res. 92 (B13) 13,832–13,840. Krien, Y., Fleitout, L., 2008. Gravity above subduction zones and forces controlling plate motions. J. Geophys. Res. 113. doi:10.1029/2007JB005270. Lallemand, S., Heuret, A., Faccenna, C., Funiciello, F., 2008. Subduction dynamics as 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. Malinverno, A., Ryan, W.B.F., 1986. Extension in the Tyrrhenian Sea and shortening in the Apennines as result of arc migration driven by sinking of the lithosphere. Tectonics 5 (2), 227–245. Mason, W., Moresi, L., Betts, P.G., Miller, M.S., 2010. Three-dimensional numerical models of the influence of a buoyant oceanic plateau on subduction zones. Tectonophys. 483, 71–79 (this issue). Moresi, L., Gurnis, M., 1996. Constraints on the lateral strength of slabs from threedimensional dynamic flow models. Earth Plan. Sci. Let. 138, 15–28. Moresi, L., Dufour, F., Mülhaus, H.B., 2003. A lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. J. Comput. Phys. 184, 476–497. Oncken, O., et al., 2006. Deformation of the Central Andean Upper Plate System — Facts, Fiction, and Constraints for Plateau Models. The Andes. Frontiers in Earth Sciences, Springer, Berlin Heidelberg. Peterson, E.T., Seno, T., 1984. Factors affecting seismic moment release rates in subduction zones. J. Geophys. Res. 89 (NB12), 233–248. Schellart, W.P., Lister, G., Toy, V.G., 2006. A Late Cretaceous and Cenozoic reconstruction of the Southwest Pacific region: tectonics controlled by subduction and slab rollback processes. Earth-Sci. Rev. 76, 191–233. 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. doi:10.1038/nature05615. Schmeling, H., et al., 2008. A benchmark comparison of spontaneous subduction models— towards a free surface. Phys. Earth Plan. Int. 171. doi:10.1016/j.pepi.2008.06.028. Scholz, C.H., Campos, J., 1995. On the mechanism of seismic decoupling and back arc spreading at subduction zones. J. Geophys. Res. 100 (b11) 22,103–22,115. Sdrolias, M., Müller, R.D., 2006. Controls on back-arc basin formation. Geoc. Geophy. Geosys. 7 (4). doi:10.1029/2005GC001090. Shemenda, A.I., 1993. Subduction of the lithosphere and back arc dynamics: insights from physical modeling. J. Geophys. Res. 98 (B9), 16167–16185. Sleep, N.H., Toksöz, M.F., 1971. Evolution of marginal basins. Nature 233 (5321), 548–550. Sobolev, S.V., Babeyko, A.Y., 2005. What drives orogeny in the Andes? Geology 33 (8), 617–620. 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. Geoc. Geophy. Geosys. 7 (3), 1–22. Stegman, D.R., Farrington, R., Capitanio, F.A., Schellart, W.P., 2010. A Regime Diagram for Subduction Styles from 3-D Numerical Models of Free Subduction. Tectonophys 483, 29–45 (this issue). Toth, J., Gurnis, M., 1998. Dynamics of subduction initiation at preexisting fault zones. J. Geophys. Res. 103 (B8), 18053–18067. Uyeda, S., Kanamori, H., 1979. Back-arc opening and the mode of subduction. J. Geophys. Res. 84 (B3), 1049–1061. van der Hilst, R.D., Seno, T., 1993. Effects of relative plate motion on the deep structure and penetration depth of slabs below the Izu–Bonin and Mariana island arcs. Earth Plan. Sci. Let. 120, 395–407. van Hunen, J., van der Berg, A.P., Vlaar, N.J., 2002. On the role of subducting oceanic plateaus in the development of shallow flat subduction. Tectonophys. 352 (3–4), 317–333. Wdowinski, S., O'Connell, R.J., England, P., 1989. A continuum model of continental deformation above subduction zones: application to the Andes and the Aegean. J. Geophys. Res. 94 (B8), 10331–10346. Zhong, S., Gurnis, M., 1995. Mantle convection with plates and mobile, faulted plate margins. Science 267, 838–843.