# Dynamics of the Ryukyu/Izu-Bonin-Marianas double subduction system

Claudio Faccenna, Adam Holt, Thorsten Becker, Serge Lallemand, Leigh Royden

► To cite this version:

Claudio Faccenna, Adam Holt, Thorsten Becker, Serge Lallemand, Leigh Royden. Dynamics of the Ryukyu/Izu-Bonin-Marianas double subduction system. Tectonophysics, 2018, 746, pp.229-238. 10.1016/j.tecto.2017.08.011. hal-01891901

HAL Id: hal-01891901

https://hal.umontpellier.fr/hal-01891901v1

Submitted on 13 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
MIT Open Access Articles

*Dynamics of the Ryukyu/Izu-Bonin-Marianas double subduction system*

The MIT Faculty has made this article openly available. **Please share**
how this access benefits you. Your story matters.

**Citation:** Faccenna, Claudio et al. "Dynamics of the Ryukyu/Izu-Bonin-Marianas double subduction system." Tectonophysics 746 (October 2018): 229-238 © 2017 Elsevier B.V.

**As Published:** http://dx.doi.org/10.1016/j.tecto.2017.08.011

**Publisher:** Elsevier BV

**Persistent URL:** https://hdl.handle.net/1721.1/125860

**Version:** Author's final manuscript: final author's manuscript post peer review, without publisher's formatting or copy editing

**Terms of use:** Creative Commons Attribution-NonCommercial-NoDerivs License
# Dynamics of the Ryukyu/Izu-Bonin-Marianas double subduction system

Claudio Faccenna<sup>a,*</sup>, Adam F. Holt<sup>b</sup>, Thorsten W. Becker<sup>c</sup>, Serge Lallemand<sup>d</sup>, Leigh H. Royden<sup>b</sup>

## ARTICLE INFO

### Keywords:

Western Pacific

Subduction zone

Trench migration

Numerical modelling

## ABSTRACT

Trench motions represent the surface expression of the interaction between subducting plates and the underlying mantle, but the inherent dynamics are not fully understood. One interesting case is the migration of the Izu-Bonin-Marianas trench (IBM) that accommodates the subduction of the Pacific beneath the Philippine Sea Plate (PSP), which is in turn subducting beneath the Eurasian plate along the Ryukyu trench. The history of the IBM trench is dominated by fast, episodic retreat from 40 to 15 Ma. However, around 10–5 Ma, the IBM trench reversed its motion from retreat to advance. The switch in trench motion occurred soon after the breakoff of the PSP slab along the Ryukyu trench and the onset of new subduction, and represents a fundamental change in the dynamics of the western Pacific subduction zones. Here, Following the modelling study of Čížková and Bina (2015), which suggested a link between IBM trench advance and Ryukyu subduction, we run 2-D numerical experiments to test the influence of a newly formed Ryukyu slab on the established IBM subduction zone, we run two-dimensional numerical experiments to test the influence of this newly formed Ryukyu slab on the IBM subduction zone. The results from our geodynamic model compare favourably with the reconstructed trend, indicating that the switch in trench motion along the IBM trench may indeed be related to the onset of a new subduction zone along the Ryukyu trench. Our analysis substantiates the idea that advancing trench motions in the western Pacific are due to the establishment of a double subduction system. Further analysis of such dynamics provides insights for the mechanisms controlling subducting plate and trench motions and mantle force transmission.

## 1. Introduction

Seafloor trenches represent the sites where oceanic lithosphere subducts into the mantle. However, motion of the trenches themselves is also an important component of the surface expression of mantle convection. Assuming negligible accretion or erosion at the trench, trench kinematics may be reconstructed by subtracting backarc deformation from upper plate motion (Dewey, 1980). Geologic reconstructions show that trenches are indeed not stationary features but migrate in space and time. Present-day migration patterns indicate that trenches migrate backward (trench rollback) or forward (advancing trench) with respect to the upper plate (Heuret and Lallemand, 2005).

The rate of trench migration depends on the absolute plate motion reference frame. The “correct” choice of reference frame is not always clear, and the way that trench migration is partitioned between advancing and retreating motion also depends on the reference frame (Funiciello et al., 2008; Schellart et al., 2008; Becker et al., 2015). However, in some systems like the northern Philippine Sea region,

where two subduction systems exist in close proximity, the sense of motion of one trench relative to one another and the sense of motion of both trenches relative to the uppermost plate (stable Eurasia) is unambiguous. On the Eastern boundary of the Philippine Plate, the Izu-Bonin Trench is presently advancing toward Eurasia, and on the Western boundary, the Ryukyu Trench is retreating, and the two trenches are converging with time.

Several geological examples illustrate that trench migration is not constant over time and reconstructions of back-arc spreading histories show that the migration of trenches may be episodic and/or short lived (e.g. Taylor and Karner, 1983). The Izu-Bonin-Marianas subduction zone is an excellent example of such a system, as geological reconstructions show that it has undergone a dramatic shift from retreating to advancing trench motion between 5 and 10 Ma.

Many factors may influence the rate and direction of trench migration in single slab subduction systems. These include plate-mantle coupling (e.g., Bellahsen et al., 2005; Faccenna et al., 2009; Funiciello et al., 2008), slab interactions with the upper-lower mantle transition

<footnote><sup>a</sup> Dipartimento Scienze, Università Roma Tre, Rome, Italy</footnote>

<footnote><sup>b</sup> Department of Earth, Atmospheric and Planetary Sciences, M.I.T., Cambridge, MA, USA</footnote>

<footnote><sup>c</sup> Institute for Geophysics and Department of Geological Sciences, Jackson School Geosciences, The University of Texas at Austin, Austin, TX, USA</footnote>

<footnote><sup>d</sup> CNRS, Laboratoire Geosciences Montpellier, Univ. Montpellier II, Montpellier, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: faccenna@uniroma3.it (C. Faccenna).</footnote>
(Kincaid and Olson, 1987; Christensen, 1996), three-dimensional (3-D) plate geometry and plate width (Funiciello et al., 2003; Stegman et al., 2006; Schellart et al., 2007), plate kinematics (Heuret et al., 2007; Lallemand et al., 2008), overriding plate properties (Capitanio et al., 2010; Yamato et al., 2009), and mantle rheology (Holt and Becker, 2017). In multi-slab systems, such as that existing in the Philippine Sea Plate region, the interaction between proximal slabs can have a major impact on trench migration rates, as explored through geologic observations and dynamic subduction models applied to various subduction regions (Jagoutz et al., 2015; Kiraly et al., 2016; Holt et al., 2017). Double slab systems are interesting in their own right but also serve to explore the parameters controlling single slab and trench dynamics, motivating study of natural examples.

Recently, Čížková and Bina (2015) used 2-D subduction models to suggest that the advancing motion of the Izu-Bonin-Mariana Trench, which accommodates west-dipping subduction along the eastern boundary of the Philippine Sea Plate, can be attributed to the presence of a second west-dipping subduction zone located along the western boundary of the Philippine Sea Plate (the Ryukyu Trench). Here, we focus on the temporal evolution of this two-slab subduction system, particularly the mechanism by which a retreating Izu-Bonin-Mariana Trench can transition to an advancing mode of trench migration. Therefore, we examine how the coupled evolution of these two slabs is reflected in the rate and direction of Izu-Bonin-Mariana Trench migration.

We begin with a tectonic reconstruction of the Izu-Bonin-Mariana system during the last 20 Ma followed by numerical tests that capture the first order dynamics of double subduction. An important aspect of this study is our exploration of the switch from trench retreat to advance at the Izu-Bonin-Mariana system, and its dependence on subduction along the Ryukyu Trench.

## 2. Evolution of the Izu-Bonin system

### 2.1. Observed geometry and kinematics

The Izu-Bonin-Marianas (IBM) and Ryukyu trenches accommodate west-dipping subduction along the eastern and western margins of the Philippine Sea Plate (PSP), respectively. The geometry of these subducting plates at depth has been illuminated by several studies of seismic tomography and seismicity (e.g. van der Hilst et al., 1991; van der Hilst and Seno, 1993; Bijwaard et al., 1998; Fukao et al., 2001; Miller et al., 2006; Li et al., 2008; Fukao and Obayashi, 2013). Along the northeast margin of the PSP, the Izu-Bonin slab is characterized by a 50° to 70° west-dipping Wadati-Benioff zone, which flattens into the transition zone as indicated by a long, high velocity anomaly that sits just above the 660-km discontinuity. Northward, below Japan, the slab shallows to a dip of about 30° whereas, southward, below the Mariana Trench, the slab becomes nearly vertical, and even overturns locally (Fig. 1). The Ryukyu Trench is currently retreating eastward, producing back-arc extension in the Okinawa Trough. The IBM Trench presently advances to the west (Fig. 1). This kinematics pattern is present in the reference frame that fixes Eurasia, which represents the upper plate of the subduction system, and in the other reference frames commonly adopted (Funiciello et al., 2008; Di Giuseppe et al., 2009). The magnitude of trench motion depends on the reference frame; higher values of advance for the IBM Trench and lower values of retreat for the Ryukyu Trench are obtained using a Pacific hotspot reference frame such as HS3 (~6 cm/yr and ~1 cm/yr, respectively; Gripp and Gordon, 2002), relative to a No-net-rotation (NNR) reference frame (~2 cm/yr and ~6 cm/yr respectively; DeMets et al., 2010) or to spreading-aligned absolute reference frame (Fig. 1, ~4 cm/yr and ~2 cm/yr, respectively; Becker et al., 2015). Accordingly, the subduction velocity at Ryukyu Trench, which depends upon the IBM and Ryukyu Trench velocity, varies from ~7 cm/yr in HS3 (Gripp and Gordon, 2002), ~2 cm/yr in NNR (DeMets et al., 2010) or to ~4 cm/yr spreading-

aligned absolute reference frame (Fig. 1, Becker et al., 2015).

The tectonic evolution of the region is here summarized using a schematic reconstruction (Fig. 2) and a table listing the main tectonic events (Table 1). The Philippine Sea Plate is surrounded by subduction zones, and has been formed by sea-floor and back-arc spreading along at least three different axes (Figs. 1, 2). The first and the oldest portion of the plate, the West Philippine Basin, formed during the Eocene from 54 to 30 Ma (Deschamps and Lallemand, 2002 and reference therein; Fig. 2a; Table 1). At ~30 Ma, spreading moved to form the Shikoku and Parece-Vela basins, and spreading occurred until ~15 Ma (Chamot-Rooke et al., 1987; Okino et al., 1994, 1998; Sdrolias et al., 2004). Three spreading stages have been detected in the Shikoku and Parece-Vela basins as the spreading direction changed from E-W to NE-SW at ~20 Ma (Okino et al., 1994, 1998; Sdrolias et al., 2004). The change in spreading direction may be related to clockwise rotation of the Philippine Sea Plate (Sdrolias et al., 2004). The half spreading rate ranged from 2 to 3 cm/yr to 5-6 cm/yr and was accompanied by a total of trench retreat of ~1000 km along the northern portion of the Izu-Bonin Trench. The third extensional episode, which continues today, occurred behind the Mariana Trough, which began rifting at ~7 Ma (Hussong and Uyeda, 1981; Fryer, 1995; Stern et al., 2003; Kato et al., 2003; Martinez et al., 1995).

The evolution of the subduction system is complex due to initiation and cessation of subduction in some places and northward migration of the Ryukyu-Izu Bonin triple junction. Subduction of the Pacific plate at the IBM Trench began at ~54 Ma (Cosca et al., 1998; Ishizuka et al., 2008, 2011; Deschamps and Lallemand, 2002; Ishizuka et al., 2011; Lallemand, 2016). Because of the northward migration of the Ryukyu-IBM triple junction, subduction in the northern side of the Ryukyu was first produced by the Pacific and later progressively replaced by the PSP (Fig. 2b). The subduction of the PSP was interrupted by the entrance at trench of the northern extent of the Gagua Ridge which is a paleo-transform plate boundary separating a Cretaceous deep lithosphere (Huatung Basin; Deschamps et al., 2000) from the shallow Eocene lithosphere (West Philippine Basin). Following Lallemand et al. (2001) this transform fault should have entered at trench producing an episode of slab break-off (Lallemand et al., 2001; Fig. 2c; Table 1). Following slab break-off, subduction of the Philippine Sea Plate was reinitiated along the Ryukyu Trench, resulting in the modern Ryukyu subduction system (Lallemand, 2016). As discussed in Section 2.2, the timing of slab break off can be reconstructed by the length of the present-day Ryukyu slab and should have occurred around ~10 Ma. Just after this episode of slab break-off, a fundamental change in kinematics of trench motion occurred around 5-8 Ma. Prior to this time, the IBM Trench underwent (eastward) trench retreat, but afterwards changed direction and began to advance westward (Carlson and Melia, 1984). The total magnitude of the IBM Trench advance is ~70 km, from ~10 Ma to the present day (Seno and Maruyama, 1984). Paleomagnetic data also provide constraints on the timing of both clockwise rotation and on the northward drift of the PSP (Hall et al., 1995), indicating trench advance over the last 5 Myr (Hall, 2002). Le Pichon and Huchon (1987) estimated that at 3 Ma the triple junction between Ryukyu, Izu, and Japan was located to the west of the present-day position by ~30 km, moving 30 km eastward, hence advancing, since that time.

The switch in the direction of trench migration, from retreating to advancing, represents a fundamental change in the dynamics of the western Pacific subduction zones (Seno and Maruyama, 1984; van der Hilst and Seno, 1993). It has been associated with a re-organization of major tectonic boundaries along the western margin of the Pacific margin and is probably associated with an increase in the speed of the Pacific Plate (Faccenna et al., 2012).

### 2.2. Reconstruction of the northern Philippine Sea Plate region

To better visualize the evolution of the subduction, we reconstruct the Ryukyu, Izu-Bonin system along a section at a latitude of 28° N,
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Present-day plate kinematic setting of the western Pacific, showing the subducting slab contours at 50 km depth intervals (RUM model of Gudmundsson and Sambridge, 1998), plate and trench velocities (plate: gray vectors, MORVEL-NNR56, Argus et al., 2011, colored: trench motions from Heuret and Lallemand, 2005), both in the spreading-aligned absolute reference frame of Becker et al. (2015). Also shown is the sea-floor age from Müller et al. (2008). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<caption>**Fig. 2.** Simplified tectonic evolution of the region in five steps from 20 Ma onward (modified from Deschamps and Lallemand, 2002; plate motion from Seton et al., 2012). The double dashed line indicates the Gagua Ridge (Fig. 2c-e from Lallemand, 2016), the single dashed line the position of the active spreading centre and the single dashed line indicates the trace of the cross section of Fig. 3. The arrow indicates the motion of the plate from Seton et al. (2012).</caption>
<caption>**Table 1**</caption>

List of the main tectonic events in and around the IBM and Ryukyu Trench. Reference number: 1. Cosca et al. (1998); 2. Ishizuka et al. (2008), 3. Ishizuka et al. (2011); 4. Deschamps and Lallemand (2002); 5. Chamot-Rooke et al. (1987); 6. Okino et al. (1994); 7. Okino et al. (1998); 8. Hussong and Uyeda (1981); 9. Fryer (1995); 10. Stern et al. (2003); 11. Kato et al. (2003); 12. Kizaki (1986); 13. Shinjo (1999); 14. Letouzey and Kimura (1985); 15. Kimura (1996); 16. Park (1996); 17. Sibuet et al. (1998); 18. Tamaki et al. (1992); 19. Kimura et al. (2005); 20. Kimura et al. (2014); 21. Lallemand et al. (1998); 22. Barrier and Angelier (1986); 23. Lu and Hsu (1992); 24. Huang et al. (1997); 25. Lallemand et al. (2001); 26. Malavieille et al. (2002); 27. Seno and Maruyama (1984); 28. Carlson and Melia (1984); 29. Huchon et al. (1986); 30. Jolivet et al. (1989); 31. Shimizu and Itaya (1993); 32. Lallemand (2016); 34. Huang et al. (2013); 35. Cao et al. (2014)

<table><thead><tr><th>Ages (read "from 20 to 22 Ma" for 20)</th><th>54</th><th>52</th><th>30</th><th>28</th><th>26</th><th>24</th><th>22</th><th>20</th><th>18</th><th>16</th><th>14</th><th>12</th><th>10</th><th>8</th><th>6</th><th>4</th><th>2</th><th>0</th><th>References</th></tr></thead><tbody><tr><td><strong>Major regional-scale geodynamic event</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>IBM arc volcanism</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1, 2, 3</td></tr><tr><td>West Philippine Basin opening</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>4</td></tr><tr><td>Shikoku &amp; Parece Vela back-arc basins opening</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>5,6,7</td></tr><tr><td>Mariana Trough opening</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>8,9,10,11</td></tr><tr><td>Central Ryukyu arc volcanism</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>12,13</td></tr><tr><td>Okinawa Trough rifting</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>14,15,16,17</td></tr><tr><td>Japan Sea opening</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>18</td></tr><tr><td>SW Japan clockwise rotation</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>19,20</td></tr><tr><td>SW Japan forearc volcanism</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>19,20</td></tr><tr><td>Cessation of PSP subduction beneath SW Japan</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>19,20</td></tr><tr><td>Subduction at Philippine Trench</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>21</td></tr><tr><td>Arc-continent collision in Taiwan</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>22,23,24,25,26</td></tr><tr><td>PSP motion change</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>27,28,29,20,31</td></tr><tr><td>Ryukyu slab break-off</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>25,26</td></tr><tr><td>Shikoku slab break-off</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>27</td></tr><tr><td>Kyushu-Palau Ridge slab tearing</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>28,29,26</td></tr></tbody></table>

adopting the plate velocity model of Sdrolias and Müller (2006), and tectonic setting based mainly on the reconstruction of Deschamps and Lallemand (2002). This reconstruction is simplified and intended to give a schematic representation of the evolution and interaction of the Izu-Bonin and Ryukyu slabs, and therefore does not include the three-dimensional complexity of the subduction setting (Wu et al., 2016). The geometry of the Izu-Bonin and Ryukyu slabs is interpreted from tomographic images (e.g., Li et al., 2008; Fig. 3a) and seismicity. In our interpretation, the length of the present-day Ryukyu slab approaches the transition zone (Fig. 3a), confirming previous interpretations (Bijwaard et al., 1998; Widiyantoro et al., 1999; Wang et al., 2008; Li and van der Hilst, 2010; Wei et al., 2012, 2015) and is detached from a flat high velocity body lying on the transition zone in front of the Ryukyu Trench.

At the latitude of our section, prior to 25 Ma, the Pacific Plate was subducting directly beneath the paleo-Ryukyu Trench. Subduction of the Pacific Plate continued here until the triple junction, where the paleo-Izu Bonin and paleo-Ryukyu trenches intersected, migrated north past 28° N, the latitude of our reconstructed section. Therefore, along this cross-section, there was first subduction of the Pacific Plate and later subduction of the Philippine Sea Plate (Fig. 3g, f). The Pacific slab at that time was probably quite long, likely with a shallow dip, as is the case for the Pacific slab segment presently subducting beneath Japan (Figs. 2a, 3g). Two to three high P-wave velocity anomalies that are present in the transition zone in front of the Ryukyu slab are probably related to this older Pacific Plate subduction phase (Li et al., 2008). At the latitude of our section, immediately after passage of the triple junction, the two slabs were probably very close to one another, then progressively separated during the opening of the Shikoku Basin (Figs. 2b, 3e). This was followed by probable break-off of the Philippine Sea slab (Figs. 2c, 3d) along the northern extent of the Gagua Ridge paleo-plate boundary (Lallemand et al., 2001). Following slab break-off, subduction of the Philippine Sea Plate was reinitiated along the Ryukyu Trench, resulting in the modern Ryukyu subduction system (Lallemand, 2016). The timing of the slab break-off can be deduced by the length of the present-day slab, which is ~400 km divided by subduction rate. Assuming that slab break-off occurred at a depth of

100 km the amount of post-break off subduction is 300 km. The present day Ryukyu subduction rate is around ~4–5 cm/yr as it is the IBM Trench advancing rate. However, during the last 10 Ma, the IBM Trench advance rate was lower than the present day one (~4–5 cm/yr), as it switched from retreat to advance, in the order of 2–3 cm/yr. Correspondingly, also the average Ryukyu subduction was lower than the present-day and of the same order. We may then estimate that timing of slab break off occurred around ~10–12 Ma.

A relevant question regarding the dynamics of the PSP system concerns the switch from trench retreat to trench advance of the IBM Trench and its relationship with the Ryukyu Trench. A range of models have been proposed. Van der Hilst and Seno (1993) suggested that the overall clockwise rotation of the Izu Bonin-Mariana Trench is related to the on-going penetration of the Mariana slab into the lower mantle. These authors proposed that the recent advance of the Izu-Bonin-Mariana Trench could be related to the northward propagation of this process, and that trench advance is about to initiate along the northern Izu-Bonin Trench region. Miller et al. (2004, 2006) attributed the change of the motion of the Izu-Bonin Trench to the subduction of the Necker-Marcus (Ogasawara Plateau) aseismic ridge, at ~8 Ma. Faccenna et al. (2009) proposed that the switch from retreat to advance was induced by the increase in seafloor age at the trench as suggested by reconstructions of global trench migration vs. age of the subducting oceanic lithosphere, and both laboratory and numerical experiments (Faccenna et al., 2007; Di Giuseppe et al., 2009).

Other studies have linked the switch in trench motion to the geometry of double subduction around the PSP. Carlson and Melia (1984) proposed that the advancing motion of the trench could be related to suction exerted by the fast retreating (Ryukyu) upper plate. Čížková and Bina (2015) proposed that the suction of the Ryukyu slab produced advancing motion of the Izu-Bonin slab. However, the reason why trench motion should have changed from retreating to advancing during the last 8–5 Ma is not clear. Building on the work of Čížková and Bina (2015), we therefore explore the interaction between the Ryukyu and IBM slabs as a cause for the change in motion of the IBM Trench through time, considering that the change from retreating to advancing motion of the IBM post-dated Gagua Ridge subduction, and possible
<figure><img src="image_3.png" /><figcaption>Fig. 3. Possible reconstruction of along a cross section at 26° N (Fig. 2) of the mantle structure. a) Tomography model from Li et al. (2008); b–g) six stage evolution of the subduction zone. The subduction zone has been restored backward in time considering the position of the trench and the subduction velocity as from plate motion of Seton et al. (2012). PSP is Philippine Sea Plate.</figcaption></figure>

slab break off and onset of modern Ryukyu subduction (Fig. 2c, d), by a few million years. We use numerical models of the single and double subduction process to investigate the feasibility of a double subduction induced switch in the trench motion direction. We additionally investigate the mechanical controls on both the onset time of the switch, and the eventual magnitude of the trench advance velocity.

### 3. Model set up

We use dynamically consistent (i.e., no external forces/velocities imposed) thermo-mechanical numerical models to examine whether the initiation of an additional subduction zone can induce a transition from a retreating trench to an advancing trench. Our approach is to run a standard single subduction zone mode, in which the slab represents that of the Izu-Bonin, to a time-step at which mature subduction is

established (for the Pacific slab at ~10 Ma), and then insert an additional subduction zone corresponding to modern Ryukyu. We then compare the subsequent trench evolution of IBM in the single slab model (i.e., without the model modification), with a model that contains an added subduction zone (i.e., the modified model). We use the finite-element code CitcomCU (Moresi and Gurnis, 1996; Zhong, 2006) to model subduction in a Cartesian domain. CitcomCU solves the equations governing convection in an incompressible viscous fluid with negligible inertia (Boussinesq and infinite Prandtl number approximations).

The domain size is 10,560 km (x) × 1320 km (z) × 13.2 km (y, i.e. the computations are pseudo-2-D), and all boundaries are free slip. Finite element dimensions vary from 6 km in the refined, plate boundary region to 14.5 km elsewhere. The refined region extends from x = 4225 km to x = 7400 km, for depths shallower than 400 km. We parameterize the 660 km phase transition with a factor of 50 viscosity increase, compatible with geoid observations (e.g. Hager, 1984; Mitrovica and Forte, 2004), but omit buoyancy effects and additional phase transitions. Both the subducting and overriding lithospheric plates have initial lengths of 4000 km, and are initially placed at a horizontal distance of 1280 km away from the side boundaries.

For the single slab model, we use a model setup similar to that presented in Holt et al. (2015), their Fig. 1, but with a high viscosity lower mantle (see Table 2 for all model parameters). In order to focus on first-order slab-slab interactions, we use a simplified subduction setup, both mechanically and rheologically. We consider an initially constant temperature (1573 K) and therefore constant viscosity (10²⁰ Pa s) sub-lithospheric mantle. Thermally, we impose constant temperature (273 K) on the upper, surface boundary, and zero heat flux on the base and side boundaries.

The initial thickness of the subducting and overriding plates is 89 km and 59 km, respectively, and these plates are initially isothermal. These lithospheric thicknesses were chosen to provide depth-integrated buoyancies, B, equivalent to those expected for the average ages of the subducting Pacific (80 Ma age gives B = 6.1 × 10⁶ kg/m²) and Philippine Sea plates (35 Ma age yields B = 4.0 × 10⁶ kg/m²),

<table><caption>**Table 2**

Geometrical, rheological and mechanical parameters of the numerical models.</caption><thead><tr><th>Symbol</th><th>Parameter name</th><th>Reference value (others tested)</th><th>Units</th></tr></thead><tbody><tr><td colspan="4"><strong>Geometrical</strong></td></tr><tr><td>h</td><td>Domain height</td><td>1320</td><td>km</td></tr><tr><td>l</td><td>Domain length</td><td>10,560</td><td>km</td></tr><tr><td>h<sub>PP</sub></td><td>Lithospheric thickness (Pacific Plate)</td><td>89</td><td>km</td></tr><tr><td>h<sub>PSP</sub></td><td>Lithospheric thickness (Philippine Plate)</td><td>59</td><td>km</td></tr><tr><td>h<sub>crust</sub></td><td>Crustal thickness</td><td>15</td><td>km</td></tr><tr><td>L<sub>SP</sub></td><td>Initial subducting plate length</td><td>4000</td><td>km</td></tr><tr><td>L<sub>OP</sub></td><td>Initial overriding plate length</td><td>4000</td><td>km</td></tr><tr><td>z<sub>notch</sub></td><td>Initial slab depth</td><td>150</td><td>km</td></tr><tr><td>R<sub>notch</sub></td><td>Initial slab radius of curvature</td><td>150</td><td>km</td></tr><tr><td>Δx</td><td>Initial trench-trench distance</td><td>1500 (1250, 1750, 2250)</td><td></td></tr><tr><td colspan="4"><strong>Mechanical and rheological</strong></td></tr><tr><td>ΔT</td><td>Plate-asthenosphere temperature contrast</td><td>1300</td><td>K</td></tr><tr><td>α</td><td>Thermal expansivity</td><td>1.6 x 10<sup>-5</sup></td><td>K<sup>-1</sup></td></tr><tr><td>ρ<sub>o</sub></td><td>Asthenosphere density</td><td>3300</td><td>kg/m<sup>3</sup></td></tr><tr><td>Δρ</td><td>Plate-asthenosphere density contrast</td><td>68</td><td>kg/m<sup>3</sup></td></tr><tr><td>η<sub>lith</sub></td><td>Lithospheric viscosity cut-off</td><td>5 x 10<sup>22</sup></td><td>Pa s</td></tr><tr><td>η<sub>mant</sub></td><td>Asthenospheric viscosity</td><td>1 x 10<sup>20</sup></td><td>Pa s</td></tr><tr><td>η<sub>crust</sub></td><td>Crustal viscosity</td><td>5 x 10<sup>19</sup> (1 x 10<sup>19</sup>, 2 x 10<sup>20</sup>)</td><td>Pa s</td></tr><tr><td>λ</td><td>Byerlee pore-pressure prefactor</td><td>0.15 (0.1)</td><td>—</td></tr><tr><td>a</td><td>Byerlee coefficient of friction</td><td>0.6</td><td>—</td></tr><tr><td>b</td><td>Byerlee cohesion</td><td>60</td><td>MPa</td></tr></tbody></table>
computed using a half-space cooling model. The initial lithosphere-asthenosphere temperature contrast is set to 1300 K. To generate the negative buoyancies stated above, we therefore use a thermal expansion of $1.59 \times 10^{-5} \text{ K}^{-1}$ and a reference mantle density of 3300 kg/m$^3$.

We adopt a lithospheric rheology similar to that described in Holt et al. (2015), where the effective viscosity is composed of the harmonic mean of a Newtonian and a plastic component. The Newtonian component of the lithospheric viscosity is controlled by a maximum viscosity cut-off and is factor 500 more viscous than the asthenosphere. For the plastic component, we use a Byerlee type yield stress, $\tau_{yield} = \lambda(a \sigma_L + b)$, where a is the coefficient of friction (0.6), b is the cohesion (60 MPa), $\sigma_L$ is the lithostatic pressure ($= \rho_0 g z$), and $\lambda$ is a “pore pressure” pre-factor (0.15). This pre-factor value is broadly consistent with that used in previous studies (e.g. Enns et al., 2005: 0.1; Di Giuseppe et al., 2008: 0.08), and is low enough to trigger plastic weakening in the bending region, yet high enough to avoid slab break-off. In Supplementary material, we also examine the trench motion systematics produced by models with either a purely viscous rheology or a visco-plastic rheology with a reduced yield stress, $\lambda = 0.1$ (Fig. S1).

To decouple the two plates, we insert a 15 km thick crustal layer within the subducting plate, which has a constant viscosity of $5 \times 10^{19} \text{ Pa s}$ (e.g. Holt et al., 2015). As the decoupling layer strength has been shown to exert a strong control on slab dynamics (e.g., Čížková and Bina, 2013), we have tested that stronger and weaker layers do not change the first order systematics (Fig. S2). To initiate subduction, we allow the subducting plate to extend below the overriding plate to a depth ($z_{init}$) of 150 km, with an initial radius of curvature ($R_{init}$) of 150 km (Fig. 4a).

In our reference model, we choose a model time of 5.4 Myr to insert Ryukyu (at a distance of 1500 km from IBM), which corresponds to the instance when the slab first impinges on the strong lower mantle (cf. Fig. 4b, c). Such a mature model subduction stage is chosen following our reconstruction (Figs. 2, 3). The Ryukyu subduction zone is added to the model viscosity and density fields as a proto-slab with the same curvature as that used to initiate IBM subduction (i.e. $R_{init} = z_{init} = 150$ km), but reduced subducting plate thickness corresponding to the younger PSP. In addition, a weak crustal, decoupling

layer is added onto the middle plate, extending from the new notch to a distance of ≈ 50 km to the right of the mature trench. This decouples the Ryukyu slab from the leftmost, Eurasian, plate. At the time-step at which Ryukyu is inserted, there is a discontinuity in the model evolution and so the trench velocity is taken to be the mean of that in the single slab model, and that in the double slab model. This modified subduction geometry is then allowed to evolve dynamically, and then the rear slab (i.e. IBM) behavior is compared with that observed in the unperturbed, single slab model (cf. Fig. 4b, c).

<figure><img src="image_4.png" /><figcaption>Fig. 5. Trench motion rate as a function of time for the reference single and double slab models (as in Fig. 4). At the insertion time-step, the trench retreat rate is taken as the mean average between the single slab model, and the modified model with an additional, incipient subduction zone.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Viscosity field plots illustrating the evolution of our reference, b) single slab, and, c) double slab models. Panel a) depicts the initial, single slab viscosity field and velocity solution. The top panels with b) and c) (i.e. $t = 5.4$ Myr) show the time at which the additional subduction zone is inserted for the double slab model. Subsequent panels show the contrasting evolution of the slab morphology for the single and double slab models, and black arrows indicate the trench motion directions. Plate velocities are shown in the bottom panels, with positive values indicating motion toward the right.</figcaption></figure>
4. Modelling results

In the reference single slab model, with a visco-plastic rheology (Fig. 4b), the subducting plate enters the mantle through a combination of subducting plate motion and trench retreat, as is typically observed in similar dynamic subduction models (e.g. Enns et al., 2005; Stegman et al., 2006). In the reference frame of the model side boundaries, the subducting plate velocity increases to a maximum rate of 102 mm/yr before the slab impinges on the lower mantle. Simultaneously, the trench retreats rapidly, increasing from an initial rate of ≈ 30 mm/yr to a maximum of ≈ 80 mm/yr (Fig. 5). As the slab interacts with the high viscosity lower mantle, the slab is folded and finally flattens while penetrating into the lower mantle slowly. This lower mantle interaction causes a reduction in both the plate and trench velocities. The subducting velocity reduces to ≈ 60 mm/yr, and the trench retreat velocity gradually reduces throughout the subsequent duration of the model, by ≈ 2.7 mm/yr per Myr.

In the model within which we insert an additional slab (Fig. 4c), the style of trench motion is very different (Fig. 5). The additional slab (Ryukyu), which itself exhibits rapid trench retreat (mean trench velocity, $V_T$ ≈ 118 mm/yr), exerts a strong influence on the behavior on the mature, rear slab (IBM). When the additional notch is initially inserted, the Izu-Bonin slab remains in the retreating regime. However, as the Ryukyu slab subducts progressively deeper into the asthenosphere, and the two slabs become closer together, the amount of slab pull acting on the middle plate (PSP) increases and the IBM slab switches from a retreating to an advancing mode. This switch occurs ~6 Myr after Ryukyu subduction initiation, and our modelled Izu-Bonin slab advances at rates of up to 23.2 mm/yr (Fig. 5). Because trench retreat tends to reduce, and advance tends to increase slab dip (e.g., Christensen, 1996), we observe near vertical rear-slab dips, in the models containing an additional slab (Fig. 4c).

In addition to inducing a switch from a retreating to an advancing trench, the initiation of the second subduction zone increases the subduction velocity of the rear plate (Pacific), in the reference frame of the model box (Fig. S3). For equivalent time periods (after the rear slab has reached the lower mantle), the mean rear plate velocity is ≈ 58 mm/yr in the single slab model, and increases to ≈ 72 mm/yr in the model with an additional, inserted subduction zone.

In Supplementary material, we show that the trench motion systems observed for the reference model are unchanged by the specific rheology assumed for the lithospheric plates. For lithospheric plates with both a viscous rheology and a reduced yield stress ($\lambda$ = 0.1), we also observe a switch from a retreating to an advancing Pacific subduction induced by Ryukyu subduction (Fig. S1). Motivated by the inherent uncertainty regarding the appropriate viscosity for parameterizing the crustal, decoupling layer (e.g., Androvičová et al., 2013), we also show that the first order trench motion trends are unchanged in models which employ either a stronger ($\eta_{crust} = 2\eta_{mantle}$) or a weaker ($\eta_{crust} = 0.1\eta_{mantle}$) crustal layer to decouple the plates (Fig. S2).

Because the PSP plate is highly 3-D, with variable trench-perpendicular distance between the two slabs/trenches as a function of along-strike location, we have conducted additional tests to examine the effect of inserting Ryukyu subduction at a range of initial, horizontal distances from the Pacific subduction zone ($\Delta x$ = 1250, 1500, 1750, 2250 km), appropriate for the PSP geometry. The largest trench-trench distance of 2250 km is approximately equivalent to the widest portion of the Pacific plate double slab system, from the centre of the Marianas Trench to the centre of the Ryukyu Trench. Fig. 6a shows the evolution of trench motion velocities for these four double slab cases. For each case, as the inserted subducting slab becomes more mature, and the two trenches become closer together (primarily due to Ryukyu Trench rollback), the trench retreat rate of the modelled IBM Trench reduces until it begins to advance. Increasing initial slab separation increases the time it takes for the trenches to become a certain distance apart, and

so delays the time at which a shift from trench retreat to trench advance occurs. Going from an initial trench separation of 1250 km to an initial trench separation of 2250 km, this trench motion transition time is delayed by ≈ 8 Myr. However, we do ultimately observe a switch from Pacific Trench retreat to advance for all of the initial trench separations. Similarly, we find that an acceleration of the modelled Pacific Plate (relative to the plate velocity observed in the single slab models) occurs for all of the initial trench separations (Fig. S3).

As the precise morphology of the Pacific slab at the time of Ryukyu slab break-off and re-initiation is not known, we also examine models within which Ryukyu subduction is initiated at different stages in the evolution of Pacific subduction. Fig. 6b shows the effect of initiating Ryukyu subduction at a time significantly before the Pacific slab tip reaches the lower mantle (model time ≈ 2.2 Myr), and a time after the slab tip has had significant interaction with the lower mantle (time ≈ 7.4 Myr), i.e. before and after the reference case shown in Figs. 4 and 5. Initiating subduction at a later stage in the Pacific subduction does not modify the trench motion systematics. However, initiating Ryukyu subduction before the Pacific slab touches the lower mantle causes the magnitude of the double slab effect on trench motion to be reduced (Fig. 6b). While Ryukyu subduction still causes both a drastic reduction in the trench retreat rate and eventual switch to trench advance, the switch is substantially delayed and the attained trench advance is at a reduced maximum rate of 12 mm/yr (Fig. 6b). In general, an important result of all of the tests described is that, despite at what trench-trench distance or relative time Ryukyu is inserted to the system (Fig. 6), the switch from a retreating to an advancing IBM Trench occurs only after the tip of the IBM slab is sufficiently supported by the stronger lower mantle.

5. Discussion and conclusions

The Izu-Bonin area is a key region for studying the dynamics of trench migration. Čížková and Bina (2015) first analyzed the influence of a double subduction system on trench migration and applied that to the present-day Izu-Bonin-Marianas setting in order to explain the advancing motion of the Izu-Bonin. Here, we combine a geological reconstruction of the region and idealized numerical models to explain why, between 10 and 5 Ma, trench motion here changes from retreating to advancing. Our tectonic reconstruction and numerical modelling suggest that the change in trench motion is related to initiation of the most recent phase of subduction along the Ryukyu Trench. In our reconstruction, this latest phase of subduction follows a period of slab break off that occurred after subduction of the Gagua Ridge. At any one point, passage of the triple junction between the paleo Ryukyu and Izu Bonin trenches is followed by ridge subduction, slab break off and, finally, re-initiation of modern Ryukyu subduction.

For an idealized model geometry, our numerical models demonstrate the ability of an incipient subduction zone (Ryukyu) to exert forces strong enough to change a nearby mature subducting slab (Izu-Bonin) from a retreating to an advancing trench mode. This occurs because the newly initiated subduction boundary causes subduction of the overriding plate of the Izu-Bonin system. Slab pull from the frontal, Ryukyu slab is transmitted through the plate separating the Philippine Sea and Pacific plates, producing a horizontal force which acts to pull the rear trench (IBM) toward the front trench (Ryukyu), producing advance of the IBM Trench. At a slab depth of 660 km, our modelled Ryukyu slab has an available slab pull force of ≈ 25 TN/m and we find that this is sufficient to ultimately trigger a transition in the trench motion direction of the mature IBM Trench. Upon initiation of Ryukyu subduction in our reference model (Fig. 4c), it takes ~6 Myr for the trench to start to advance. This delay occurs because the Ryukyu subduction must be mature enough to be able to transmit a large slab pull force (proportional to slab depth), and the slab pull of the rear IBM slab must be strongly supported by the stronger lower mantle. This modelled delay time is broadly consistent with the reconstructed time delay
<figure><img src="image_6.png" /><figcaption>Fig. 6. Trench motion rate as a function of time for models with various, a) initial trench separations, and, b) double slab (Ryukyu) initiation times. Additionally indicated in panel b is the period during which the rear (Pacific) is sinking through the lower viscosity upper mantle.</figcaption></figure>

between spreading ridge subduction and Izu-Bonin Trench advance.

We also examine the effect of both the initial trench-trench distance (as this is highly variable for the present Philippine plate geometry) and the stage of IBM subduction at which Ryukyu subduction is initiated (Fig. 6). For a range of initial trench separations up to 2250 km (the maximum present day Ryukyu-Izu Bonin separation), all models exhibit a switch to Izu-Bonin trench advance. However, the time lag between the initiation of Ryukyu subduction and the switch to IBM Trench advance is greater for larger initial trench separations. Because the middle PSP is gradually consumed during subduction, the slabs become closer through time, and so the effectiveness of Ryukyu slab pull transmission through the middle Philippine Sea Plate increases. This time lag for larger initial trench-trench distances can therefore be understood by considering that the two slabs/trenches must be at or below a threshold distance apart (~600 km for our model setup) for the required force transmission. In other words, it is the time delay at which the trench motion switch occurs that is a strong function of the initial slab separation, rather than whether advancing trench motion ultimately occurs.

We find that Ryukyu subduction initiated at an earlier stage of Pacific Plate subduction ($t \approx 2.2$ Myr) induces a switch from a retreating to an advancing IBM Trench at a much later time of $\approx 20$ Myr after subduction initiation (Fig. 6b). However, Ryukyu subduction that is initiated during, or after, slab folding atop the stronger lower mantle results in strong advance of the Izu-Bonin Trench 5-7 Myr after subduction is initiated (Fig. 6b). This short delay is consistent that inferred from our tectonic reconstruction, and a mature/deep Pacific slab geometry is also consistent Ryukyu re-initiation at ~10 Ma (given an Izu-Bonin-Marianas initiation at ~50 Ma and average convergence rate of > 50 mm/yr). In the model in which Ryukyu is inserted at an earlier stage, we find that the IBM slab tip impinges on the lower mantle viscosity increase with an increased dip angle (due to a reduced trench retreat rate after Ryukyu insertion), which causes significant slab folding at depth (e.g. Christensen, 1996). This slab folding is apparent as an oscillation in the trench migration rate evolution (Fig. 6b), and appears to result in the time delay of the switch between IBM Trench retreat and advance. In addition to the model requirement that the IBM slab is in contact with the strong lower mantle, this further implies that a switch from a retreating to an advancing trench is not solely a function of transmitted slab pull, but also strongly depends on other factors

influencing the IBM force balance such as slab interaction with a higher viscosity lower mantle or a phase transition.

While our 2-D model is able to produce the first order Izu-Bonin Trench motion trends, we suggest that the three-dimensional geometry of the system may significantly complicate the evolution. As the numerical models are 2-D, mantle flow through the lower mantle is required to accommodate the motion of the two slabs toward each other as the middle plate is subducted. However, in a 3-D setting, toroidal flow around slab edges may be a more energetically favourable way of evacuating asthenosphere from this central region, and this variability in mantle flow style may influence the plate kinematics (e.g., Funiciello et al., 2003; Schellart, 2004). Indeed, single slab subduction models have demonstrated that 3-D mantle flow can play a large role in dictating plate kinematics with, for example, narrower plates tending to produce more rapid rates of trench retreat (e.g., Stegman et al., 2006; Schellart et al., 2007). While similar 3-D effects are likely important in a Philippine Plate style geometry, particularly at the northern terminus of the Ryukyu/Nankai slab, the complex slab geometry of this region (e.g., highly variable trench separation) complicates the direct application of these studies.

While we therefore choose to focus on the less computationally demanding 2-D models, previous work has shown that the same first order plate kinematics are present in a 3-D double slab model, in which rollback-driven mantle flow is almost exclusively toroidal (cf. Fig. 5b in Holt et al., 2017). In the 3-D models, as the two trenches/slabs converge, there is an increase in both the extensional stress in the middle plate and the mantle pressure between the two slabs. These first order patterns are also present in our 2-D models (e.g. Fig. S4), albeit with somewhat different stress magnitudes owing to the differing model setups (including 2-D vs. 3-D geometry). Therefore, even in our simplified 2-D picture, the trench motions obtained by the geodynamic model can, to first order, be compared favourably with the reconstructed trend. Our analysis substantiates the idea that advancing trench motions in the western Pacific are due to the establishment of a double subduction system. Further joint analysis of geological constraints and geodynamic models of subduction may elucidate the dynamics of both isolated and regionally interacting slab systems.
Acknowledgement

Evgenii Burov dedicated several papers to the problem of migration of trenches in the presence of a complex lithospheric rheology. His vision about these and other geodynamical processes profoundly influenced our understanding on the way lithosphere and mantle interacts over geological time scale. We dedicate this paper to him, and thank the editor, associate editor, and reviewers for their helpful comments that improved a previous version of this manuscript. TWB and AFH were partially supported by NSF EAR 1215720.

Appendix A. Supplementary data

Supplementary data to this article can be found online at http://dx.doi.org/10.1016/j.tecto.2017.08.011.

References

Androvičová, A., Čížková, H., van den Berg, A.P., 2013. The effects of rheological decoupling on slab deformation in the Earth's upper mantle. Stud. Geophys. Geod. 57, 460-481.

Argus, D.F., Gordon, R.G., DeMets, C., 2011. Geologically current motion of 56 plates relative to the no-net rotation reference frame. Geochem. Geophys. Geosyst. 12, Q11001. http://dx.doi.org/10.1029/2011GC003751.

Barrier, E., Angelier, J., 1986. Active collision in eastern Taiwan: the Coastal Range. Tectonophysics 125 (1-3), 39-72.

Becker, T.W., Schaeffer, A.J., Lebedev, S., Conrad, C.P., 2015. Toward a generalized plate motion reference frame. Geophys. Res. Lett. 42, 3,188-3,196.

Bellaھs, N., Faccenna, C., Funiciello, F., 2005. Dynamics of subduction and plate motion in laboratory experiments: insights into the "plate tectonics" behavior of the Earth. J. Geophys. Res. 110. http://dx.doi.org/10.1029/2004JB002999.

Bijwaard, H., Spakman, W., Engdahl, R., 1998. Closing the gap between regional and global travel time tomography. J. Geophys. Res. 103, 30055-30078.

Cao, L., Wang, Z., Wu, S., Gao, X., 2014. A new model of slab tear of the subducting Philippine Sea Plate associated with Kyushu-Palau Ridge subduction. Tectonophysics 636, 158-169. http://dx.doi.org/10.1016/j.tecto.2014.08.012.

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

Carlson, R.L., Melia, P.J., 1984. Subduction hinge migration. Tectonophysics 102, 1-16.

Chamot-Rooke, N., Renard, V., Le Pichon, X., 1987. Magnetic anomalies in the Shikoku basin: a new interpretation. Earth Planet. Sci. Lett. 83, 214-228.

Christensen, U., 1996. The influence of trench migration on slab penetration into the lower mantle. Earth Planet. Sci. Lett. 140, 27-39.

Čížková, H., Bina, C., 2013. Effects of mantle and subduction-interface rheologies on slab stagnation and trench rollback. Earth Planet. Sci. Lett. 379, 95-103.

Čížková, H., Bina, C., 2015. Geodynamics of trench advance: insights from aphilippine-sea-style geometry. Earth Planet. Sci. Lett. 430, 408-415.

Cosca, M.A., Arculus, R.J., Pearce, J.A., Mitchell, J.G., 1998. 40Ar/39Ar and K-Ar geochronological age constraints for the inception and early evolution of the Izu-Bonin-Mariana arc system. Island Arc 7, 579-595. http://dx.doi.org/10.1111/j.1440-1738.1998.00211.x.

DeMets, C., Gordon, R.G., Argus, D.F., 2010. Geologically current plate motions. Geophys. J. Int. 181, 1-80.

Deschamps, A., Lallemand, S., 2002. The West Philippine Basin: an Eocene to early Oligocene back arc basin opened between two opposed subduction zones. J. Geophys. Res. 107 (B12), 2322. http://dx.doi.org/10.1029/2001JB001706.

Deschamps, A.E., Monié, P., Lallemand, S.E., Hsu, S.-K., Yeh, K.Y., 2000. Evidence for Early Cretaceous oceanic crust trapped in the Philippine Sea Plate. Earth Planet. Sci. Lett. 179, 503-516.

Dewey, J.F., 1980. Episodicity, sequence and style at convergent plate boundaries. In: Strangway, D.W. (Ed.), The Continental Crust and Its Mineral Deposits. Geological Association of Canada, Special Paper 20pp. 553-573 pp. 804.

Di Giuseppe, E., van Hunen, J., Funiciello, F., Faccenna, C., Giardini, D., 2008. Slab stiffness control of trench motion: insights from numerical models. Geochem. Geophys. Geosyst. 9. http://dx.doi.org/10.1029/2007GC001776.

Di Giuseppe, E., Faccenna, C., Funiciello, F., van Hunen, J., Giardini, D., 2009. On the relation between trench migration, seafloor age, and the strength of the subducting lithosphere. Lithosphere 1 (2), 121-128.

Enns, A., Becker, T.W., Schmeling, H., 2005. The dynamics of subduction and trench migration for viscosity stratification. Geophys. J. Int. 160, 761-775.

Faccenna, C., Heuret, A., Funiciello, F., Lallemand, S., Becker, W.T., 2007. Predicting trench and platemotion from the dynamics of a strong slab. Earth Planet. Sci. Lett. 257, 29-36.

Faccenna, C., Di Giuseppe, E., Funiciello, F., Lallemand, S., van Hunen, J., 2009. Control of seafloor aging on the migration of the Izu-Bonin-Mariana trench. Earth Planet. Sci. Lett. 288 (3-4), 386-398. http://dx.doi.org/10.1016/j.epsl.2009.09.042.

Faccenna, C., Becker, T.W., Lallemand, S., Steinberger, B., 2012. On the role of slab pull in the cenozoic motion of the Pacific. Geophys. Res. Lett. 39, L03305. http://dx.doi.org/10.1029/2011GL050155.

Fryer, P., 1995. Geology of the Mariana Trough. In: Taylor, B. (Ed.), Backarc Basins: Tectonics and Magmatism. Plenum Press, New York, NY, pp. 237-279.

Fukao, Y., Obayashi, M., 2013. Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity. J. Geophys. Res. 118, 5920-5938. http://dx.doi.org/10.1002/2013JB010466.

Fukao, Y., Widiyantoro, S., Obayashi, M., 2001. Stagnant slabs in the upper and lower mantle transition region. Rev. Geophys. 39, 291-323.

Funiciello, F., Faccenna, C., Giardini, D., Regenauer-Lieb, K., 2003. Dynamics of retraining slabs (part 2): insights from 3D laboratory experiments. J. Geophys. Res. 108 (B4). http://dx.doi.org/10.1029/2001JB000896.

Funiciello, F., Faccenna, C., Heuret, A., Lallemand, S., Di Giuseppe, E., Becker, T.W., 2008. Trench migration, net rotation and slab-mantle coupling. Earth Planet. Sci. Lett. 271, 233-240.

Gripp, A.E., Gordon, R.G., 2002. Young tracks of hotspots and current plate velocities. Geophys. J. Int. 150, 321-361.

Gudmundsson, O., Sambridge, M., 1998. A regionalized upper mantle (RUM) seismic model. J. Geophys. Res. 103, 7121-7136.

Hager, B.H., 1984. Subducted slabs and the geoid: constraints on mantle rheology and flow. J. Geophys. Res. 89 (B7), 6003-6015.

Hall, R., 2002. Cenozoic geological and plate tectonic evolution of SE Asia and the SW Pacific: computer-based reconstructions, model and animations. J. Asian Earth Sci. 20, 353-431.

Hall, R., Ali, J.R., Anderson, C.D., Baker, S.J., 1995. Origin and motion history of the Philippine Sea Plate. Tectonophysics 251, 229-250.

Heuret, A., Lallemand, S., 2005. Plate motions, slab dynamics and back-arc deformation. Phys. Earth Planet. Inter. 149, 31-51.

Heuret, A., Funiciello, F., Faccenna, C., Lallemand, S., 2007. Plate kinematics, slab shape and back-arc stress: a comparison between laboratory models and current subduction zones. Earth Planet. Sci. Lett. 256, 473-483. http://dx.doi.org/10.1016/j.epsl.2007.02.004.

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 Planet. Sci. Lett. 120, 375-407.

van der Hilst, R.D., Engdahl, E.R., Spakman, W., Nolet, G., 1991. Tomographic imaging of subduction lithosphere below northwest Pacific island arcs. Nature 353, 37-42.

Holt, A.F., Becker, T.W., 2017. The effect of a power-law mantle viscosity on trench retreat rate. Geophys. J. Int. 208, 491-507.

Holt, A.F., Buffett, B.A., Becker, T.W., 2015. Overriding plate thickness control on subducting plate curvature. Geophys. Res. Lett. 42, 3802-3810.

Holt, A.F., Royden, L.H., Becker, T.W., 2017. The dynamics of double slab subduction. Geophys. J. Int. http://dx.doi.org/10.1093/gji/ggw496.

Huang, C.Y., Wu, W.Y., Chang, C.P., Tsao, S., Yuan, P.B., Lin, C.W., Xia, K.Y., 1997.

Tectonic evolution of accretionary prism in the arc-continent collision terrane of Taiwan. Tectonophysics 281 (1-2), 31-51.

Huang, Z., Zhao, D., Hasegawa, A., Umino, N., Park, J.H., Kang, I.B., 2013. Aseismic deep subduction of the Philippine Sea Plate and slab window. J. Asian Earth Sci. 75, 83-94.

Huchon, P., Barrier, E., de Bremaecker, J.C., Angelier, J., 1986. Collision and stress trajectories in Taiwan: a finite element model. Tectonophysics 125 (1-3), 179-191.

Hussong, D.M., Uyeda, S., 1981. Tectonic processes and the history of the Mariana Arc: a synthesis of the results of Deep Sea Drilling Project Leg 60. Initial Rep. Deep Sea Drill. Proj. 60, 909-929.

Ishizuka, O., Yuasa, M., Sakamoto, I., Kanayama, K., Taylor, R.N., Umino, S., Tani, K., Ohara, Y., 2008. Earliest Izu-Bonin arc volcanism found on the submarine Bonin Ridge. Eos. Trans. AGU 89 (53) (Fall Meet. Suppl., Abstract V31A-2106).

Ishizuka, O., Tani, K., Reagan, M.K., Kanayama, K., Umino, S., Harigane, Y., Sakamoto, I., Miyajima, Y., Yuasa, M., Dunkley, D.J., 2011. The timescales of subduction initiation and subsequent evolution of an oceanic island arc. Earth Planet. Sci. Lett. 306, 229-240. http://dx.doi.org/10.1016/j.epsl.2011.04.006.

Jagoutz, O., Royden, L.H., Holt, A.F., Becker, T.W., 2015. Anomalously fast convergence of India and Eurasia caused by double subduction. Nat. Geosci. 8, 475-478.

Jolivet, L., Huchon, P., Rangin, C., 1989. Tectonic setting of Western Pacific marginal basins. Tectonophysics 160, 23-47.

Kato, T., Beavan, J., Matsushima, T., Kotake, Y.J., Camacho, T., Nakao, S., 2003. Geodetic evidence of back-arc spreading in the Mariana Trough. Geophys. Res. Lett. 30 (12), 1625. http://dx.doi.org/10.1029/2002GL016757.

Kimura, G., 1996. Collision orogeny at arc-arc junctions in the Japanese Islands. Island Arc 5 (3), 262-275.

Kimura, J.-I., Stern, R.J., Yoshida, T., 2005. Reinitiation of subduction and magmatic responses in SW Japan during Neogene time. Geol. Soc. Am. Bull. 117 (7-8), 969-986.

Kimura, G., Hashimoto, Y., Kitamura, Y., Yamaguchi, A., Koge, H., 2014. Middle Miocene swift migration of the TTT triple junction and rapid crustal growth in southwest Japan: a review. Tectonics 33, 1219-1238. http://dx.doi.org/10.1002/2014TC003531.

Kincaid, C., Olson, P., 1987. An experimental study of subduction and slab migration. J. Geophys. Res. 92, 13832-13840.

Kiraly, A., Capitanio, F., Funiciello, F., Faccenna, C., 2016. Subduction zone interaction: controls on arcuate belts. Geology 44 (8). http://dx.doi.org/10.1130/G37912.1.

Kizaki, K., 1986. Geology and tectonics of the Ryukyu Islands. Tectonophysics 125 (1), 193-207.

Lallemand, S., 2016. Philippine Sea Plate inception, evolution, and consumption with special emphasis on the early stages of Izu-Bonin-Mariana subduction. Prog. Earth Planet. Sci. 3 (15). http://dx.doi.org/10.1186/s40645-016-0085-6.

Lallemand, S., Popoff, M., Cadet, J.-P., Deffontaines, B., Bader, A.-G., Pubellier, M., Rangin, C., 1998. The junction between the central and the southern Philippine
Trench. J. Geophys. Res. 103 (B1), 933-950.

Lallemand, S., Font, Y., Bijwaard, H., Kao, H., 2001. New insights on 3-D plates inter-
action near Taiwan from tomography and tectonic implications. Tectonophysics 335
(3-4), 229-253.

Lallemand, S., Heuret, A., Faccenna, C., Funiciello, F., 2008. Subduction dynamics as
revealed by trench migration. Tectonics 27, Te3014. http://dx.doi.org/10.1029/
2007tc002212.

Le Pichon, X., Huchon, P., 1987. Central Japan triple junction revisited. Tectonics 6,
35-45.

Letouzey, J., Kimura, M., 1985. Okinawa Trough genesis: structure and evolution of a
backarc basin developed in a continent. Mar. Pet. Geol. 2 (2), 111-130. http://dx.doi.
org/10.1016/0264-8172(85)90002-9.

Li, C., van der Hilst, R.D., 2010. Structure of the upper mantle and transition zone be-
neath southeast Asia from traveltime tomography. J. Geophys. Res. 115. http://dx.
doi.org/10.1029/2009JB006882.

Li, C., Van der Hilst, R.D., Engdahl, E.R., Burdick, S., 2008. A new global model for P-
wavespeed variations in Earth's mantle. Geochem. Geophys. Geosyst. 9, Q05018.
http://dx.doi.org/10.1029/2007GC001806.

Lu, C.Y., Hsu, K.J., 1992. Tectonic evolution of the Taiwan mountain belt. Petrol. Geol. Taiwan 27, 21-46.

Malavieille, J., Lallemand, S.E., Dominguez, S., Deschamps, A., Lu, C.-Y., Liu, C.-S., Schnürle, P., the ACT scientific crew, 2002. Arc-continent collision in Taiwan: new marine observations and tectonic evolution. Geol. Soc. Am. Spec. Pap. 358, 187-211.

Martinez, F., Fryer, P., Baker, N.A., Yamazaki, T., 1995. Evolution of back-arc rifting: Mariana Trough, 20-24N. J. Geophys. Res. 100, 3807-3827.

Miller, M.S., Kennett, B.L.N., Lister, G.S., 2004. Imaging changes in morphology, geo-
metry, and physical properties of the subducting Pacific plate along the Izu-Bonin-
Mariana arc. Earth Planet. Sci. Lett. 224 (3-4), 363-370.

Miller, M.S., Kennett, B.L.N., Toy, V.G., 2006. Spatial and temporal evolution of the subducting Pacific plate structure along the western Pacific margin. J. Geophys. Res. 111, B02401. http://dx.doi.org/10.1029/2005JB003705.

Mitrovica, J.X., Forte, A.M., 2004. A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data. Earth Planet. Sci. Lett. 225, 177-189.

Moresi, L.N., Gurnis, M., 1996. Constraints on the lateral strength of slabs from three-
dimensional dynamic flow models. Earth Planet. Sci. Lett. 138, 15-28.

Müller, R.D, Sdrolias, M., Gaina, C., Roest, W.R., 2008. Age, spreading rates and
spreading symmetry of the world's ocean crust. Geochem. Geophys. Geosyst. 9http://
dx.doi.org/10.1029/2007GC001743. Q04006.

Okino, K., Shimakawa, Y., Nagaoka, S., 1994. Evolution of the Shikoku basin. J. Geomagn. Geoelectr. 46, 463-479.

Okino, K., Kasuga, S., Ohara, Y., 1998. A new scenario of the Parece Vela basin genesis.
Mar. Geophys. Res. 20 (1), 21-40.

Park, J.O., 1996. Seismic Stratigraphy and Tectonic Evolution of the Southern Ryukyu Island Arc. Ocean Res. Inst. Univ. of Tokyo (Ph. D. thesis).

Schellart, W.P., 2004. Kinematics of subduction and subduction-induced flow in the upper mantle. J. Geophys. Res. 109, B07401. http://dx.doi.org/10.1029/2004JB002970.

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

Schellart, W.P., Stegman, D.R., Freeman, J., 2008. Global trench migration velocities and slab migration induced upper mantle volume fluxes: constraints to find an Earth reference frame based on minimizing viscous dissipation. Earth-Sci. Rev. 88,

118-144.

Sdrolias, M., Müller, R.D., 2006. Controls on back-arc basin formation. Geochem. Geophys. Geosyst. 7, Q04016. http://dx.doi.org/10.1029/2005GC001090.

Sdrolias, M., Roest, W.R., Müller, R.D., 2004. An expression of Philippine Sea plate rotation: the Parece Vela and Shikoku Basins. Tectonophysics 394, 69-86. http://dx.doi.org/10.1016/j.tecto.2004.07.061.

Seno, T., Maruyama, S., 1984. Paleogeographic reconstruction and origin of the Philippine Sea. Tectonophysics 102, 53-84.

Seton, M., Muller, R.D., Zahirovic, S., Gaina, C., Torsvik, T.H., Shephard, G., Talsma, A., Gurnis, M., Turner, M., Maus, S., Chandler, M., 2012. Global continental and ocean basin reconstructions since 200 Ma. Earth Sci. Rev. 113 (34), 212-270. http://dx.doi.org/10.1016/j.earscirev.2012.03.002. (ISSN 0012-8252).

Shimizu, S., Itaya, T., 1993. Plio-Pleistocene arc magmatism controlled by two over-
lapping subducted plates, central Japan. Tectonophysics 225 (3), 139-154.

Shinjo, R., 1999. Geochemistry of high Mg andesites and the tectonic evolution of the
Okinawa Trough Ryukyu arc system. Chem. Geol. 157 (1), 69-88.

Sibuet, J.C., Deffontaines, B., Hsu, S.K., Thareau, N., Le Formal, J.P., Liu, C.S., Party, A.C.T., 1998. Okinawa trough backarc basin: early tectonic and magmatic evolution. J. Geophys. Res. 103 (B12), 30245-30267.

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. Geochem. Geophys. Geosyst. 7, Q03012. http://dx.doi.org/10.1029/2005GC001056.

Stern, R.J., Fouch, M.J., Klemperer, S.L., 2003. An overview of the Izu-Bonin-Mariana subduction factory. In: Eiler, J. (Ed.), Inside the Subduction Factory. Geophys Monogr 138 American Geophysical Union, Washington, D.C., pp. 175-222.

Tamaki, K., Suyehiro, K., Allan, J., Ingle, J.C., Pisciotto, K., 1992. Tectonic synthesis and implications of Japan Sea ODP drilling. Proc. Ocean Drill. Program Sci. Results 1 (27-128), 1333-1350.

Taylor, B., Karner, G.D., 1983. On the evolution of marginal basins. Rev. Geophys. Space Phys. 21, 1727-1741.

Wang, Z., Huang, R., Huang, J., He, Z., 2008. P-wave velocity and gradient images be-
neath the Okinawa Trough. Tectonophysics 455, 1-13. http://dx.doi.org/10.1016/j.
tecto. 2008.03.004.

Wei, W., Xu, J., Zhao, D., Shi, Y., 2012. East Asia mantle tomography: new insight into
plate subduction and intraplate volcanism. J. Asian Earth Sci. 60, 88-103.

Wei, W., Zhao, D., Xu, J., Wei, F., Liu, G., 2015. P and S wave tomography and anisotropy in Northwest Pacific and East Asia: constraints on stagnant slab and intraplate volcanism. J. Geophys. Res. 120, 1642-1666. http://dx.doi.org/10.1002/
2014JB011254.

Widiyoranto, S., Kennett, B.L.N., van Der Hilst, R.D., 1999. Seismic tomography with P and S data reveals lateral variations in the rigidity of deep slabs. Earth Planet. Sci. Lett. 173, 91-100.

Wu, J., Suppe, J., Lu, R., Kanda, R., 2016. Philippine Sea and East Asian plate tectonics since 52 Ma constrained by new subducted slab reconstruction methods. J. Geophys. Res. 121, 4670-4741. http://dx.doi.org/10.1002/2016JB012923.

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

Zhong, S., 2006. Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature and upper mantle temperature. J. Geophys. Res. 111. http://dx.doi.org/10.1029/2005JB003972.