Figure 1 . Structure-tectonic maps of the Celebes Sea and surrounding area. (a) Location of the study area—Celebes Sea (CS)—within the wider Southeast Asian region, including the Sulu Sea (SS), the South China Sea (SCS), the north Banda Sea (NBS) and the south Banda Sea (SBS). Thick lines indicate subduction zones; triangles on the lines indicate the down-dip direction. Plate motions are indicated by the arrows (Mustafar et al., 2017; Kreemer et al., 2014). Red triangles indicate volcanoes (Global Volcanism Program, 2013). (b) Regional topographic map. The yellow hollow circles are heat flow measurements (Fuchs et al., 2021). The pink dashed lines show the major magnetic anomalies (Weissel, 1980). The irregular shape with transparent colour gradient represents the subducted plate; colours indicate depth. The dotted lines with triangles indicate trench retreat. Events in the EHB catalogue (Engdahl et al., 2020) are represented by black circles.
Model setup
We used the geodynamic code ASPECT (Advanced Solver for Problems in Earth’s ConvecTion, version 2.3.0-pre) to setup a series of 2D subduction models. The initial model consists of oceanic plate and continental plate on either side (Figure 2a). According to the evolutionary history of NSSZ, the age of oceanic plate is taken to be 40 Ma. The oceanic lithosphere geotherm corresponds to a half-space cooling model for a seafloor age of 40 Ma. The thickness of oceanic plate is therefore set to 82 km. The oceanic crust consists of two layers and is 10 km thick in total. The composition for the mantle part of the oceanic lithosphere and for the asthenosphere are taken to be identical. The density of them is both set to 3300 kg/m3. The continental plate consists of upper crust, lower crust and lithospheric mantle. The continental crust is set to 30 km based on the current crustal thickness of Sulawesi (Fauzi et al., 2021). Because the initial lithospheric thickness, density of lithospheric mantle and ocean-continent transition angle have never been definitely settled, we vary the values of these three parameters to investigate their effect on the spontaneous subduction (Table 1). Temperature structure of the continental plate is controlled by the lithospheric thickness and is defined by a steady state profile from 0 ℃ at the surface to 1300℃ at the base of lithosphere (Figure 2a).
Spontaneous and induced subductions are modelled by controlling the boundary conditions and density differences between oceanic plate and asthenosphere. The spontaneous models (Experiment 1~24) have all free slip boundaries (Figure 2a) and the Gabbro-Eclogite phase transition for the oceanic crust was implemented. Density of oceanic crust increases from 3000 kg/m3 to 3450 kg/m3 at pressures of ~1.9 GPa (~60 km depth). In the induced model (Experiment 25~30), horizontal far field forces are represented by specifying a convergence velocity (V in) on one side of the overriding plate; the Gabbro-Eclogite phase transition for the oceanic crust was not executed.