Factors on the Mesozoic Transition From Flat to Steep Subduction of the Paleo‐Pacific Plate Beneath South China: Thickened Oceanic Crust and Subduction Rate

In the Mesozoic, the age of Paleo‐Pacific plate that flatly subducted beneath the South China was ∼80 Ma, resulting in a cold and dense oceanic slab. The transition process from flat to steep subduction of this old slab and factors that affect it are still unclear. Using 2‐D thermomechanical numerical models, we investigate the dependence of the transition mode on the size of TOC (length and thickness) and the subduction rate. The results indicate that under low subduction rates (e.g., 1 cm/yr), even a 100‐km long eclogitized TOC can induce slab delamination beneath continental lithosphere. The maximum subduction rate for slab delamination to occur increases with an increasing length or thickness of TOC. However, when the subduction rate is over a certain critical value (e.g., ≥4 cm/yr), oceanic slab rolls back. The presence of the TOC delays the initial rollback time, as compared to cases with no TOC. In addition, with a TOC of small size (e.g., 100‐km length and 20‐km thickness) under a certain range of subduction rate, an 80‐Ma old slab can maintain a long‐lasting flat subduction. The comparisons between the numerical results and global plate motion models indicate that the northwestward motion of South China Block since ca. 195 ± 5 Ma induced the delamination of oceanic slab. The arrival of the TOC in the Paleo‐Pacific Plate to the continental arc should be later than the previously estimated ca. 250 Ma. Determination of the TOC size requires three‐dimensional simulations and more detailed geological records.


10.1029/2022GC010523
3 of 26 simulations suggest that the 1,300-km-wide intracontinental orogeny and widespread magmatism in the Mesozoic South China are affected by the flat subduction of the PPP and the subsequent steepening processes (Chu et al., 2012;Dai et al., 2020;Z. Li & Li, 2007; Z. F. Y. Wang et al., 2011). From ca. 265 Ma to ca. 190 Ma, the orogenic front propagated cratonward from the coastal region in the SCB, suggesting the cratonward motion of the flatly subducted PPP (Z. . The arrival of TOC to the continental arc in the ca. 250 Ma is speculated that induced the flat subduction of the PPP (Z. . From ca. 190 Ma to ca. 150 Ma, the large-scale magmatic activities spread from the center of the orogen to the surrounding regions, showing a northwestward younging trend to the continental interior and a southeastward younging trend to the coastline. It is assumed to be induced by slab delamination in the middle of subducted slab (Dai et al., 2020;Z. Li & Li, 2007). After ca. 150 Ma, the postorogenic intrusions and volcanics in southeast China show a younging trend toward the coastline (Figure 2). This trend is induced by the southeastward rollback of the subducted slab (K. T. Sun, 2006;Zhou & Li, 2000).
In the transition process from flat-to steep-slab subduction, the size (length and thickness) of the TOC and the subduction rate have been proposed as two important factors (Arrial & Billen, 2013;Dai et al., 2020;Huangfu et al., 2016). The thickened oceanic crusts, including seamounts, oceanic plateaux and aseismic ridges, apply opposing effects before and after the eclogitization phase change. In the initial subduction stage, TOC promotes flat subduction, owing to its lower density (3,000 kg/m 3 ) compared with the lithospheric mantle (3,300 kg/m 3 ). However, as temperature and pressure increase in the flat subduction process, the basalt-to-eclogite phase transformation substantially increases the density and the negative buoyancy of the crust (Arrial & Billen, 2013). Whether the eclogitization of TOC (oceanic plateau) is more favorable for slab delamination (Dai et al., 2020;Z. Li & Li, 2007) or slab rollback (Arrial & Billen, 2013;S. Liu & Currie, 2016) is controversial. Besides the TOC, the subduction rate also controls the transition mode of a flat oceanic slab. Models with 300 and 500 km long TOC and subduction rates varying at an interval of 1 cm/yr have been investigated by Dai et al. (2020). Their results show the effect of the length of TOC on the deformation of overriding continental lithosphere. Moreover, their results suggest that a rapid subduction rate leads to slab rollback, whereas a slow subduction rate leads to slab delamination. Nevertheless, the effect of TOC size on slab evolution is still unclear. For example, under the same subduction rate, how does slab evolution vary with the size of TOC? And how does the range of subduction rate that can lead to delamination vary with TOC size? Since the TOC size remains unclear and the subduction rate of the PPP is somewhat controversial (Koppers et al., 2001;W. Sun et al., 2007), we adopt a wider range of TOC size and subduction rates than the previous study of Dai et al. (2020) to answer the above questions. Moreover, results of the reconstructed plate motion models are incorporated in the present study to constrain the evolution of plate motions in our numerical simulation.
In the present study, we use two dimensional thermal-mechanical numerical models to investigate the effects of TOC size and the subduction rate on the slab transition process from flat-to steep-slab subduction for cases with a complete subduction of TOC. In the present models, the size of TOC varies from 0 to 400 km in length and 15-25 km in thickness. The subduction rate varies from 0 to 5.0 cm/yr. Based on the numerical results, the geophysical records and geological data, we discuss the flat subduction and following steepening process of the PPP beneath the SCB in the Mesozoic, including the arrival time to the continental arc and the size of the TOC.

Geological Setting of the Study Area
The SCB locates south of the Qinling-Dabie Orogen, east of the Tibetan Plateau and west of the Pacific Plate. It is composed of the Yangtze and the Cathaysia blocks (Figure 2), which collided in the Neoproterozoic and formed the Jiang-Shao fault zone (Charvet et al., 1996;J. F. Chen et al., 1991). The southeast China is defined here as the east of the SCB, including most of the Cathaysia Block and the eastern Yangtze Block (denoted by the red rectangular box in Figure 2). The Mesozoic southeast China gradually transformed from the Tethyan tectonic domain to the Pacific tectonic domain since the Late Permian (J. H. Li et al., 2017;Song et al., 2016;Y. Q. Zhang et al., 2012). This transformation was significantly affected by the flat subduction and the following steepening of the PPP beneath the SCB (X. Li et al., 2009;Z. Li & Li, 2007), resulting in the Mesozoic tectonic evolution and the widespread magmatism distribution in southeast China (Dai et al., 2020).
The widespread Mesozoic igneous rocks in southeast China are mostly composed of granitoids and the equivalent volcanic rocks (L. Liu, Xu, & Xia, 2014;Zhou et al., 2006). The Triassic (250-200 Ma) igneous rocks in southeast China, mainly composed of the peraluminous S-type granite (T. Sun et al., 2005), are scattered in Hunan Province, the Wuyi Mountain Range and the southeast coastal region. This distribution was affected by the collisions of the SCB with the other plates, including the North China Block, the Indochina Block and the PPP (L. Liu, Xu, & Xia, 2014). The cratonward flat subduction of the PPP beneath the SCB between ca. 265 Ma and ca. 190 Ma resulted in the cratonward migration of the foreland fold-and-thrust belt (Z. . From early to late Jurassic (200 Ma to 160 Ma), the magmatic activities mostly occurred in the Nanling Range and the Wuyi Mountain Range, generating intrusive rocks such as the calcareous granite, the I-type granite (X. H. Li et al., 2003) and effusive rocks mainly composed of basalt and rhyolite with a few andesite (S. Z. Li et al., 2018).   Wang et al., 2003, Y. J. Wang et al., 2005Ye et al., 2013) suggested that crustal partial melting has been induced by mantle-derived basaltic underplating (L. Liu, Qiu, et al., 2014;Zhou & Li, 2000). Between ca. 180 Ma and ca. 152 Ma, the back-arc environment transformed from extrusion to tension with the participation of mantle-derived magma (X. Li et al., 2009;M. L. Zhang et al., 2018), indicating the transition from flat-to steep-slab subduction. The ca. 190 Ma Keshubei A-type granite in southern Jiangxi has been attributed to slab foundering and a following slab delamination (Dai et al., 2020;Z. Li & Li, 2007). From late Jurassic to late Cretaceous (150 Ma to 80 Ma), magmatic activities are widely distributed in the coastal regions and the middle to lower reaches of the Yangtze River, mostly composed of I-type granite (M. L. Zhang et al., 2018) and A-type granite (C. H. Chen et al., 2016;Mao et al., 2009). Over this period, the magmatism distribution showed a coastward younging trend, which suggests slab rolled back southeastward.

Numerical Model Design
To study the dynamic transition process from flat-to steep-slab subduction in the study region, the i2vis code (Gerya & Yuen, 2003a, 2003b, based on the finite difference method and the marker-in-cell technique, is used for two dimensional (2D) thermal-mechanical simulations. The numerical methodology of the code is detailed in Supporting Information S1. A series of large-scale models (4,000 × 670 km) are constructed with a nonuniform 699 × 134 rectangular grid. The grid resolution varies from 2 × 2 km in the subduction zone to 30 × 30 km in the boundary zone. The lithological structure in the model is represented by 7.48 million Lagrangian markers (Gerya & Yuen, 2003a;Ten et al., 1999). For the initial compositional layer distribution, the lengths of the oceanic plate and the continental plate are both 2,000 km. The age of the oceanic lithosphere is set as 80 Ma to simulate the subduction of the PPP as an old oceanic lithosphere (Müller et al., 2016). The thickness of oceanic lithosphere is set based on the half space cooling model (Turcotte & Schubert, 2002) with a 2-km-thick sediment and an 8-km-thick oceanic crust. The frontal part of TOC locates 200 km away from the trench. The length of the TOC varies from 100 to 400 km in different models with the thickness varying from 8 to 25 km. The 120-km-thick continental lithosphere (An & Shi, 2006) is composed of a 6-km-thick sediment, a 14-km-thick upper continental crust, a 15-km-thick lower continental crust and an 85-km-thick lithospheric mantle. Previous studies (Huangfu et al., 2016) suggest that a small initial subduction angle facilitates the formation of flat subduction. According to these studies and our numerical tests, a value of 12° is set for all the cases. This value leads to flat subduction at the initial stage for all the cases. A low density and viscosity fluid (air/water) layer is applied above the sediment layer at the top of the model, which allows the calculation of topography evolution (Gerya & Yuen, 2003a;Matsumoto & Tomoda, 1983).
For the initial distribution of compositional layers, a weak zone with low viscosity is set between the continental and the oceanic domains for initial subduction process ( Figure 3). While the transition of basalt-to-eclogite starts at a temperature of 400°C and a pressure of 1.5 GPa (Hacker et al., 2003), the basalt and gabbro in oceanic crust can be metastable with low water activity (Ahrens & Schubert, 1975;Hacker, 1996;Spohn & Neugebauer, 1978). As a result, the significant density increase of oceanic crust is delayed before the ideal conditions for eclogitization is satisfied, which induces a later occurrence of steepening process (Dai et al., 2020;Hacker et al., 2003). To simplify, the temperature and pressure for the transition of the oceanic crust to the eclogite in our numerical models is set as 650°C and 2.5 GPa, respectively, similar to van Hunen et al. (2002avan Hunen et al. ( , 2002b. Once the pressure and the temperature conditions for eclogitization are both satisfied, the density of oceanic crust is instantly increased to 3,600 kg/m 3 (Dai et al., 2020;Hacker, 1996;van Hunen et al., 2002b). The rock properties and viscous flow laws used in our numerical simulations are shown in Tables 1 and 2, respectively.
The initial temperature distribution uniformly varies from 0°C at the top to 1300°C at the bottom of the lithospheric mantle, with a downward temperature gradient of 0.5°C/km in the asthenospheric mantle (Turcotte & Schubert, 2002). The thermal boundary condition includes two vertical boundaries with zero horizontal heat flux, a constant temperature condition (273 K) at the top and an infinity-like external boundary at the bottom (Burg & Gerya, 2005;Gerya et al., 2008;Ueda et al., 2008). A constant temperature condition is specified at a parallel boundary located at the 1,670 km depth with y = − external ΔL , where external is the temperature at 1,670 km depth (2045°C) and Δ is the distance (1,000 km) from the bottom to the external parallel boundary.

Table 1 Material Properties in Numerical Models
in the computational domain (Turcotte & Schubert, 2002). The external free-slip condition for the bottom boundary of the model follows Δ + = 0 , = 0 , where Δ is the vertical distance (1,000 km) from the bottom boundary to the external boundary, and are horizontal and vertical velocity at the boundary. Both the mantle flow and the vertical heat flux transfer through the bottom permeable boundary.
The prescribed subduction velocity of oceanic lithosphere represents the convergence rate, is dominated by an internal domain adjacent to the outer boundary of oceanic lithosphere (Z. Li et al., 2011). In this paper, we mainly discuss a decreasing or a constant subduction rate. The largest subduction rate in the models is set as 5 cm/yr. The model adopts a subduction rate varying with time and the transition between different subduction rates is set as instantaneous. Due to the differences in the size of the TOC and the initial distance away from the trench, the time of the initial subduction of the TOC and its initial eclogitization vary in different models. Different from the model setting of Dai et al. (2020), the entire subduction process, from the initial unsubducted passive margin to the eventual steep subduction state, is separated into four stages by three signature events: (a) the reach of the front of TOC beneath the trench (b) the initial eclogitization of TOC (c) the complete eclogitization of TOC. Accordingly, the four different stages are: (a) Before the TOC subducts beneath trench; (b) After the TOC starts to subduct beneath the trench and before the TOC starts eclogitization; (c) After the initial eclogitization of TOC and before its complete eclogitization; (d) After the complete eclogitization of the TOC.

Modeling Results
In most numerical cases, flat subduction transforms into steep subduction through either rollback or delamination. These two transform modes are shown in Sections 4.1 and 4.2 respectively, demonstrating the effect of subduction rate on the transform mode. Meanwhile, the effect of TOC size is shown in Sections 4.3 and 4.4. TOC size is shown not only to affects the sagging position in case of delamination (Section 4.3), but can also determine the fate of subduction. In Section 4.4, it is shown that with a relatively small (∼100 km long and 20 km thick) TOC, a long-lasting (more than 80 Myr) flat subduction can be maintained as long as the subduction rate lies within a certain range (e.g., 1.8-2.0 cm/yr).

A Low Subduction Rate Leads to the Delamination of a Slab With TOC
Our results suggest that a low subduction rate is favorable for the occurrence of slab delamination. As will be shown in Section 5.2, the range of subduction rate that can induce delamination expands with increased TOC size. For a 100 km long and 20 km thick TOC, our results show that delamination occurs if the subduction rate is smaller than 1.4 cm/yr (see Movie S1).
For a certain size of TOC, the subduction rate directly determines the location of the subducted TOC, which affects the location of slab sagging and the process of delamination. A sufficiently low subduction rate keeps the subducted TOC away from the frontal part of the flat slab, and therefore prevents slab rollback. A comparison Note. η 0 is the reference viscosity, which is calculated: η 0 = (1/A D ) × 10 6n . φ eff is the effective internal frictional angle implemented for plastic rheology. Other parameters are illustrated in Table 1. The molten felsic are used for partial molten sediment and continental crust. References: Kirby (1983), Kirby and Kronenberg (1987), Ranalli (1995), and Ranalli and Murphy (1987). between numerical results of M3-3 ( Figure 4) and M3-4 ( Figure 5) demonstrates the effect of subduction rate on the evolution process with a 400 km long TOC. In both cases, the TOC has subducted to the trench within 7 Myr under a subduction rate of 5 cm/yr in stage 1. In M3-3, it takes about 21 Myr for the initial eclogitization to occur in the TOC (Figure 4a). With a decreased subduction rate of 1 cm/yr since stage 3 (after the initial eclogitization of the TOC), the flat-slab subduction gradually transforms to steep-slab subduction through slab delamination (Figure 4). In stage 3, the initial eclogitization begins at the front of the TOC, which has subducted for the longest time beneath the continental lithosphere ( Figure 4a). As a result of the strong interplate coupling, the slab keeps subducting flatly with an increasing cross-sectional area of eclogitized oceanic crust enveloped by the upper continental and the lower oceanic lithospheric mantle (Figure 4b).
After most part of the TOC gets eclogitized, the eclogitized oceanic crust accumulates and sinks in the middle of the subducted slab (Figure 4b), driving the sinking of the surrounding oceanic lithospheric mantle. The sagging eclogitized oceanic crust eventually sinks into the asthenospheric mantle (Figure 4c), and this process increases the subduction angle in stage 4. The further inland eclogitized oceanic crust gets enclosed by the surrounding lithospheric mantle during the sagging process (Figures 4c and 4d). Owing to the initial sagging of the slab close to the trench, the further inland residual eclogitized oceanic crust loses the forward momentum to break the frontal enclosure of surrounding lithospheric mantle (Figure 4c). Affected by its high density, this residual eclogitized oceanic crust delaminates into the lower asthenospheric mantle after the initial sagging of the slab, inducing mantle upwelling from the inland area to the trench (Figures 4e and 4f). As the lithospheric mantle retreats and sinks, the surrounding hot asthenospheric mantle takes up its original position and provides a higher thermal influence on the bottom of the continental lithosphere (Figures 4e and 4f).
In the cases of slab delamination, the numerical results of models M3-3 to M3-8 also indicate that the time history of the subduction rate affects the location of slab sagging. Compared with the subduction rate in model M3-3, the subduction rate in model M3-4 is smaller (3 cm/yr) in stage 2 and higher (2 cm/yr) since stage 3. In model M3-4, the initial eclogitization in the TOC takes about 22 Myr (Figure 5a), which is similar to that in model M3-3 ( Figure 4a). With a smaller subduction rate in stage 2, the TOC moves to a shorter distance landward (Figures 4a and 5a). Therefore, the front of the TOC is closer to the trench at the beginning of stage 3. However, the higher subduction rate in stage 3 (2 cm/year in Figure 5 as compared with 1 cm/ year in Figure 4) pushes the eclogitized TOC further inland. The eclogitized oceanic crust accumulates near the front of the slab (Figure 5c), instead of the middle of the slab (Figure 4c). Therefore, slab sagging occurs closer to the front of the slab (Figure 5d). The following retreat from the front of the slab is accompanied with upwelling of the asthenospheric mantle (Figures 5e  and 5f). The break-off of the dense slab front during stage 3 (Figure 5b) is caused by the thermal relaxation with the temperature increase, which decreases the strength of the slab interior and activates thermomechanical necking in dislocation creep region (Gerya et al., 2004;Zlotnik et al., 2008).
For an 80 Ma oceanic slab, a series of numerical models (M2-1 to M2-5) indicates that a normal-thick oceanic crust will only induce rollback, even under a low subduction rate of 1 cm/yr. Our numerical results show that under the subduction rate of 1 cm/yr since stage 3, the minimum length of the TOC for slab delamination to occur is approximately 100 km.

A High Subduction Rate Leads to Slab Rollback Regardless of the Size of the TOC
Under a sufficiently high subduction rate (e.g., 3 cm/yr), the slab rolls back once the TOC subducts to the frontal part of the flat slab. In model M3-2, with a subduction rate of 3 cm/yr since stage 3, instead of being enclosed between the continental and the oceanic lithospheric mantle in model M3-3 (Figures 4b and 4c), the eclogitized oceanic crust breaks through the enclosure of lithospheric mantle and subducts to the front of flat slab (Figures 6a  and 6b). As the slab rolls back, the eclogitized oceanic crust peels away from the overlying continental crust (Figure 6b). This process is accompanied with the upwelling of asthenospheric mantle and the thinning of the continental plate (Figure 6b). Compared with M3-3 and M3-4, no significant accumulation of eclogitized oceanic crust is shown in the process of slab rollback. In the process of rolling back, the slab breaks off at the boundary between the eclogitized and the uneclogitized TOC (Figure 6c), mainly owing to the density difference between eclogite (3,600 kg/m 3 ) and basalt (3,000 kg/m 3 ). Subsequently, the residual TOC gets eclogitized and induces the subsequent rollback (Figure 6d). In the process of rolling back, the slab transforms from flat to steep subduction in approximately 5 Myr, much faster than the slab delamination process which takes ∼40 Myr.
Comparing the numerical results of models M2-2 and M3-2, the presence of TOC does not significantly accelerate or slow down the slab rollback process. The model M2-2 with normal-thick oceanic crust demonstrates a transitional process from flat subduction to slab rollback similar to that in M3-2, except for the starting time of slab rollback and the absence of slab break-off (Figure 7). Without the TOC, the slab rolls back slightly earlier due to the less buoyancy force from oceanic crust against the driving force from the slab dense front. The dense slab front mainly induces the rollback process ( Figure 7b). In addition, the normal-thick oceanic crust provides insufficient density contrast for the occurrence of slab break-off under the subduction rate shown in M2-2 (Figures 7c  and 7d). The high subduction rate leads to a continuous subduction process and a relative low slab inner temperature, which helps to prevent the break-off of the slab dense front.

The Effect of TOC Size on the Evolution of Flat-Slab Subduction
In the case of slab delamination, apart from the subduction rate, the size of the TOC can also affect the location of slab sagging and the scale of delamination. We test a number of models with various sizes of the TOC, from 100 to 400 km in length and from 15 to 25 km in thickness, under different subduction rates (Table 3). The numerical results show that the size of TOC affects the location of slab sagging.
Under the same subduction rate, slab sagging occurs further inland as TOC decreases in length. Figure 8 demonstrates the various sagging locations with the length of the TOC ranging from 100 to 400 km under a subduction rate of 3 cm/yr in stage 2 and 2 cm/yr since stage 3 at ∼53 Myr. In model M4-1, the slab with a 100 km long TOC keeps flat subduction for more than 80 Myr after the initial subduction ( Figure 8a). The eclogitized oceanic crust is trapped in the flat subduction channel with an increasing temperature, resulting in a decreasing viscosity and density ( Figure S1 in Supporting Information S1). No evident sagging or delamination is observed at ∼53 Myr in this case. In model M4-2, the slab with a 200 km long TOC shows evident sagging ∼600 km landward away from the trench at the time of 51.6 Myr (Figure 8b). A further increase of the length of the TOC to 300 km or 400 km leads to slab sagging ∼400 km landward away from the trench (Figures 8c and 8d). The similar sagging locations in models M4-3 and M4-4 with 300 and 400 km long TOC indicate that the effect of the increase in the TOC length on the sagging position of the slab is limited.
The slabs with different lengths but the same cross-sectional area of TOC show similar subduction mode between flat subduction and steep subduction (Arrial & Billen, 2013). Our numerical results further indicate that the slab with same cross-sectional area of TOC shows similar transition process. For example, models M4-2 and M6-3 with same cross-sectional area of 4,000 km 2 but different thicknesses of 15-km and 20-km TOC show similar sagging location and start time of delamination as shown in Figure 8b. However, in the case of a 25-km thick TOC, the slab quickly evolves into steep subduction, regardless of the cross-sectional area of the TOC. Under a high subduction rate (e.g., 5 cm/yr), the slab rolls back from the front since 21 Myr and breaks off during the process of rolling back, which is similar to model M3-2 in Figure 6. In model M6-6, under the same subduction rate, the slab with a 25-km thick TOC sags (Figures 9b-9d) about 12 Myr earlier than the results of model M3-4 with a similar cross-sectional area but thinner TOC (Figures 5b and 5c). In addition, the location of slab sagging in M6-6 ( Figure 9d) is also much closer to the trench than that in M3-4 ( Figure 5c). Moreover, the thicker TOC in model M6-6 (25 km thick) also leads to an earlier (∼5 Myr) occurrence of eclogitization in the middle of the   Note. Stage 1: Before the TOC subducts beneath trench; Stage 2: After the subduction of the TOC subducts beneath trench and before the TOC starts eclogitization; Stage 3: After the initial eclogitization of TOC and before its complete eclogitization; Stage 4: After the complete eclogitization of the TOC.

Table 3
Description of the Numerical Models  TOC than models with a 20-km-thick TOC (Figure 9b). This is caused by the vertical temperature gradient and the deformation of the TOC during the flat subduction process, resulting in a faster increase of the temperature in the middle part than the front part. According to the vertical temperature gradients in the oceanic lithosphere, a thicker TOC indicates a higher average temperature of the TOC. Compared with a 20 km thick TOC, the temperature at the bottom of the 25 km thick TOC is easier to reach beyond the eclogitization condition. Together with the uplift in the middle of the TOC, the TOC gets eclogitized as shown in Figure 9b.

The Formation of a Long-Lasting Flat Subduction
Among the numerical results, some cases indicate that the oceanic slab with a 100 long TOC can maintain a long-lasting flat subduction under a certain range of subduction rate. For example, delamination occurs in model M5-1 with a TOC of 100 km length and 20 km thickness under a low subduction rate of 1.5 cm/yr (Figure 10a). In model M5-4, with a higher subduction rate of 2.2 cm/yr, the eclogitized TOC subducts further inland and results in the occurrence (48.1 Myr) of slab rollback (Figure 10b). However, in models M4-1 and M5-3, once the eclogitized TOC subducts ∼500 km landward under a relatively low subduction rate of 1.8-2.0 cm/yr since stage 3, the slab keeps a long-term flat subduction for more than 80 Myr with episodic break-offs at the front of the slab (Figures 10c and 10d). Instead of simply subducting landward with the oceanic lithospheric mantle, the eclogitized TOC maintains beneath the continental lithosphere (see Movie S2).
However, the flat slab with a 200 km or longer TOC cannot maintain such a long-lasting flat subduction under the range of the subduction rate considered in our numerical models (≤5 cm/yr). These results suggest that a long-lasting flat subduction requires both a relatively small cross-sectional area of the TOC and a subduction rate within an appropriate range. Compared with the results for slabs with TOC of a normal-thick oceanic crust, the enclosure of flat subduction channel at the frontal part of the slab with a 100 km long TOC stops the advance of oceanic crust and sediment. The sagging at the middle part of the flat slab leads to the uplift at the frontal part of the slab, which stops the subduction of materials enclosed in the flat subduction channel (Figure 10c). Besides, the accumulation of the eclogitized oceanic crust owing to its increased density leaves the upper part of the flat subduction channel for the partial molten sediment with a relatively low density and viscosity. When the partial molten sediment (dark and light yellow layers) reaches the front end of flat subduction channel, it occupies the upper part of the flat subduction channel (Figures 10c and 10d). Therefore, the eclogitized oceanic crust cannot break through the enclosure of the lithospheric mantle and partial molten sediment.
In addition, the dense frontal slab breaks off episodically from the flat slab, which significantly decreases the driving forces on the front of the slab. Without breaking off, this dense frontal part of the slab will cause rollback. Meanwhile, this long-lasting flat subduction requires an appropriate range of subduction rates, which is not sufficiently high to induce rollback and not sufficiently low to induce delamination. The appropriate subduction rate and coupling between partial molten sediment and eclogitized oceanic crust provide enough resisting forces against the additional driving force of the 100 km long eclogitized oceanic crust that tends to drag the slab downwards. When the subduction rate is beyond the range which maintaining flat subduction, the slab evolves into steep subduction. After the slab has formed a long-lasting flat subduction stage, even a change in the subduction rate within a certain range will not induce the transition to the steep subduction. For example, when the subduction rate changes within 2.4 cm/yr to 1.5 cm/yr after 63 Myr in model M5-3 (Figure 10a), the slab can still keep flat subduction.

Discussion
Our numerical results show that the flat-to steep-slab subduction transformation process is affected by both the subduction rate and the size of TOC. A high subduction rate is favorable for slab rollback, whereas a low subduction rate together with a large-scale TOC is more likely to result in slab delamination. These results provide the bases for us to analyze the flat subduction and steepening processes of the PPP beneath the SCB in the Mesozoic.

How the TOC Affects the Fate of Flat Subduction?
The presence of the TOC exerts opposing effects on the slab before and after the eclogitization phase change. The uneclogitized oceanic crust is beneficial to flat subduction (Dai et al., 2020;Huangfu et al., 2016), since it provides an additional positive buoyancy force through its lower density. Compared with the slab with a normal-thick oceanic crust, the slab with a TOC shows a delayed rollback starting time (Figures 6 and 7). In contrast, after eclogitization, the TOC promotes interplate decoupling with its density increasing to 3,600 kg/m 3 . Flat subduction begins to transform to steep subduction when the negative buoyancy force associated with eclogitization breaks the balance between the driving forces and the resisting forces. The transition to steep subduction occurs in the form of the slab delamination or rolling back, depending on the location where the negative buoyancy force is applied. When the negative buoyancy force provided by eclogitized oceanic crust or the dense frontal part of the slab is applied on the front of the slab (Figure 6a), it promotes slab rollback (Figure 6b). In contrast, when the negative buoyancy force of eclogitized oceanic crust locates within the flat subduction channel, delamination occurs (Figures 4c-4e and 5c-5e). This additional negative buoyancy force is strongly affected by the size of the TOC, including the length and the thickness.
The slabs with different size of TOC as a seamount (100 km), seamount group or oceanic plateau can all induce the sagging and delamination. Under the same subduction rate, the size of TOC determines the location of slab sagging and even the transition mode. When the additional driving force locates in the middle of flat subduction channel, it dominates the delamination process and affects the location of slab sagging. The additional driving force from gravity of a longer or thicker eclogitized TOC is more likely to break the force balance against the additional driving force earlier during the flat subduction stage (Figure 1). As a result, the location of slab sagging is closer to the trench. On the other hand, once the additional driving force locates at the frontal part of the flat slab, it promotes the slab rollback. For example, under a 3.5 cm/yr subduction rate in stage 2 and a 2 cm/yr subduction rate since stage 3, the 100-km long eclogitized TOC leads to slab rollback and the 400-km long TOC leads to slab delamination (red circles in Figure 11).

The Effects of Subduction Rate on the Steepening Process
The effects of subduction rate are different between cases with and without a TOC. For the M2-2 with a normal-thick oceanic crust under a 3 cm/yr subduction rate, slab rolls back in about 10 Myr once the slab subducts beyond the innermost location of the flat subduction channel (Figure 7). The higher the subduction rate, the earlier the slab rolls back. A decrease of the subduction rate cannot turn rollback into delamination, but only delays the starting time of the rollback process.
For the cases with a TOC, the location of the TOC at a certain time is determined by the subduction rate. If the eclogitized TOC has subducted beyond the innermost location of the flat subduction channel, the balance between the driving forces and the resisting forces is broken irreversibly. The oceanic lithosphere rolls back even if the subduction rate decreases to a low value of 0.5 cm/yr at a later stage. A relative high subduction rate (>2.4 cm/yr) since the initial subduction of the TOC (≥100 km long and 20 km thick) leads to slab rollback. Whether the slab break-off occurs in the rollback process depends on both the subduction rate and the size of TOC. The warming up of slab is mainly controlled by the thermal conductivity, the surrounding temperature and the duration of conduction. A higher subduction rate results in a shorter duration of warming up and thus a cooler slab. A cooler temperature helps to maintain the strength of the slab and prevent the slab from breaking off. In addition, the size of TOC also plays a role in slab break off. The larger the size, the larger the driving forces induced by eclogitized TOC is, which increases the likelihood of slab break-off.
With a sufficiently low subduction rate (≤2.4 cm/yr) since the initial subduction of the TOC (≥100 km long and 20 km thick), especially after the initial eclogitization in the TOC (the beginning of stage 3), instead of purely rolling back, the slab delaminates ( Figure 11). A continuous low subduction rate since stage 2 results in the ponding of eclogitized TOC in the flat subduction channel. This eclogitized TOC imposes a sufficient negative buoyancy force on the slab, and therefore the slab delaminates. Within the range of subduction rate at which delamination occurs, a higher subduction rate since the initial subduction of the TOC pushes the TOC further inland, resulting in the sagging position of the TOC further away from the trench. This is evident by comparing M3-3 ( Figure 4c) and M3-4 ( Figure 5c).
In summary, a low subduction rate after the initial subduction of TOC is necessary to induce sagging in the middle of the slab and therefore the following slab delamination ( Figure 11). The range of the subduction rate that can induce the sagging in the middle of the slab varies with the length of the TOC. The range shrinks to lower values as the length of TOC decreases ( Figure 11).

The Dynamic Factors for the PPP Delamination Beneath the SCB in the Mesozoic
Previous studies based on the geological records (Z. Z. Li et al., 2012; and numerical simulations (Dai et al., 2020) suggest that the orogen and magmatism in the Mesozoic SCB is strongly affected by the PPP delamination process. To figure out the triggering factor for the PPP delamination, we consider the kinematic conditions in the Mesozoic SCB. The Mesozoic movements of PPP and SCB are mainly constrained by the global plate motion model (Müller et al., 2016) and also incorporate findings from other paleogeographic and paleomagnetic researches (Bao et al., 2013;Isozaki et al., 2010;S. Liu et al., 2017;Maruyama et al., 1997;Yang et al., 2016). The plate movement reconstruction (see Movie S3) shows the animation of principal motion for the oceanic and the continental lithospheres ( Figure 12). For the oceanic lithosphere (PPP), it kept a ∼7 cm/yr subduction rate along the NW direction in the Mesozoic before ca. 150 Ma. After that, the PPP moved NNW-ward with a gradually increasing subduction rate, resulting in a 14 ± 5 cm/yr subduction rate along the NW direction after ca. 125 ± 5 Ma. For the continental lithosphere (SCB), it moved along the SEE direction with a rate of 4.5 cm/yr before ca. 195 ± 5 Ma. Between ca. 195 ± 5 Ma and ca. 150 ± 10 Ma the overall motion of SCB gradually turned to the northwest direction. After that, the SCB shows a general southward motion with a rate of 2.5 ± 0.5 cm/yr.
Before ca. 195 ± 5 Ma, as a result of the SE-ward motion of the SCB, the relative subduction rate of the PPP beneath the SCB is higher than 5 cm/yr. The previous study interprets that the arrival of an oceanic plateau at the continental arc induced the flat subduction, which caused the cartonward migration of the orogen from ca. 250 Ma (Z. . However, our numerical results indicate that the flat subduction of PPP can also be induced with a normal-thick oceanic crust (Figure 7a). In addition, if the TOC in the PPP subducted beneath the SCB since ca. 250 Ma, at the ca. 190 Ma the TOC would have subducted 800 km inland from the trench and induce slab rollback rather than slab delamination under the high subduction rate during that period. As a result, the arrival of the TOC at the continental arc is expected to be later than the previous estimation.
It is speculated that the trenchward motion of the SCB before ca. 195 ± 5 Ma increased the pressures on the PPP upper surface (Figure 13a). Together with the strong interplate viscous coupling and additional buoyancy from TOC, these factors promoted the flat subduction of the PPP (Huangfu et al., 2016). The strong interplate coupling between PPP and SCB resulted in the compressive deformation in the SCB as a fold-and-thrust belt (Figure 13a). We also notice that the crucial motion change of the SCB and the PPP corresponds well to the starting time of the slab delamination. The retreat time of the SCB in ca. 195 ± 5 Ma (Figures 12b and 13b) is close to the sagging time of the PPP (Dai et al., 2020). The retreat of the continental lithosphere since ca. 195 ± 5 Ma decreased the relative subduction rate of the PPP beneath the SCB (Müller et al., 2016) and remove the tectonic pressure from continental lithosphere on the PPP upper surface (Figure 13b). According to our numerical results, this motion provided a favorable kinematic condition for the occurrence of delamination. The decreased subduction rate and eclogitization of TOC induced slab sagging and delamination.
The sagging and delamination of the PPP is accompanied by the mantle upwelling (Figures 13b and 13c), which induced the large-scale anorogenic magmatism evolved since ca. 190 Ma (H. Z. Li & Li, 2007;Z. Li et al., 2012). The rollback of the PPP since ca. 150 Ma induced the SE-ward younging trend of magmatism in the coastal area of the SCB (Figure 13d). The range of the subduction rate between ca. 195 ± 5 Ma and ca. 150 Ma is inferred from numerical results and plate motion reconstruction. Our numerical results indicate that with a 400-km long and 20-km thick TOC, the maximum subduction rate for delamination to occur is ∼2.4 cm/yr (Figure 11). With a smaller TOC, the maximum subduction rate for delamination to occur decreases ( Figure 11). Otherwise, under a higher subduction rate, the PPP would rollback directly without delamination.
Based on the relative motions of the SCB and the PPP (Müller et al., 2016), the estimated minimum subduction rate after ca. 195 ± 5 Ma is 1.6 cm/yr. This minimum value is estimated from the extreme case with an 8 cm/yr northward motion of PPP and a 4 cm/yr NW-ward motion of SCB (Figure 12c), leading to a 1.6 cm/yr velocity along the NW direction. From our numerical results and the recorded plate motions, we infer that the relative subduction rate is 2.0 ± 0.4 cm/yr between ca. 195 ± 5 Ma and ca. 150 Ma.   slabs with different TOC lengths result in a similar size of the concave-up flexural extent with the same sagging location (Figure 8). This extent explains the shallow marine subsidence in the SCB during Early Jurassic (Z. . As a result, the length of TOC cannot be simply speculated from the numerical results. To obtain the precise size of TOC, three-dimensional numerical simulations and a detailed plate motion model consistent with geological records are required.

Model Limitations
One of the limitations of the present study is that the effects of water on material properties and on the eclogitization process are not taken into account. The water content in the rock, the dehydration process and the transfer of the water can all affect the material properties (including density and viscosity), the eclogitization and the melting process (Z. Li, 2020;. The water content in the oceanic crust determines the temperature and the pressure conditions for eclogitization to occur. The time that eclogitization occurs varies with water contents, and this affects the transition process from flat to steep subduction. In addition, the dehydrated water from the slab can weaken the overriding continental lithosphere and change the ratio of partial melting. To improve accuracy, quantitative simulations on the aqueous fluid activity are required in further studies. Another limitation of our models is the setting of plate motion. The subduction rate in the present study is simplified to change instantly at the time of the three signature events and remains constant over each stage, whereas a more smooth and continuous variation of the subduction rate is more reasonable. In addition, the three-dimensional motions of the SCB and the PPP were more complex than what we estimated (Jiang et al., 2013;Koppers et al., 2001;Lü et al., 2005;W. Sun et al., 2007). Therefore, three-dimensional numerical simulations with a more smooth variation of the subduction rate are worth pursuing to study the impact of complex plate motions on the transition of the subduction process in future studies.

Conclusion
Based on the numerical simulations, we study the effects of the TOC size and the subduction rate on the fate of flat subduction. In terms of flat subduction, the results reveal that with a low initial subduction angle of 12°, an 80 Ma old oceanic slab with a normal-thick oceanic crust can subduct flatly beneath the 120-km thick continental lithosphere. The dependence of the transition mode on the TOC size and the subduction rate suggests that a higher subduction rate or a smaller cross-sectional area of TOC is favorable for slab rollback, while slab delamination requires both a low subduction rate and the presence of TOC. In the process of slab delamination, with the same location of slab sagging, the continental lithosphere shows similar concave-up flexural extent. Parameters, including the size of TOC and the subduction rate, can both affect the sagging location. In addition, under a certain range of subduction rate, an 80 Ma old slab with a small size area of TOC (e.g., 100 km long and 20 km thick) can maintain a long-lasting flat subduction. However, this flat subduction stage is unstable, a slight variation of the subduction rate will lead to the transition to steep subduction, through either delamination or direct rollback.
From the comparison between the numerical results and global plate motion models, we suggest that the arrival time of the TOC to the continental arc of the SCB should be later than ca. 250 Ma. The precise arrival time to the continental arc depends on the thickness of TOC and the precise subduction rate. After that, the motion change of the SCB from SE-ward to NW-ward since ca. 195 ± 5 Ma decreased the relative subduction rate to 2.0 ± 0.4 cm/ yr between ca. 195 ± 5 Ma and ca. 150 Ma, which induced the PPP delamination. To obtain the more accurate size of TOC, three-dimensional numerical simulations and a detailed plate motion model are required.

Data Availability Statement
The figures of numerical models are produced by the Matlab of MathWorks. All related data are provided in Zenodo (https://doi.org/10.5281/zenodo.7047321).