3 Model Results

3.1 Parameters that affect spontaneous subduction

Since the subduction initiation of NSSZ occurred at 5–9 Ma and the maximum depth of the subducted slab is about 260 km at present. The simulation ends when the subducted slab reaches a depth of 260 km or the time reaches 10 myr. If the oceanic slab does not exceed a depth of 150 km within 10 myr, we denot it by “No” (Figur 2b and Table 1). If the maximum depth of the subducted slab is deeper than 150 km but shallower than 260 km (Figure 2b), we list the depth. For the experiments where the slab reached 260 km depth within 10 myr, we list this time (Figure 2b).

Experiments show transition from fast to slow subduction with increasing thickness of the continental lithosphere. Spontaneous subduction initiation is almost impossible to occur with 70 km or more thick continental plate (Experiment 9–12, 21–24). The effect of density contrast is opposite to that of the thickness of the continental plate. A larger density contrast promotes spontaneous subduction initiation. The ocean-continent boundary angle has a nonlinear effect on the spontaneous subduction initiation. Either too large (60°) or too small (15°) angle are all adverse to spontaneous subduction initiation. The thickness and density of the continental lithospheric mantle and the ocean-continent boundary angle all affect whether subduction initiation occurs spontaneously. It suggests that the force due to the density difference between the continental and oceanic plate is the main force leading to the spontaneous subduction initiation.

3.2 Trench retreat velocity in spontaneous and induced model experiments

Trench location and speed are inferred from the location and speed of the edge of the continental crust that is the farthest forward in the model. In all spontaneous and induced model experiments, trench retreat speed increases and then decreases during the first 1–2 myr (Figure 2c). At the beginning of the simulation period, the density difference between the continental and oceanic crusts gives rise to lateral stresses that push the continental crust towards the ocean; as the difference between the stresses decreases, trench retreat speed decreases gradually. This process is active in approximately 2 myr.

After the first 2 myr in the spontaneous model (Figure 2c), subduction depth increases, the trench retreat velocity increases correspondingly. The increasing begins slowly, then more and more rapidly. It is because as the oceanic plate sinks, there are more and more eclogite-face oceanic crust to increase the negative buoyancy of the subducted slab. In the induced model, the trench retreat speed increasing slowly, then becoming almost constant. The final speed is slightly greater than the specified velocity V _in.

3.3 Deformation and heat flow of the overriding plate

The deformation of the overriding plate varies in the different subduction models. In the spontaneous subduction model, the rollback of the subducted plate obviously causes stretching of the overriding plate. The crustal thickness decreased from 30 km to about 23 km (eg. Experiment 2, 3). In the induced subduction model, the overriding plate thrusts over the subducted plate, and the overriding plate is slightly thinned (Experiment 25–30). The latest published geological profile of the North Arm of Sulawesi indicates lithospheric significant extension between the Late Miocene and Pliocene (Nugraha et al., 2022, Advokaat et al., 2017), which is consistent with the results of the spontaneous subduction model.

The heat flow of the overriding plate is closely related to the lithospheric thickness either in spontaneous or induced model experiments. It is opposite to the continental lithospheric thickness. The values of heat flow listed in Table 1 is the calculated surface heat flow which lies at a distance of 100 km above the top of the subducted slab. It corresponds to the location of the maximum heat flow (68 mW/m²) on the North Arm of Sulawesi (Figrue 1b). In all the experiments, the model with 75 km thick continental lithosphere has the closest value to that observation. The heat flow values in models have a 70 km or thinner lithosphere are all greater than the maximum observed values.

4 Discussion

4.1 Implications for transition of passive margins into subduction zones

The mechanism of subduction initiation at passive margins have been studied by several numerical simulations (Faccenna et al., 1999; Goren et al., 2008; Marques et al., 2014; Turcotte et al., 1977; Zhou et al., 2020). Our intention is not repeat other studies and prove the feasibility of spontaneous subduction initiation at passive margins, but to estimate the subduction mechanism through the constraints of the geological case–NSSZ. The spontaneous and induced subduction models represent the extreme cases of subduction driven by horizontal far-field forces and of subduction driven by the negative buoyancy of the down-going plate itself. Examples of completely spontaneous or induced mechanisms are difficult to reproduce the evolution of NSSZ—subduction depth of 260 km after approximately 5–9 myr, 68 mW/m² heat flow of the overriding plate and obvious extended continental crust. Because the mainly force that cause the oceanic plate to bend and sink in spontaneous subduction initiation is lateral stress, which is caused by the density difference between the continental plate and the oceanic plate (Nikolaeva et al., 2010). A thick crust and hot lithosphere are necessary to produce the sufficiently large lateral stress (Nikolaeva et al., 2010). According to the present crustal structure and heat flow of North Arm of Sulawesi, spontaneous subduction initiation is unlike to occur in the NSSZ. Therefore, an external force is necessary to push the continental plate towards the oceanic plate in the early stage of subduction initiation. For the NSSZ, the external force may originate mainly from the collision of the surrounding plates. Because since the Cenozoic, Sulawesi has been at the center of the convergence of surrounding plates (Hall, 2012). This induced mechanism dominates the early stage of subduction initiation, which takes about 2 myr. Once the oceanic crust reaches the depths of ~60 km where eclogite transformation can occur, the main force controlling the subduction begin to change from the horizontal external force to the vertical negative buoyancy of the subducted slab. The subsequent subduction, trench retreat and extension of continental crust are all mainly depend on the negative buoyancy of the subducted slab, as exhibitions in spontaneous model experiments (Table1 and Figure 2). In addition, adjacent areas in the same tectonic setting can also provide insight into the subduction process.

The Sula Deep, Tolo Trough and Cotobato Trench located near the North Sulawesi Trench (Figure 3a, d) are considered to examples of three stages in the subduction initiation. Sula Deep and Tolo Trough is suggested to record the earliest stages and second stage in subduction development, respectively (Hall, 2019). Cotobato Trench is interpreted as the third stage where a seismically defined slab can be recognized (Hall, 2019). The North Sulawesi Trench represents the final stage, which indicates that subduction initiation is end and a mature subduction zone has been formed. From the topography of these trenches (Figure 3), we can see that each trench has a different landward slope. The areas (Sula Deep and Tolo Trough) in early stage of subduction initiation, has a very steep landward slope (Figure 3b, c). As the subduction development (Cotobato Trench and North Sulawesi Trench), the landward slope changes from steep to gentle (Figure 3e, f).

There are obvious differences in landward slope angles vary in different mechanisms models. In the spontaneous numerical model, the slope gradually smooth with the subduction development (Figure 3g). While the slope in induced model is steeper than that in spontaneous model and the slope can be maintained at a steeper state under the push of far field forces (Figure 3h). Combining the steeper slope in the Sula Deep and Tolo Trough and gentle slope in Cotobato Trench and North Sulawesi Trench, we could also suggest the subduction slab controlled the retreat of the trench.