Numerical models of trench migration in continental collision zones

Numerical models of trench migration in continental collision zones V. Magni, J. van Hunen, F. Funiciello, and C. Faccenna Università “Roma Tre”, Dipartimento di Scienze Geologiche, L.go S. Leonardo Murialdo 1, 00146, Rome, Italy “Alma Mater Studiorum” Università di Bologna, Dipartimento di Geofisica, V. le Berti Pichat 8, 40127, Bologna, Italy Durham University, Department of Earth Sciences, Durham DH1 3LE, UK


Introduction
Slab pull caused by the negative buoyancy of subducting lithosphere represents a main driving contribution for plate motion (e.g.Elsasser, 1969;Forsyth and Uyeda, 1975).When continental lithosphere enters the subduction zone, its positive buoyancy acts as a resisting force to sinking and, therefore, leads to a slowdown, and eventually to the end of subduction (e.g.Cloos, 1993).After collision three possible scenarios may occur: (1) continental crustal material arrives at the trench and subducts together with the lithospheric mantle (e.g.Ranalli et al., 2000;Regard et al., 2003;Toussaint et al., Figures Back Close Full (2) part of the continental buoyant crust is delaminated from the lithospheric mantle and accreted onto the overriding plate, allowing continuation of plate convergence (delamination) (e.g.Cloos, 1993;Chemenda et al., 1996;Kerr and Tarney, 2005;Vos et al., 2007) or (3) the slab breaks off as a consequences of tensile stress from the pull of the sinking oceanic lithosphere in the mantle connected to the continental part of the plate (e.g.Wortel and Spakman, 1992;Davies and von Blanckenburg, 1995).
Which scenario takes place depends on several factors including convergence rate, thickness, rheological and physical composition and thermal structure of the continental lithosphere (van den Beukel and Wortel, 1987).De Franco et al. (2008) suggested that the geometry of the incoming continental plate, the tectonic setting, the rheology and the length of the sinking oceanic plate controls whether the incoming crust subducts entirely, separates partially or entirely from the lithospheric mantle or blocks the trench, likely leading to slab break-off.
An important aspect of continental collision is how plates deform.Many laboratory models of continental collision show that a subducting plate will deform internally only if the force distribution varies laterally along the subduction zone, e.g. by the asymmetrical entrance of continental material along the trench (Bellahsen et al., 2003, Regard et al., 2005, Faccenna et al., 2006).The along-trench variation of the subducting plate plays an important role in convergent margin deformation, and is often invoked to explain the indenter-like geometry of plate boundaries (e.g. the Arabia-Eurasia and Introduction

Conclusions References
Tables Figures

Back Close
Full India-Eurasia plate boundaries).Hence, plate deformation and trench migration play a critical role in a continental collision system.Although continental collision eventually leads to the cessation of subduction, the motion of the trench during continental subduction has been poorly investigated.The trench migration has been largely studied only for purely oceanic subducting plates.Heuret and Lallemand (2005) suggested that in these scenarios trench migration appears partitioned between advancing and retreating.Laboratory (Funiciello et al., 2003a;Bellahsen et al., 2005) and numerical (Funiciello et al., 2003b;Stegman et al., 2006;Di Giuseppe et al., 2008;Di Giuseppe et al., 2009;Capitanio et al., 2010b) models of oceanic subduction showed that trench migration is controlled by geometrical (thickness and width of the plate) and rheological (density contrast, viscosity ratio) parameters of the subducting lithosphere.In the continental subduction scenario, Royden (1993) proposed the same two end-member classes: the advancing style (e.g.Alps and Himalayas) and the retreating style (e.g.Appennines and Carpathians).In the retreating subduction system the lithospheric mantle detaches from the buoyant continental crust and it keeps subducting into the mantle (Cloos, 1993;Kerr and Tarney, 2005;Brun and Faccenna, 2008;Bialas et al., 2010).This crustal delamination process is possible if there is a layer of weak material within the continental crust, for example, the lower crust (Ranalli et al., 2000;Toussaint et al., 2004).In contrast, in advancing subduction systems, hundreds of kilometres of continental lithosphere can subduct without significant crustal delamination (Royden, 1993).The overturned morphology of the subducted lithosphere inferred from tomographic images provides evidence of trench advancing in continental collision zones like the India-Eurasia system (Replumaz et al., 2004(Replumaz et al., , 2010) ) and Arabia-Eurasia (Hafkenscheid et al., 2006).For the advancing case, the driving forces are debated.The advancing motion is often explained with far field stresses related to global plate motion that push the plates and force the trench to migrate towards the overriding plate.For example, to explain Indian indentation rates the ridge push of Indian ocean (Chemenda et al., 2000), the presence of the R éunion plume (van Hinsbergen at al., 2011;Becker and Faccenna, 2011)  slabs (Li et al., 2008) has been invoked.Capitanio et al. (2010a) suggested that in the India-Eurasia collisional system an imbalance between ridge push and slab pull can develop and cause trench advance and indentation.To explain the advancing of India into the Arabia-Eurasia system Becker and Faccenna (2011) proposed that mantle drag exerted on the base of the lithosphere by a large-scale mantle flow component with an active upwelling component is likely the main cause for the ongoing indentation.
In this study, we investigate the dynamic behaviour of the trench during continental collision without any external forces, i.e. driven entirely by the changes in the force balance due to the collision process.We explore the robustness of these results with a parametric study where geometrical and rheological constraints of the system are systematically changed in order to define their effect on the kinematics of continental collision.

Numerical methods and governing equations
We use Citcom, a finite element code that solves thermal convection problems in a Cartesian geometry (Moresi and Gurnis, 1996;Zhong et al., 2000).The nondimensional governing equations for conservation of mass, momentum, energy and composition, assume incompressible flow and adopt the Boussinesq approximations:

Conclusions References
Tables Figures

Back Close
Full with symbols defined in Table 1.The equations are solved with an iterative conjugate gradient solver.The thermal and compositional Rayleigh numbers in Eq. ( 2) are defined, respectively, as: where δp c = 600 kg m −3 is the density difference between the continental crust (with composition C=1) and the mantle (with C=0) and ∆T is the temperature drop over model domain.The composition function C is advected through a particle tracking technique (Ballmer et al., 2007;Di Giuseppe et al., 2008), in which tracers move with the flow velocity and give the compositional information to discriminate the buoyant continental crust from the lithospheric and mantle material.The initial tracer distribution is random with a minimum density of 40 tracers per element.
To decouple the converging plates, we impose a narrow low viscosity zone (Fig. 1) (about 20 km wide with a viscosity of 10 20 Pa s) between the incoming and the overriding plate for the top 50 km and a 200 km wide mantle wedge with same viscosity up to 100 km depth (van Hunen and Allen, 2011).For the aim of this study is essential that the trench is allowed to freely migrate.In our models the weak zone moves with the velocity of the overriding plate: each timestep the motion of a sub-vertical series of ∼100 tracers very close to the trench inside the overriding plate is monitored and the position of the trench is updated by moving the weak zone between the plates accordingly.Differential motion of the tracers at different depth allows for a change in slab dip angle, and allows the underlying mantle wedge to change its shape.
The mesh element size varies from 5 km to 15 km in each dimension, where the better resolution is used to resolve the narrow weak zone area.Introduction

Conclusions References
Tables Figures

Back Close
Full

Model setup
To study the effect of continental collision on the dynamics of trench migration we model subduction in a 2-D rectangular domain.The bottom of the domain corresponds to the upper-lower mantle discontinuity at 660 km, and all models have a 1:5 aspect ratio (Fig. 1).
The closure of an ocean basin and subsequent continental collision is modelled.In the initial setup the incoming plate consists of a subducting oceanic lithosphere with a continental block initially about 500 km from the trench (Fig. 1).The overriding plate is entirely continental.To initiate subduction without imposing any external forces, the initial oceanic slab extends to ∼300 km depth.The initial temperature field for the oceanic lithosphere is calculated following the half-space cooling solution for a 50 Myr old plate (Turcotte and Schubert, 2002).For the continental lithosphere, temperature extends linearly from 0 • C at the surface to the mantle reference temperature at 150 km depth.A 40 km deep continental crust is characterized by positive buoyancy, which creates a resisting force to subduction when the continental block inside the subducting plate arrives at the trench.The oceanic lithosphere is modelled without oceanic crust, since it does not affect the subduction dynamics significantly (Cloos, 1993).
The top boundary has a fixed temperature of 0 • C, whereas the other boundaries have a fixed mantle temperature T m = 1350 • C. Velocity boundary conditions are freeslip on all but the bottom boundary, where a no-slip condition is applied (Fig. 1).

Rheology
The viscosity of the system is temperature and stress dependent.The strain rate is accommodated by both diffusion and dislocation creep (Hirth and Kohlstedt, 2003;Karato and Wu, 1993;Korenaga and Karato, 2008).The effective viscosity η for each mechanism is defined as:

Conclusions References
Tables Figures

Back Close
Full with symbols as defined in Table 1.For diffusion creep n = 1, whereas we assume n = 3.5 for dry dislocation creep (Karato and Wu, 1993).The activation energy E * = 360 kJ mol −1 is an average of published values, that vary from 250 to 540 kJ mol −1 (Karato and Wu, 1993).A brittle yielding behaviour is calculated close to the surface to reduce the strength of the lithosphere (Byerlee, 1978).The viscosity for this rheology is given by: with τ y is the yield stress described as where τ max is the maximum yield stress and τ 0 + µp 0 is the Byerlee's law, with τ 0 is the yield stress at the surface, µ is the friction coefficient, and p 0 is the lithostatic pressure.At any point the effective viscosity is the minimum of viscosity values derived from each mechanism described above.To account for other weakening mechanisms at low temperatures, we impose a maximum viscosity value, that is varied between 10 23 and 10 24 Pa s.
The mantle wedge with constant, low viscosity is a simplification to simulate the presence of a fluid-rich area above the slab due to the dehydration of the subducting plate.This method allows a good decoupling between the plates during the subduction process.The plates may be allowed to move freely using a 30 km wide low viscosity areas at both the left and right top corners of the model domain.The subducting plate is free to move in all models, whereas, the overriding plate is either free or fixed at the right hand side (Fig. 1).

Numerical results
We performed two sets of model calculations.In the first set, we investigate the effect of the presence of buoyant continental lithosphere on trench motion.All models have Introduction

Conclusions References
Tables Figures

Back Close
Full the same initial geometry setup.Given the sensitivity of results to plate strength, differences among the models concentrate on two rheological parameters: the maximum viscosity of the lithosphere (η max ) and the friction coefficient mu.η max varies between 10 23 and 10 24 Pa s, whereas µ varies between 0 and 0.1.In models with µ = 0 the maximum yield stress τ max and the surface yield stress τ 0 are both 200 MPa, whereas, in models with µ = 0.1, τ 0 = 40 MPa and τ max = 400 MPa.All of these parameters have a strong influence on the subduction dynamics because they control the strength of the lithosphere and, therefore, its capability to deform and bend.We compare these models with a similar set of models with a purely oceanic subducting plate, i.e. by replacing the continental block with oceanic lithosphere, which affects the average density and thermal thickness.
The reference continental model C1 (Fig. 2) has a maximum viscosity value of 10 23 Pa s and depth-independent yield strength (i.e.µ = 0).The model dynamics can be subdivided in 4 phases: (1) the sinking of the slab into the upper mantle, (2) the unbending of the slab due to its interaction with the 660 km discontinuity, (3) the arrival of the continental plateau at the subduction zone and (4) the necking and break-off of the slab at depth.
In the first phase the trench retreats (Figs. 2 and 3) and its velocity progressively decreases while the slab approaches the 660 km discontinuity.Afterwards, during the second phase, because of the unbending of the slab at depth, the retreating velocity increases in order to accommodate the deep deformation (Figs. 2 and 3).After about 3.5 Myr from the beginning of the process, the continental plate reaches the zone where the slab needs to bend to enter the trench (phase 3; Fig. 3b).At this stage the trench velocity decreases until, at about 5 Myr, the continent arrives at the subduction zone, and the trench starts to advance (Fig. 3).At this point the buoyancy of the continental crust chokes the subduction process, which dramatically slows it down.After about 7 Myr since the arrival of the continent at the trench (see Fig. and the buoyancy of the shallow continental crust (phase 4).During this last phase the trench motion slightly increases its velocity for a short period and then it completely stops (Fig. 3).The break-off occurs after about 13 Myr from the collision and the continent reaches its maximum depth of ∼160 km.The amount of trench advancing is about 100 km from 5 to 22 Myr (Fig. 3).
Results from the continental model (C4) with an increased maximum viscosity value of 10 24 Pa s and µ = 0.1 are shown in Figs. 4 and 5. Although, in comparison to model C1, the cold parts of the slabs becomes stronger, material near the surface has a reduced strength due to the depth-dependent yield strength.Compared to C1 the amount of retreat in the first phases of subduction is larger.This is mostly due to the lower yield strength at the surface, which reduces the effective viscosity of the plate in the bending zone where stresses are higher.The collision occurs after about 7 Myr since subduction started and the trench begins to advance (Figs. 4 and 5).The slab starts necking at about 16 Myr and it breaks after 25 Myr.The trench advancing is clearly accelerated during this phase (Fig. 5).The continental crust arrives at a depth of almost 200 km (Fig. 4).Compared to C1, it takes longer for the slab to break because the higher maximum viscosity makes the slab stronger.Moreover, the continental crust arrives about 40 km deeper compared to C1, mostly because the lower yield strength at the surface allows for easier bending and, therefore, favours the subduction.In this model the amount of advancing is about 220 km in 20 Myr (Fig. 5).
When the reference model C1 is compared to model O1, with a purely oceanic lithosphere but otherwise the same rheological parameters of C1, significant differences are evident (Fig. 6).In the first two phases of subduction (sinking of the slab and unbending) the dynamics of the system is the same.Figure 6d-e shows that trench velocities and positions are almost identical until about 4 Myr.In model C1, 4 Myr corresponds to the arrival of the continental plateau to the bending zone close to the trench, whereas, in model O1 the system reaches the steady-state condition of retreating.
Comparing O1 with all the other oceanic subducting plate models (Fig. 7) we observe that the amount of retreat and trench velocity slightly change because of the difference Introduction

Conclusions References
Tables Figures

Back Close
Full show qualitatively the same advancing behaviour when collision occurs (Fig. 7).The far end of the overriding plate is not free to move, but still, the trench is able to retreat and advance over about 40 km in response to the dynamics of the system.This is due to internal deformation of the overriding plate.This scenario does not occur in the models with a stronger lithosphere with η max = 10 24 Pa s (C7 and C8), because the overriding plate is too strong to deform.In that case, the trench is thus almost stationary during every phase of subduction (Fig. 7b-d).

Discussion
Despite different rheological parameters, all model calculations of continental subduction show the same kinematic behaviour (Fig. 8b), in which the trench starts to advance once the continent arrives at the subduction zone.The advancing mode in continental collision scenarios is at least partly driven by an intrinsic feature of the system (Fig. 9).
The locking of the subduction zone on arrival of the continent at the trench induces the deepest part of the slab to move towards backward, and the slab to progressively steepen.This, in turn, triggers return flow around the slab creating two small-scale cells (Fig. 9).This flow drives the upper plate towards the overriding plate and results in an advancing trench.This process pursues until the slab reaches a vertical position.
The lower the slab angle is before the continental plateau arrives at the trench, the more advancing is obtained for the slab to reach its vertical position.The total amount of trench advancing depends on the dip of the slab before the collision.e.g. in model C4 the slab is flatter than in C1, and therefore the post-collisional trench advance is larger (Figs. 2 and 4).
After the vertical position is reached, the advancing velocity decreases.Later, we observe a new peak in the trench velocity plots due to slab detachment (Figs.5b).Over time, the slab weakens thermally, starts stretching, and eventually leads to slab necking and break-off (see van Hunen and Allen, 2011, for a more detailed discussion of slab break-off conditions).After slab detachment, the trench motion stops few million years later.During stretching of the slab and, eventually, the break-off, the deepest part of the slab is allowed to sink and this activates again the small-scale flow circuit (Fig. 9).The induced mantle circulation can exert shear tractions at the base of nearby plates (Conrad and Lithgow-Bertelloni, 2002).This leads to a new "pulse" of trench advancing (Figs.3b and 5b) and an overturning of the slab (Fig. 9).The resulting geometry of our models during slab detachment resembles that of the experimental results in (Regard et al., 2008) where the return flow during the break-off leads to slab overturning.Furthermore, during the stretching and subsequent slab break-off, the buoyant, previously subducted, continental material starts rising to the surface (Fig. 9).This exhumed continental material between the plates pushes them apart, thereby promoting trench migration (Brun and Faccenna, 2008).The weaker the lithosphere is, the more the exhumation is accommodated by internal plate deformation, whereas, for stronger lithosphere the exhumation leads to a larger divergence of the overriding plates.Indeed, the increase of velocity in trench motion is more evident in model C4 than in the other continental models (Fig. 8).
Models with the friction coefficient µ = 0.1, i.e. a depth-dependent yield strength, show a larger amount of both retreating and advancing (Fig. 8a) than the ones with a constant (high) yield strength.This is because the lithosphere can more easily deform near the weak surface.Moreover, models with a fixed overriding plate and η max = 10 24 Pa s show no trench motion because the strength of the lithosphere is too high to deform it.Since internal lithospheric deformation is common in collisional settings (as expressed in e.g.back-arc basins and orogeny), this suggests that the effective viscosity of Earth's lithosphere is less than 10 24 Pa s.
Although our numerical models are not intended to reproduce the detailed behaviour of any particular subduction zone, obtained results are able to highlight some key Introduction

Conclusions References
Tables Figures

Back Close
Full features of continental collision.In natural cases, evidences of trench advancing have been recognized in continental collision zones as India-Eurasia and Arabia-Eurasia (Replumaz et al., 2010;Hatzfeld and Molnar, 2010).The values of advancing we obtain in our model (50-220 km; Fig. 8) are consistent with the amount of trench migration recorded in the Arabia-Eurasia system (50-150 km; Hatzfeld and Molnar, 2010).On the other hand, our kinematic results cannot explain all of the advancing motion in the India-Eurasia collision (more than 1000 km in 50 Ma) (Matte et al, 1997;Replumaz et al., 2004;Guillot et al., 2003).Hence it is likely that external forces play a role, such as the pull of neighbouring slabs (Li et al., 2008), or the drag exerted on the base of the lithosphere by a large-scale, convective "conveyor mantle belt" (Becker and Faccenna, 2011;Cande and Stegman, 2011).Trench advancing has been observed also in the Melanesian arc system where the Ontong Java plateau, the largest igneous province on Earth, collided with the arc between 25 Ma (Hall, 2002;Knesel et al., 2008) and 10 Ma (Mann and Taira, 2004).Buoyancy analyses (Cloos, 1993) show that oceanic plateaus with abnormally thick basaltic crust resist subduction and can cause collisional orogenesis.Many authors proposed that after the collision, the anomalously thick lithosphere of the Ontong Java Plateau jammed the trench that led to the break-off of the Pacific plate (Hall and Spakman, 2002) and ultimately to a reversal of subduction polarity (Petterson et al., 1997;Knesel et al., 2008).

Conclusions
We performed a parametric numerical study to investigate the intrinsic effects of continental subducting lithosphere on trench migration.Our model results show how trenches can start advancing when continental collision occurs.These results are consistent with natural collisional zones as Arabia-Eurasia and India-Eurasia, where trench advancing is recorded in geological features.We propose that the advancing mode is induced by return flow triggered by the lock of the subduction zone and the Introduction

Conclusions References
Tables Figures

Back Close
Full following steepening of the slab generating a small-scale convection cell.The same process activates again when the deepest oceanic part of the slab sinks, thanks to the stretching and, then, the detachment of the slab, resulting in a last pulse of trench advancing.Moreover, we found that the total amount of advancing depends on the dip of the slab before the collision: the lower the slab angle is before the continental plateau arrives at the trench, the more the trench advances.Advancing behaviour is observed without imposing any external forces, and over distances similar to those observed for the Arabia-Eurasia collision.These results show that there is not necessarily a need to invoke external plate driving forces to explain trench advancing motion during continental collision.

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | or the pull of adjacent Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2 panel 12.6 Myr), the slab starts showing the first indications of necking and eventually it breaks because of the interaction of the two opposite forces: the pull of the oceanic slab at the depth Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |in the slab strength due to the variations of the rheological parameters, but the overall subduction dynamics remains unchanged.Continental models with a fixed overriding plate and η max = 10 23 Pa s (C5 and C6) Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |