Numerical models of slab migration in continental collision zones

Continental collision is an intrinsic feature of plate tectonics. The closure of an oceanic basin leads to the onset of subduction of buoyant continental material, which slows down and eventually stops the subduction process. In natural cases, evidence of advancing margins has been recognized in continental collision zones such as India-Eurasia and Arabia-Eurasia. We perform a parametric study of the geometrical and rheological influence on subduction dynamics during the subduction of continental lithosphere. In our 2-D numerical models of a free subduction system with temperature and stress-dependent rheology, the trench and the overriding plate move self-consistently as a function of the dynamics of the system (i.e. no external forces are imposed). This setup enables to study how continental subduction influences the trench migration. We found that in all models the slab starts to advance once the continent enters the subduction zone and continues to migrate until few million years after the ultimate slab detachment. Our results support the idea that the advancing mode is favoured and, in part, provided by the intrinsic force balance of continental collision. We suggest that the advance is first induced by the locking of the subduction zone and the subsequent steepening of the slab, and next by the sinking of the deepest oceanic part of the slab, during stretching and break-off of the slab. These processes are responsible for the migration of the subduction zone by triggering small-scale convection cells in the mantle that, in turn, drag the plates. The amount of advance ranges from 40 to 220 km and depends on the dip angle of the slab before the onset of collision.


Introduction
Slab pull caused by the negative buoyancy of subducting oceanic 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.McKenzie, 1969;Cloos, 1993).After continental collision, different 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., 2004); (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 consequence 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; Wong-A-Ton and Wortel, 1997;Duretz et al., 2011;van Hunen and Allen, 2011).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 control whether the incoming crust subducts entirely, separates partially or entirely from the lithospheric mantle or blocks the trench, likely leading to slab break-off.

V. Magni et al.: Numerical models of slab migration
In many areas where continental collision occurs, seismic tomography has revealed the presence of gaps in the subducting lithosphere, interpreted as evidence of detached lithosphere within the slab (e.g.Mediterranean, Carpathians, eastern Anatolia) (Carminati et al., 1998;Wortel and Spakman, 2000;Faccenna et al., 2004;Piromallo and Faccenna, 2004;Faccenna et al., 2006;Hafkenscheid et al., 2006;Lei and Zhao, 2007).Other suggestions of slab detachment have been provided by geochemical studies of magmatism in collision zones (Keskin, 2003;Ferrari, 2004;Qin et al., 2008).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).
Once the continental collision occurs, the position of the trench becomes a zone that marks a suture between the two plates.The mechanisms that drive the motion of this zone during the evolution of collision are still poorly understood.The combination of geological observations, such as structural and volcanological data, with interpretation of tomographic images may provide insights to infer the tectonic evolution of the collisional system.
The trench migration has been largely studied only for purely oceanic subducting plates.Many authors suggested that in these scenarios trenches can both advance and retreat (e.g.Carlson and Melia, 1984;Jarrard, 1986;Garfunkel et al., 1986;Heuret and Lallemand, 2005).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.Apennines 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;Gogus et al., 2011).In contrast, in advancing subduction systems, hundreds of kilometres of continental lithosphere can subduct without significant crustal delamination (Royden, 1993).
Evidence of advancing margins has been recognized in continental collision zones (e.g.Replumaz et al., 2010;Hatzfeld and Molnar, 2010 and reference therein).Trench advance has been observed also in the Melanesian arc system where the Ontong Java Plateau (Hall, 2002), the largest igneous province on Earth, collided with the arc between 25 Ma (Knesel et al., 2008) and 10 Ma (Mann and Taira, 2004).The overturned morphology of the subducted lithosphere inferred from tomographic images in continental collision zones like the India-Eurasia (Replumaz et al., 2004;Replumaz et al., 2010) and the Arabia-Eurasia (Faccenna et al., 2006;Hafkenscheid et al., 2006) systems suggests the occurrence of an advancing migration of the plate margin.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 slab 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) or the pull of adjacent slabs (Li et al., 2008) has been invoked.Capitanio et al. (2010a) suggested that the dense Indian continental slab provides a significant driving force for the current India-Asia motion.To explain the advance 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 slab and the margin between the plates 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 non-dimensional governing equations for conservation of mass, momentum, energy and composition assume incompressible flow and adopt the Boussinesq approximations: 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 ρ 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.
The mesh elements are quadrilateral with a linear velocity and constant pressure, and their size varies from 5 km to 15 km in each dimension, where the better resolution is used in the trench region.

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 upperlower 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 and has an initial curvature radius of 500 km.The initial thermal structure of the oceanic lithosphere is determined by a half-space cooling model 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).
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 150 km depth that simulates the weakening of the area above the slab due to the slab dehydration and the mantle hydration (Billen and Gurnis, 2001;van Hunen and Allen, 2011).In our models, we refer to the trench position as the location of the shallow weak zone.In oceanic subduction, this refers to the area where the subducting slab begins to descend beneath the overriding plate, whereas, in continental collision systems, this is a rather wide deformation zone between the pro-and the retro-wedge (Willet et al., 1993;Beaumont et al., 1996).For the aim of this study, it is essential that the trench is allowed to freely migrate.In our models, the weak zone moves horizontally with the velocity of the overriding plate: at each time step, the horizontal velocity of a point within the overriding plate and close (about 30 km) to the trench is taken from the global model velocity field.This velocity is used to calculate the new position of the trench.The shallow weak zone between the plates has a fixed shape, whereas the shape of the weak mantle wedge changes to follow the variations of the dip of the slab during the model evolution.The position of a set of a hundred tracers within the slab at a depth interval of 50-150 km is used to re-shape the mantle wedge according to the slab dip.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 free-slip on all but the bottom boundary, where a no-slip condition is applied to model the effect of the high viscosity lower mantle acting as a rigid boundary (Fig. 1).The plates may be allowed to move freely using 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 on the right-hand side (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 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, which vary from 250 to 540 kJ mol −1 (Karato and Wu, 1993).In addition, the viscosity for the yielding mechanisms is calculated to reduce the strength of the lithosphere: τ y is the yield stress described as where τ max is the maximum yield stress and τ 0 + µp 0 is the Byerlee's law (Byerlee, 1978), τ 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, which is varied between 10 23 and 10 24 Pa s.The same rheology is assumed for mantle and crustal material without any internal rheological layering in the lithosphere.

Numerical results
We performed two sets of model calculations (Table 2).In the first set, we investigate the effect of the presence of buoyant continental lithosphere on trench motion.All models have 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 µ. η 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 into 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 continent at the subduction zone and (4) the necking and break-off of the slab at depth.
In the first phase, the trench retreats and its velocity progressively decreases while the slab approaches the 660 km discontinuity (Fig. 2).The decrease of the retreating velocity is linked to the fact that the slab steepens.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 (Fig. 2c).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. 2a).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. 2b-c).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. 2a 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 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. 2c).The break-off occurs after about 13 Myr from the collision and the continent reaches its maximum depth of ∼160 km.The amount of advance is about 100 km from 5 to 22 Myr (Fig. 2b).
Results from the continental model (C4) with an increased maximum viscosity value of 10 24 Pa s and µ = 0.1 are shown in Fig. 3. Although, in comparison to model C1, the cold parts of the slabs become stronger, and 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 (Fig. 3).The slab starts necking at about 16 Myr, and it breaks after 25 Myr.The trench advance is clearly accelerated during this phase (Fig. 3b-c).The continental crust arrives at a depth of almost 200 km (Fig. 3a).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 advance is about 220 km in 20 Myr (Fig. 3b).
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. 4).In the first two phases of subduction (sinking of the slab and unbending), the dynamics of the system is the same.Figure 4d-e shows that trench positions and velocities 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.At this stage, the direction of trench migration in the two models is opposite: the trench starts to advance in models C1, whereas, in model O1, it keeps retreating.Comparing O1 with all the other oceanic subducting plate models (Fig. 5), we observe that the amount of retreat and trench velocity slightly change because of the difference 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) show qualitatively the same advancing behaviour when collision occurs (Fig. 5a).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. 5b).

Dynamics of continental slab migration
Despite different rheological parameters, all model calculations of continental subduction show the same kinematic behaviour (Fig. 6b), 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. 7).The locking of the subduction zone on arrival of the continent at the trench induces the deepest part of the slab to move backwards, and the slab to progressively steepen.This, in turn, triggers a return flow below the slab creating a small-scale cell that drives the subducting plate towards the overriding plate and results in an advancing margin (Fig. 7).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 advance is obtained for the slab to reach its vertical position.The total amount of advance depends on the dip of the slab before the collision.For example, in model C4 the slab is flatter than in C1, and therefore the post-collisional advance is larger (Figs. 2 and 3).
After the vertical position is reached, the advance velocity decreases.Later, we observe a new peak in the trench velocity plots due to slab detachment (Figs.2c and 3c).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 a 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.In particular, one vigorous convection cell forms below the slab and a weaker one above the slab (Fig. 7c).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 advance (Figs.2c and  3c) and an overturning of the slab (Fig. 7).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. 7).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 the trench motion is more evident in model C4 than in the other continental models (Fig. 6).
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 features of continental collision.Tomographic images beneath continental collisions show that slabs are steep and often detached (e.g.India, Replumaz et al 2004;Northern Apennines, Lucente et al., 1999;Faccenna et al., 2001;Piromallo and Morelli, 2003;Carpathians, Wortel and Spakman, 2000;Alps, Piromallo and Faccenna, 2004;Spakman et al., 2004).Furthermore, the values of advance we obtained in our model (50-220 km; Fig. 6a) are consistent with the amount of migration of the deformed belt 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).

Model limitations
Models with the friction coefficient µ = 0.1, i.e. a depthdependent yield strength, show a larger amount of both retreating and advance (Fig. 6a) than the ones with a constant (high) yield strength.This is because the lithosphere can more easily deform near the weak surface.In our models, the lithosphere is not free to move vertically since we use a free slip boundary condition for the surface.Therefore, a weaker lithosphere (with µ = 0.1) is probably more appropriate to overcome this simplification (see Enns et al., 2005;Di Giuseppe et al., 2008 for more details).We observe that in all models the amount of retreat during the oceanic subduction and the advance during continental subduction are approximately the same.This is because, when the slab arrives at the bottom boundary, it is almost vertical and its dip decreases (slightly in model C1 and dramatically in model C4) (Figs. 2 and 3) during the retreating phase.Then, the advancing phase continues until the vertical position is reached again (and a little after that, especially in model C4) and, therefore, approximately at the position of the trench in the early stage of the model.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.However, internal lithospheric deformation is common in collisional settings (as expressed in e.g.back-arc basins and orogeny).Modelling results are difficult to directly export to a natural system, given the several assumptions made in the models (such as absence of possible external forces).Keeping in mind these assumptions, however, we propose that a viscosity of 10 24 Pa s as too high to model the effective viscosity of Earth's lithosphere.
The mantle flow patterns that lead to the trench advance in our models are two-dimensional.To test the effect of the toroidal flow around the slab edges, we run some preliminary 3-D models of continental subduction and we found the same advancing behaviour.Furthermore, in our models we do not change the length of the oceanic slab, whereas in natural cases continental collision might occur after a long-lived oceanic subduction.However, we expect that the dynamics of the system would not change with a longer oceanic slab.Though, this might have an effect on the amount of trench migration, since the dip of the slab once the continent arrives at the trench might be different.Finally, in our models we assume that the discontinuity between upper and lower mantle is an impermeable barrier, whereas in nature slabs can penetrate and might become anchored in the lower mantle.This is likely to have an effect on the dynamics of subduction and, therefore, on trench migration.However, even taking into account the presence of the lower mantle, the advancing mode is still expected since the mechanisms described above would still occur.

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 such as Arabia-Eurasia and India-Eurasia, where the advancing style 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 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 advance.Moreover, we found that the total amount of advance 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 the migration of the subduction zone towards the overriding plate during continental collision.

Fig. 1 .
Fig. 1.Initial model setup with dimensions and mechanical and thermal boundary conditions.The model domain is 660 km deep, and it has an aspect ratio of 1:5.The continental parts are modelled with a thick lithosphere (yellow) and a 40 km buoyant crust (blue).Areas with an imposed low viscosity (i.e. the reference mantle viscosity) are outlined in red.

Fig. 2 .
Fig. 2. Results of the reference model C1 (η max = 10 23 Pa s and µ = 0): (a) viscosity plots of the model evolution; arrows show the velocity field.Thick arrows indicate the trench position at each time step; the dotted line shows the trench position at the beginning of the model.(b) Trench position and (c) trench velocity during model evolution.Dashed lines correspond to panels of the viscosity plots in (a).Negative velocity indicates trench retreating.Different colours and symbols indicate the four different subduction phases: black triangles (phase 1), magenta squares (phase 2), grey diamonds (phase 3) and blue dots (phase 4).

Fig. 3 .
Fig. 3. Results of the model C4 (η max =10 24 Pa s and µ = 0.1): (a) viscosity plots of the model evolution; arrows show the velocity field.Thick arrows indicate the trench position at each time step; the dotted line shows the trench position at the beginning of the model.(b) Trench position and (c) trench velocity during model evolution.Dashed lines correspond to panels of the viscosity plots in (a).Negative velocity indicates trench retreating.Different colours and symbols indicate the four different subduction phases: black triangles (phase 1), magenta squares (phase 2), grey diamonds (phase 3) and blue dots (phase 4).

Fig. 4 .
Fig. 4. (a-b-c) Viscosity plots of purely oceanic subducting plate model O1; arrows show the velocity field.Thick arrows indicate the trench position at each time step.Comparison between model C1 (red dots) and O1 (blue squares): (d) Trench migration and (e) trench velocity (negative velocity indicates trench retreating) during model evolution.

Fig. 6 .
Fig. 6.(a) Amount of advance for all continental models grouped by the changing parameters: free op = free overriding plate, fix op = fixed overriding plate.η max is the maximum viscosity and µ is the friction coefficient.(b) is trench migration of all continental models.

Fig. 7 .
Fig. 7. Schematic evolution of the continental models C1 (left column) and its viscosity plot in the column on the right (the continental material is outlined in white): (a) during oceanic subduction the trench retreats; (b) the continental collision causes the lock of the subduction zone and a consequent steepening of the slab; (c) the formation of small-scale convection cells and the exhumation of the buoyant continental material cause the advance of the trench after continental collision and during the slab detachment.

Table 1 .
Symbols, units and default model parameters

Table 2 .
Summary of the models used in this parametric study.