Terrane Collision‐Induced Subduction Initiation: Mode Selection and Implications for Western Pacific Subduction System

In the regime of plate tectonics, the subduction of an oceanic plate generally terminates with the collision and accretion of continental terranes. Then, a new subduction zone may form in the neighboring oceanic plates, which is defined as the terrane collision‐induced subduction initiation (SI). Based on the analyses of the western Pacific subduction system in the Cenozoic, three types of collision‐induced SI have been observed: subduction polarity reversal, subduction transference and far‐field subduction. However, the dynamics and controlling factors of SI mode selection after terrane collision are not clear. In this study, a multi‐terrane collision model has been conducted with variable rheological strength of continental terranes and different convergence velocities. The model results indicate that the relative strength of the terranes controls the SI mode selection, with the new subduction zone tending to form beneath weaker terranes. In addition, the higher convergence velocity can facilitate the collision‐induced SI. An analytical study of force balance has been further conducted, which provides a mechanical explanation for the numerical prediction of a weak overriding terrane as a favorable SI site. The numerical models and force balance analyses are further compared with the natural cases in the western Pacific subduction system. This indicates that subduction polarity reversal is the most favorable mode after terrane collision in the western Pacific, possibly due to the weakness of overriding plate during the preceding subduction‐induced fluid/melt activity. This comprehensive study provides systematic constraints for the dynamics of collision‐induced subduction jump, especially for the western Pacific subduction zones in the Cenozoic.

There are two classical types of SI, that is, spontaneous or induced, which are distinguished by the sources of driving forces (Stern & Gerya, 2018).The spontaneous SI refers to the collapse of oceanic plate along transform faults or fracture zones, mainly driven by the local density contrast between plates.Previous statistical and numerical studies indicate that spontaneous SI is not easy and requires more critical conditions (Arcay et al., 2020;Crameri et al., 2020;Lallemand & Arcay, 2021;Zhong & Li, 2021).The induced SI is driven by, or at least accompanied with, far-field forces, such as terrane collision (Stern, 2004; G. X. Yang, 2022;Zhong & Li, 2020), mantle plume (Baes et al., 2020;Gerya et al., 2015;van Hinsbergen et al., 2021), large-scale mantle flow traction and mid-ocean ridge push (Baes & Sobolev, 2017;Lu et al., 2015;Zhong & Li, 2019) or neighboring slab pull (Govers & Wortel, 2005;Zhou et al., 2020).Among these induced SI cases, terrane collision indicates the closure of an old subduction zone; the continued convergence may lead to strain localization and finally formation of a new subduction zone (Lallemand & Arcay, 2021;Liu et al., 2021;Sun et al., 2021;Yan et al., 2021;Zhong & Li, 2020).Terrane collision-induced SI occurs widely (e.g., Li et al., 2023;G. X. Yang, 2022; and references therein), for example, during the long-term evolution of the Tethyan system as well as the western Pacific subduction system.
Based on the analyses of western Pacific subduction zones (Figure 1), we find that a number of SI cases are related to the terrane collision, and several different types of SI can be resulting after the terrane collision, for example, subduction polarity reversal, subduction transference, and far-field subduction (Figure 1).The most favorable SI case after terrane collision in the western Pacific is the subduction polarity reversal, in which the SI occurs in the far-end of the overriding terrane and has an opposite subduction direction.The natural cases of subduction polarity reversal may include the Aleutians, Bowers Ridge, Proto-Kamchatka, Kuril, Philippines, Palawan, N-Sulawesi, Flores, Wetar, New Guinea-Pocklington, North Melanesian-Solomons-Vitiaz, New-Britain-San Cristobal-New Hebrides, Ralik-Gilbert-Samoa, and so on, in the western Pacific region (Lallemand & Arcay, 2021;Vaes et al., 2019).On the other hand, the subduction transference may also be resulting after terrane collision, which generally occurs behind the subducting low-density terrane and has the same polarity as the preexisting subduction zone.The natural cases of subduction transference include Kamchatka, Sulu, and Trobriand in the western Pacific region (Schellart et al., 2006;Wu et al., 2016), and this mode is considered the dominant case during Tethyan evolution (Li et al., 2023;Zhong & Li, 2020, 2022).As a rare case, terrane collision may also lead to a far-field SI in the far-end of the neighboring oceanic lithosphere in the lower plate and having an opposite subduction direction, which may be the case of Sangihe and Halmahera with forming a divergent double-sided subduction of an oceanic plate.In contrast to these three modes of SI, the accumulated force during terrane collision may not be enough to conquer the resistance and trigger the SI, which will lead to a mode with no subduction jump (Tao et al., 2020;Vogt & Gerya, 2014;S. H. Yang et al., 2018;Zhong & Li, 2020).
In the previous studies, the different modes of collision-induced subduction jump have generally been investigated separately by different types of geodynamic models.These studies indicate that the favorable conditions of subduction polarity reversal include the positively buoyant and rheologically strong subducting plate, young and weak overriding plate, weak mantle wedge, and so on (Almeida, Riel, Rosas, Duarte, & Kaus, 2022;Almeida, Riel, Rosas, Duarte, & Schellart, 2022;Baes et al., 2011;Dong et al., 2022;Riel et al., 2023;Sun et al., 2021;Wang et al., 2022; S. X. Zhang & Leng, 2021).Alternatively, the subduction transference can be resulting if the overriding plate is rheologically strong or has continental affinity, with other favorable conditions including young and moderate-sized oceanic plateau with a strongly depleted mantle root, existence of weakness at the farend passive margin of the subducting terrane, the forced convergence, and so on (Liu et al., 2021;Tao et al., 2020;Yan et al., 2021;S. H. Yang et al., 2018;Zhong & Li, 2020).Far-field SI after terrane collision has never been studied, possibly due to its rare occurrence in nature as well as its dynamic difficulty.
Based on the natural cases in the western Pacific subduction zones and the previous numerical studies, the key scientific question is that why the terrane collision can result in such different modes of SI; or what is(are) the controlling factor(s) of the SI mode selection after terrane collision?In order to answer this question, a large-scale multi-terrane numerical model is conducted, with variable rheological strength of terranes and variable convergence velocity.Furthermore, an analytical study of force balance is conducted, which provides a better understanding of the numerical model.The numerical and analytical results are finally compared with the geological records of the western Pacific subduction system, with providing a comprehensive understanding of the terrane collision-induced SI and shedding a new light on the SI dynamics and plate tectonics.

Numerical Model Configuration
The numerical simulations are conducted with the original code I2VIS (Gerya, 2010).The processes of hydration, partial melting and phase transitions at 410 and 660 km are included in the numerical models.The specific version  (Wessel et al., 2013), with the data from Tozer et al. (2019).
and detailed numerical methods are described in Li et al. (2019) and also summarized in Supporting Information S1.The two-dimensional numerical model is configured within a domain of 5,000 km in length and 1,400 km in depth, as shown in Figure 2a and Figure S1b in Supporting Information S1 with a different colorbar.The spatial resolution of the model is 5 km in the horizontal direction, while that in the vertical direction is 1 km from 0 to 150 km and gradually reduces to 10 km downward to the bottom.The model includes three continental terranes and three oceanic plates.Terrane 1 is located in the left part of the model and has a length of 942 km.Terranes 2 and 3 are both 700 km in the reference models, the effect of which has been tested.The crust of these continental terranes has the same thickness of 35 km, including an upper crust of 20 km and a lower crust of 15 km.The The snapshot represents a model stage after convergence for 10 Myr, with a total amount of convergence of 500 km.The inset-enlarged figures show the deformed OCT after convergence.The dark line within the lithospheric mantle indicates the lithosphere cooling and growth, that is, the newly formed lithospheric material below this line.The colorbar indicates the composition field of the model: 1-sticky air; 2-sea water; 3, 4-sediment; 5-continental upper crust; 6-continental lower crust; 7-oceanic upper crust; 8-oceanic lower crust; 9-oceanic lithospheric mantle; 10-asthenosphere; 11-hydrated mantle; 12-hydrated mantle with serpentinization; 13-lithospheric mantle of Terrane 1; 14-lithospheric mantle of Terrane 2; 15-lithospheric mantle of Terrane 3; 16-partially molten sediments; 17-partially molten continental crust; 18-partially molten oceanic crust; 19-partially molten mantle.
thickness of continental lithosphere is 140 km for all the three terranes.On the other hand, Oceanic Plates 1 and 2 are both set horizontally and have the same length of 1,150 km, whereas the downgoing slab is composed of two parts, that is, a horizontal segment of 358 km in length and a bending segment to the depth of 250 km with a dip angle of about 30°.The oceanic lithosphere is composed of a 3-km-thick basalt layer as the upper crust and a 5km-thick gabbro layer as the lower crust.The thickness of the oceanic lithosphere is controlled by the age, which is initially set to be 20-Myr-old in the reference models.A weak zone is set between the subducting slab and Terrane 3 to ensure the location of the first subduction and collision.A "sticky air" layer is emplaced on the surface of the model to accommodate the crustal deformation, gravity isostasy and topography (Crameri et al., 2012;Schmeling et al., 2008), which is 10-km-thick above the continental terrane and 12-km-thick above the oceanic plate.
For the temperature field configuration, the top and bottom boundaries of the model are set to be 273 and 2248 K, respectively.The thermal structure of oceanic plate is emplaced according to the half-space cooling model (Turcotte & Schubert, 2002).The continental lithosphere is controlled by surface temperature of 273 K and an initial temperature of 1623 K at the bottom of the lithosphere, as well as a linear gradient in between.The "sticky air" layer remains the constant temperature of 273 K.An initial thermal gradient of 0.5 K/km is applied for the sub-lithospheric mantle.The left and right boundaries of the model are adiabatic with no horizontal heat flux.
For the mechanical boundary condition, free-slip is satisfied for all the boundaries.Within the model domain, Terrane 1 is pushed with a constant velocity, and the velocity patch is located in a constant position relative to Terrane 1 during the model evolution, as shown by the yellow arrow in Figure 2.With a continuous convergence of 5 cm/yr for about 10 Myrs, the downgoing oceanic plate has been totally subducted and Terrane 2 starts subduction and collision with Terrane 3 (Figure 2b).This stage of terrane collision will be employed as the initial model, with the oceanic lithosphere thickened after cooling for another 10 Myrs (cf.Figures 2a and 2b), in order to further investigate the mode selection of SI/jump.

Reference Models
First, we are going to show four typical model results with one mode of "no subduction jump," as well as three different modes of collision-induced SI, that is, "subduction polarity reversal," "subduction transference," and "far-field subduction."These different modes are mainly controlled by the lithospheric mantle strength, represented by the plastic friction coefficient λ • sin (φ) of Terrane 1, Terrane 2, and Terrane 3, which varies from 0.6 to 0.001 according to previous investigations of fluid/melt weakening activity, as well as the heterogeneity of mantle rocks (Gerya et al., 2015;Li et al., 2016Li et al., , 2019)).The relationship between the composite rheological strength and plastic friction is shown in Text S1 in Supporting Information S1, while the yield stress profiles with different plastic friction coefficients are shown in Figure S1c in Supporting Information S1.In this set of models, a constant convergence velocity of 5 cm/yr is applied.
During the preceding model evolution, the oceanic subduction-induced hydration weakens the subduction channel.After the collision between Terranes 2 and 3, the external push and local slab pull drive the Terrane 2 into the channel.However, the low-density of continental terrane resists its subduction and results in slab break-off (Figure 3a).Consequently, the slab pull loses and the hard collision between Terranes 2 and 3 leads to strong resistance for plate convergence.However, because all the blocks (continental terranes and oceanic plates) are rheologically strong, the SI is not favorable.Thus, the plate convergence has to be accommodated by the underthrusting of continental Terrane 2 (Figure 3).During this process, larger and larger resistance has been built, which means that the driving force of plate convergence is increasing in order to maintain the constant kinematic convergence rate (Zhong & Li, 2019).In this model, the SI is not predicted even with the fully subduction of Terrane 2, indicating this is not a favorable case for subduction jump.In nature, the plate convergence may cease before the fully subduction of Terrane 2 and the new subduction zone will not form.
The beginning of model evolution is similar (cf.Figures 3 and 4), with slab break-off after terrane collision.During the continuous plate convergence, the onset of SI occurs in the overriding plate (Figure 4b), because the overriding Terrane 3 is rheologically weak and thus favors strain localization.Subsequently, Oceanic Plate 2 subducts continuously beneath the low-density Terrane 3 in the opposite direction to the previous subduction, which is defined as the mode of subduction polarity reversal.
During the mature stage of Oceanic Plate 2 subduction, the trench retreat results because the right-end of Oceanic Plate 2 is fastened to the rigid boundary wall.The fast trench retreat drives the back-arc extension and further relieves the collision between Terranes 2 and 3 (Figure 4).It indicates that the formation of a new subduction zone can release the previous convergent force, or even lead to extension.

Subduction Transference
In this model (Figure 5), the Terranes 1 and 3 are rheologically strong; instead, the lithospheric mantle of Terrane 2 is weaker, with λ • sin (φ) = 0.1 initially and decreasing linearly to 0.01 with the accumulated plastic strain to ε plas = 1.Other parameters are identical as the above cases.
The initial collision process is similar to that of subduction polarity reversal (cf.Figures 4 and 5).However, after the slab break-off, the passive margin between Oceanic Plate 1 and Terrane 2 is the favorable site for strain localization (Figure 5b).Consequently, a new subduction zone is formed in Oceanic Plate 1, with the same polarity as the previous subduction, which is defined as the mode of subduction transference (Figure 5c).Although Terrane 1 is still pushed rightward, the slab pull from the subducted Oceanic Plate 1 dominates and drives the back-arc extension and the collision termination between Terranes 2 and 3 (Figure 5d).

Far-Field Subduction
In this model (Figure 6), the Terranes 2 and 3 are rheologically strong; the lithospheric mantle of Terrane 1 is weaker, with λ • sin (φ) = 0.1 initially and decreasing linearly to 0.01 with the accumulated plastic strain to ε plas = 1.Other parameters are identical as the above cases.
After the hard collision between the strong Terranes 2 and 3, the downgoing Terrane 2 is forced to bend, with the collision zone thickening, surface uplift and orogeny.Then, the continuous convergence leads to strain localization in the passive margin between the Oceanic Plate 1 and the relatively weak Terrane 1 (Figure 6b).A mode of far-field subduction is resulting (Figures 6c and 6d).It is interesting that the subducting Oceanic Plate 1 drives the far-field stretching of the present lower plate and also decouples the Terranes 2 and 3 in the previous collision zone (Figure 6d), the dynamics of which has been systematically studied by S. Yang et al. (2021).
In this case, the rheological strength of Terrane 1 cannot be too low; otherwise, the first subduction between Terranes 2 and 3 will lead to the far-field stretching of the weak Terrane 1.Thus, the strain localization and breakup between Terrane 1 and Oceanic Plate 1 will be resulting before the downgoing slab break-off.Consequently, the far-field subduction may be induced by stretching, rather than terrane collision as expected.Therefore, a medium mantle strength is applied for Terrane 1 (Figure 6) rather than the weakest one.

Mode Selection With Variable Rheological Strength of Terranes
The above reference models show that the lithospheric strength of continental terranes can strongly influence the modes of collision-induced SI.In this section, we try to build a phase diagram to clearly demonstrate the rules for such mode selection (Figure 7).In the 3-D phase diagram, the three axes stand for the lithospheric strength of the three continental terranes, respectively, which are further represented by the plastic friction coefficient (λ • sin (φ)) of the lithospheric mantle.The two values for each point (e.g., 0.6, 0.1) indicate the initial friction coefficient and that after the plastic strain weakening, respectively.For example, λ • sin(φ) = (0.6, 0.1) means the plastic friction coefficient of lithospheric mantle is 0.6 initially and decreases linearly to 0.1 with the accumulated plastic strain to ε plas = 1, while λ • sin(φ) = (0.1, 0.1) means the plastic friction coefficient is 0.1 initially and remains constant with the variation of ε plas .All the other conditions are identical, including the constant convergence velocity of 5 cm/yr.
In the diagram A (Figure 7), the lithospheric strength of Terrane 3 is kept constant with a high value; thus, the mode section is dependent on the relative strength of Terranes 1 and 2. If Terrane 1 is weaker than Terrane 2, the mode of far-field subduction beneath Terrane 1 is preferred, as marked in red (Figure 7, A).In contrast, when Terrane 2 is weaker than Terrane 1, the mode of subduction transference is favored, as marked in orange (Figure 7, A).As a special case, when the lithospheres of Terranes 1, 2, and 3 are identically strong, no subduction jump is predicted as marked in gray (Figure 7, A), the evolution of which is shown in Figure 3.
Diagram B in Figure 7 describes the mode selection when the lithospheric strength of Terrane 3 is medium.Similar to diagram A, the mode selection between far-field subduction and subduction transference is still controlled by the relative strength of Terranes 1 and 2. However, when the Terranes 1 and 2 are both rheologically stronger than Terrane 3 with an intermediate strength, subduction polarity reversal is favored in which the weak Terrane 3 is the overriding plate, as marked in blue in the diagram B.
In the diagram C (Figure 7), the weakest Terrane 3 dominates the mode selection of collision-induced SI.The mode of subduction polarity reversal is the most favorable case, which occurs when Terranes 1 and 2 are not as weak as Terrane 3. The mode of far-field subduction is predicted when Terrane 1 is rheologically weak with λ • sin (φ) ≤ 0.01.Alternatively, the mode of subduction transference is predicted when Terrane 2 is rather weak.
From diagram A to C in Figure 7, the lithospheric strength of Terrane 3 gradually decreases and the mode of subduction polarity reversal gradually dominates.In nature, the Terrane 3, as the overriding plate of the preceding subduction zone, is more likely to be rheologically weak due to the hydration-induced lithospheric weakening  3; Circle for the model with subduction polarity reversal in Figure 4; Square for the model with subduction transference in Figure 5; Triangle for the model with far-field subduction in Figure 6.(Li, 2020(Li, , 2022)).Thus, the mode of subduction polarity reversal is generally favored.The systematic model results indicate that the collision-induced SI prefers to occur with a relatively weak overriding plate.

Effect of Convergence Velocity
The constant convergence velocity of 5 cm/yr has been applied in all the aforementioned models.In order to better understand the effect of convergence velocity on the mode selection of terrane collision-induced SI, comparable models are conducted with either a lower convergence velocity of 2 cm/yr or a higher value of 10 cm/yr.
The phase diagram with a low convergence velocity of 2 cm/yr is exhibited in Figure 8a, which demonstrates a similar trend of mode selection as in the diagram with 5 cm/yr (Figure 7).This indicates that the relative lithospheric strength of continental terranes dominates the mode selection and the SI prefers to occur with a weak overriding plate.The mode with no subduction jump is also predicted when all the three terranes are relatively strong; in addition, this mode occupies a wider parameter range, indicating that slower convergence may resist the collision-induced SI.On the one hand, a lower convergence velocity generally represents a lower boundary force, which thus resists the strain location and SI formation (Faccenna et al., 1999;Gurnis et al., 2004;Zhong & Li, 2019).On the other hand, a lower convergence velocity means that the oceanic plate cools and thickens faster than that with a higher convergence velocity; thus, the resulting thicker and stronger oceanic lithosphere also resists SI (McKenzie, 1977;Mueller & Phillips, 1991).
The phase diagram with a high convergence velocity of 10 cm/yr is shown in Figure 8b.The results demonstrate similarly that the collision-induced SI favors occurrence with a weak overriding terrane.In this regime, the mode with no subduction jump has never been predicted, even with all the strong continental terranes.This confirms that the collision-induced SI is more likely to be induced with a higher convergence velocity or larger force.
In the regime with the high convergence velocity of 10 cm/yr, multiple modes may be selected simultaneously or one after another, depending on the competition of strain localization in different domains.For example, when the lithospheric strength of Terranes 2 and 3 are both medium, but Terrane 1 is strong, the modes of subduction transference and polarity reversal occur in the same model, as marked by the yellow ellipsoid symbol in Figure 8b.This model evolution is further shown in Figures 9a-9d, which demonstrates that the subduction polarity reversal occurs shortly after the collision between Terranes 2 and 3 (Figures 9a and 9b).Then, the subduction transference occurs due to the fast convergence that cannot be accommodated by the subduction of Oceanic Plate 2 (Figures 9c  and 9d).This model evolution is dependent on the competition of subduction choices beneath Terranes 2 and 3 with similar rheological strength.Another case of subduction competition is marked with a yellow diamond in Figure 8b, and its detailed evolution is shown in Figures 9a'-9d'.After the terrane collision and slab break-off, the inception of subduction transference occurs (Figure 9b'), which is replaced later by subduction polarity reversal (Figure 9c').Finally, the subduction transference ceases and the mode of subduction polarity reversal dominates (Figure 9d').
Such kind of complex and multiple mode selection is related to the high convergence velocity as well as the similarly weak terranes.When an oceanic plate attempts to subduct at a favorable site, it requires a certain amount of time for the plate bending, collapse and formation of a new subduction channel.During this period and under a fast convergence, the strain localization may be achieved in another site.Thus, multiple subduction zones may be predicted and they will dynamically compete with each other (Figure 9).In contrast, if the convergence velocity is not high enough, the force accumulated during the formation of the first SI is insufficient for a new subduction jump (Figures 7 and 8a).This explains why the multiple mode competition is predicted only in the regime with a high convergence velocity (Figure 8b).

Sensitivity of Other Parameters
Several additional parameters may affect the terrane collision-induced SI mode selection, for example, the initial age of oceanic plate, the density of the terrane lithospheric mantle, the width of terranes, the shape of oceancontinent transition zone and the absence of external push.The parameters and results of all these sensitivity test models are listed in Table S1.
The age of oceanic plate controls its thermal structure and thus the rheological resistance to SI.In the reference model, the initial oceanic ages are set to be 20 Ma, with a lithospheric thickness of about 50 km, which thickens rapidly during the initial model evolution before terrane collision.As a sensitivity test, additional models are conducted with an older initial oceanic plate of 70 Ma, as shown in Figure S6 in Supporting Information S1.All the other parameters are identical to those in Figure 7.The SI mode distribution in Figure S7 in Supporting Information S1 with older oceanic plates is generally consistent with that in Figure 7 with younger ones, which confirms the conclusion that the collision-induced SI prefers to occur with a relatively weak overriding plate.It is worth noting that the terrane subduction is generally deeper in the models with older oceanic plates (cf. Figure S8b in Supporting Information S1 and Figure 4d), because larger driving force is required for older and thicker oceanic plates to yield and initiate a new subduction zone (Zhong & Li, 2019, 2020).
The reference density of the terrane lithospheric mantle in the reference model is set to be 3,200 kg/m 3 , assuming a refractory mantle composition with a depletion of about 3%.In order for comparison, a series of additional models have been conducted with fertile terrane lithospheric mantle (i.e., reference density = 3,300 kg/m 3 ).The model results are summarized in Figure S9 in Supporting Information S1, which exhibits a generally identical mode selection pattern as in Figure 7.This indicates that the depletion degree and related density variation of the terrane lithospheric mantle do not have a critical impact on the mode selection of collision-induced SI.
The width of collisional terranes may affect the SI mode selection.In order to test the possible influences, additional models are conducted with narrower Terranes 2 and 3, the widths of which are reduced from 700 km (Figure 7) in the reference models to 300 km.The model configuration is shown in Figure S11 in Supporting Information S1, and the phase diagram of terrane collision-induced SI is shown in Figure S12 in Supporting Information S1.The mode with no subduction jump occupies a relatively wider parameter range in the new regime with narrower terranes (Figure S12 in Supporting Information S1) than that with wider terranes (Figure 7).However, the general trend of collision-induced SI mode selection is consistent (cf. Figure S12 in Supporting Information S1 and Figure 7).
The ocean-continental transition (OCT) zone is the key region for strain localization and SI occurrence; thus, the geometry of OCTs may play a certain role in SI mode selection.In order to test this effect, a set of additional models are conducted with a new geometry of OCTs (Figure S14 in Supporting Information S1).The model results are summarized in Figure S15 in Supporting Information S1.The comparison between the model results in Figure 7 and Figure S15 in Supporting Information S1 indicates that the geometry of OCTs does not significantly affect the mode selection of collision-induced SI.
In all the previous models, the collision-induced SI is driven by both the local collision and external push.In addition, a series of models without the external push after collision have been further conducted.In these models, the convergence velocity is applied only for 10 Myr until the beginning of collision between Terranes 2 and 3. Subsequently, the models evolve freely under the slab pull of the downgoing plate.Other parameters are the same as the models with older initial oceanic plates of 70 Ma (Figure S7 in Supporting Information S1).The results indicate that collision-induced SI rarely occurs without an external push unless the overriding or far-field terrane is extremely weak.However, such SI cases are more like spontaneous SI under certain disturbances, rather than collision-induced SI.Meanwhile, subduction transference cannot be induced without the external push, as shown in Figure S17 in Supporting Information S1.

Force Balance Analysis
The numerical models reveal that the oceanic plate prefers subducting beneath the relatively weak overriding terranes, which means that the relative rheological strength of Terranes 1, 2, and 3 plays an important role in the mode selection of SI.In order to better understand the dynamics of terrane collision-induced subduction jump, a corresponding force balance analysis is conducted (Figure 10).
The plate convergence is driven by the prescribed kinematic velocity and the subducting slab pull.However, the slab break-off generally occurs after the terrane collision, which means the slab pull is absent in the following process of SI.In this case, we consider a steady state when the convergence force is balanced by the positive buoyancy of the subducted continental Terrane 2 as well as some other minor forces, for example, the possible mantle flow traction induced by the detached sinking slab.In this situation, the convergence force in the left Oceanic Plate 1 (F c1 ) should be comparable to that in the right Oceanic Plate 2 (F c2 ), as shown in Figure 10. (1) SI usually occurs at the ocean-continent transition (OCT) zone in such a tectonic regime.Consequently, the OCT zone will be taken as the representative area for analysis in each system.First, the Oceanic Plate 2 domain will be discussed.The effective viscosity of the OCT can be approximated from the neighboring Terrane 3 and Oceanic Plate 2, as η 3 = 2 • η t3 • η oc2 / (η t3 + η oc2 ).In addition, the mean thickness of the OCT can be calculated as h 3 = (h t3 + h oc2 )/ 2. The strain rate at the OCT in the overriding plate is Oceanic Plate 1 has two OCTs.In a similar way as the overriding plate domain (Equation 2), the strain rates at the left and right OCTs of Oceanic Plate 1 are as follows, respectively: Figure 10.Schematic illustration of force balance analysis.η t1 , η t2 , and η t3 are the viscosity of Terranes 1, 2, and 3, respectively, with h t1 , h t2 , and h t3 as the corresponding plate thickness.η oc1 and η oc2 are the viscosity of Oceanic Plates 1 and 2, respectively, with h oc1 and h oc2 as the corresponding plate thickness.η 1 , η 2 , and η 3 denote the viscosity of the three ocean-continental transitions, respectively.F c1 is the convergence force in Oceanic Plate 1, while F c2 is a similar one but in Oceanic Plate 2.
In the initial numerical model, the thicknesses of Terranes 1, 2, and 3 are the same, while the Oceanic Plates 1 and 2 have the same age (20 Ma) and thus the same thickness due to the half-space cooling model.The effective viscosities of Oceanic Plates 1 and 2 are also the same in the original models.Therefore, h t1 = h t2 = h t3 , h oc1 = h oc2 , and η oc1 = η oc2 .In order to compare the effects of the rheological strength of Terranes 1, 2, and 3, Equations 2-4 can be reformulated to: (5) For Equations 5-7, all the other parameters are identical except the effective viscosities of the continental terranes (η t1 , η t2 , η t3 ) which thus determine the strain rate of the OCTs, for example, a weaker terrane leading to a higher strain rate at the corresponding OCT.Due to the strain (rate) weakening, the lithospheric strength with a higher strain rate will decrease intensively, which may promote the strain localization and SI.Thus, the mode of subduction polarity reversal is favored with the weakening of Terrane 3, whereas the mode of subduction transference is favored with the weakening of Terrane 2. The mode of far-field subduction follows the same rule.If all the three terranes are similarly strong, the strain rates are rather low at the three OCTs, which explains the mode with no subduction jump (Figure 3).In addition, higher convergence velocity with similarly low strength of all terranes will equally increase the strain rate of all OCTs, indicating that there is more than one proper site for SI.Then, the competition of SI mode selection is resulting, similar to the cases shown in Figure 9.
It is worth noting that subduction nearly always initiates in the far-field when all the three terranes have the same strength.In an ideal state, the tectonic plates are considered rigid, so there is no loss of convergence force during the transmission process.However, the plates in this model undergo ductile deformation due to compression, resulting that the father from the convergence velocity patch (shown as the yellow arrow in Figure 2a), the lower convergence force is presented.In other words, the convergence force on Terrane 3 is a bit lower than that on Terrane 2, and similarly, the force on Terrane 2 is a bit lower than that on Terrane 1.That explains the dominant case of far-field SI when the mantle strengths of Terranes 1-3 are the same.However, this point is ignored in the force balance analysis, because the force loss is very low due to the strain in the plate interior is rather low.
The force analysis confirms and mechanically explains the numerical results of mode selection of terrane collision-induced SI, which prefers to occur beneath a weaker terrane.In addition, a higher convergence velocity will promote the SI due to the corresponding higher convergence force.Subduction mode competition will be resulting in the regime with similarly low terrane strength and high convergence velocity.

Comparisons With Previous Models
The terrane collision-induced subduction jump has been investigated by both laboratory and numerical modeling; however, the previous studies generally focus on one specific mode, or to say that, the different modes of subduction polarity reversal, transference, and far-field subduction are generally studied separately.
The dynamics of subduction polarity reversal has been widely investigated.The early laboratory models of Chemenda et al. (1997) reveal that the overriding plate may fail and produce a new subduction zone with a reversed polarity when the distance between the pre-existing trench and island arc exceeds a critical value.Further on, the laboratory models of Regard et al. (2008) demonstrate that the slab break-off occurs after the entrance of terranes with positive buoyancy at the trench, which may cause a slowdown of pre-existing subduction and trigger a subduction polarity reversal.Both of the above laboratory modeling studies indicate that the strong pushing on the overriding plate contributes to the formation of a new subduction zone with opposite polarity.Recently, a number of numerical modeling studies have been conducted to better understand the controlling factors of collision-induced subduction polarity reversal.S. X. Zhang and Leng (2021) confirms the possibility of both spontaneous and induced subduction polarity reversal by modifying the plastic strength of overriding plate, which highlights the importance of weak overriding lithosphere.The effect of oceanic plateau collision with continental terranes has also been studied, which may lead to the break-off of the subducted slab and result in subduction polarity reversal (Sun et al., 2021).Further on, Wang et al. (2022) and Dong et al. (2022) combines the effect of arc rheological strength and oceanic plateau collision, proposing that the overriding arc strength plays a more important role in the subduction polarity reversal and the effect of plateau collision is secondary.Similarly, Almeida, Riel, Rosas, Duarte, and Kaus (2022) and Almeida, Riel, Rosas, Duarte, and Schellart (2022) have conduct 3-D and 2-D numerical models, respectively, suggesting that the favorable conditions of subduction polarity reversal include the older (stronger and thicker) downgoing plates and relatively younger (weaker, thinner, and narrower) overriding plates.The above-mentioned laboratory and numerical modeling studies indicate that the relatively weak overriding plate plays an important role in the general formation of subduction polarity reversal.A consistent rule has been obtained in the current numerical and analytical studies, which confirm the critical role of overriding terrane weakening on the subduction polarity reversal (Figures 7 and 8).
Besides the mode of subduction polarity reversal, the terrane collision may also induce other modes of SI, for example, subduction transference and far-field subduction; however, these modes have not been considered in the above modeling studies.
For the terrane collision-induced subduction transference, this mode has been less studied comparing to the polarity reversal.Tetreault and Buiter (2012) simulate terrane accretion during collision, indicating that the accretion of crustal units occurs when the crust of future allochthonous terrane is relatively weak.In this case, the mature subduction transference is not predicted.Similar modes without subduction jump have been obtained in Vogt and Gerya (2014) and S. H. Yang et al. (2018).Further on, other major factors for oceanic plateau collision and accretion have been studied, revealing that a thinned continental margin of the overriding plate, weak layers in the oceanic lithosphere, and a young oceanic plateau are favorable conditions for subduction transference (Tao et al., 2020).In addition, the dynamics of subduction transference have been systematically studied in Zhong and Li (2020) by the new numerical model with a constant convergence force, which reveals that a large force is generally required to trigger and sustain subduction transference, and this required force will decrease with the young and weak oceanic lithosphere.The above studies indicate that the rheological strength of the incoming terrane or oceanic plateau plays important roles in the subduction transference, and the forced convergence is generally required.
For the far-field subduction, the laboratory and numerical modeling has rarely been conducted to investigate the dynamics of such kind of SI during terrane collisions.Holt et al. (2017) and Q. W. Zhang et al. (2017) have presented 3-D numerical models to reproduce the geodynamic processes of divergent double-sided subduction.However, pre-existing subducted slabs are prescribed in the original model setup of both studies, which means they are focusing on the geodynamic effects of double subduction rather than the initiation of the far-field subduction zone.
The current comprehensive numerical study combines all the above modes with a multi-terrane-collision model, and more importantly, obtains the controlling factor of the mode selection among them.The relative lithospheric strength of collisional terranes plays a critical role in determining the site and mode of subduction jump.The SI prefers to occur beneath the weakest overriding terrane.

Tectonic Background of Lithospheric Weakness
From the numerical models and force balance analyses, it indicates that the weaker terrane controls the site of SI occurrence.Thus, the tectonic background and mechanism of lithospheric weakening are discussed here.Laboratory experiments show that the presence of water in the middle and lower continental crust of the island arc may significantly reduce its strength (Kohlstedt et al., 1995).In addition, the magmatism of the active island arc will lead to a higher thermal gradient, which will also decrease its rheological strength.Thus, the (relic) island arc is generally considered as one of the widely distributed weak terranes on the Earth (Almeida Zhang & Leng, 2021).Similarly, the back-arc basins are generally composed of relatively weaker plates/terranes, due to the hydration from the subducting slab, relatively higher strain rate and lower lithospheric age (Arcay et al., 2006;Billen & Gurnis, 2001).Moreover, the transform faults or oceanic fracture zones are vertical weak belts cutting through the whole plate, leading to lithospheric weakening, which are widely considered to facilitate the SI (Baes et al., 2011;Boutelier & Beckett, 2018;Hall et al., 2003;Mueller & Phillips, 1991;Shuck et al., 2021;Sutherland et al., 2020;Toth & Gurnis, 1998;Zhong & Li, 2023;Zhou et al., 2020).As a summary, the lithospheric weakness may come from variable tectonic processes; this study does not aim to distinguish them, but instead applies variable degrees of terrane weakening to test their effects on the mode selection of terrane collision-induced SI.

Implications for Western Pacific Subduction System
In the western Pacific region, a number of terrane collision-induced SI cases have been observed/proposed (Figure 1), which cover all the possible modes as shown in this study (Figures 7 and 8).Thus, these natural cases are further compared with the numerical models to obtain the more general and robust rules of collision-induced SI dynamics.Because the western Pacific region is very large and has many subduction zones, it is further divided into three domains, that is, Northwest Pacific, Southeast Asian Pacific, and Southwest Pacific, which are discussed separately.It is worth noting that the plate motion histories in some regions are not clearly determined but have to be deduced from the plausible plate reconstructions.Thus, we are focusing on the more general trend of terrane collision and SI occurrence rather than deciphering the detailed evolution of specific subduction zones.

Northwest Pacific Domain
The Northwest Pacific has experienced continuous southeast-northwest convergence since the Cenozoic (Figure 11), with complete subduction of the Kula Plate, Kronotsky Plate, and Olyutorsky Plate, as well as the formation and accretion of Olyutorsky arc and Kronotsky arc (Lallemand & Arcay, 2021;Vaes et al., 2019;G. X. Yang, 2022;and references therein).The complex evolution is accompanied by multiple stages of terrane collision-induced SI, which provides a good opportunity for studying the collision-induced SI mode selection.
In the northwestern region, the Eurasian Plate subducts beneath the Olytorsky Plate, generating the Olyutorsky arc which breaks into two segments later (Figure 11a) (Domeier et al., 2017;Müller et al., 2016;Shapiro et al., 2008;Vaes et al., 2019).During 55-50 Ma, the north Olyutorsky arc collides with the west Kamchatka, which terminates the southeastward subduction of the Eurasia Plate (Stage-A in Figure 11b).Then, a new northwestward Proto-Kamchatka subduction zone forms beneath the weak Olyutorsky arc (Stage-A' in Figure 11b), as a typical case of subduction polarity reversal (G.X. Yang, 2022).On the other hand, the south Olyutorsky arc collides with Sakhalin at around 45 Ma, inducing the initiation of northwestward Kuril subduction beneath the weak Olyutorsky arc (Stage-A' in Figure 11c), which is a similar case of subduction polarity reversal (Vaes et al., 2019).These two stages of subduction polarity reversal form a long subduction zone, as shown by A'-A' in Figure 11c.
In the southeastern region, the Kula Plate has subducted beneath the North America and Kronotsky Plate since 60 Ma, generating the Kronotsky arc (Figures 11a and 11b).Then, the Kronotsky arc moves northwestwards and firstly collides with the Central Aleutian subduction zone (C-& W-Aleutians), which results in a southwest dipping subduction zone at the Bowers Ridge (Stage-B' in Figure 11c), which can be regarded as another case of subduction polarity reversal.However, the Proto-Kamchatka subduction drags the westward motion of Kronotsky arc.The subduction at the Bowers Ridge does not maintain for a long period, which starts at about 34 Ma and demises at 26 Ma, as revealed by the volcanic arc records.Afterward, the C-&W-Aleutians recovers subduction (Figure 11d) (Vaes et al., 2019;Wanke et al., 2012).The Kronotsky arc keeps the northwestward moving until it collides with the accreted Olyutorsky arc, blocks the Proto-Kamchatka subduction zone, and induces a new Kamchatka subduction at around 13 Ma (Lallemand & Arcay, 2021;Vaes et al., 2019), as a typical case of subduction transference (Stage-A'' in Figure 11e).
As a summary, the reconstruction of evolution history and terrane collision-induced SI in this region mainly contains two series, that is, A-A'-A'' and B-B' (Figure 11).First, the collision and accretion of the Olyutorsky arc to the Eurasian plate induces two cases of subduction polarity reversal, which are the Proto-Kamchatka and Kuril subduction zones (Stage-A' in Figure 11).Subsequently, the collision of Kronotsky arc with the Olyutorsky arc induces subduction transference, forming the Kamchatka subduction zone (Stage-A'' in Figure 11).As an additional case, the Kronotsky arc collides temporarily with the Central Aleutian subduction zone during its northwestward movement, which results in a subduction polarity reversal at the Bowers Ridge (Stage-B' in Figure 11).In all these collision-induced SI cases, the oceanic plate favors subduction beneath a relatively weak terrane, such as the Proto-Kamchatka and Kuril subduction beneath the relatively weak Olyutorsky arc, the Kamchatka subduction beneath the relatively weak Kronotsky arc, and the subduction at the Bowers Ridge beneath the weak mantle wedge induced by the pre-existing C-&W-Aleutians subduction.These natural cases confirm the rules obtained from the numerical and analytical models in this study.

Southeast Asian Pacific Domain
Southeast Asia is one of the most complicated tectonic regions on the Earth (Figure 12) and is surrounded by several major subduction zones, including the Izu-Bonin-Mariana subduction in the northeast, the Sumatra-Java-Banda subduction in the south, and the South China Sea plate subduction in the northwest.The oceanic basins within this region are generally small and have younger and weaker lithospheres.The oceanic plates are separated by many island arcs and continental terranes, which thus provide a good opportunity for studying the dynamics of collision-induced SI.
The northwestward motion and clockwise rotation of the Philippine Sea Plate are accompanied by the late Eocene-Oligocene collision with the Caroline Plate (Wu et al., 2016).The north-south-trending eastward subduction of Manila, Negros and Cotobato (Figure 12) may develop contemporaneously around the middle Miocene from transform boundaries due to the change in the Philippine Sea Plate motion (Lallemand & Arcay, 2021;Pubellier & Morley, 2014).Then, the Palawan and Sulu continental terranes are dragged eastward and collided with the weak Philippine archipelago terranes in the late Miocene.The terrane collision leads to a subduction polarity reversal with forming the southwestward subduction of the Philippine Sea Plate at around 5 Ma (Stage-A' in Figure 12) (Lallemand & Arcay, 2021).During this process, the Palawan continental block collided with the Palawan arc, which promotes the subduction transference from Palawan to the young Sulu basin beneath the weak Sulu ridge at 15-12 Ma (Stage-B-B' in Figure 12) (Lallemand & Arcay, 2021;Pubellier & Morley, 2014;Wu et al., 2016).
In the center of this domain, the Molucca Sea plate was once a relatively wide oceanic basin with almost 1,600 km in width (Hall & Spakman, 2015;Wu et al., 2016).At around 30-25 Ma, the Molucca Sea plate starts to subduct westward; but later, a new subduction zone is formed on the other side of the Molucca Sea plate at about 10 Ma (Stage-C' in Figure 12).The present-day structures indicate that the westward dipping slab has already reached ∼900 km, whereas the eastward dipping slab has a depth of about 375 km (Wu et al., 2016).Such kind of divergent double subduction has rarely been observed in the global subduction system, which may be related to the far-field subduction after a terrane/plateau collision on the other side of the oceanic plate.
In the south of this domain, the inception of the newly formed Flores and Wetar subduction zone has been regarded as another case of subduction polarity reversal (Stage-D' in Figure 12).With continuous northward convergence, the Australia Plate collides with the Banda arc at 9.8-9.5 Ma (Keep & Haig, 2010) or 7.8 Ma (Tate et al., 2015).The collision leads to stress accumulation and shortening of Wetar segment and further results in the inception of subduction polarity reversal at about 4.5 Ma (Tate et al., 2015).Meanwhile, the subduction of Roo Rise as a massif of seamounts in the Java trench provides external force for the shortening and possible SI at the west Flores since 2 Ma (Silver et al., 1983).Then, the east Wetar and west Flores link together as the early stage of subduction polarity reversal with southward dipping beneath the weak arc belt (Stage-D' in Figure 12).
In this complex tectonic region, a large number of subduction zones have been formed in the Cenozoic, including all the three modes of terrane collision-induced SI.These natural cases confirm that the oceanic plate favors SI beneath the weaker overriding plate, such as the Philippine plate subduction beneath the Philippine archipelago terranes, the Sulu plate subduction beneath the possibly weak Mindanao, the inception of Flores and Wetar beneath the weak island arc belt generated by the Sumatra-Java-Banda trench as well as the far-field subduction of Molucca Sea Plate beneath Halmahera island.

Southwest Pacific Domain
In this domain, the Ontong Java Plateau (OJP) plays a key role in the collision-induced SI (Figure 13).The OJP is one of the largest igneous provinces on the present-day Earth, with a crustal thickness of 30-35 km and an area of 1.86 × 10 6 km 2 (Miura et al., 2004).The Pacific Plate subducts along the North Melanesian-Solomons-Vitiaz trench (Figures 13a and 13b).The collision between the OJP and the North Melanesian Arc System at ∼15 Ma leads to the southwestward SI of Trobiand as a subduction transference from the North Melanesian to the Solomon Sea (Stage-A' in Figure 13b) (Schellart et al., 2006).Then, the OJP collides with the Vanuatu Arc first in the east at 16-10 Ma (Mann & Taira, 2004), inducing the subduction polarity reversal across the weak Vanuatu Arc, forming the New Hebrides subduction zone (Stage-A' in Figure 13c) at around 10-8 Ma (Auzende et al., 1988).With continuous convergence, the collision of OJP and Vanuatu Arc propagates westward gradually, reaching the Solomons at around 10 Ma.It leads to the initiation of San Cristobal and New Britain subduction beneath the weak North Melanesian Arc at 10-5 Ma, as a typical subduction polarity reversal (Stage-A' in Figure 13d) (Mann & Taira, 2004).Under the competition with the New Britain-San Cristobal subduction zone, the Trobriand fails and ceases at about 0.5 Ma (Mann & Taira, 2004).With the increasing resistance of northward New Britain-San Cristobal-New Hebrides subduction, a new subduction zone may be initiated along the Samoa-Gilbert-Ralik (Kroenke & Walker, 1986;Lallemand & Arcay, 2021;Okal et al., 1986).This new SI (Stage-A'' in Figure 13e), as a subduction polarity reversal, may play an important role in the future convergence between the Australian and Pacific Plate after the cessation of New Britain-San Cristobal-New Hebrides subduction.This subduction termination may be due to the entering of buoyant blocks, such as the Louisiade Plateau, Renell Ridge, d'Entrecasteaux Zone, or Loyalty Ridge, into the trench (Lallemand & Arcay, 2021).
As a summary, the large-scale convergence between Pacific and Australian Plates since 20 Ma has mainly experienced three stages of terrane collision-induced SI.The first stage is the subduction transference from the southward North-Melanesian-Solomons-Vitiaz subduction to the southward Trobriand subduction.It is followed by the subduction polarity reversal of the northward New Britain-San Cristobal-New Hebrides as the second stage.Then, the newly formed southward Samoa-Gilbert-Ralik subduction may be a new inception of subduction polarity reversal as the third stage.The North Melanesian Arc and Vanuatu Arc, formed during the previous subduction-induced magmatism, act as the weak overriding plate in the SI.

Competition During Natural SI Mode Selection
During the tectonic evolution of the western Pacific, an interesting phenomenon is the competition between different SI modes, which finally results in a single "winner," or both.For example, in the Southwest Pacific, the subduction transference mode of Trobriand competes with the subduction polarity reversal mode of the New Britain and San Cristobal subduction zones (Figure 13).Finally, the latter wins with the cessation of the earlier Trobriand subduction.This case is similar to the Model B in Figure 9.
In the Southeast Asian Pacific (Figure 12), the northwest dipping Celebes Sea plate subduction is induced by the collision of continental terrane and Palawan arc as a case of subduction polarity reversal from the Palawan trench (Lai et al., 2021).Later on, the southward Sulu plate SI wins the competition and dominates the convergence as a case of subduction transference with the cessation of the earlier Celebes Sea subduction.
The SI mode competition can also be observed in the Northwest Pacific (Figure 11), for example, the complex competition between Bowers Ridge and C-&W-Aleutians subduction zones.After the arriving of Kronotsky arc, the pre-existing C-&W-Aleutians subduction ceases, leading to the formation of a new subduction zone at Bowers Ridge; however, when the Kronotsky arc moves away, the previous C-&W-Aleutians subduction restarts again (Lallemand & Arcay, 2021;Vaes et al., 2019).
The numerical models reveal that the SI mode competition prefers high convergence velocity as well as similarly weak overriding terranes (Figures 7 and 8).Thus, the reason for the above-mentioned mode competition may be that the terranes in the Aleutians, Sulu and Trobriand regions are similarly weak, respectively, due to the magmatism of previous subduction.Meanwhile, the high convergence velocity in these regions may also be a key factor, with a convergence rate of about 10 cm/yr between the Australian and Pacific Plates (Hall, 2002;Schellart et al., 2008;Tregoning et al., 1998), about 9 cm/yr between the Philippine Sea Plate and South Asia (Wu et al., 2016), and about 10 cm/yr of Pacific Plate beneath Aleutians (Schellart et al., 2008).

Conclusions
Systematic numerical models and analytical studies have been conducted to investigate the dynamics of terrane collision-induced SI, focusing on the Cenozoic cases in the Western Pacific.The main conclusions of this study include: 1.The rheological strength of overriding plate controls the mode selection of collision-induced SI; a new subduction zone prefers to form beneath a weaker overriding terrane.2. The convergence velocity also plays an important role in the terrane collision-induced SI; faster convergence promotes the SI, whereas slower convergence fails to generate enough strain localization for the SI. 3. Subduction polarity reversal is the most widely distributed SI mode in the Western Pacific, possibly due to the weakness of the overriding plate during the preceding subduction-induced fluid/melt activity, including the Olyutorsky arc, Kronostky arc, Philippine archipelago, Banda arc, North Melanesian arc, Vanuatu arc and so on.4. Subduction transference is resulting after a collision when the overriding plate is rheologically strong or has continental affinity, for example, the Kamchatka subduction zone. 5. Collisional-induced far-field subduction has rarely been observed, because it is generally difficult for the collision-induced compressional stress to be transferred across the whole lower plate to the far-end.The Sangihe and Halmahera case may provide as a single candidate for this mode in the western Pacific.6.The competition between SI modes is predicted when the convergence velocity is high and the collision terranes are similarly weak.This may explain the complex evolution histories of Bowers Ridge and Aleutians, Sulu trench and northwest Celebes subduction, and Trobriand and New Britain-San Cristobal-New Hebrides.

Figure 1 .
Figure 1.Distribution of subduction zones in the western Pacific region (modified after Lallemand & Arcay, 2021; and references therein).The colors of subduction zones indicate different types of collision-induced subduction initiation, with the legends shown at the bottom-left corner of the map and graphically illustrated below the map.The background topographic map is plotted using Generic Mapping Tools (Wessel et al., 2013), with the data from Tozer et al. (2019).

Figure 2 .
Figure 2. (a) Original model configuration with three continental terranes marked as Terrane 1, Terrane 2, and Terrane 3, two oceanic plates marked as Oceanic Plate 1 and Oceanic Plate 2, as well as a pre-existing subduction zone marked as Downgoing Slab.A weak zone is applied initially between the Downgoing Slab and Terrane 3 to ensure the location of the first subduction.The model is pushed by a constant velocity block marked by the yellow arrow in Terrane 1.The yellow dashed lines indicate the phase transitions at 410 and 660 km in depth.The inset enlarged figures represent the location of ocean-continent transitions (OCT).(b) The beginning of terrane collision as the initial model in the following study.The snapshot represents a model stage after convergence for 10 Myr, with a total amount of convergence of 500 km.The inset-enlarged figures show the deformed OCT after convergence.The dark line within the lithospheric mantle indicates the lithosphere cooling and growth, that is, the newly formed lithospheric material below this line.The colorbar indicates the composition field of the model: 1-sticky air; 2-sea water; 3, 4-sediment; 5-continental upper crust; 6-continental lower crust; 7-oceanic upper crust; 8-oceanic lower crust; 9-oceanic lithospheric mantle; 10-asthenosphere; 11-hydrated mantle; 12-hydrated mantle with serpentinization; 13-lithospheric mantle of Terrane 1; 14-lithospheric mantle of Terrane 2; 15-lithospheric mantle of Terrane 3; 16-partially molten sediments; 17-partially molten continental crust; 18-partially molten oceanic crust; 19-partially molten mantle.

Figure 3 .
Figure 3. Model evolution with no subduction jump.The lithospheric mantles of Terranes 1, 2, and 3 are identically strong, with the plastic friction coefficient λ • sin (φ) = 0.6 initially and decreasing linearly to 0.1 with the accumulated plastic strain to ε plas = 1.The colors indicate rock types as in Figure 2. The detailed evolution of the viscosity field is shown in Figure S2 in Supporting Information S1.

Figure 4 .
Figure 4. Model evolution with subduction polarity reversal.The lithospheric mantles of Terranes 1 and 2 are identically strong, with the plastic friction coefficient λ • sin (φ) = 0.6 initially and decreasing linearly to 0.1 with the accumulated plastic strain to ε plas = 1.However, the lithospheric mantle of Terrane 3 is weaker, with λ • sin (φ) = 0.1 initially and decreasing linearly to 0.01 with ε plas = 1.The colors indicate rock types as in Figure2.The detailed evolution of the viscosity field is shown in FigureS3in Supporting Information S1.

Figure 5 .
Figure 5. Model evolution with subduction transference.The lithospheric mantles of Terranes 1 and 3 are identically strong, with the plastic friction coefficient λ • sin (φ) = 0.6 initially and decreasing linearly to 0.1 with the accumulated plastic strain to ε plas = 1.However, the lithospheric mantle of Terrane 2 is weaker, with λ • sin (φ) = 0.1 initially and decreasing linearly to 0.01 with ε plas = 1.The colors indicate rock types as in Figure 2. The detailed evolution of the viscosity field is shown in Figure S4 in Supporting Information S1.

Figure 6 .
Figure6.Model evolution with far-field subduction.The lithospheric mantles of Terranes 2 and 3 are identically strong, with the plastic friction coefficient λ • sin (φ) = 0.6 initially and decreasing linearly to 0.1 with the accumulated plastic strain to ε plas = 1.However, the lithospheric mantle of Terrane 1 is weaker, with λ • sin (φ) = 0.1 initially and decreasing linearly to 0.01 with ε plas = 1.The colors indicate rock types as in Figure2.The detailed evolution of the viscosity field is shown in FigureS5in Supporting Information S1.

Figure 7 .
Figure 7. Phase diagram of terrane collision-induced subduction initiation (SI) with a convergence velocity of 5 cm/yr.The three axes represent the lithospheric mantle strength of Terrane 1, Terrane 2, and Terrane 3, respectively, with the values for the plastic friction coefficient λ • sin (φ) .The colors in the diagram are used for distinguishing the SI modes, as graphically illustrated at the bottom.The reference models are marked in the diagram with different symbols: Star for the model with no subduction jump in Figure3; Circle for the model with subduction polarity reversal in Figure4; Square for the model with subduction transference in Figure5; Triangle for the model with far-field subduction in Figure6.

Figure 8 .
Figure 8. Phase diagram of terrane collision-induced subduction initiation with different convergence velocity (V x ).(a) Phase diagram with a relatively low convergence velocity of 2 cm/yr.(b) Phase diagram with a relatively high convergence velocity of 10 cm/yr.The blue-orange mixed cases represent the occurrence of subduction transference and subduction polarity reversal together during the model evolution, as shown in the cartoon at the bottom-right of the figure.The symbols of ellipse and diamond represent the special cases with a competition between contrasting modes, which will be discussed later in Figure 9.Other details are the same as in Figure 7.

Figure 9 .
Figure 9. Special model cases with a competition between contrasting modes.Both subduction transference and polarity reversal are predicted in Model A, whereas the onset of subduction transference is further replaced by the polarity reversal in Model B. The cases and parameters of Models A and B are marked in Figure 8b with the symbols of ellipse and diamond, respectively.The colors represent the effective viscosity with the colorbar shown at the bottom.

Figure 11 .
Figure 11.The evolution history of the Northwest Pacific during Cenozoic (modified after Vaes et al., 2019; and references therein).The capital characters (A and B) indicate the different tectonic systems, while the prime and double primes (e.g., A', A'') denote the different subduction stages.The detailed meanings of symbols are shown in the legend.The blue arrow indicates the main direction of Pacific plate motion.

Figure 12 .
Figure 12.The tectonic structures and subduction zones in the Southeast Asian Pacific (modified after Lallemand & Arcay, 2021; and references therein).The capital characters (A, B, C, D) indicate the different tectonic series, while the prime (e.g., A') denotes the different subduction stages.The detailed meanings of subduction modes are shown at the bottom and other properties are the same as in Figure 11.

Figure 13 .
Figure 13.The reconstruction of Southwest Pacific evolution since 20 Ma (modified after Mann & Taira, 2004; and references therein).The detailed meanings of subduction modes are shown in the legend and other properties are the same as in Figure 11.