Imaging of the Upper Mantle Beneath Southeast Asia: Constrained by Teleseismic P-Wave Tomography

Abstract

It is of great significance to construct a three-dimensional underground velocity model for the study of geodynamics and tectonic evolution. Southeast Asia has attracted much attention due to its complex structural features. In this paper, we collected relative travel time residuals data for 394 stations distributed in Southeast Asia from 2006 to 2019, and 14,011 seismic events were obtained. Then, teleseismic tomography was applied by using relative travel time residuals data to invert the velocity where the fast marching method (FMM) and subspace method were used for every iteration. A novel 3D P-wave velocity model beneath Southeast Asia down to 720 km was obtained using this approach. The tomographic results suggest that the southeastern Tibetan Plateau, the Philippines, Sumatra, and Java, and the deep part of Borneo exhibit high velocity anomalies, while low velocity anomalies were found in the deep part of the South China Sea (SCS) basin and in the shallow part of Borneo and areas near the subduction zone. High velocity anomalies can be correlated to subduction plates and stable land masses, while low velocity anomalies can be correlated to island arcs and upwelling of mantle material caused by subduction plates. We found a southward subducting high velocity body in the Nansha Trough, which was presumed to be a remnant of the subduction of the Dangerous Grounds into Borneo. It is further inferred that the Nansha Trough and the Dangerous Grounds belong to the same tectonic unit. According to the tomographic images, a high velocity body is located in the deep underground of Indochina–Natuna Island–Borneo–Palawan, depth range from 240 km to 660 km. The location of the high velocity body is consistent with the distribution range of the ophiolite belt, so we speculate that the high velocity body is the remnant of thee Proto-South China Sea (PSCS) and Paleo-Tethys. This paper conjectures that the PSCS was the southern branch of Paleo-Tethys and the gateway between Paleo-Tethys and the Paleo-Pacific Ocean. Due to the squeeze of the Australian plate, PSCS closed from west to east in a scissor style, and was eventually extinct under Borneo.

1. Introduction

The Proto-South China Sea (PSCS) refers to the ocean that was once located in the southeastern part of China in the area where the South China Sea (SCS) now stands [1,2,3,4]. In a narrower sense, it was believed that the PSCS was a subduction of the SCS region at the Eurasian continental margin from Izonaki to Eurasia, and that the formation of the PSCS was due to the rollback of the Izonaki subduction plate to form a spreading environment and the pulling apart of the post-arc basin [5,6,7]. Numerous scholars have carried out extensive research in the PSCS remnant zone. A set of ophiolite remained from the PSCS has been identified using geological data, which is distributed in an arc along Natuna Island–Kalimantan Island–Palawan Island [8]; others have applied the tomography method to the region, and it was inferred that the PSCS remained in the deep mantle of the southern Philippines region [9], and it has also been suggested that the remnant of the PSCS presents in the mantle of the SCS basin [10,11]. Taylor and Hayes suggested that the PSCS existed in the Middle Triassic to Late Cretaceous Tertiary based on magnetic anomalies, paleo-magnetic, and paleontological evidence [3]. While Liu defined the PSCS as an eastern extension of the Middle Tethys based on regional geological and geophysical data and inferred a Late Permian date for its emergence [1]; Wu inferred the earliest existence of the PSCS in the Triassic period based on Mesozoic stratigraphic and fossil characteristics in the southern SCS [12]; Zhou et al. inferred that the PSCS has existed since the Late Jurassic–Early Cretaceous based on deep-sea radiolarians in the serpentine suite of the Natuna–Kalimantan–Palawan distribution, and that it may be a back-arc sea formed by subduction of the Pacific plate [13]; Metcalfe argued that the PSCS was formed when the Cretaceous East Asian continental margins pulled apart by studying the tectonic shelf of the Epiphytic Sundaland [14]; Lu studied the influence of Tethys and Paleo-Pacific on the PSCS through regional geological, paleontological, gravity, and magnetic anomalies [15]; numerous scholars have also determined the existence of the PSCS from 70 Ma to 10 Ma through plate reconstruction methods [7,16].

Considering the location of the PSCS, which is the junction of the Paleo-Pacific Ocean and Tethys, the Tethys tectonic domain has been extensively studied by previous researchers [17,18,19,20,21], and some scholars had also commented on the Paleo-Pacific of tectonic evolution for reconstruction [22,23,24,25,26]. However, due to the limitations of data and research methods as well as the complexity of the tectonic evolution of the Paleo-Pacific and Tethys tectonic domains, the relationship between the PSCS and Tethys tectonic domains and the Paleo-Pacific has not been clarified clearly.

Seismic tomography has the advantage of a large detection depth and has been widely applied to study the deep structure and tectonic evolution [27]. Extensive excellent work has been carried out by previous researchers in SE Asia [28,29,30,31,32]. Based on previous studies, this paper used the latest seismic data that can be collected in Southeast Asia to study its velocity structure. Relative travel time residuals data from 2006 to 2019 were collected for 394 stations in the study area. The teleseismic tomography method was used to invert the P-wave velocity structure of the study area. Finally, a three-dimensional P-wave velocity model was obtained, which was beneath SE Asia. Combined with the results of previous studies, the distribution range of the PSCS, the location of subduction extinction, and its relationship with the Paleo-Pacific and Tethys tectonic domain can be summarized as follows. Figure 1 shows a sketch of the tectonic and topography in Southeast Asia.

2. Data

The seismic data used in this paper were primarily derived from the Incorporated Research Institutions for Seismology (IRIS) [33] and International Seismic Center (ISC) [34]. The research range is (15° S–30° N, 90° E–130° E). To satisfy the requirement of the subsequent inversion calculation, data were collected according to the following criteria: (1) The seismic events used for the inversion have magnitudes greater than 5.0; (2) The distances from the station to hypocenter were between 30° and 90°, minimizing the impact of the complex structure of the lower mantle and core mantle boundary [35]; 3. Seismic events used for inversion must be recorded by at least five stations simultaneously. The final collection was 14,011 seismic events recorded from 394 seismic stations in the study area, and 180,225 useful P-wave relative travel time residuals were obtained for inversion. Compared to previous work, this study added data from one new station at the Dangerous Grounds. Therefore, the results of this study are more reliable than those of the previous work. The distribution of stations in the study area and the distribution of seismic sources used for inversion calculation are shown in Figure 2a,b, respectively.

To obtain the relative travel time residual data for inversion calculations, it is necessary to pick up the P-wave first arrival time of the collected raw seismic waveform data. Manually picking each waveform data will be a very time-consuming work. Seismologists have proposed a number of proven methods for picking up the P-wave first arrival time based on the principle of waveform correlation [36,37,38]. A graphical interface program for picking up relative travel time residuals was developed by Zhang [38], which combined two methods of picking up the relative travel time residuals together. This program greatly improved the efficiency of picking up relative travel time residuals. The work of picking up relative travel time residuals in this paper benefited from this program.

The process of picking up relative travel time residuals based on the above method is shown in Figure 3 [39]. First, the appropriate recoding of teleseismic events were downloaded and filtered according to the study acquisition requirements. Second, the raw data was band-pass filtered before the pickup work. The frequency range of 0.01–1 Hz was chosen to avoid filtering out the relatively high-frequency information of the teleseismic events. The first arrival time of the P-wave can be better preserved. You can also adjust the range based on the sampling rate of the data. The next step is to select the theoretical travel time provided by the ak135 model to align the waveform data with each other, and perform the inter-correlation process on its theoretical travel time. We used the adaptive stacking (astack) method [40] and multi-channel cross correlation (MCC) method [36] to superimpose the waveform data, and then compared the difference between the residuals obtained by the two picking methods. Waveforms with a P-wave first arrival error greater than 0.2 s and a shape obviously different will be deleted. Repeat the above steps until all the picked results are within the error range. This process of picking up relative travel time residuals guarantees accuracy as much as possible [41]. Figure 4 shows the waveform pickup window for a magnitude 6.1 earthquake on 22 June 2008.

3. Methodology

The principle of the teleseismic tomography imaging method is shown in Figure 5. First, we need to dissect the subsurface structure beneath the receiver array in the study area into a gridded model. This model will be placed into a globally symmetric reference model [42,43,44]. Figure 5 illustrates the principle of teleseismic tomography. As shown in Figure 5, a local gridded model is located beneath the receiver array, which contains 3D variations in velocity. Regions outside the local gridded model are assumed to be spherically symmetric. The advantage of this pattern is that the travel time from the seismic source to the local gridded model boundary can be predicted quickly. The travel time inside the local gridded model is calculated by ray-tracing techniques for laterally heterogeneous media [45]. The premise of the above principle is that the lateral variations outside the local gridded model will not significantly affect the relative travel time residuals. In other words, all relative travel time residuals at the stations are then assumed to be only caused by the local anomalies.

The teleseismic tomography process can be divided into model parameterization, forward step, and inversion step. The following description of the inversion scheme is described according to these steps.

3.1. Model Parameterization

Model parameterization refers to the process of using a series of constant or changing nodes or blocks to characterize velocity variations. The following form is used to formulate the slowness field within a 3D spherically symmetric grid of nodes [46]:

$$s\left(r,\theta ,\varphi \right)={\sum }_{i,j,k=1}^4{S}_{ijk}{C}_i\left(R\right)C_j\left(\Theta \right)C_k\left(\Phi \right)$$

(1)

where $S_{ijk}$ represents the slowness at the nodes of the $4\times 4\times 4$ grid surrounding the point $\left(r,\theta ,\varphi \right)$. $C_i\left(R\right)$, $C_j\left(\Theta \right)$, and $C_k\left(\Phi \right)$ are known as the cardinal splines, and $R$, $\Theta$, and $\Phi$ are local coordinates of $r$, $\theta$, and $\varphi$ [46].

Conventional cubic spline interpolation has been widely applied in travel time tomography [47,48,49,50]. Spline passing through node values is essential for traditional cubic spline interpolation. The form of interpolation function cubic B-splines used in this paper was similar to Equation (1). This was locally supported and did not necessarily pass through the node values.

3.2. Forward Step

The next step of teleseismic tomography scheme is the forward step. The aim of the forward step in the inverse problem is calculating the theoretical travel times of waves.

In travel time inversion, the conventional method of calculating source-receiver travel time was the ray-tracing technique. More recently, the eikonal equation has been widely applied in travel time inversion [47]. The following eikonal equation governed the path of seismic waves in the local gridded model:

$$\left|{\nabla }_{x}T\right|=s\left(x\right)$$

(2)

In Equation (2), $T$ represents travel time and $s\left(x\right)$ represents slowness as a function of position $x$. The core of the FMM method is to solve Equation (2) with an upwind entropy that satisfies finite differences. In this paper, we used the following form, which has been employed by other authors [48,49]:

$${\left[\begin{array}{c}max{(D_a^{-r}T,-D_b^{+r}T,0)}^2\\ max{(D_c^{-\theta }T,-D_d^{+\theta }T,0)}^2\\ max{(D_e^{-\varphi }T,-D_f^{+\varphi }T,0)}^2\end{array}\right]}_{i,j,k}^{\frac{1}{2}}=S_{i,j,k}$$

(3)

In Equation (3), $\left(i,j,k\right)$ denote grid increment variables in $\left(r,\theta ,\varphi \right)$. The integer variable a, b, c, d, e, f is used to determine the order of accuracy of the upwind finite-difference operator used in each of the six portions. If $\left(r,\theta ,\varphi \right)$ represents a spherical coordinate system, the first two gradient operators in each of the three directions can be described by the following form:

$$\begin{array}{l}{D}_1^{-r}{T}_{i,j,k}=\frac{T_{i,j,k}-T_{i-1,j,k}}{\delta r}\\ D_2^{-r}{T}_{i,j,k}=\frac{3T_{i,j,k}-4T_{i-1,j,k}+T_{i-2,j,k}}{2\delta r}\\ D_1^{-\theta }{T}_{i,j,k}=\frac{T_{i,j,k}-T_{i,j-1,k}}{r\delta \theta }\\ D_2^{-\theta }{T}_{i,j,k}=\frac{3T_{i,j,k}-4T_{i,j-1,k}+T_{i,j-2,k}}{2r\delta \theta }\\ D_1^{-\varphi }{T}_{i,j,k}=\frac{T_{i,j,k}-T_{i,j,k-1}}{rcos\theta \delta \varphi }\\ D_2^{-\varphi }{T}_{i,j,k}=\frac{3T_{i,j,k}-4T_{i,j,k-1}+T_{i,j,k-2}}{2rcos\theta \delta \varphi }\end{array}$$

(4)

The operator used in Equation (4) depends on the availability of upwind travel times and the max order allowed. By default, we used a mixed order form which used empirical $D_2$ operators but reverted to $D_1$ if $T_{i-2,j,k}$, $T_{i,j-2,k}$, or $T_{i,j,k-2}$ was unavailable [47]. Equation (3) describes the finite difference form for calculating the travel time associated with a specified grid node. However, the order of updating grid nodes have not been specified. To satisfy causality, the values of T should change from small to large. To satisfy the above rule, FMM systematically constructs the travel time in a downwind manner based on known values in the upwind direction. This is called the narrowband approach. Figure 6 illustrates the principle of the multi-stage FMM approach. A detailed explanation of this method can be found in [49].

3.3. Inversion Step

The essence of seismic tomography inversion is to solve the optimization problem. The subspace method can be minimized simultaneously along the search direction, and across a subspace of the model space, which is illustrated in Figure 7. The following form is the most common objection function used in inversion:

$$S(m)=(g(m)-d_{obs}{)}^T{C}_d^{-1}(g(m)-d_{obs})+\epsilon {(m-m_0)}^T{C}_m^{-1}(m-m_0)$$

(5)

There is a basic assumption before using gradient methods that when $S\left(m\right)$ is smooth enough, the following form can be obtained:

$$S\left(m+\delta m\right)\approx S\left(m\right)+\widehat{\gamma }\delta m+\frac{1}{2}\delta m^T\widehat{H}\delta m$$

(6)

In the above form, $\delta m$ indicates a disturbance to the current model, and $\widehat{\gamma }=\frac{\partial S}{\partial m}$ indicates the gradient vector; $\widehat{H}=\frac{{\partial }^{2}S}{\partial m^2}$ indicates the Hessian matrix. Combining the above partial differences gives the following forms [50]:

$$\widehat{\gamma }=G^T{C}_d^{-1}\left[g\left(m\right)-d_{obs}\right]+\epsilon C_m^{-1}\left(m-m_0\right)$$

(7)

$$\widehat{H}=G^T{C}_d^{-11}G+{\nabla }_m{G}^T{C}_d^{-1}\left[g\left(m\right)-d_{obs}\right]+\epsilon C_m^{-1}$$

(8)

In Equations (7) and (8) $G=\frac{\partial g}{\partial m}$ denotes the $F{r}^{\prime }echet$ matrix of partial derivatives. For the case of constant slowness blocks $G=\left[l_{ij}\right]$, $l_{ij}$ indicates the ray segment length of the $i^{th}$ ray in the $j^{th}$ block.

The general theory for the normal subspace inversion method has been explained by some scholars [51,52,53]. The subspace method limits the minimization of the quadratic approximation of S(m) to an n-dimensional subspace of model space, so that the disturbance $\delta m$ occurs in the space spanned by a set of n M-dimensional basis vectors $\left\{a^j\right\}$:

$$\delta m=\sum _{j=1}^n{\mu }_j{a}^j=A\mu$$

(9)

In Equation (9), ${\mu }_j$ is the component defining the length of the corresponding vector $a^j$ that minimizes the quadratic form of $S\left(m\right)$ in the space. Equation (6) takes the following form:

$$S\left(m+\delta m\right)=S\left(m\right)+\sum _{j=1}^n{\mu }_j{\widehat{\gamma }}^T{a}^j+\frac{1}{2}\sum _{j=1}^n\sum _{k=1}^n{\mu }_j{\mu }_k[a^k{]}^T\widehat{H}[a^j]$$

(10)

and searches for the minimum values of $S$ considering $\mu$:

$$\frac{\partial S\left(m\right)}{\partial {\mu }_q}={\widehat{\gamma }}^T{a}^q+\sum _{k=1}^n{\mu }_k[a^k{]}^T\widehat{H}[a^j]=0$$

(11)

When $q=1,\dots ,n$. the following form can be obtained from Equation (9):

$$\mu =-{[A^T\widehat{H}A]}^{-1}{A}^T\widehat{\gamma }$$

(12)

where $\delta m=A\mu$, and $\widehat{H}=G^T{C}_d^{-1}G+\epsilon C_m^{-1}$, and the following form obtained:

$$\delta m=-A{[A^T\left(G^T{C}_d^{-1}G+\epsilon C_m^{-1}\right)A]}^{-1}{A}^T\widehat{\gamma }$$

(13)

The quantities $A$, $\widehat{\gamma }$, and $G$ are re-evaluated between successive iterations. The construction of the basis vector $\left\{a^j\right\}$ for the subspace approach is typically done with the steepest upward vector and rate of change in the model space [51,52,53,54].

The subspace method has been applied to invert the velocity, interface structure, and hypocenter location using local earthquake and artificial source travel times [53,54,55,56]. The inversion scheme has been successfully applied in northern Tasmania and the Murray Basin [57,58]

4. Crust Correction

When using teleseismic tomography to study the velocity structure of the upper mantle, the complex crustal structure may affect the tomographic results of the upper mantle [60,61,62,63,64]. To ensure the accuracy of the upper mantle velocity imaging results, crust correction is needed to eliminate the influence of the crust. Many scholars have carried out research in this area [65,66,67]. In this paper, the method proposed by Jiang [67] is used for crustal correction of teleseismic imaging. The idea of this method is to obtain the relative travel time residuals inside the crust before the inversion calculation, and then process the relative travel time residuals of the entire ray.

The main steps of crustal correction are:

(1) Select the appropriate 1-D crustal model and 3-D crustal model. The 1-D crustal model in this article was obtained from AK135 [68], and the 3-D crustal model was the CRUST1.0 [69];

(2) Calculate the travel time of the ray in the 1-D model and the 3-D model of the crust respectively, and then calculate the relative travel time residual in the crust, the calculation formula, takes the following form:

$$\delta T_{crust}=T_{3D}-T_{1D}={(\sum _l\frac{h_l}{\mathit{cos}{\theta }_l\times V_l})}_{3D}-{(\sum _k\frac{h_k}{\mathit{cos}{\theta }_k\times V_k})}_{1D}$$

(14)

where $\delta T_{crust}$ denotes the relative travel time residuals in the crust. $T_{3D}$ and $T_{1D}$ represent the travel time in the crust, respectively. $h_l$ and $h_k$ represent the thickness of each layer in 1-D and 3-D models, respectively. ${\theta }_l$ and ${\theta }_k$ represent the incident angle of the ray at each interface in 1-D and 3-D models, respectively.

(3) Subtract the relative travel time residual in the crust from the relative travel time residual of the ray, and use it as the inversion data of teleseismic tomography.

$$t=\left(T_{obs}-T_{cal}\right)-\delta T_{crust}$$

(15)

where $t$ indicates the relative travel time residuals of rays. $T_{obs}$ and $T_{cal}$ represent the observed travel time and theoretical travel time, respectively.

Figure 8 shows the distribution of the average relative residuals at each station before and after crust correction and the change of average relative travel time residuals. From Figure 8a,b, the difference between the average of each station before and after crustal correction to travel time residuals is very small, ranging from −2.4 s to 2.4 s after crustal correction. Figure 8c reveals that the variation of the average travel time residual of each station ranged from −0.09 s to 0.09 s.

5. Resolution Test

Before displaying the tomographic imaging, we used the checkerboard test method to check the stability and resolution of the model [70]. The specific steps of the checkerboard test are as follows: (1) Establish an initial model 1 and a model 2 based on the initial model 1 with positive and negative disturbance velocity anomalies added; (2) Use model 2 as the initial model to calculate the theoretical relative travel time residuals of each ray of the observation system consisting of stations and events in the study area; (3) Use model 1 as the initial model and use the theoretical travel time obtained in step (2) to calculate residuals used in the inversion step, and model 3 was obtained through inversion calculations. Comparing the velocity anomaly images of model 2 and model 3, the smallest anomaly scale that could be recovered in model 3 was the resolution of the model [41].

During the checkerboard test in this study, a plus or minus 0.3 km/s perturbation was added to the theoretical model. In the vertical direction, the resolutions of the 40 km, 60 km, and 90 km grids were tested, and the results suggested that when the vertical grid was 90 km, the model had the best recovery effect. In the horizontal direction, the checkerboard test results revealed that the most optimal grid space was 1.65° × 1.65°.

Figure 9 shows the checkerboard test results in the horizontal, latitude, and longitude directions. Figure 9a(I) is a horizontal slice of the input model at 100 km, and Figure 9a(II), Figure 9a(IV), and Figure 9a(VI) are the slices of the output model at 100 km, 400 km, and 600 km, respectively. According to the results of the checkerboard test, when the depth was 100 km, the recovery effect of the output model was good in most areas of the study area. Only the southwestern area of the study area and the central area of the SCS basin showed poor recovery effects. The main reason is that these two areas are both marine areas so there are fewer stations and sparse rays. Figure 9a(III), Figure 9a(V), and Figure 9a(VII) show the distribution of ray paths at 100 km, 400 km, and 600 km, respectively. When the depth increased, the recovery effect of the South China Sea basin gradually improved. Figure 9b show the results of the checkerboard test in the latitude direction. Figure 9b(I) is a slice of the input checkerboard model along the 0° direction, and Figure 9b(II) and Figure 9b(III) are slices of the output model along the 0° and 12° N, respectively. It can be seen that the recovery effect of the slice in the 0° direction was very good. The recovery effect of the velocity anomalies in the area of 110°E–115°E on the slice along 12°N was not very ideal, while the recovery effect of the velocity anomalies in other areas was very good. Figure 9c shows the results of the checkerboard test in the longitude direction. Figure 9c(I) shows the slice of the input model along 104°E. Figure 9c(II) and Figure 9c(III) are the slices of the output model along 104°E and 114°E, respectively. It can be seen that the velocity anomalies of the slice along 104°E recovered well. The overall recovery effect of thee slice along 112° E was very good. The shallow recovery effect of 10°–15° was slightly worse, which was caused by the scarcity of stations in the ocean area. In general, the effect of the chessboard test is relatively satisfactory. The velocity model obtained in this study had a large scale, so it was mainly used to study the large-scale structural evolution. According to the results of the checkerboard test, the model used in this study had a denser ray distribution in the upper mantle and mantle transition zone, so the inversion results in this depth range are more reliable.

In order to apply the above-mentioned tomographic method to obtain the 3D velocity model beneath SE Asia, the local model beneath the receiver array was first parameterized. The structure under 394 stations in the study area was divided into 27,225 nodes. During the inversion, the vertical grid space was 90 km, and the grid space in the longitude and latitude directions was 1.65° (~160 km). The above parameters were selected based on the checkerboard test results. The depth of the local gridded model in the vertical direction was down to 720 km, the range in the longitude direction was 40°, and the range in the latitude direction was 45°. The reference velocity model was given by ak135. After six iterations of inversion, an inversion result was obtained. The observed relative travel time residuals and the theoretical relative travel time residuals were highly matched. Figure 10 shows the distribution of the observed relative travel time residuals and the post-inversion relative travel time residuals, respectively. Observed data and the inversion results were normally distributed, indicating that the initial model evolved along the direction of the fitting anomaly after the inversion calculation. This further verifies the reliability of the inversion results, indicating that the inversion results are reasonable.

6. Results

Map views of the P-wave tomography at nine depths are shown in Figure 11. Overall, the velocity anomalies at different depths had good correspondence with the associated tectonic units. When the depth was smaller than 200 km, it can be seen that the southeastern Tibet Plateau, the South China Plate, the Philippines, the Indochina Peninsula and Java showed obvious high-velocity anomalies. The Dangerous Grounds and the western region of Borneo also exhibited high-velocity anomalies, while low-velocity anomalies were mainly located in areas such as the Malay Peninsula, Sumatra, and Borneo. As the depth increased, the anomaly distribution characteristics of the study area changed. When the depth was between 250 km and 630 km, the Tibet Plateau and the South China Plate and the Philippines still exhibited the characteristics of high-velocity anomalies. The area from the Andaman Sea along Sumatra to Java changed from low-velocity to high-velocity, just like a high-velocity anomaly covered with a low-velocity anomaly. This phenomenon was also found in the previous study [33]. It is worth noting that the high-velocity anomaly area of the Indochina Peninsula had become larger and showed a trend of spreading from north to south. The area outside the western part of Borneo also changed from a low-velocity anomaly in the shallow to a high-velocity anomaly. The high-velocity anomalies in these two areas tended to be integrated. The overall route is shown in Figure 11, from the depth of 200 km–630 km with red arrows. This area was the focus of this paper. In addition, the velocity anomalies in the SCS basin also changed from high-velocity anomalies to low-velocity anomalies. The characteristics of velocity anomalies below 700 km were like those at 630 km. High-velocity anomalies mainly occur in the subduction zone, shallow ocean basin, and continental area, while low-velocity anomalies occur in the deep ocean basin and associate locations with high-velocity anomalies. Previous studies have suggested that the low-velocity area may reflect the upwelling of mantle material, and was mainly caused by plate subduction and lower mantle collapse [71,72].

Figure 12 and Figure 13 show the slices along the east–west and north–south directions, respectively. As can be seen in the east–west slices (Figure 12), the shallow part of the Andaman Sea and the SCS basin exhibited high-velocity anomalies. The deep areas of the Philippines, Sumatra, and Java exhibited high-velocity anomalies. The low-velocity anomalies in the shallow areas covered the high-velocity anomalies in the deep area. As can be seen in the north–south slices (Figure 13), the Indochina, the South China Plate, and the shallower part of the SCS basin exhibited high-velocity anomalies, while in the subduction zone such as Sumatra and Java, it is characterized by shallow low-velocity anomalies covering deep high-velocity anomalies. The slices results are more supportive of the possibility that the low-velocity region near the subduction zone was due to material upwelling caused by plate subduction into the mantle [71,72].

We then compared our 3D P-wave velocity model in the study area with previous studies on seismic tomography imaging conducted in the same region [73,74,75,76,77] or on a global scale [78,79,80,81]. Our results were basically consistent with the previous results: high-velocity anomalies were located below the plate boundary and island arcs such as the southeastern Tibet Plateau, South China Plate, the Philippines, Sumatra, and Java and low-velocity anomalies were located in the deep part of the SCS basin and the shallow part of Borneo. What is worth noting in the results of this study is the high-velocity anomaly discovered on the Indochina–Natuna Island–Borneo–Palawan, which may indicate the location of the sutures of the PSCS, which will be discussed in detail in the next section.

7. Discussion

The three-dimensional velocity model obtained in this paper can provide more evidence for studying the tectonic evolution of the PSCS. According to previous studies on the PSCS, some scholars believed that the stitching position of PSCS was in central Borneo [8,82]. Some authors have found that the magnetic anomaly of the satellite shows that the lower part of Natuna Island exhibited the same low magnetic anomaly as the ocean basin area, and then inferred that it may be the location of the PSCS suture in the area [83]. While agreeing with the above conclusion, this paper argues that the PSCS suture should extend northeastward to Palawan, and extend northwestward to Indochina. The ophiolites on the Lupar line in the southern SCS are distributed along the line from Naturna Island–Northern Borneo–Southern Palawan [84,85]. Ophiolite, as a representative of ancient lithosphere fragments, indicates that the above path may be a suture zone of the PSCS. According to the characteristics of the velocity anomalies at different depths, during the demise of the PSCS, some plates may remain in the mantle, which provides guidance for us to restore the position of the PSCS (Figure 10).

According to the position of the PSCS sutures delineated in the horizontal slices, the velocity anomalies on the east–west slices through the sutures were analyzed. In the slice along 16°N (Figure 12a), a high-velocity body was located in the area of 96°E–106°E and a depth range of 240 km–660 km (red solid area in Figure 12a). The high-velocity body was covered by the low velocity body, which was also revealed in a previous tomographic model [31,86]. We inferred the high-velocity body to a remnant slab of Paleo-Tethys. In the slice along 12°N, in the area of 104°–106°, in the depth range of 500 km–720 km, a similar high-velocity body was also found (red solid area in Figure 12b), where we inferred that the high-velocity body was also the remnant slab of Paleo-Tethys. In the 8°N slice, only a small area of high-velocity body was located below 600 km near 105° (Figure 12c red dotted area), where the high-velocity body may also be the remains of Paleo-Tethys. The smaller scale may be caused by the upwelling and melting of the hot material caused by the subduction of the plate. Coincidentally, in the slice along 4°N, a high-velocity body subduction eastward was found within the range of 108°E–117°E (red solid area in Figure 12d), and the deepest point reached below 660 km. The low-velocity anomalies above Borneo were caused by the upwelling of the material due to the subduction of the plate into the mantle. A high-velocity body subduction eastward was also found along the 0°N slice, in the range 106°E–116°E (Figure 12e), with a low-velocity anomaly caused by the upwelling of the mantle material above it. Some scholars concluded that the PSCS subducted mainly southward beneath Borneo during 45 Ma and 16 Ma [6,87]. Previous researchers have also found a high-velocity anomaly dipping S–SE below northern Borneo [30,88]. Therefore, we infer that the high-velocity body is a remnant slab of the PSCS subducting below Borneo. However, we believe that the slab-pull of the PSCS was not the main factor for the opening of the SCS. A high-velocity body subduction eastward (black solid area in Figure 12) was located on the westernmost side of the east–west slice. According to its position distribution, it should indicate the subduction zone of New-Tethys. It is worth noting that on the 4°N section, another high-velocity body (yellow solid area in Figure 12d) was located on the west side of the high-velocity belt representing the subduction of New-Tethys, and it is speculated that this high-velocity body may be the subduction of Mid-Tethys. Based on the above results, we have revised the positions of Mid-Tethys and Paleo-Tethys proposed by our predecessors (Figure 15).

Similarly, based on the location of the PSCS suture delineated in the horizontal slices, the velocity on the north–south slice passed through the suture location anomalies for analysis. In the slice along 104°E (Figure 13a), a high-velocity body was located within the depth range of 240 km–660 km in the area of 15°N–20°N (red solid area in Figure 13a). This high-velocity body coincides with the position of the high-velocity body found in Figure 12a. It can be inferred that the high-velocity body is a remnant plate of Paleo-Tethys, and the high-velocity body is covered by the low velocity body. This phenomenon has been revealed by previous study [31,86]. In the slice along 108°E, due to less rays in this region only in the range 10°N–13°N, the high-velocity bodies with suspected Paleo-Tethys remnant slabs were located below 600 km (the area circled by the red dashed line in Figure 13b). On the 112°E slice, in the range of 4°N–2°S, a high-velocity body was found diving southward to below Borneo (Figure 13c red dashed line circled the area) up to 600 km deep, and the high-velocity body was presumed to be a remnant plate of the PSCS. On the slice along 116°E, within the range of 7°N–4°N, a high-velocity body subduction southward (Figure 13d black solid line circled the area) was found within the range, and the deepest point reached below 220 km. Some scholars believed that the slab-pull of the PSCS will cause the southward movement of the Dangerous Grounds and the thinning of the crust [89]. Therefore, we speculate that the high-velocity body is a remnant plate subduction from the Nansha Trough below Borneo. Based on the velocity anomaly signature, we inferred that the Nansha Trough and the Dangerous Grounds were the same tectonic unit. The northward high-velocity body was located at the southernmost end of the north–south slices, and it can be inferred that this high-velocity zone indicated the New-Tethys subduction. It is worth noting that the shallower Borneo exhibited low velocity anomalies that were caused by material upwelling due to the penetration of the Mid-Tethys and the PSCS to the mantle (Figure 14).

Based on the above viewpoint, combined with the velocity anomalies characteristic of the Indochina–Natuna Island–Borneo-Palawan, we concluded that the PSCS from Indochina is in a nearly north–south direction, passing northeast of Natuna Island to western Borneo and turning to a nearly east–west direction. It finally turned north–eastward in northern Borneo, reaching Palawan and connecting with the suture belt of the Paleo-Pacific Ocean. The exact path is shown by the solid line with red arrows in Figure 11.

In summary, our results suggest that the PSCS is both the southern branch of Paleo-Tethys and a channel between Paleo-Tethys and the Paleo-Pacific Ocean. The location of the PSCS suture should be distributed along the Indochina–Natuna Island–Borneo–Palawan. The specific route is shown in the solid black line with arrows in Figure 14 and Figure 15. We speculate that the northeast end of the PSCS is connected with the subduction zone of the Paleo-Pacific near Palawan Island. Due to the northward backlog of the Australian plate, the PSCS finally closed from west to east in a scissor-like fashion, and died out under Borneo. As shown in Figure 15, the blue dotted line indicates the location of the sutures of Paleo-Tethys–PSCS–Paleo-Pacific.

8. Conclusions

We presented here a three-dimensional P-wave velocity beneath SE Asia using extensive relative travel time residual data obtained from ISC (International Seismic Center) [34] and IRIS (Incorporated Research Institutions for Seismology) [33]. Although the stations are mainly located in regions such as land and island arcs, the abundance of seismic events ensured the reliability of the results of this study.

Our results revealed the velocity structure of the upper mantle beneath the SE Asia. It was found that high-velocity anomalies were mainly located in regions such as the southeastern Tibetan Plateau, the Philippines, Sumatra, Java, and deep Borneo, and low velocity anomalies were mainly located in the deep parts of the SCS basin and in the shallower areas near the subduction zone such as Borneo. High- velocity anomalies mainly suggest subduction plates and stable continental blocks, while low- velocity anomalies suggest island arcs and upwelling of mantle material caused by subduction plates. We found a southward subducting high-velocity body in the Nansha Trough, which is presumed to be a remnant of the subduction of the Dangerous Grounds into Borneo. We further inferred that the Nansha Trough and the Dangerous Grounds belong to the same tectonic unit.

A high-velocity body was located in the deep beneath the Indochina–Natuna Island–Borneo–Palawan. The location of the high-velocity body is consistent with the distribution range of the ophiolite belt, so it is speculated that the high-velocity body indicating the PSCS and the Paleo-Tethys. Residue. PSCS was the southern branch of Mid-Tethys, and it was also the channel connecting Paleo-Tethys and the Paleo-Pacific Ocean. Due to the squeeze of the Australian plate, PSCS closed from west to east in a scissor-style, and eventually died out under Borneo. We believe that the slab-pull of the PSCS was not the controlling factor for the opening of the SCS.