Subduction Initiation During Collision‐Induced Subduction Transference: Numerical Modeling and Implications for the Tethyan Evolution

The collision‐induced subduction transference is a composite dynamic process including both the terrane collision/accretion and the subduction initiation (SI) at the neighboring passive margin. This process occurred repeatedly during the evolution of Tethyan systems, with multiple ribbon‐like continents or microcontinents drifting from Gondwana in the southern hemisphere and accreting to the Eurasian continent since Paleozoic. In the previous numerical studies, the dynamics of terrane collision and induced SI are investigated individually, which however need to be integrated to study the controlling factors and time scales of collision‐induced subduction transference. Systematic numerical models are conducted with variable properties of converging plates and different boundary conditions. The model results indicate that the forced convergence, rather than pure free subduction, is required to trigger and sustain the SI at the neighboring passive margin after terrane collision. In addition, a weak passive margin can significantly promote the occurrence of SI, by decreasing the required boundary force to reasonable value of plate tectonics. The lengths of subducted oceanic slab and accreting terrane play secondary roles in the occurrence of SI after collision. Under the favorable conditions of collision‐induced subduction transference, the time required for SI after collision is generally short within 10 Myr, which is consistent with the general geological records of Cimmerian collision and the following Neo‐Tethyan SI. In contrast, the stable Indian passive margin and absence of SI in the present Indian Ocean may be due to the low convergent force and/or the lack of proper weak zones, which remains an open question.


Introduction
In the plate tectonics and the Wilson cycle, the plate convergence starts with the subduction initiation (SI) and terminates with the continental collision and orogeny. Because of the low density of continental crust, the continental subduction-collision process cannot last for very long time as the oceanic subduction. Thereby, the plate convergence will generally stop after the orogeny with a time scale of several tens of million years, although the special case of India-Asia collision is still ongoing for ≥50 Myr (DeCelles et al., 2014;Najman et al., 2010;Yin & Harrison, 2000;Zhu et al., 2015). On the other hand, with the collision of a buoyant block, a new subduction zone may form in the neighboring plate to accommodate the continuous convergence, which is defined as the "collision-induced subduction transference" (Figures 1a and 1b; Niu et al., 2003;Stern, 2004;Zhu et al., 2011;. It may play a crucial role in the accretion and amalgamation of continents (Stern, 2017).
The collision-induced subduction transference, although without clear Cenozoic examples, is quite popular during the evolution of Tethyan system, which experienced multiple subduction-collisionsubduction processes, as well as the assembly or accretion of terranes (Figure 1c; Stampfli & Borel, 2002;Şengör et al., 1988;Metcalfe, 2013;Wan et al., 2019;Zhu et al., 2011Zhu et al., , 2013. The ribbon-like continents or microcontinents continuously drifted from Gondwana in the southern hemisphere and accreted to the Eurasian continent throughout the evolution of Tethyan system since Paleozoic (Figure 1c). The Tethyan system began with the breakup of the Galatian super terrane from the Gondwana, resulting in the opening of Paleo-Tethyan Ocean. The collision between the Galatian terrane and Laurasia finally closed the Rheic Ocean, with the subduction transferred to the Paleo-Tethys (Stampfli et al., 2013). Then, the Cimmerian terranes broke up from the Gondwana in late Paleozoic, leading to the opening of Neo-Tethyan Ocean. The collision of the Cimmerian terranes with Laurasia closed the Paleo-Tethys, after which the subduction was transferred to the Neo-Tethys in Mesozoic (Stampfli & Borel, 2002;Wan et al., 2019;Zhu et al., 2013). Finally, the Indian plate broke up from Gondwana in Cretaceous and collided with the Eurasian plate in the early Cenozoic (DeCelles et al., 2014;Najman et al., 2010;Searle, 2019;Yin & Harrison, 2000;Zhu et al., 2015). Although the collision has lasted for a long time, there is still no clear sign for the SI in the present-day Indian Ocean (Stern, 2004;. Thus, it is still a mystery and challenge that whether and when the Indian oceanic plate will begin subduction. From the summarized Tethys evolution history, a crucial issue for the continuous plate convergence is the SI after the terrane collision and accretion (Wan et al., 2019), that is, the collision-induced subduction transference (Figures 1a and 1b). Its dynamics is rarely investigated and thus still not resolved. The previous correlated numerical studies are generally focusing on two points. One is the subduction of anomalous blocks within an oceanic plate, including the island arc, oceanic plateaus, or continental fragments, which is generally focusing on the resulting flat subduction (Mason et al., 2010;van Hunen et al., 2002van Hunen et al., , 2004 or the contrasting modes of accretion (e.g., Selzer et al., 2008;Tao et al., 2020;Tetreault & Buiter, 2012;Vogt & Gerya, 2014;Yang et al., 2018). Notably, a "trench jump" after collision is predicted in some of the models (Tetreault & Buiter, 2012;Vogt & Gerya, 2014;Yang et al., 2018), which is however caused by the detachment of weak and buoyant crust of the accreting terrane, rather than initiating a new subduction zone in the neighboring oceanic plate. Another type of model focuses on the SI at passive continental margins, which only deal with two plates and a transition between them (e.g., Baes & Sobolev, 2017;Nikolaeva et al., 2011;Rey et al., 2014;Toth & Gurnis, 1998;Ulvrova et al., 2019;Zhong & Li, 2019). These models indicate that the SI at passive continental margin is not easy, which generally requires special conditions, for example, (1) a thin, weak, and very buoyant continental lithosphere (e.g., Nikolaeva et al., 2011;Marques et al., 2013Marques et al., , 2014Rey et al., 2014); (2) a prescribed weak transition zone between the continental and oceanic plates (e.g., Baes et al., 2011;Toth & Gurnis, 1998); (3) driven by downward mantle flow (e.g., Baes & Sobolev, 2017); or (4) driven by a boundary stress/force (e.g., Zhong & Li, 2019). On the other hand, the natural examples of Atlantic and Indian passive margins, neighboring to relatively old oceanic lithospheres, are generally stable and difficult for SI (Cloetingh et al., 1989;Mueller & Phillips, 1991;Niu et al., 2003;Zhong & Li, 2019). Thus, it indicates that the collision-induced subduction transference during Tethyan evolution should not be so easy, because the multiple Tethyan Oceans generally have an old and thick lithosphere at the passive margin (Müller et al., 2008;Stampfli & Borel, 2002;Stampfli et al., 2013). In this study, we aim to combine the models with terrane collision/accretion and SI at the passive margins, in order to understand the dynamics of collision-induced subduction transference as well as the key point of Tethyan evolution.
Another important issue for the collision-induced subduction transference is the time span between the terrane collision and SI in the neighboring plate. In the Tethyan system, for example, the collision of the Qiangtang-Lhasa terrane, which occurred in Early Cretaceous (~140 Ma), may trigger the northward SI (~137 Ma) of the Indus-Yarlung Zangbo Tethyan Ocean, that is, a branch of the Neo-Tethys (Zhu et al., 2013). The time for the collision between the Iranian terrane and Eurasia was at about 228-201 Ma, whereas the neighboring SI at about 187 Ma (Wan et al., 2019), which indicates the SI occurring at about 20 Myr after the collision. Notably, the times and time spans between collision and neighboring SI in Tethyan system are not well constrained due to the vague geological records. Specially, the Indian oceanic plate does not begin subduction after the India-Eurasia collision for ≥50 Myr (DeCelles et al., 2014;Najman et al., 2010;Searle, 2019;Yin & Harrison, 2000;Zhu et al., 2015), which makes it even more difficult to answer the question from geological observations (Niu et al., 2003;Stern, 2004;.
The required conditions and time scales of collision-induced subduction transference are rarely studied and thus remain unresolved. In order to study these problems, systematical numerical models are conducted with variable properties of converging plates and different boundary conditions. The model results are further compared with the geological records of Tethys and shed new lights on the Tethyan dynamic evolution as well as the plate tectonic theory.

Original Model Setup
The two-dimensional original model is 4,000 km in length and 1,400 km in depth, as shown in Figure 2a. The model includes two oceanic plates and two continental plates. A spreading oceanic plate (1,000 km) is configured on the left side of the model, which is neighbored by a drifting continental terrane (1,000 km). Another oceanic plate (600 km) is set between the drifting continental terrane and the overriding continental plate (1,400 km). A weak zone is initially applied between the subducting oceanic plate and the overriding continental plate. Wet olivine rheology is applied for the weak zone, whereas dry olivine rheology is used for normal mantle rocks (Karato & Wu, 1993). The crust of continental plate is set to be 35 km, including a 20km-thick upper crust and a 15-km-thick lower crust. For the drifting continental plate, the lithospheric thickness is set to be 100 km, whereas the overriding continental lithosphere is 140 km thick. The oceanic lithosphere is composed of a 3-km-thick basalt layer as the upper crust, a 5-km-thick gabbro layer as the lower crust, and a lithospheric mantle layer with the thickness dependent on the age (Turcotte & Schubert, 2002). The subducting oceanic plate is 100 km in thickness and 180 Ma in age, whereas the age of spreading oceanic plate is varied in the numerical studies. A "stick air" layer is applied on the top of the model to allow the deformation of the crustal surface (Crameri et al., 2012;Schmeling et al., 2008). In Figure 2. (a) Original model configuration with an overriding continental plate, a subducting oceanic plate, and an initial weak zone in between. Wet olivine rheology is applied for the weak zone, whereas dry olivine rheology is used for normal mantle rocks (Karato & Wu, 1993). A drifting continental terrane is incorporated in the oceanic plate. The subducting plate is pushed at two different positions as marked by the orange rectangles and arrows, with the same and constant velocity (v x = 10 cm/year). The inset figures show two different configurations of the neighboring passive margin, that is, with a weak zone or not. (b) The beginning of continental collision after 6 Myr with total convergence of 600 km. This snapshot is used as the initial model for the following study, which is either pure free subduction without any pushing, or pushed by a constant force at the left tip of the spreading oceanic plate. The yellow dashed lines denote the 410and 660-km discontinuities. Colors indicate the rock types, specified by the color grid: 1, stick air; 2, sea water; 3 and 4, sediments; 5 and 6, continental upper and lower crust, respectively; 7 and 8, oceanic upper and lower crust, respectively; 9, lithospheric mantle; 10, asthenosphere; 11, hydrated mantle; 12, serpentinized mantle; 13, partially molten sediments; 14, partially molten continental upper crust; 15, partially molten continental lower crust; 16, partially molten oceanic crust; 17, partially molten mantle. The zoomed-in subduction zone shows the hydrated and partially molten rock types. the model, a stick air layer of 10 km thick is set above the continental plate and 12 km above the oceanic plate, in order for the gravity isostasy.
The temperature conditions of the top and bottom boundaries are fixed, with 273 and 2,248 K, respectively. The left and right boundaries are adiabatic, with no heat flux across them. For the thermal structure of oceanic lithosphere, half-space cooling model is applied (e.g., Turcotte & Schubert, 2002). For the continental lithosphere, a linear temperature gradient is applied with 273 K on the surface and 1623 K at the bottom of lithosphere. An adiabatic temperature gradient of 0.5 K/km is applied for the sublithospheric mantle. The phase transitions at 410 and 660 km are applied , with the Clapeyron slopes of 2.5 and −1.0 MPa/K, respectively.
All the velocity boundary conditions of the model are set to be free slip. The subducting plate is pushed with a constant velocity (v x = 10 cm/year) in the left side of the spreading oceanic plate and the middle of drifting continental terrane, as shown in Figure 2a with two orange rectangles and arrows. The same constant convergent velocities are applied for 6 Myr in the original model in order to avoid deformation in the passive margin between the spreading oceanic plate and drifting continental plate.

Model Setup and Boundary Conditions From Collision
The original model in Figure 2a is pushed by a constant velocity of 10 cm/year for 6 Myr, with a total convergence of 600 km. It leads to a snapshot with the beginning of continental collision, as shown in Figure 2b, after which the previous boundary conditions with constant convergent velocity are cancelled. Thus, the original model driven by boundary velocity is only used for the first subduction and provides the slab pull force for the further evolution of the model.
In the following study, the initial collision snapshot (Figure 2b) is employed as the initial model. Two different types of boundary conditions are further applied and compared. One is the self-consistent free subduction driven purely by the slab pull, that is, no pushing anywhere. Another one is pushed by a constant force at the left tip of the spreading oceanic plate (Figure 2b), the effect of which is combined with the slab pull from the subducted slab. The conditions for the SI of spreading oceanic plate after the continental collision will be systematically investigated. In this study, the SI is defined as the slab bending smoothly and reaching the depth of 150-200 km. Thus, the time scale of SI after collision is defined by the time span between collision and slab reaching~150 km.

Free Subduction After Collision
In this set of models, no prescribed convergent boundary conditions are applied after the initial collision. In order to study the collision-induced subduction transference and SI at the passive margin, two different models with contrasting rheological strengths of passive margin are applied. In Model-1, no weak zone is prescribed (Figure 3a). In Model-2, an initial weak zone with wet olivine rheology (Karato & Wu, 1993) is applied in the ocean-continental transition (OCT) zone (Figure 3b).
The evolutions of these two models are similar before the slab break-off (cf. Figures 3a and 3b). The dehydration occurs during the previous subduction of oceanic crust, which results in a weak subduction channel. The drifting continental lithosphere subducts under the fixed overriding continental plate, driven by the slab pull from the subducted oceanic plate. However, due to the low density, a part of the continental crust is detached and exhumed to the surface or the crustal level. Another part of continental crust could be dragged into the mantle. The buoyancy of subducted continental crust competes with the slab pull from the sinking oceanic plate. Finally, the slab break-off occurs around the OCT. Consequently, the subducted continental plate migrates upward with the loss of slab pull (i.e., eduction; Duretz et al., 2012), which pushes the neighboring passive continental margin and modifies the stress and strain rate fields as shown in Figures 4a and 4b. In Model-1 without any prescribed weak zone, the eduction-induced deformation is rather limited that no clear sign for SI can be observed (Figure 3a). There is no significant deformation or stress localization in the neighboring passive margin after the break-off (Figure 4a). In Model-2 with a prescribed weak zone, an incipient subduction initiates at the beginning (Figure 3b), which stops a short time later due to the limited pushing from eduction. The deformation and stress localization emerge in the neighboring passive margin corresponding to the inception of SI, which however quickly fades with time ( Figure 4b). Thus, the collision-induced subduction transference is not applicable in such cases with free subduction after collision (Figure 3).

Forced Subduction After Collision
In nature, a horizontal plate without connecting to any subducting slab could also be dragged to move by the pushing from mid-ocean ridge and shearing from the mantle convection (Mahatsente, 2017;Niu et al., 2003;Sun, 2019;Turcotte & Schubert, 1982). In this section, we conduct systematic numerical models with a constant convergent force employed on the spreading oceanic plate (Figure 2b) to investigate the conditions and time scale of collision-induced subduction transference. The models are classified into two groups, that is, without or with a prescribed weak zone at the passive margin.

No Prescribed Weak Zone at the Neighboring Passive Margin
Systematic numerical models are conducted with variable convergent boundary forces and variable ages of the spreading oceanic plate ( Figure 5). The model results indicate that the collision-induced SI at the neighboring passive margin requires a relatively high boundary force, that is, >8.0 × 10 12 N/m, which slightly increases with older oceanic lithosphere. If the boundary force is larger enough for a specific passive margin, the SI generally occurs within 10 Myr after the initial collision ( Figure 5). This time scale of SI is also dependent on the boundary force, that is, the larger boundary force resulting in earlier SI occurrence.
The detailed evolutions of two representative models are further demonstrated, that is, Model-3 with oceanic lithosphere of 80 Ma and boundary force of 9.0 × 10 12 N/m, as well as Model-4 with oceanic lithosphere of 80 Ma and boundary force of 8.5 × 10 12 N/m (Figures 6 and 7). In Model-3, the SI occurs at the neighboring  (Figure 7a) demonstrates a significant stress localization at X = 1,600-2,000 km, that is, around the OCT in the neighboring passive margin (Figure 6a). This process lasts for about 2 Myr and finally leads to SI with high strain rate and release of stress (Figure 7a).
As a comparison, the relatively low boundary force in Model-4 prevents the collision-induced SI at the neighboring passive margin (Figure 6b). In this regime, the general stress building occurs slowly and in a wider region, which leads to weak strain localization in the neighboring passive margin (Figure 7b). Thus, the SI is not predicted. In more details, there are two peak stress localization processes around the passive margin. The first one occurs at about 10 Myr after collision, which is corresponding to the detachment and exhumation of the subducted continental crust. However, it just lasts for a short time and does not result in SI. The later stress localization occurs at around 30 Myr (Figure 7b), corresponding to the slab break-off, which again does not lead to abrupt failure of the passive margin. Finally, although the extensive continental collision occurs for >30 Myr, no SI is predicted at the neighboring passive margin.

Existence of a Prescribed Weak Zone at the Neighboring Passive Margin
The existence of a prescribed weak zone at the passive margin dramatically reduces the required convergent force for SI after collision. For example, the boundary force required for SI at the neighboring passive margin with 60-Ma oceanic lithosphere is~4.0 × 10 12 N/m (Figure 8), which is about half of that in the models without weak zone. Two modes of SI are obtained according to systematic numerical studies, which are the earlier SI after initial collision and the later SI after slab break-off (Figure 8). In the models with younger oceanic plate and relatively larger boundary force, the SI results immediately after the initial collision. However, with older oceanic plate and relatively lower boundary force, the SI occurs much later after the slab break-off. Additionally, it is worth noting that models with very low convergent forces may encounter serious problems in the solver convergence, which thus prevents modeling with convergent forces lower than 3.0 × 10 12 N/m. The evolution of two contrasting models are demonstrated in Figures 9 and 10, which are either SI after collision (Model-5 with oceanic lithosphere of 80 Ma and boundary force of 4.0 × 10 12 N/m) or SI after slab break-off (Model-6 with oceanic lithosphere of 80 Ma and boundary force of 3.0 × 10 12 N/m).   . Two contrasting modes of forced subduction initiation (SI) at the neighboring passive margin with a prescribed weak zone, that is, earlier SI after initial collision versus later SI after slab break-off, dependent on the age of oceanic lithosphere and the convergent boundary force.

Journal of Geophysical Research: Solid Earth
In Model-5, the SI occurs at the neighboring passive margin soon after the continental collision (Figure 9a). The significant stress localization in the passive margin leads to the failure of prescribed weak zone and the occurrence of SI (Figure 10a). The forced convergence is accommodated by both subduction zones, that is, the old one on the right and the newly formed one on the left (Figure 9a). However, through the model evolution, the strain rate in the old subduction zone decreases, whereas that in the new subduction zone increases (Figure 10a). It indicates that the dominance of convergence switches gradually from the old subduction zone to the new one.
In Model-6, the relatively low convergent force is not enough to trigger SI immediately after collision (Figure 9b). Although SI does not occur, there is still deformation localized in the weak zone after collision, which is clearly demonstrated by the stress building and resulting higher strain rates (Figure 10b). However,

Journal of Geophysical Research: Solid Earth
the deformation in the passive margin is too slow for the occurrence of SI. At about 18.5 Myr after collision, the slab break-off occurs with eduction of continental plate, which results in an additional push on the neighboring passive margin. Finally, the SI is induced, which is further sustained by the continuous convergent boundary force.

Effect of the Length of Accreted Continental Terrane
In the previous models, the length of drifting continental terrane within the subducting oceanic plate is constant of 1,000 km. During the Tethyan evolution, the multiple terranes accreted to the Eurasian plate may be wider (e.g., Indian plate) or narrower (e.g., Lhasa terrane). However, it is worth noting that the Lhasa terrane is included in the Cimmerian block, the length of which is hard to constrain. In order to understand the effect of the terrane's length on the SI at the neighboring passive margin, comparable experiments are further conducted with shorter terrane of 500 km.
Model-7 and Model-8 are comparable to Model-3 ( Figure 6a) and Model-4 (Figure 6b), respectively. The only difference is the application of a shorter drifting continental terrane (500 km). The evolutions of these two additional models are similar to Model-3 (Figure 6a), in which the SI occurs at the neighboring passive margin after the continental deep subduction (Figure 11). The comparisons between Model-8 (SI in Figure 11b) and Model-4 (no SI in Figure 6b) indicate that the shorter drifting continental terrane favors the collisioninduced subduction transference.
In order to better understand the effects of continental terrane's length, we further conduct a set of numerical models with shorter continental terrane and compare them to the previous models with longer

10.1029/2019JB019288
Journal of Geophysical Research: Solid Earth continental terrane (Figure 12). The results indicate that the length of drifting continental terrane does not play significant roles in the models with higher boundary forces. However, in the regime with lower boundary forces, the collision of a shorter continental terrane leads to SI a bit earlier (by 0-2 Myr) and easier (by decreasing the required boundary force of 0.6 × 10 12 N/m) than a longer one (Figure 12).

Effect of the Length of Firstly Subducted Oceanic Slab
In the previous models, the length of the firstly subducted oceanic slab is constant of 600 km, which leads to the slab arriving at around 660-km discontinuity, that is, the boundary between the upper and lower mantles. Further varying the length of subducted slab may play a role in modifying the force of slab pull. In order

Journal of Geophysical Research: Solid Earth
to test its effects on the SI at the neighboring passive margin, additional experiments are conducted with longer subducted oceanic slab of 1,000 km (Figures 13 and 14).
In the regime with relatively lower convergent force of 8.5 × 10 12 N/m, the evolution of Model-10 is very similar to Model-4 although with different lengths of subducted slab (cf. Figures 13b and 6b). No SI is predicted at the neighboring passive margin of both models. Alternatively, in the regime with relatively higher convergent force of 9.0 × 10 12 N/m, the collision-induced SI at the passive margin occurs a bit later for~2 Myr in Model-9 with a longer subducted oceanic slab of 1,000 km than that in Model-3 with a shorter subducted slab of 600 km (cf. Figures 13a and 6a). Figure 14 further summarizes and illustrates the effects of subducted slab length on the time of SI at the neighboring passive margin. It indicates that the time of SI after collision is generally later (by~2 Myr) in the models with a longer subducted slab of 1,000 km than those of 600 km, because the larger slab pull delays the stress building and lithospheric collapse. However, the length of subducted slab does not affect the threshold of convergent force required to trigger SI at the neighboring passive margin. Thus, it plays secondary roles in the collision-induced subduction transference.

Controlling Factors of Collision-Induced Subduction Transference
The collision-induced subduction transference is a complex dynamic process, which is affected by not only the present subduction/collision itself but also the properties of the neighboring passive margin and the boundary conditions. The current numerical models indicate that the resistance of convergence from continental collision leads to stress building in the neighboring plate and passive margin (e.g., Figures 4 and 7). On the other hand, the slab break-off and the resulting eduction can trigger the incipient SI at the neighboring passive margin (e.g., Figures 3b and 9b). However, the steady SI and the final development of mature subduction require additional forces, that is, the boundary convergence force. In Model-2 without boundary force (Figure 3b), the incipient SI does not evolve into a self-sustained subduction zone because of the short and limited convergence from slab break-off. In contrast, the similarly induced SI is sustained by the Figure 12. Comparisons among models with shorter (500 km) and longer (1,000 km) drifting continental terranes. The circles show the results of numerical models as in Figure 5 with spreading oceanic plate of 80 Ma and longer continental terrane of 1,000 km, whereas the triangles for the comparable models with shorter continental terrane of 500 km. All the other parameters are identical.  Figure 9b). Thus, the convergent force on the drifting oceanic plate is a critical factor for the occurrence of collision-induced subduction transference. In addition, the increasing of convergent force could extensively reduce the time for SI after collision, as shown in Figure 5.
The rheological strength of the passive margin is controlled by the age of neighboring oceanic lithosphere as well as the possible presence of weak zones. The numerical models indicate that the collision-induced SI at the passive margin without weak zone generally requires a high convergent force of above 8.0 × 10 12 N/m ( Figure 5). However, the force of ridge push is estimated to be around 3.0 × 10 12 N/m (Ghosh et al., 2006;Harper, 1975;Mahatsente, 2017;Turcotte & Schubert, 1982), which is not large enough for the collisioninduced SI. The existence of a weak zone with wet olivine rheology (Karato & Wu, 1993) at the OCT can significantly reduce the required convergent force for SI after collision (Figure 8), that is, no higher than 3.0-5.0 × 10 12 N/m. On the other hand, the age of the neighboring oceanic plate plays a second-order role in the collision-induced subduction transference, comparing to the effects of weak zone and boundary force, although the required force for triggering SI increases slightly with the oceanic lithospheric age in both models with and without weak zones (Figures 5 and 8). The effects of oceanic age are two-folded, that is, on viscosity and density, respectively. In the viscosity aspect, the old oceanic lithosphere increases the thickness and rheological strength of the passive margin, which leads to difficulty in its collapse and further SI. In the density aspect, the negatively buoyancy of oceanic plate increases with age and thus contributes to the gravity instability and further the SI. However, the increase of strength and resistance to SI dominates, which finally results in the positive correlation between the oceanic age and the required convergent force for SI at the neighboring passive margin (Figures 5 and 8).

Implications for the Dynamics of Tethyan Evolution
The Tethyan system is a large-scale and long-living orogenic system on the Earth since the Paleozoic. The most special and intriguing character of Tethys is the multiple terranes or microcontinents drifting from Gondwana, and accreting to the Eurasia (or Laurasia) continent (Figure 1c; Wan et al., 2019). The drifting and collision processes of the Galatian supercontinent and the following SI of the Paleo-Tethys in the Early Paleozoic are much far away from present and consequently not well constrained from the Figure 14. Comparisons among models with shorter (600 km) or longer (1,000 km) subducted oceanic plate. The circles show the results of numerical models as in Figure 5 with spreading oceanic plate of 80 Ma and subducted oceanic plate of 600 km, whereas the stars for the comparable models with longer subducted oceanic plate of 1,000 km. All the other parameters are identical.

10.1029/2019JB019288
Journal of Geophysical Research: Solid Earth geological records. Thus, we focus on the collision of the Cimmerian terranes with Eurasia and the following SI of the Neo-Tethys. More geological records can be obtained to constrain this collision-induced subduction transference.
The Neo-Tethyan tectonic belt is very long for~10,000 km from the western Europe to the southeastern Asia. Here we focus on the regions with relatively well constrained geochronologic data, as shown in Figure 15. In the middle Tethys, the collision of Iranian block with the Eurasia is recorded by the deformation of early Mesozoic strata in Fariman basin, which indicates that the time of collision is around 228-201 Ma (Wan et al., 2019;Zanchi et al., 2016). The subduction of the neighboring Neo-Tethyan oceanic plate in Iranian region started at about 187-200 Ma (Chiu et al., 2013;Davoudian et al., 2016;Wilmsen et al., 2009) (Metcalfe, 2013). As a summary, the SI of Neo-Tethys is generally following the collision of Cimmerian plates with Eurasia, with a time delay of <20 Myr as shown in Figure 15. It is worth noting that various geological records are previously used as the proxy for the start of collision, that is, the exhumation of eclogite, the strata deformation, and the metamorphic rocks. Thus, the timing of collision is a bit varied among different studies (Aitchison et al., 2007;DeCelles et al., 2014;Najman et al., 2010;Searle, 2019;Yin & Harrison, 2000;Zhu et al., 2015). On the other hand, the geological records for SI are even more complex (Guilmette et al., 2018;Hall, 2018;Pandey et al., 2019;Parlak et al., 2019;Patriat et al., 2019;Searle, 2019;Stern, 2004;Stern et al., 2012).

10.1029/2019JB019288
Journal of Geophysical Research: Solid Earth according to the magmatic records. The current numerical models indicate that the SI occurs in the neighboring plate within 10 Myr after collision, under the favorable conditions of collision-induced subduction transference ( Figure 5). Considering all the uncertainties as discussed above, the numerical models are consistent with the observations in the middle-eastern Tethys (Figure 15).
Another issue for discussion is whether and when will the SI occur in the Indian Ocean as a response to the India-Asia collision, which is a challenging problem. The collision has already occurred for more than 50 Myr (Lippert et al., 2014;Najman et al., 2010;Searle, 2019;Yin & Harrison, 2000;Zhu et al., 2015); however, there is still no clear sign for the SI to the south of Indian continent, although some guesses have been suggested (Niu et al., 2003;Stern, 2004;. Several reasons have been proposed for the absence of SI in the Indian Ocean, although none has been confirmed, for example, the continuous shortening of overriding Tibetan plateau, the large width or the triangular shape of the Indian plate, and the strength of an old and stable OCT (Cloetingh et al., 1989;Stern, 2004;. Based on our numerical models, the difficulty for the SI in the Indian Ocean may be due to the low convergent force and/or the lack of proper weak zones at the passive margin. The convergent force between the Indian plate and the Eurasian continent is estimated by the GPE to be about 3.0 × 10 12 N/m (Ghosh et al., 2006;Schmalholz et al., 2014), which is lower than the required force for triggering SI for the neighboring Indian Oceanic plate with over 100-Ma lithosphere, if no weak zone is present ( Figure 5). In this case, the convergence between plates is mainly accommodated by the deformation in the Himalayan collision zone, which is consistent with Model-4 ( Figure 6b). On the other hand, the absence of SI in the Indian Ocean may also be a result of lack of proper weak zone according to our numerical models. The existence of weak zone is a common case on the Earth, which could be the faults, hydration zones, and highly deformed/fracture zones (Arcay et al., 2020;Gurnis et al., 2004;Leng & Gurnis, 2011;Zhou et al., 2018). These weak zones may not collapse into subduction zone under the typical tectonic force of ridge push of about 3.0 × 10 12 N/m (Mahatsente, 2017;Turcotte & Schubert, 1982) as shown in Figure 8. However, a short but strong impulse (e.g., slab break-off) can result in an incipient subduction zone, which is then sustained by a low convergent force (e.g., 3.0 × 10 12 N/m) in Model-6 ( Figure 9b). Obviously, the Indian plate has experienced the collision and slab break-off (Aitchison et al., 2007;Kohn & Parkinson, 2002;Najman et al., 2010;Yin & Harrison, 2000;Zhu et al., 2015). The absence of SI may indicate the lack of proper weak zones in the Indian oceanic plate to localize the incipient subduction.

Conclusions
The collision-induced subduction transference includes two aspects, namely, the terrane collision/accretion and the SI at the neighboring passive margin, the dynamics of which are generally investigated individually. In this study, we combine these two regimes into an integrated model, and conduct systematic numerical experiments. The main conclusions include the following: 1. The boundary force-driven convergence is required to trigger and sustain the SI at the neighboring passive margin after terrane collision and accretion. In contrast, the self-consistent force variations in the existing subduction/collision system (e.g., induced by continental subduction and/or slab break-off) are not enough for the subduction transference, although they can indeed trigger the incipient SI. 2. The existence of weak zone at the passive margin can significantly promote the occurrence of SI, by decreasing the required boundary force to the reasonable value of plate tectonics. In addition, the age of oceanic lithosphere also plays a certain role by affecting the strength of passive margin, that is, easier SI for younger oceanic plate under the same boundary force. 3. The length of subducted oceanic slab or the accreting continental terrane plays secondary roles in the occurrence of SI after collision. 4. Under the favorable conditions of collision-induced subduction transference, the time required for SI after collision is generally short within 10 Myr, which is strongly dependent on the convergent boundary force, but weakly on the age of oceanic lithosphere. 5. The SI of the Neo-Tethyan oceanic plate generally occurred shortly after the collision between the Cimmerian terranes and Eurasian continent, which may indicate the relatively large convergent force and/or weakened passive margin. In contrast, the stable Indian passive margin and absence of SI in the Indian Ocean may due to the low convergent force and/or the lack of proper weak zones, which still requires further studies.

10.1029/2019JB019288
Journal of Geophysical Research: Solid Earth