Oceanic lithosphere thickens and strengthens as it grows older. Worldwide databases reveal that the age of the slab to a certain extent controls the subduction style. Old and thick (and consequently strong) slabs show a trench “advance,” while younger, thin (weak) slabs are migrating in retreating style (trench “rollback”). We performed numerical models to show that this configuration could be the result of the dynamic equilibrium between gravity and resisting forces. In particular we show that energy dissipation caused by bending and unbending of the slab, although less important than other resisting forces, could be a primary control on trench migration. Our results fit well with global compilations of kinematic data from modern subduction zones in two reference frames with different amounts of net rotations. Based on the age at which the transition from retreating to advancing style occurs, we propose an effective lithosphere/mantle viscosity of ~200 during bending of the lithosphere into the subduction zone.


Lithosphere bends and subducts into the mantle at trenches. Geological data show that backarc regions undergo episodic deformation, possibly related to the intermittent migration of trenches (Taylor and Karner, 1983). But trench motion cannot be directly measured and can only be inferred by removing the deformation rate of the arc from the upper plate motion (Dewey, 1980). Molnar and Atwater (1978) first advanced the concept that the gravitational forces acting on old and dense subducting lithosphere produced the backward migration of trenches. The relationships between trench-plate kinematics and lithospheric age at trenches have been explored and suggest that differences in tectonic styles and geologic development in the overriding plate may be favored by the slab age (Carlson et al., 1983; England and Wortel, 1980; Jarrard, 1986; Heuret and Lallemand, 2005; Buffett and Rowley, 2006). However, Carlson and Melia (1984) showed that the Izu Bonin–Mariana subduction zone, the oldest subducting lithosphere on Earth, is presently advancing in the hotspot reference frame toward the Philippine overriding plate. Advancing slab is characterized by a trench motion directed toward the overriding plate, while slab retreat (or rollback) occurs when the trench migrates seaward with respect to the underlying mantle. Heuret and Lallemand (2005), using a global compilation of the absolute motion of upper plates and trenches, and backarc deformation rates and slab ages, found an inverse relationship between trench migration and age of the lithosphere at the trench, where old subducting plates advance toward the upper plate and young ones retreat. This compilation then illustrates that Western Pacific trenches migrate toward the upper plate, whereas the Eastern Pacific ones retreat. This relationship contradicts the common idea that old slabs roll back (Molnar and Atwater, 1978; Dewey, 1980; Garfunkel et al., 1986) and poses new questions on the dynamics of subduction and its relationships with trench kinematics. More recently, laboratory models (Bellahsen et al., 2005; Faccenna et al., 2007; Funiciello et al., 2008), numerical experiments (Di Giuseppe et al., 2008) along with conceptual models (Lallemand et al., 2008) suggested that the strength of the slab, in particular its resistance to bend, could exert a primary control on the style of trench migration. Alternatively, from a global view, it has been proposed that the asymmetric plate configuration across the Pacific and the tectonic forcing caused by upper plate motion could force the trench migration (Husson et al., 2008; Nagel et al., 2008).

Here, we quantify how lithospheric age influences the kinematics of the subduction process with numerical simulations. Our results show that thermal aging increases the lithospheric strength so that old ocean (>80 Ma) forces the trench to advance while young ocean retreats. In particular, numerical calculations fit natural data if effective lithosphere/mantle viscosity ratios in the bending zone are ~200. Our results imply that the age discrepancy in the Pacific, by causing differences in advancing-retreating styles, may cause the overall westward motion of the Pacific plates.


We use the updated worldwide database published by Heuret and Lallemand (2005) and Lallemand et al. (2005), which consists of transects every 2 degrees of latitude-longitude covering nearly 36,000 km of trenches. We analyzed 146 transects for a total of 36 trenches extracted from 180 cross sections of oceanic subduction zones that are “not perturbed” by subduction of ridges or oceanic plateaus. We also excluded certain subduction zones (Tonga, Sandwich, Philippine, Yap, Luzon, and New Hebrides), where lateral vigorous toroidal mantle flow produces along-strike variation of the trench velocity (Heuret and Lallemand, 2005).

Because trench migration rates are usually quite low and of the same order of magnitude as the net rotation (NR), the choice of the reference frame can be relevant (Funiciello et al., 2008; Becker and Faccenna, 2009; Schellart et al., 2008; Schellart, 2008). Therefore, we present our results in two different hotspot reference frames: HS3 (Gripp and Gordon, 2002) and SB04 (Steinberger et al., 2004). These two frames differ on the choice of the selected hotspot tracks (Pacific for HS3 and Indo-Atlantic for SB04), resulting in a different amount of NR of the lithosphere with respect to the mantle (HS3: 0.44°/Ma; SB04: 0.17°/Ma; Becker and Faccenna, 2009). Figure 1A shows the variation of the trench velocities, Vt, on a seafloor-age map (Müller et al., 1997) for HS3 and SB04. The signs of most trench motions (i.e., advancing or retreating) are similar for the two reference frames, but rates and directions differ. Note that only for the case of the Marianas–NE Japan–Kuriles–Kamchatka trench, the style of migration changes as a function of the reference frame, as all the Eastern Pacific trenches retreat and most of the Western Pacific ones advance. Depicting the normal-to-trench component of Vt versus age of the subducting seafloor in the HS3 (Fig. 1B) and SB04 (Fig. 1D) reference frames shows that in both cases young slabs retreat (Vt> 0), whereas old ones advance (Vt< 0). However, the linear correlation is stronger for HS3 than for SB04. This relationship flattens out in a no-net-reference frame, but a net rotation naturally appears in a model Earth with lateral viscosity variations (Ricard et al., 1991), and a net rotation close to the SB04 one properly matches azimuthal anisotropy measurements (Becker, 2008). As in previous works (e.g., Heuret and Lallemand, 2005; Goes et al., 2008; Schellart, 2008) the age at the trench is taken as representative of the age of the slab at subduction. The threshold age between the two trench migration styles occurred ca. 70 and 80 Ma for HS3 and SB04, respectively. The correlation between plate velocity, Vp, and lithospheric age (Figs. 1D and E) shows somewhat less evidently that old plates move faster than young plates.

Of course, other parameters could play important roles in the kinematics of the subduction system. The fluid budget, for example, could influence trench migration by weakening the upper plate (Gerya et al., 2008; Faccenda et al., 2008). The width of the slab appears as the main controlling factor as larger slabs stir larger volumes of mantle material (Funiciello et al., 2003; Bellahsen et al., 2005; Schellart et al., 2007). Another important parameter is represented by the length of the plate or the distance between the trench and the ridge, influencing the potential energy budget at the trench (Nagel et al., 2008). Among others, the upper plate motion shows significant correlation with trench migration, plate motion, and back deformation (Otsuki, 1989; Uyeda and Kanamori, 1979; Heuret and Lallemand, 2005). In the case of the Eastern Pacific trenches, for example, it has been demonstrated that the motion of the upper plate and the resulting backarc deformation could indeed produce a decrease in the plate motion (Iaffaldano et al., 2006) and, at global scale, can produce an overall reorganization of the Pacific mantle (Husson et al., 2008). However, if this relationship appears significant on the eastern side, it is less intuitive for application on the western side, where the upper plates are either stationary or move away from the trenches. Under these conditions, subducting and overriding plates are decoupled and, as demonstrated by numerical models (Pacanovsky et al., 1999) and analogue models (Heuret et al., 2007), the role of the suction force is minor in the overall equilibrium. Finally, the possibility that the Pacific asymmetry could influence the trench migration and plate motion would result in a bimodal distribution between Vt and Vp or age, whereas Figure 1 shows a somewhat linear relationship.

Although we are aware of the possible influence of other external forcing, our modeling simulations are focused on testing the role of the lithospheric thermal strengthening on the overall kinematics of the subduction system.


We constructed three-dimensional (3-D) numerical models of viscous flow with compositional buoyancy by using the parallel finite element code CITCOM (Moresi and Gurnis, 1996; Zhong et al., 2000). The code solves the continuity equation, the momentum equation, and the equation for chemical advection. The mantle is treated as an incompressible viscous medium, with an infinite Prandtl number and Boussinesq approximations. The composition function is a switch between two end-member components representative of the oceanic lithosphere and the upper mantle, respectively. Compositional information is carried by a large number of particles, more than a hundred per finite element, in order to ensure numerical accuracy. Temperature is not explicitly solved; therefore, thermal effects and phase changes are neglected. In addition, density variations owing to thermal expansion are replaced by an equivalent compositional density variation (Enns et al., 2005; Stegman et al., 2006; Di Giuseppe et al., 2008). The 3-D computational domain has an aspect ratio of 7 × 7 × 1 in x-, y-, and z-directions, respectively, representing a square box that is ~4600 km wide/large and 660 km deep. In all, 32 elements are in the vertical direction and 192 in both horizontal directions.

External forces are not imposed: In this closed system, gravity forces arising from density variations drive the entire dynamics.

Lithosphere has a finite strength that can be affected by several factors, such as temperature, pressure, presence of water, and applied stresses. When tectonic stress exceeds the material strength, the rock fails by brittle deformation (Kohlstedt et al., 1995). In order to model such brittle failure in the lithosphere, we assume a visco-plastic rheology such that the material fails when a certain yield stress is exceeded, while the viscosity is purely Newtonian elsewhere. The upper mantle viscosity is fixed at 1020 Pa·s. Brittle failure occurs when the second invariant of the stress tensor exceeds a yield value τyield:
with μf the friction coefficient and P the lithostatic pressure (Byerlee, 1968). The friction coefficient μf has been fixed at 0.08 (Di Giuseppe et al., 2008). Owing to large bending stresses, the strength of the plate in the bending zone is reduced, resulting in a lower effective viscosity. Thus, the weakening can be expressed by an effective viscosity ηeff:
where ηlo is the imposed plate viscosity (varying between 4 × 1022 and 2 × 1023 Pa·s), and ηyield is equal to τyield/2ε̇II, with ε̇II the second invariant of the strain rate.

The boundary conditions are described by a free-slip top and by a no-slip bottom and side-walls. The 660-km discontinuity is simulated as an impermeable barrier. Even if slabs can penetrate into the lower mantle (e.g., van der Voo et al., 1999), this assumption is validated by previous work, which showed that penetration through the discontinuity is temporarily inhibited if the viscosity contrast between the upper and lower mantle is larger than one order of magnitude (e.g., Guillou-Frottier et al., 1995; Funiciello et al., 2003), by the effects of the endothermic phase change (e.g., Ringwood and Irifune, 1988), and if the time scale of the analyzed process is limited (e.g., Davies, 1995). In fact, at the present day, most subduction appears to be driven only by upper-mantle slab buoyancy (Faccenna et al., 2007; Goes et al., 2008).

The overriding plate is not modeled. The presence of the overriding plate and its interaction with other plates, the backarc, and the mantle wedge are known to be important in the slab dynamics (Billen and Gurnis, 2001; van Keken, 2003; Becker and Faccenna, 2009; Heuret et al., 2007), but this is beyond the scope of our investigation. Thereby, the setup is as simple as possible, and the stress regime on the overriding plate cannot be inferred.


In order to simulate the influence of the lithospheric age and strength on trench motion, Vt, a set of 22 models has been constructed. The imposed viscosity contrast between lithosphere and upper mantle, ηloum, ranges from 400 to 2000. The age variation is obtained by changing the lithospheric thickness, h, from 58 km to 114 km. Plate thickness is converted to lithospheric age, t, assuming the half-space cooling model according to which the lithospheric thickness is proportional to the square root of its age (Turcotte and Schubert, 1982; Carlson et al., 1983): h = 2(κt)1/2 [m], with t in [s], and κ~10−6 m2s−1 being the thermal diffusivity, so that h refers to the depth of the isotherm, T = 0.84Tm (Davies, 1999).

Two different styles of subduction are observed when varying the geometrical and rheological model parameters: (1) a retreating style (or slab rollback) with Vt > 0, and (2) an advancing style with Vt < 0. Sometimes a slab begins its migration in advancing style and ends it in retreating style after folding on the discontinuity (Bellahsen et al., 2005; Di Giuseppe et al., 2008). A weak and/or dense slab creates a retreating trench (slab rollback) and a low Vp. Advancing trenches result from subduction of strong and/or light slabs, producing a backward-reclined slab geometry and a high Vp (Fig. 2).

The results of the numerical models show that a stronger slab has difficulty to unbend during its descent into the upper mantle: The trench accommodates this backward reclined configuration by advancing toward the upper plate. A weak slab, conversely, easily unbends and folds at depth, thereby favoring slab rollback.

Figure 3A shows a positive trench velocity (trench rollback), Vt, for young slabs, which flips to a negative velocity for old slabs. The critical lithospheric age for the transition between the advancing and retreating trench is variable and scales inversely with the slab-mantle viscosity contrast. These model results fit the following mathematical expression well:
where Vto is the amplitude of the trench velocity, to marks the threshold age at which the transition from retreating style (Vt > 0) to advancing style (Vt < 0) occurs, and Δt is ca. 1 Ma, being the reference interval of time. We fixed this parameter to a reasonably fitting value of 1 Ma, although the data density is not sufficient to allow a more thorough parameter investigation. Both the threshold age, to, and the amplitude of the fitting curves, Vto, progressively increase with a decreasing viscosity ratio ηloum (Fig. 3); to is ca. 41 Ma, and Vto ~27.6 mm/a for a viscosity ratio of 2000; to is ca. 82 Ma, and Vto is ~37.2 mm/a for ηloum = 400. Vt ≅ 0 refers to a stable trench. In models such as ours, where the 660-km discontinuity is modeled as an impermeable barrier, Vt ≅ 0 would imply slab accumulation, which, over long periods of time, is not a stable solution. This explains why our resulting Vt values are rarely close to zero and why equation 3 does not fit the observations that show a clear linear trend (Fig. 1). On Earth, slabs can enter the lower mantle and may even stay fixed by anchoring in the more viscous lower mantle.

Figure 3B describes the behavior of plate velocity versus lithospheric age. Old slabs appear characterized by fast plate motion, and young slabs by slow plate motion. The trend is given by an inverse linear proportionality between the age of the plate and its velocity. The velocity of a plate of the same thickness increases with an increase in the viscosity contrast ηloum.


Lithosphere cools, thickens, and contracts in time, inducing changes in the strength of the lithosphere and its buoyancy. In the numerical model setup presented in this paper we tested the effects of lithospheric aging on trench migration. Our results show that the style of trench migration is strongly controlled by the lithospheric aging. Young slabs, and consequently thin, weak, and light slabs, enhance the retreating style of subduction characterized by trench rollback, low rates of incoming plate velocity, and shallow slab dip (Fig. 2). Old slabs, and consequently thick, strong, and heavy slabs, enhance the advancing style of subduction characterized by trench advancement, high rates of incoming plate velocity, and deep slab dip (Fig. 2). In particular the lithosphere is hardly deformed after being bent and subducted at the trench, and it reaches the 660-km discontinuity maintaining a backward-reclined shape. In nature, most of the tomographic images show subducting plates horizontally deflected or flattened in the upper or lower mantle transition zone or, less frequently, penetrating into the lower mantle (Fukao et al., 1992, 2001), whereas full backward-reclined slab shapes are not observed (except perhaps for the case of India: van der Voo et al., 1999; and below the Taiwan subduction zone, Lallemand et al., 2001). This can be explained by the fact that the 660-km discontinuity is not fully impermeable, as in our simplified models, or because the present-day trench advancement style could represent only a fraction of a longer and more complex history of trench migration (van der Hilst and Seno, 1993).

The effect of thermal thickening of the oceanic lithosphere is twofold: It increases the negative buoyancy of the subducting plate (Carlson et al., 1983; Molnar and Atwater, 1978; England and Wortel, 1980), and also its strength. Plate strength may be evaluated by its resistance to bending: Thick slabs show high resistance to deformation and bending, whereas thin plates are more easily deformed (Conrad and Hager, 1999; Capitanio et al., 2007; Di Giuseppe et al., 2008; Lallemand et al., 2008). Our models confirm that the lithospheric yield strength has a stronger effect on subduction dynamics than the (somewhat small) differences in slab density owing to plate age (Billen and Hirth, 2007; Billen, 2008). Despite the fact that on Earth the plate density increases with lithospheric aging, the effects of thermal contraction on density contrast can be assumed constant in time (Korenaga, 2007). We find, by fitting our numerical results to a scaling parameter, that the lithospheric stiffness, S, is a function of the viscosity ratio, thickness, and density contrast (Di Giuseppe et al., 2008). The competition between the slab strength, which is strictly related to the bending resistance, and the slab pull-force controls the trench behavior. A stronger slab has difficulty unbending inside the mantle, and so attains a backward reclined configuration, forcing the trench to advance (Di Giuseppe et al., 2008).

This result is confirmed by Lallemand et al. (2008), who recognize the resistance of the slab to bending as the main factor that influences subducting plate dynamics, because the plate stiffness increases with age faster than the slab pull for plates older than 80 Ma. Lithospheric strength has also been shown to play a dominant role in slab behavior according to the variations of the present-day potential temperatures (e.g., Tetzlaff and Schmeling, 2000; van Hunen et al., 2000; Billen and Hirth, 2007) and on the viability of plate tectonics during the Precambrian (van Hunen and van den Berg, 2008). In the Precambrian Earth, the balance of the temperature and the melting dehydration led to subduction, slab break-off, or no subduction at all, if the plate was strong or too weak, respectively.

A key mechanism affecting the lithospheric strength in our models is the presence of pseudo-plasticity at the bending zone. Because of pseudo-plasticity the lithospheric viscosity undergoes a significant reduction relative to the imposed initial value, ηlo (Enns et al., 2005; Stegman et al., 2006; Di Giuseppe et al., 2008), and the trench is weakened differently according to model parameters—i.e., the viscosity ratio and thickness (Fig. 4). In young (thin) plates the relatively shallow weakening meets the asthenosphere, affecting the slab until it reaches great depths, (i.e., green color at the bending region, Fig. 4A and C), whereas in old (thick) slabs the asthenosphere is too deep, and the plates stay strong there (i.e., yellow color, Fig. 4B and D). At the same time, plates with lower viscosity ratios (e.g., ηloum = 5 × 102) are on the whole weaker than more viscous slabs (light colors, Fig. 4C and D), and therefore are more deformable. But is a similar mechanism operating on the Earth? Subduction implies a bending lithosphere at the trench. How this mechanism weakens the lithosphere while bending occurs is still poorly understood, but a systematic decrease of shallow lithospheric strength by brittle and ductile failure is in any case expected from rock rheology (e.g., Brace and Kohlstedt, 1980; Ranalli, 1992). In particular, the presence of bending-related faulting that cuts across the crust, penetrating deeply into the mantle is confirmed by multibeam bathymetry and seismic reflection images (Ranero et al., 2005), high-resolution tomography (Zhao et al., 2000; Calò et al., 2009), and studies of lithospheric flexure (McAdoo et al., 1985; Billen and Gurnis, 2005). Faults promote the hydration of the cold crust, weakening lithospheric strength by possible serpentinization (Ranero et al., 2003; Calò et al., 2009; Faccenda et al., 2008) and magma escape (Hirano et al., 2006). These data illustrate how the empirically calibrated weakening mechanism operating in our models is likely appropriate also for nature.

Modeling results can be difficult to export to a natural system because of their dependency on a simplified applied starting condition. Hence, other mechanisms not modeled here may influence the subduction style. Among those parameters is the width of the plate, which may be relevant, as pointed out in other modeling results (Bellahsen et al., 2005; Stegman et al., 2006; Schellart et al., 2007; Di Giuseppe et al., 2008). In general, wider plates stir a larger amount of mantle material. Thereby, different amounts of the return mantle flow could modify the slab migration (Funiciello et al., 2006; Piromallo et al., 2006) as well as the evolution of the trench shape (Morra et al., 2006). However, in our models (Di Giuseppe et al., 2008) the strength of the plate does not change by varying the width of the plate. Therefore, this parameter has not been considered in this work. In addition, the presence of the lower mantle and phase transitions–viscosity increase can strongly deform the slab and potentially modify its behavior at depth (e.g., Bunge et al., 1997; Lithgow-Bertelloni and Richards, 1998; Tetzlaff and Schmeling, 2000; Cizkova et al., 2002; Billen and Hirth, 2007). Finally, as introduced before, the kinematics of the overriding plate may affect trench migration and overriding plate deformation (Heuret et al., 2007). Husson et al. (2008) highlight the effects of the upper plates on the net rotation of the Pacific—i.e., the forces exerted by the South American plates. The decrease of the convergence rate between the Nazca plate and South America may have been caused by the growth of the Andes (Iaffaldano et al., 2006). Despite these oversimplifications, however, our numerical models provide an efficient, systematic means for supplying insightful information about the role played by plate strength in influencing the style of subduction. In this sense the capability of the slab's advance-retreat can be expressed in terms of an effective viscosity, ηeff. Our models show that ηeff at the trench inversely scales with plate age (Fig. 5). Similar behavior was recognized in nature by Zhao et al. (2000), who, in analyzing their tomographic data, recognized how existing fractures on the seafloor cut young, thin slabs and old, thick slabs differently. A larger lithospheric weakening occurs at the bending zone for younger slabs, whereas there is a smaller weakening for older plates.

It is now tempting to examine where in Figure 5 natural trench data are distributed, speculating further on the possibility of ηeffum matching the natural value. The critical age that marks the transition between old, advancing trenches in the Western Pacific and young, retreating trenches in the Eastern Pacific occurred ca. 70 and 80 Ma for the HS3 and the SB04 reference frames, respectively (Fig. 1B and D). This age also is relevant, as it represents the threshold over which the thermal thickening of the oceanic lithosphere flattens (Parsons and Sclater, 1977) and marks the transition from compressional to extensional regimes behind subduction zones (England and Wortel, 1980). Using our modeling results we find that this age corresponds to ηeff / ηum, ranging between 100 and 200 (Fig. 5).

The different trench-plate velocity distributions shown in Figure 1 depend on the adopted reference frame and its NR with respect to the deep mantle. We argue that the amount of NR scales linearly with the ηeffum viscosity ratio. For example, for high ηeffum, the critical age for the variation in trench migration style decreases (i.e., ca. 40 Ma and ca. 68 Ma, for averaged effective viscosity ratios of 500 and 200, respectively). This idea of a relatively weak slab is consistent with studies on plate strength and flexural rigidity by using topography-gravity admittance (Billen and Gurnis, 2005) and arguments about subductability of purely viscous slabs (Conrad and Hager, 1999; Becker et al., 1999; Funiciello et al., 2008) and global plate mobility (Wu et al., 2008) as constrained by numerical and laboratory modeling results integrated with natural data. However, there is a discrepancy with respect to the lower value and the extremely weak slab viscosity (less than 1022) predicted by geopotential, dynamic topography and in-slab force transmission (Hager, 1984; Moresi and Gurnis, 1996; Moresi et al., 1996; Zhong and Davies, 1999; Billen et al., 2003). This difference can be explained as being due to results that came from time dependent and instantaneous models, respectively. In the former models the yield stress brings variations in slab viscosity along the slab. The latter, forced to use a viscosity cutoff and a uniform slab viscosity, appear more sensitive to the presence of the weakest trench region and give lower range values of slab viscosity and, in turn, ηeffum. Hence, the yield stress likely plays a fundamental role in controlling the degree of slab coupling with the surrounding mantle (Billen, 2008).

Our result confirms previous laboratory modeling results (Faccenna et al., 2007) that identify the pull of the upper mantle slab as the driving force that directly pulls the plates. However, the two modeling approaches differ in the partitioning of plate and trench motion, which is modulated by the resisting forces—i.e., the bending resistance of the slab at the trench and shear drag between the slab and mantle. Because the latter is almost constant for upper mantle slabs, the bending resistance, although subordinate (Di Giuseppe et al., 2008), is here playing a leading role in the trench kinematics and, in turn, on plate motion. This and previous (Faccenna et al., 2007) regional models then differ from global model results, in which bending resistance degrades the plate motion fit (Wu et al., 2008). This difference is probably related to the different approach between regional and global models (Becker and Faccenna, 2009). Regional models have the advantage of incorporating trench migration, but at the same time neglecting the large-scale contribution given by the surrounding trenches and slabs and by the rigid plate.

If this model is correct, it could play an important role in the global pattern of plate motion, because advancing trenches on the Western Pacific plate drive plates faster than retreating ones on the Eastern Pacific. As a large fraction of net rotation is carried out by the high westward speed of the Pacific plate, we speculate that the regional advancing style of trench migration could contribute, possibly in conjunction with other possible mechanisms (Zhong and Davies, 1999; Becker, 2008), to global net rotation.


The results of our numerical models indicate that the switch from a retreating to advancing subduction style, as observed along the Western Pacific subduction zone (i.e., Heuret and Lallemand, 2005; Sdrolias and Müller, 2006), may have been caused by the aging of the lithosphere at the trench. In particular, we highlight that the style of subduction and the capability of the trench to advance or retreat is controlled by the slab strength, which depends mainly on lithospheric aging. The stiffness of the plate is proportional to its age: Old plates are thick and therefore strong, whereas young plates are thin and consequently weak. Natural cases confirm, in general, the dichotomy found with numerical models: Young, weak subduction zones (e.g., Mexico, Costa Rica, Chile) retreat, whereas old, strong subduction zones (e.g., Marianas, Izu-Bonin, Antilles, Kermadec) advance. In addition, our results confirm the important role played by the effective relative viscosity between the plate and the upper mantle in controlling trench migration. Comparison of natural and numerical data suggests an effective viscosity ratio ranging from 100 to 200 in order to obtain the variety in subduction style recognized on Earth.

We thank two anonymous reviewers and Editor R.M. Russo for their careful reviews, which improved an earlier version of this manuscript. We thank Arnauld Heuret and Serge Lallemand for providing us with the updated worldwide data set and discussions. We also thank Saskia Goes for the review of an early version of this manuscript. The research of EDG and FF has been partly supported by the EURYI (European Young Investigators) Awards Scheme (Eurohorcs/ESF, responsible F.F.) with funds from the National Research Council of Italy and other National Funding Agencies participating in the Third Memorandum of Understanding, as well as from the EC Sixth Framework Programme.