The 3D Density Structure of the South China Sea Based on Wavelet Multi-Scale Analysis of Gravity Data and Its Tectonic Implications

Abstract

Due to its unique geographical location and complex geological evolution processes, the South China Sea has been a focus of extensive research. Previous studies on the density structure of the South China Sea mostly focused on 2D density structures, with relatively limited research on 3D density structures. A comprehensive study is still needed to refine the expansion mechanism and tectonic evolution of the South China Sea. In this study, we utilized wavelet multi-scale analysis of gravity data to obtain a 3D density model of the South China Sea and discussed its tectonic evolution from the pattern of density anomalies. The inversion results show that (1) the expansion of the South China Sea caused the typical thin oceanic crust and parts of the continent–ocean transition zone may fracture due to the expansion; (2) the low-density anomaly in the upper mantle of Luzon Island may indicate partial melting or the upwelling of asthenosphere materials; and (3) the expansion of the South China Sea is influenced by multiple plate forces and uneven forces affect the distribution of high-density anomalies in the upper mantle.

1. Introduction

The South China Sea (SCS), one of China’s three marginal seas, is located at the intersection of the Eurasian Plate, the Philippine Sea Plate (PSP), the Pacific Plate, and the Indo-Australian Plate, making it the largest marginal sea in the western Pacific Ocean. Due to its unique geographical location, the SCS has developed landforms such as “trench–arc–backarc basin systems” and become a highly popular research area. Its topographic map is shown in Figure 1 [1]. However, precisely because the SCS is situated at the intersection of these plates, determining its origin and evolution has always been a complex challenge. Clarifying the expansion mechanism and tectonic evolution processes of the SCS will help us better understand the dynamic mechanisms of long-term tectonic complexity.

Previous researchers have proposed multiple models regarding the driving force of the opening of the SCS, such as the collision-induced lateral extrusion and escape model, the mantle plume model, and the back-arc spreading model [4,5,6,7,8,9,10]. Zhang et al. [11] divided the basement of the SCS into seven tectonic units and concluded that the formation of a unified basement in the SCS region was primarily attributed to two causes: the collision and amalgamation of the Cathaysian Block with the Indochina Block during the Triassic period and the collision and amalgamation of the Cathaysian Block with the Nansha Block during the Early Cretaceous. By analyzing various rift basins in the SCS, Wang et al. [12] suggested that research on the SCS should not solely focus on the sea itself but should also encompass its surrounding areas, considering them as an integrated system. These studies have provided valuable insights into our understanding of the tectonic evolution of the SCS.

Although the tectonic evolution of the SCS has been extensively explored, numerous controversies remain, such as the expansion mechanism of the SCS and the tectonic evolution of adjacent regions like the MT and Luzon Island, which require further research to support existing models or propose new ones. Some scholars have utilized gravity inversions to obtain the Moho topography model of the SCS [13,14], thereby analyzing the crust–mantle structure of the SCS. However, the deformation mechanism of the Moho topography still requires further exploration. By combining water depth data and sediment thickness information, Li et al. [15] conducted a 3D gravity inversion in the northern SCS to derive its density anomaly model, providing a different perspective for analyzing the tectonics of the SCS. Yet, geological structures exhibit not only strong lateral heterogeneity but also complex vertical characteristics due to varying tectonic sequences and the interlaced appearance of various complex geological units. Thus, the analysis of deep-seated structures is equally essential.

The vertical resolution of gravity signals is weak, and the inversion solution is not unique. Without enough constraints, “skin effects” may arise. Some studies employ depth weighting functions to impose constraints [16], but this can lead to issues that the inversion results are overly influenced by the weightings [2,17]. Some scholars use seismic tomography to obtain the 2D density structure of some profiles [18]. Some scholars proposed a 3D gravity inversion method to obtain a 3D density model [19], but this 3D inversion method relies on a priori models. Compared with these methods, the wavelet multi-scale analysis method does not rely on a priori models and weighting functions and can effectively improve the vertical resolution and robustness of inversion model.

In this study, we first utilized the wavelet multi-scale analysis method to separate the Bouguer gravity anomaly signal into wavelet approximation and details. By employing the radial logarithmic power spectrum approach, we calculated the average field source depth. Then, we performed layered density inversion based on the depth values and wavelet details. We analyzed the patterns of density anomalies in the SCS and discussed the expansion and tectonic evolution of the SCS.

2. Methods

2.1. Improved Parker Method

The Bouguer gravity anomaly is affected by signal perturbations from uneven thickness and sedimentary layers, which have a significant impact on density inversion. Therefore, it is necessary to deduct the gravity anomaly of the sedimentary layer. Since both deep and shallow field sources generate gravity signals across the entire frequency range, it is difficult to effectively separate the gravity anomaly of the sedimentary layer solely through filtering [20]. Using the improved Parker method, we calculate the gravity anomaly effect of the sedimentary layer [21,22]. The formulas for marine sediment thickness and sediment density adopt the quadratic function fitted by Yu et al. [23], as shown below:

$$∆\rho \left(z\right)=\left\{\begin{array}{c}-1.0395+0.4280\times t-0.0546\times t^2, 0<t\le 3\\ -0.4076+0.0626\times t-0.0029\times t^2, t>3\end{array}\right.$$

(1)

where $∆\rho \left(z\right)$ represents the density of marine sediment and $t=\tilde{z}-{\tilde{z}}_1$ stands for the thickness of the sedimentary layer. The expression for $∆\rho \left(z\right)$ is given:

$$∆{\rho }_s\left(\tilde{x},\tilde{y},\tilde{z}\right)=\sum _{m=0}^M{a}_m\left(\tilde{x},\tilde{y}\right)·{\left(\tilde{z}-{\tilde{z}}_1\right)}^m=\sum _m^M{b}_m\left(\tilde{x},\tilde{y}\right)·{\tilde{z}}^m$$

(2)

where $a_m\left(\tilde{x},\tilde{y}\right)$ and $b_m\left(\tilde{x},\tilde{y}\right)$ represent the coefficients of an m-degree polynomial. The forward formula is as follows [22]:

$$F\left[g^{s,c,t}\left(x,y,z_0\right)\right]=2\pi G{e}^{-\left|k\right|z_0}\sum _{n=0}^N{\displaystyle \frac{{\left(-\left|k\right|\right)}^{n-1}}{n!}}F\left[∆{\rho }_{s,c,t}\left(\tilde{x},\tilde{y},\tilde{z}\right)·\left({\tilde{z}}_2^n\left(\tilde{x},\tilde{y}\right)-{\tilde{z}}_1^n\left(\tilde{x},\tilde{y}\right)\right)\right]$$

(3)

where the positive direction of the z-axis is downward, $F$ represents the 2D Fourier transform, $z_0$ denotes the height of the computation plane, $g^{s,c,t}$ stands for the gravity value at a specific point ($x,y,z_0$), k is the wavenumber, $N$ represents the order of computation, and $∆{\rho }_{s,c,t}$ signifies the anomalous density of the integral element ($\tilde{x},\tilde{y},\tilde{z}$). The depths of the upper and lower interfaces at ($\tilde{x},\tilde{y}$) are indicated by ${\tilde{z}}_1$ and ${\tilde{z}}_2$, respectively. By substituting Equation (2) into Equation (3), the gravity disturbance caused by the sedimentary layer can be obtained:

$$F\left[g^{s,c,t}\left(x,y,z_0\right)\right]=2\pi G{e}^{-\left|k\right|z_0}\sum _{n=0}^N{\displaystyle \frac{{\left(-\left|k\right|\right)}^n}{n!}}F\left[\sum _{m=0}^M{b}_m·{\displaystyle \frac{{\tilde{z}}_2^{n+m+1}\left(\tilde{x},\tilde{y}\right)-{\tilde{z}}_1^{n+m+1}\left(\tilde{x},\tilde{y}\right)}{n+m+1}}\right]$$

(4)

2.2. Wavelet Multi-Scale Separation

Gravity anomalies are a comprehensive reflection of internal material signals. To analyze the material from different depths, it is necessary to separate and process the gravity anomaly signals. Currently, commonly used separation methods are analytic continuation, trend analysis, wavelet separation, etc. Among them, wavelet multi-scale separation is one of the most effective methods due to its advantages such as not relying on depth-weighted functions [2]. Therefore, after sedimentary layer correction processing, this study adopted the 2D wavelet multi-scale analysis method [24] to separate the gravity anomalies $\Delta \mathrm{g}\left(\phi ,\lambda \right)$ in the SCS and adjacent regions, obtaining low-frequency signals and high-frequency signals of different orders, as shown in Equation (5):

$$\Delta \mathrm{g}\left(\phi ,\lambda \right)=U_J\left(\phi ,\lambda \right)+\sum _{j=1}^J{D}_j\left(\phi ,\lambda \right)$$

(5)

where $U_j\left(\phi ,\lambda \right)$ represents the wavelet approximation, which signifies the low-frequency signal caused by deep-seated geological bodies; $D_j\left(\phi ,\lambda \right)$ represents the wavelet detail, which signifies the high-frequency signal caused by shallow geological bodies.

2.3. Estimation of the Average Field Source Depth

Wavelet detail and wavelet approximation in different frequency bands are caused by geological bodies from different depths. Therefore, after decomposing gravity anomalies, the key issue is to determine the depth of the gravity anomaly.

The wavelet multi-scale separation result $D_j(j=\mathrm{1,2},\dots ,J)$ can be expressed as follows [25]:

$$D_j\left(\phi ,\lambda \right)=\sum _{\phi }\sum _{\lambda }{A}_k{e}^{i2\pi \left(k_{\phi }\phi +k_{\lambda }\lambda \right)}{e}^{2\pi kH}$$

(6)

where $A_k$ is the amplitude, $k=\sqrt{k_{\phi }^2+k_{\lambda }^2}$ is the wave number, and $H$ is the average field source depth of $D_j$. $A_k$ is given as follows:

$$A_k=\sum _{\phi }\sum _{\lambda }{D_j\left(\phi ,\lambda \right)e}^{-i2\pi \left(k_{\phi }\phi +k_{\lambda }\lambda \right)}{e}^{\pm 2\pi kH}$$

(7)

when $H=0$,

$${{(A}_k)}_0=\sum _{\phi }\sum _{\lambda }{D_j\left(\phi ,\lambda \right)e}^{-i2\pi \left(k_{\phi }\phi +k_{\lambda }\lambda \right)}$$

(8)

We substitute Equation (8) into Equation (7) and simplify to obtain Equation (9):

$$A_k={{(A}_k)}_0{e}^{\pm 2\pi kH}$$

(9)

Since the power spectrum is the square of the amplitude, we square both sides of the equation and take their logarithms; then, we obtain the following formula:

$$P_k={{(P}_k)}_0{e}^{\pm 4\pi kH}$$

(10)

$$\ln P_k=\ln {{(P}_k)}_0\pm 4\pi kH$$

(11)

Hence, the average depth $H$ of $D_j$ can be estimated using the formula below:

$$H={\displaystyle \frac{∆\ln P_k}{4\pi ∆k}}$$

(12)

2.4. Layered Density Inversion

Based on the average field source depth H, the internal structure of the SCS is stratified and each layer is gridded using Tesseroid elements. According to the available research [26], the relationship between the gravity anomaly $D_j\left(\phi , \lambda \right)$ and the density anomaly $\Delta {\rho }_j\left({\phi }_0,{\lambda }_0\right)$ of each Tesseroid element is expressed as follows:

$$D_j\left(\phi , \lambda \right)=\sum _{{\phi }_0}\sum _{{\lambda }_0}G\Delta {\rho }_j\left({\phi }_0,{\lambda }_0\right)\Delta r_j\Delta \phi \Delta \lambda \left[L_{000}+{\displaystyle \frac{1}{24}}\left(L_{200}\Delta r^2+L_{020}\Delta {\phi }^2+L_{002}\Delta {\lambda }^2\right)\right]$$

(13)

where ${\phi }_0$ and ${\lambda }_0$ represent the central latitude and longitude of the Tesseroid element, respectively; $\Delta r_j$ is the thickness; $\Delta \phi$ and $\Delta \lambda$ are the intervals of latitude and longitude; and $L_{000}$, $L_{200}$, $L_{020}$, and $L_{002}$ are the Taylor expansion coefficients. Their expressions are as follows:

$$L_{000}={\displaystyle \frac{r_0^2\left(R-r_0\cos {\psi }_0\right)\cos {\phi }_0}{l_0^3}}$$

(14)

$$\begin{array}{c}{L}_{200}={\displaystyle \frac{\mathrm{Rcos} {\phi }_0}{l_0^3}}\{2-{\displaystyle \frac{3r_0}{l_0^2}}[5r_0-(2R+3r_0\mathrm{c}\mathrm{o}\mathrm{s} {\Psi }_0)\mathrm{c}\mathrm{o}\mathrm{s} {\Psi }_0]\\ +{\displaystyle \frac{15r_0^3}{l_0^4}}{\mathrm{s}\mathrm{i}\mathrm{n}}^2 {\psi }_0(r_0-R\mathrm{c}\mathrm{o}\mathrm{s} {\psi }_0)\}\end{array}$$

(15)

$$\begin{array}{c}{L}_{020}={\left({\displaystyle \frac{r_0}{l_0}}\right)}^3\mathrm{c}\mathrm{o}\mathrm{s} \phi (1-2{\mathrm{s}\mathrm{i}\mathrm{n}}^2 {\phi }_0)\mathrm{c}\mathrm{o}\mathrm{s} \delta \lambda +{\displaystyle \frac{r_0^2}{l_0^5}}\{-R(R^2+r_0^2)\mathrm{c}\mathrm{o}\mathrm{s} {\phi }_0\\ +r_0\mathrm{s}\mathrm{i}\mathrm{n} \phi [-R{r}_0(\mathrm{s}\mathrm{i}\mathrm{n} \phi \mathrm{c}\mathrm{o}\mathrm{s} {\phi }_0-\mathrm{c}\mathrm{o}\mathrm{s} \phi \mathrm{s}\mathrm{i}\mathrm{n} {\phi }_0\mathrm{c}\mathrm{o}\mathrm{s} \delta \lambda )+\\ \mathrm{s}\mathrm{i}\mathrm{n} {\phi }_0\mathrm{c}\mathrm{o}\mathrm{s} {\phi }_0(2R^2+4r_0^2-3R{r}_0\mathrm{s}\mathrm{i}\mathrm{n} \phi \mathrm{s}\mathrm{i}\mathrm{n} {\phi }_0)]+\\ r_0^2\cos \phi \cos \delta \lambda \left(1-2{\sin }^2 {\phi }_0\right)\times \left[r_0+\mathrm{Rcos} \phi \cos {\phi }_0\cos \delta \lambda \right]\\ +R{r}_0^2\cos \phi \sin {\phi }_0\cos {\phi }_0\cos \delta \lambda \\ \times [3\mathrm{s}\mathrm{i}\mathrm{n} \phi \mathrm{c}\mathrm{o}\mathrm{s} {\phi }_0-4\mathrm{c}\mathrm{o}\mathrm{s} \phi \mathrm{s}\mathrm{i}\mathrm{n} {\phi }_0\mathrm{c}\mathrm{o}\mathrm{s} \delta \lambda ]\}\\ +{\displaystyle \frac{5R{r}_0^3}{l_0^7}}\{-R(R^2+r_0^2)\mathrm{s}\mathrm{i}\mathrm{n} {\phi }_0+r_0^2\mathrm{c}\mathrm{o}\mathrm{s} \phi \mathrm{s}\mathrm{i}\mathrm{n} {\phi }_0\mathrm{c}\mathrm{o}\mathrm{s} {\phi }_0\mathrm{c}\mathrm{o}\mathrm{s} \delta \lambda \times \\ \left(r_0+\mathrm{Rcos} \phi \cos {\phi }_0\cos \delta \lambda \right)\\ +r_0\sin \phi \left[2R^2-r_0^2-R{r}_0\cos {\psi }_0+{\sin }^2 {\phi }_0\right]\\ \times (R^2+2r_0^2-R{r}_0\mathrm{s}\mathrm{i}\mathrm{n} \phi \mathrm{s}\mathrm{i}\mathrm{n} {\phi }_0)\}\\ \times \left(\sin \phi \cos {\phi }_0-\cos \phi \sin {\phi }_0\cos \delta \lambda \right) \end{array}$$

(16)

$$\begin{array}{c}{L}_{002}={\left({\displaystyle \frac{r_0}{l_0}}\right)}^3\cos \phi {\cos }^2 {\phi }_0\times \{\mathrm{c}\mathrm{o}\mathrm{s} \delta \lambda -{\displaystyle \frac{3R}{l_0^2}}[2r_0\cos \phi \cos {\phi }_0{\sin }^2 \delta \lambda \\ +(R-r_0\mathrm{c}\mathrm{o}\mathrm{s} {\psi }_0)\mathrm{c}\mathrm{o}\mathrm{s} \delta \lambda ]\\ +{\displaystyle \frac{15R^2{r}_0}{l_0^4}}\mathrm{c}\mathrm{o}\mathrm{s} \phi \mathrm{c}\mathrm{o}\mathrm{s} {\phi }_0(R-r_0\mathrm{c}\mathrm{o}\mathrm{s} {\psi }_0){\mathrm{s}\mathrm{i}\mathrm{n}}^2 \delta \lambda \}\end{array}$$

(17)

where $r_0=R-H_j$, $\delta \lambda ={\lambda }_0-\lambda$, and $l_0=\sqrt{R^2+r_0^2-2R{r}_0\mathrm{c}\mathrm{o}\mathrm{s} {\psi }_0}$, simplifying Equation (13) to obtain

$$D_j\left(\phi ,\lambda \right)=\mathit{B}\Delta {\rho }_j\left({\phi }_0,{\lambda }_0\right)$$

(18)

where $\mathit{B}$ represents the correlation coefficient matrix between gravity anomalies and density anomalies. Using the Tikhonov regularization method [27], the density anomaly value $\Delta {\rho }_j$ can be solved:

$$\Delta {\rho }_j\left({\phi }_0,{\lambda }_0\right)={\left({\mathit{B}}^T\mathit{B}+\mathit{\alpha }I\right)}^{-1}{\mathit{B}}^T{D}_j\left(\phi , \lambda \right)$$

(19)

where $I$ is the identity matrix and $\mathit{\alpha }$ is the regularization factor.

3. Data

3.1. Bouguer Gravity Anomaly

The study area ranges approximately between longitudes 108° to 123°E and latitudes 5° to 23°N. The Bouguer gravity anomaly is from the WGM2012 model [28] (Figure 2). The range of gravity anomaly values is 10–565 mGal. The gravity anomaly value is the highest in the SCSB, gradually decreases towards surrounding trenches and island arcs, and becomes the lowest in the continental area. This phenomenon may be related to the thin oceanic crust in the SCS [29,30].

3.2. Sediment-Free Gravity Anomaly

Figure 3a shows the sediment thickness in the SCS [32], while Figure 3b illustrates the sediment gravity disturbances calculated by the Parker method, which range from −40 to 0 mGal. The absolute values of gravity anomalies exhibit a characteristic pattern centered on the sea basin, gradually increasing outward, with extreme values primarily distributed in subduction zones such as oceanic trenches and troughs. According to Figure 3a, the SCS region has many sedimentary basins with varying thicknesses. The YGHB has a thickness larger than 15 km. A few basins, such as the ZMB, have sediment layers with thicknesses of 3–7 km, while most basins have sediment layers with thicknesses of 1–2 km. By deducting the gravity effect caused by sediments, we obtained the gravity anomaly without sediments (Figure 4). After eliminating the gravitational effects of sedimentary layers, the regional characteristics of the seafloor spreading ridge become more pronounced.

4. Results

4.1. Wavelet Detail Gravity Anomalies

Based on previous experiments on the selection of optimal wavelet basis [33], we chose the “Coif3” wavelet basis and performed 2D wavelet multi-scale separation on the Bouguer gravity anomaly. The first eight-order wavelet approximation and wavelet detail gravity anomalies are shown in Figure 5 and Figure 6. Subsequently, the average field source depths corresponding to the wavelet detail signals were calculated using the radial logarithmic power spectrum method (Figure 7). The segment with the most intense amplitude corresponds to the power of the abnormal signal. We selected the point near the peak as the starting point and chose the turning point where the slope of the curve changes significantly as the end point. The points between starting point and end point were then linearly fitted to obtain the slope. Finally, we calculated the depth. The depth of each layer is presented in

Table 1. Depth range and thickness of each layer.

Number of Layers	Depth Range (km)	Average Field Source Depth (km)	Thickness (km)
1	0.0~1.6	0.8	1.6
2	1.6~3.2	2.4	1.6
3	3.2~7.0	5.1	3.8
4	7.0~12.0	9.5	5.0
5	12.0~22.8	17.4	10.8
6	22.8~33.4	28.1	10.6
7	33.4~51.2	42.3	17.8
8	51.2~70.4	60.8	19.2

.

As shown in Figure 6, the first- and second-order wavelet detail signals are very weak and are considered as noise. The gravity anomaly for D3 is approximately −40 to 40 mGal. The gravity anomaly signals are larger in the areas around the SCSB, such as the south of the SWSB, the west of Luzon Island, and the west of the SCSB. These regions are the continent–ocean transition (COT) zone. Meanwhile, the gravity anomaly signals in the SCSB are obvious near the seafloor spreading ridge. The gravity anomaly for D4 is approximately −50 to 80 mGal. Compared to D3, there is a slight increase in both the values and volumes of gravity anomaly signals in D4, with more pronounced changes observed in the MT region. The gravity anomaly for D5 is approximately −100 to 90 mGal, with an average field source depth of 18.42 km. There are many pronounced gravity anomaly signals, indicating that the tectonic composition at this depth is very complex and there is a boundary characteristic. The gravity anomaly for D6 is approximately −130 to 120 mGal. Both the values and the volume of gravity anomaly signals continue increasing. The Zhongnan–Liyue fault zone separates two negative gravity anomaly signals on the east and west sides [30]. The signals in the SCSB increase compared to previous wavelet details. The gravity anomaly for D7 is approximately −90 to 110 mGal. Signals continue increasing in the Luzon Island area, the SCSB, and the COT, while signals in the continental region begin to decrease, indicating that the structures in the continental region are relatively stable in this depth range. The gravity anomaly for D8 is also approximately −130 to 120 mGal. Compared to D7, there is no significant fluctuation in the gravity anomaly, suggesting that the structures at this depth are stable. There is a positive gravity anomaly signal in the northeast and southwest of the SCSB. Additionally, there is a pronounced negative gravity anomaly signal in the Luzon Island area. Since Luzon Island was formed after the SCS stopped opening and was subducted to the PSP [34], this negative signal may be associated with plate subduction.

In summary, the range and volume of gravity anomaly signals increase with depth. In D3 and D4, the anomalous signals in the COT and near the MT are relatively prominent. The anomalous gravity signals in D5 mainly originate from the Moho topography undulations, and the anomalous gravity signals in the sea basin area begin to strengthen. In D6, the negative anomaly signals on the two sides of the Zhongnan–Liyue fault zone have distribution characteristics similar to the terrain. In D7 and D8, there are obvious positive gravity anomaly signals in the northeast and southwest of the SCSB, which are presumed to be related to the opening of the SCS. In D7, there are positive gravity anomaly signals in the eastern part of the SCSB, while in D8, there are prominent negative gravity anomaly signals in the south of Luzon Island. These two gravity anomaly signals are adjacent to each other and located at the plate boundary. Combined with the formation of Luzon Island, it can be inferred that the anomalous signals in these two places are related to plate subduction.

4.2. Layered Density Anomaly

Each layer was modeled by a Tesseroid body with a spacing of 0.1° × 0.1°. Subsequently, the density anomaly was calculated using Equation (19). Simultaneously, based on the STW105 model [35], we calculated the average density of each layer (

Table 2. The average density value of each layer.

Number of layers	3	4	5	6	7	8
Average density value (g/cm3)	2.60	2.66	2.90	3.06	3.30	3.31

) and superimposed the average density onto the corresponding density anomaly. The results are presented in Figure 8. As first- and second-order signals are weak and considered as noise, the corresponding density anomaly results are not shown.

According to Figure 8, the density in D3 is approximately 1.2 to 3.5 g/cm³. There is no significant density variation in the SCSB. In contrast, the density variation is very pronounced in the COT. In the islands south of the SWSB, positive and negative density anomalies are alternately distributed. The density in D4 ranges from approximately 1.8 to 3.8 g/cm³. Compared with that in D3, the volume of density anomalies in D4 is larger. Obvious density anomalies begin to appear at the seafloor spreading ridge and are most likely caused by the expansion of the SCS around 15Ma ago [36]. The range of density values in D5 is approximately 1.8 to 3.9 g/cm³. The volume of density anomalies continues to increase, and significant density variations are observed in all regions except the continental area. These significant density variations exhibit strong lateral heterogeneity and indicate highly active tectonic movements, consistent with the Moho topography [2]. Notably, the plate boundary at 16°N acts as a “dividing line” for the MT, with the density variation in the north of the boundary line being significantly weaker than that to the south. Previous studies have indicated significant structural differences between the northern and southern segments of the MT, with the northern segment characterized as an accretionary margin and the southern segment as an erosional margin [37]. This difference is highly consistent with the distribution of density anomalies, indicating a strong correlation between them. The range of density values in D6 is approximately 2.5 to 3.6 g/cm³. Compared to that in D5, the amplitude of density in D6 decreases, but the volume of density anomalies continues to increase. The southern Zhongnan–Liyue fault zone separates two low-density anomalies. The density in D7 is approximately 2.7 to 3.8 g/cm³. There are high-density anomalies in the SWSB, the SWTB, and the southern segment of the MT. Among them, the high-density anomalies in the SWSB and the SWTB are related to the expansion of the SCS [34], while the high-density anomaly in the southern segment of the MT may be due to the subduction of the SCS plate towards Luzon Island [37]. The range of density variation in D8 is approximately 2.9 to 3.6 g/cm³. There are high-density anomalies in the SWSB-NWB area and to the north of the ESB, shaped like elongated ellipses. A low-density anomaly is present in southern Luzon Island, with the MT to its west and the plate boundary line to its north.

In summary, the volume of density anomalies continues to grow with increasing depth. The range of density variations reaches its maximum at D5, with a difference of more than 2.1 g/cm³ between the highest and lowest values. This occurs because D5 is located near the Moho topography, where there are significant density changes. In D5, the distribution of density anomalies in the northern and southern segments of the MT differs significantly. In D3, D4, and D5, the density of the COT exhibits strong lateral heterogeneity. In D6, the Zhongnan–Liyue fault zone appears to separate the two low-density anomalies. In D7 and D8, high-density anomalies are present in both the northeastern and southwestern parts of the SCSB. The possible causes of these phenomena are discussed in the following sections.

4.3. Gravity Anomalies of the Crust and Upper Mantle

To better analyze the major signal sources of the South China Sea, we conducted a radial logarithmic power spectrum calculation on the sediment-free gravity anomaly (Figure 4), and the results are shown in Figure 9. Two abnormal signals were identified in the gravity anomaly spectrum, corresponding to average field source depths of 10.30 km and 48.57 km, respectively, and regarded as shallow and deep anomalies. The details of these anomalies are presented in

Table 3. Depth range and thickness of two anomaly signals.

	Depth Range (km)	Average Field Source Depth (km)	Thickness (km)
Shallow anomaly signal	0.0~19.2	9.6	19.2
Deep anomaly signal	19.2~73.6	46.4	54.4

.

Since the depth range of shallow anomalies roughly aligns with that of the crust, we classified shallow anomalies as crustal gravity anomalies and deep anomalies as upper-mantle gravity anomalies. Combining the depth of each layer in Table 1, we concluded that the sum of the third-, fourth- and fifth-order wavelet details can represent crustal anomalies, while the sum of the sixth-, seventh-, and eighth-order wavelet details can represent upper-mantle anomalies, as shown in Figure 10.

The range of the crustal anomaly is approximately −110 to 110 mGal. Along the seafloor spreading ridge, there is a strip-shaped negative gravity anomaly. The gravity anomaly signals in the COT area and near the MT are very obvious. The range of the upper-mantle anomaly is approximately −160 to 160 mGal, with prominent gravity anomalies in the SCSB and the Luzon Island region.

4.4. Density Anomalies of the Crust and Upper Mantle

The depth information of the crustal and upper-mantle anomaly is presented in

Table 4. Depth of crustal anomaly and upper-mantle anomaly.

	Depth Range (km)	Average Field Source Depth (km)	Thickness (km)
Crust anomaly	3.2~22.8	13.0	19.6
Upper-mantle anomaly	22.8~70.4	46.6	47.6

. Similarly, we obtained the average density based on the STW105 model (

Table 5. The average density of crustal anomaly and upper-mantle anomaly.

	Crust Anomaly	Upper-Mantle Anomaly
Average density value (g/cm3)	2.75	3.23

). By combining the depth from Table 4 and the average density from Table 5, we performed layered density inversion to obtain the crust and upper-mantle density anomalies, as shown in Figure 11.

According to Figure 11, the range of the crustal density anomaly is approximately 1.6 to 3.6 g/cm³, corresponding to a depth of 13.0 km. High-density anomalies are widespread in the study area. Significant density anomalies are present in the MT and the Luzon region, as well as in the COT region. For the upper-mantle density anomaly, the range of density is approximately 2.9 to 3.6 g/cm³, corresponding to a depth of 46.6 km. The SCSB area is dominated by high-density anomalies, while the Luzon region primarily exhibits low-density anomalies. Compared to those in the crust, the density anomalies in the upper mantle are smoother. It is evident that the density variations in the shallow depth of the SCS are more pronounced in areas such as the COT and the MT, rather than in the SCSB. In contrast, the density variations in the SCSB are more significant at greater depth.

5. Discussion

5.1. Comparison with Previous Models

Yu et al. [2] also conducted a wavelet multi-scale analysis of gravity data in the SCS. During the wavelet multi-scale decomposition process, Yu et al. [2] just used the diagonal component to reconstruct the wavelet detail signals. However, we found that it is not enough. In this study, both the diagonal, vertical, and horizontal components were used for the reconstruction of wavelet detail signals. Therefore, our wavelet details are more reasonable and have a better horizontal resolution. There are high-density anomalies in the upper mantle of the SEHB and NWSB (Figure 11b). Scholars believe that the materials in this region are similar, primarily originating from the upper mantle of the Indo-China Plate. Driven by the collision between the Indo-China Plate and the Eurasian Plate, these materials were extruded eastward from the Indo-China Plate [2]. However, due to the rotation and northward migration of the PSP, the Shikoku Basin is separated from the SCS, aligning with the plate boundary in the SEHB-NWSB region [38]. The prominent high-density anomaly of the upper mantle beneath the SCSB is consistent with the high-velocity anomaly in this area [3].

5.2. Density Variations in the COT

According to Figure 8 and Figure 11, the COT exhibits pronounced density anomalies in the shallow layers. The Moho depth in the COT is approximately 15–25 km (Figure 12). Compared to the sea basin and continental shelf regions, the Moho depth variations in the COT are more dramatic (Figure 13). It is suggested that the SCSB experiences significant mantle uplift [39]. According to Figure 11, the crustal density variation in the COT is more obvious than that in the SCSB. In contrast, the upper-mantle density variation in the SCSB is more significant. We inferred that the SCSB has a typical thin oceanic crust due to the expansion of the SCS. Parts of the COT crust may fracture under tension [40]. These movements cause significant density anomalies in the shallow depth of the COT. Therefore, the crustal density structure in the COT is unstable. Meanwhile, the high-density anomalies in the shallow depth gradually expanded due to magmatic activity [41].

5.3. The Manila Trench and Luzon Island

The Manila Trench (MT) and Luzon Island (LZ) exhibit distinct density anomalies in both shallow and deep layers (Figure 8 and Figure 11). Some scholars suggested that the formation of Luzon Island resulted from the subduction of the PSP after the cessation of the SCS expansion [34]. Meanwhile, other scholars proposed that due to the interaction between the SCS and the PSP, Luzon Island underwent clockwise rotation, overriding the younger SCSB, leading the transform fault on the western side of Luzon Island to become a subduction zone, thus forming the MT [43]. Our density model supports this viewpoint. The density anomalies of the MT and LZ vary significantly with depth (Figure 13). It is suggested that the density difference between oceanic and continental areas causes lateral dragging between the plates, exerting an upward force on the plate boundaries [44]. However, due to significant density contrast across regions, the forces applied to different parts of the MT differ, resulting in notable structural disparities between the northern and the southern segments of the MT. There is an obvious low-density anomaly in the upper mantle of southern LZ (Figure 13), indicating partial melting or the upwelling of asthenosphere materials.

5.4. The Expansion of the SCS

The SCS experienced two seafloor spreading events [45]. The first seafloor spreading took place during the Late Eocene to Early Oligocene, with a NW-SE spreading direction, resulting in the formation of the NWSB and SWSB. The second seafloor spreading occurred during the Late Oligocene to Early Miocene, with an N-S spreading direction, leading to the formation of the ESB. Due to the clockwise rotation and northwest migration of the Philippine Sea Plate, the spreading direction of the SCS shifts from NW-SE to N-S. The direction of Zhongnan–Liyue fault zone is approximately N-S in the northern segment and changes to NW-SE in the southern segment [46]. We suggest that the Zhongnan–Liyue fault zone was formed during different stages of the expansion of the SCS. According to Figure 14, there are high-density anomalies in the NWSB, SWSB, SWTB, and the southern MT. These high-density anomalies may have been formed during the expansion of the SCS. Besides, the expansion of the SCS was not uniform. Scholars believe that neither the western nor the southern dynamic sources alone can adequately explain the expansion of the SCS [34,38]. Our findings strongly support this view. In addition to the traction from the PSP, other plates also exerted forces in different directions on the SCS, resulting in uneven driven forces and influencing the uneven distribution of density anomalies.

6. Conclusions

In this study, we used wavelet multi-scale analysis of gravity anomaly data and obtained a 3D density model of the SCS. Our model is more reliable than previous models. The expansion of the SCS affects both the crust and the upper-mantle structure. According to the density model, we suggest that the expansion of the SCS has caused the typical thin oceanic crust of the SCSB. Parts of the COT crust may fracture under tension. The low-density anomaly in the upper mantle of LZ may indicate partial melting or the upwelling of asthenosphere materials. Multiple upper-mantle high-density anomalies were observed in the SCS, which were formed during the expansion of the SCS. Our findings support this view that the expansion of the SCS is driven by forces from multiple directions. This study provides a reference to understand the expansion mechanism and tectonic evolution of the SCS.