3D Magnetic Inversion with Mixed L_p Norm Regularization for Unveiling the Spatial Distribution of Mesozoic Volcanic Arcs in the Northern South China Sea

1. Introduction

The South China Sea (SCS), one of the largest marginal seas in the Western Pacific, is situated at the junction of the Eurasian, Pacific, and Indochina plates (Figure 1). Its periphery exhibits three distinct types of continental margins: extensional margins to the north and south, a shear margin to the west, and a subduction margin to the east. Since the Mesozoic, the SCS has transitioned from an active to a passive margin, with rifting likely initiated on pre-existing subduction systems, followed by continental margin rifting, basin expansion, oceanic crust formation, subduction cessation, and continental collision during the Cenozoic[1]. Understanding the pre-Cenozoic tectonic architecture is critical for elucidating the evolution of the northern continental margin, but questions regarding Mesozoic subduction zones and pre-Cenozoic structures remain unresolved due to insufficient surface and subsurface evidence.

In recent years, numerous studies integrating gravity, magnetic, seismic, and drilling data have revealed the complexity of the crust in the northeastern South China Sea (SCS). A long-standing debate persists regarding whether an ancient oceanic crust exists in the region or if it represents a highly stretched continental crust[1,2,3,4,5,6,7,8,9,10,11]. A high magnetic anomaly belt, extending NE from the northwestern sub-basin edge through the Dongsha Rise and Taiwan Shoal to the western Taiwan margin, exhibits peak values of 60-200 nT. Broadband reflection seismic profiles further reveal high-velocity layers beneath the anomalies, and deep-reflection seismic profiles indicate Moho discontinuities[7]^. These geophysical features are commonly associated with Mesozoic subduction-related volcanic activity and mantle intrusion during continental rifting. Moreover, petrochemical analyses of rock and mineral samples have validated the existence of large-scale Mesozoic tectonic movements in the northern SCS, indicating a subduction-related hydrothermal context for rock and mineral formation[13,14,15]. Thus, investigating volcanic arc (subduction zone) distribution and morphology using high magnetic anomaly belt is crucial for understanding the pre-Cenozoic tectonic evolution.

Figure 1. Tectonic Map of the South China Sea and Adjacent Regions[17,18,19,20,21,22,23,24].①:Zhaoshu Liu et al., 1997; ②: Qiuju Huang et al., 2005; ③&④: Jinwei Gao et al., 2015; ⑤: Yan Pin et al., 2014; ⑥&⑦: Jinwei Gao, 2019;⑧: Chiyue Huang et al., 2001; ⑨, ⑩, ⑪, ⑫: Ryan Lester et al., 2014;⑬: Kirk McIntosh et al., 2014;⑭, ⑮, ⑯,⑰: Ryan Lester et al., 2011.

Figure 1. Tectonic Map of the South China Sea and Adjacent Regions[17,18,19,20,21,22,23,24].①:Zhaoshu Liu et al., 1997; ②: Qiuju Huang et al., 2005; ③&④: Jinwei Gao et al., 2015; ⑤: Yan Pin et al., 2014; ⑥&⑦: Jinwei Gao, 2019;⑧: Chiyue Huang et al., 2001; ⑨, ⑩, ⑪, ⑫: Ryan Lester et al., 2014;⑬: Kirk McIntosh et al., 2014;⑭, ⑮, ⑯,⑰: Ryan Lester et al., 2011.

The origin of the high magnetic anomaly belt in the northern South China Sea remains controversial, but most scholars believe it is a volcanic arc formed by magmatic processes related to Late Mesozoic subduction[11]. Di Zhou et al. (2006)[7], Hui Min et al. (2010)[10], and Zhaocai Wu et al. (2011)[11] utilized combined magnetic anomaly and gravity data to delineate the possible planar location of the volcanic arc (and ancient subduction zones) in the northern South China Sea. They suggested that the ancient subduction zone in the northern South China Sea was caused by the subduction of the ancient South China Sea beneath the South China continent. The subduction zone extends along the northern continental slope of the South China Sea, from the western edge of Taiwan Island in the southwest, passing through the southwestern depression of Taiwan to the southern part of the Chaoshan depression, and further extending southwestward to Yitong Ansha[10]; or it may be due to the subduction of the ancient Pacific beneath the South China continent, with a subduction zone extending from the southwestern Japan-Ryukyu Islands-Taiwan region southward to the Borneo and Palawan regions, and from the Meratus belt in southeastern Borneo to the Woyla belt in Sumatra. The subduction zone in the northern South China Sea originally extended from the Baiyun Depression northward, extending northeastward to the western part of Taiwan, including the Dongsha Rise[7]. Denke Hu et al. (2008)[12] used the multichannel reflection profile DPSR2002 crossing the Baiyun Depression (Figure 1, profile ②) and performed joint inversion of magnetic and gravity data. They proposed that the high magnetic anomaly corresponds to volcanic rocks at depths of 2.5 to 10 km, with a density of 2.65 g/cm³ and a magnetization strength of 0.7 A/m, which are products of a concealed oceanic crust subduction zone. Fucheng Li et al. (2018) [13] modeled the magnetic profile along the OBS survey line that crosses the Dongsha Rise, concluding that the high magnetic anomaly is associated with the root of a volcanic arc located at a depth of 18-24 km within the lower crust. However, these studies mainly rely on regional magnetic anomaly characteristics and 2D magnetic models for local sections, and lack a comprehensive exploration of the deep 3D magnetic distribution.

Magnetic 3D inversion is a primary technique for revealing the 3D magnetic structures at depth and has been widely applied in fields such as deep structural studies and resource exploration[14,15,16]. The magnetic 3D mixed Lp norm regularization inversion method (hereinafter referred to as mixed L_p norm inversion) is a newly developed technique in recent years[17,18]. Compared to single-norm regularization inversion methods, this technique offers the advantages of better delineating source boundaries, focusing the source morphology, and providing physical property values that are closer to the true characteristics of the source. We first investigate the magnetic 3D mixed L_p norm regularization inversion method. Then, using high-resolution aeromagnetic anomaly data from the northern South China Sea and deep seismic reflection data as reference constraints, we construct a reliable 3D magnetic susceptibility model of the crust through magnetic 3D mixed L_p norm inversion. The study proceeds to analyze the spatial distribution of the Mesozoic volcanic arc, providing deep insights for pre-Cenozoic tectonic studies in the South China Sea.

2. Mixed L_p Norm Regularization for 3D Magnetic Inversion

2.1. Methodology

Fournier and Oldenburg (2019)[17] proposed a mixed L_p norm regularization algorithm for geophysical inversion, which replaces the regularization term in the traditional Tikhonov regularization (L₂ norm) objective function with norms $p_{x,y,z,s}$∈ {0, 1, 2}. These norms control the regularization term and play a crucial role in determining the final model's features. The objective function for the traditional Tikhonov regularization inversion is as follows:

$$\underset{m}{\min }\varphi \left(m\right)={\varphi }_d+\beta {\varphi }_m$$

(1)

Among them, Φₘ is the regularization term of the objective function, which plays a key role in determining the characteristics of the final model. In the mixed L_p norm inversion, $p_{x,y,z,s}$ are independent values of 0, 1, or 2. The regularization term is as follows:

$${\Phi }_m={{\displaystyle \sum }}_{j=s,x,y,z}{\alpha }_j{{\displaystyle \int }}_V{\left|f_j\left(m\right)\right|}^{p_j}dV$$

(2)

In the formula:

$$f_s=m,f_x={\displaystyle \frac{𝜕m}{𝜕x}},f_y={\displaystyle \frac{𝜕m}{𝜕y}},f_z={\displaystyle \frac{𝜕m}{𝜕z}}$$

(3)

The discrete form of the formula is as follows:

$${\Phi }_m={\alpha }_s{\Vert V_{s}m\Vert }_{p_s}^{p_s}+{{\displaystyle \sum }}_{j=x,y,z}{\alpha }_j{\Vert V_j{G}_{j}m\Vert }_{p_j}^{p_j}$$

(4)

By employing various mixed norms for model regularization, the mixed L_p norm inversion method generates a range of solutions. This approach highlights the inherent multiplicity of solutions in the inversion process and mitigates the risk of over-interpretation, which can arise from considering only a single solution.

The data experiments with the methods mentioned above indicate that when sparse norms (L₀, L₁) are applied to the model vector $m$, an issue of dependence on the upper limit of physical properties arises. Conversely, when sparse norms are applied to the model gradient ∇m, a "staircase" effect is observed. The former results in the physical property values of the anomaly body in the inversion approaching the upper limit, while the latter causes the anomaly body to exhibit horizontal or vertical boundaries. To address these issues, Sun and Wei (2021)[17] further developed this approach based on the work of Fournier and Oldenburg (2019). The regularization term in the objective function of the mixed Lp norm inversion consists of four components, as shown below:

$${\Phi }_m={\alpha }_s\int {\left|f_s\left(m\right)\right|}^{p}dv+{\displaystyle \sum }_{j=x,y,z}{\alpha }_j\int {\left|f_j\left(m\right)\right|}^{q}dv$$

(5)

The discrete form of the formula is as follows:

$${\Phi }_m={\alpha }_s{\Vert W_{s}m\Vert }_p^p+{{\displaystyle \sum }}_{j=x,y,z}{\alpha }_j\Vert W_j{\nabla }_j{m}_j{\Vert }_q^q$$

(6)

Where the f is defined as in equation (3).

The parameters ${\alpha }_s$ and ${\alpha }_j$ are weighting coefficients that control sparsity and smoothness, while $p$ and $q$ are the norm exponents that regulate the smoothness of the model's small-scale features and gradients, respectively. $Ws$ and $Wj$ are weight matrices used to adjust the model's depth and other structural characteristics.

The method utilizes two norms, $p$ and $q$, to control the model$m$ and its gradient $\nabla m$, respectively. When $p=q=1 or 2,$ it simplifies to conventional L₁ and L₂ norm regularized inversion. The algorithm operates in two stages: first, it solves for the L₂ norm solution, then switches to the desired L_p norm suite, using the L₂ norm solution as the initial model to solve for the L_p norm solution. To handle the large variations in the values of the regularization function when the value of p changes, and to ensure that each component of the regularization contributes to the final solution, the Lawson algorithm is used to approximate the L_p norm. The iterative results are continuously readjusted in the Gauss-Newton solution. According to experimental results, when using the L₀₂ norm, the mixed L_p norm inversion does not encounter issues of bound dependence and staircase effects. Compared to conventional L₁ and L₂ norm inversions, the mixed L_p norm inversion provides more focused results with clearer boundaries. However, the selection of the optimal mixed norm requires empirical testing to determine the best choice in practical applications.

2.2. Theoretical Model Experiment

The theoretical model consists of two combined cubes. Cube 1 has a side length of 400 m, with its center at X = 1000 m, Y = 1000 m, Z = 500 m, and a magnetic susceptibility of 0.04 (SI). Cube 2 has a side length of 200 m, with its center at X = 2000 m, Y = 2000 m, Z = 300 m, with a magnetic susceptibility of 0.08 (SI). The Earth's magnetic field strength is 55,000 nT, with a dip angle of 90° and a declination of 0°. Both cubes have a magnetic inclination of 90° and a declination of 0°. The magnetic anomalies generated by the combined cube model at the observation plane of Z = 0 m are shown in Figure 2, which includes 2% Gaussian noise.

The white dashed lines represent the positions of the cubes, while the black dashed lines denote the locations of the profiles shown in Figure 3 and Figure 4

Figure 3 and Figure 4 show the inversion results using conventional L₁ and L₂ norms, as well as the mixed L_p norm regularization. It can be observed that the L₁ norm provides stronger boundary constraints compared to the L₂ norm, with the inversion results better restoring the model boundaries both laterally and vertically. The L₂ norm inversion results exhibit divergence in the vertical boundaries of the anomaly body. In contrast, the mixed L_p norm not only better restores the model boundaries but also produces physical property values closer to the true values. Among the mixed norms, L₀₂ provides better results compared to L₁₂.

The black box represents the location of the true model, and the black dashed lines indicate the locations of the profiles shown in Figure 5. Figure (a) and (b) show the inversion results at Z = 350 m and along the black dashed line in Figure 2 using the L1 norm regularization algorithm; Figure(c) and (d) show the inversion results at Z = 350 m and along the black dashed line in Figure 2 using the L2 norm regularization algorithm.

The black box represents the location of the theoretical model, and the black dashed lines indicate the locations of the profiles shown in Figure 5. Figure (a) and (b) show the inversion results at Z = 350 m and along the black dashed line in Figure 2 using the L₀₂ norm regularization algorithm; Figure (c) and (d) show the inversion results at Z = 350 m and along the black dashed line in Figure 2 using the L₁₂ norm regularization algorithm.

Figure 5 shows the extracted profiles from each slice along X = 1000m. The Z = 350 m horizontal slice profile results (Figure 5(a)) indicate that the horizontal boundary delineation of the anomaly body is similar for both conventional norm inversion and mixed L_p norm regularization. However, the restoration of physical properties is notably better with the mixed norm, with peak values following the order: L₁₂ > L₀₂ > L₁ > L₂. The profile extraction along the black dashed line in Figure 2 (Figure 5(b)) reveals significant differences in the delineation of the vertical boundary of the anomaly body across the various norms. Among the different norms, L₀₂ provides the most focused result, while L₁ and L₁₂ are similar, and L₂ exhibits the most divergence, with a distinct "tailing" effect along the Z direction. The restoration of physical properties follows a similar trend to the Z = 350 m horizontal slice, with peak values L₁₂ > L₀₂ > L₁ > L₂. The inversion result with L₀₂ shows a higher effective magnetization, meaning that within the theoretical model range (Z range from 300 to 700 m), the magnetization curve delineates a larger area. Overall, from the perspective of boundary delineation and physical property restoration, the mixed L_p norm regularization inversion yields better results compared to conventional L₁ and L₂ norm regularization, with L₀₂ outperforming L₀₁.

3. Aeromagnetic Data and Physical Property Data

3.1 Airborne Magnetic Anomaly Data

We obtained the aeromagnetic data for the northern South China Sea through vectorization from the 2021 China and Adjacent Sea Areas Aeromagnetic Series Map (scale 1:5,000,000) compiled by the China Aero-Geophysical Survey and Remote Sensing Center for Natural Resources. This series consolidates a substantial amount of aeromagnetic survey data collected up to 2021, along with the CCOP East Asia Sea anomaly map data (i.e., data from the CCOP GSJ AeroMAG1966 East Asia Sea anomaly map), selecting the most recent and highest precision data. Using digital mapping and cartographic techniques, the final data grid spacing is 1 km × 1 km.

Given that the northern South China Sea is located in a low-latitude region, we applied the stable pseudo-inclination correction technique along with its latitude-dependent algorithm[31], in order to obtain the corrected magnetic anomaly data for the study area (Figure 6(a)). From the data, it can be observed that the corrected magnetic anomaly values in the northern South China Sea predominantly range from -100 to 100 nT, with positive anomalies being the most common. Low-value anomalies are mainly concentrated in the depressions and sedimentary areas. Additionally, a NE-trending high-value positive anomaly belt is visible across the region, along with several sporadic high-value anomalies.

The high magnetic anomaly belt is located to the north of the Chaoshan Depression, extending along the Dongsha Rise to Taiwan Island, and is believed to be related to the volcanic arc generated by the subduction zone (Di Zhou et al., 2006; Hui Min et al., 2010) [7,10]. The sporadic high-value anomalies may be caused by isolated volcanoes (Fan et al., 2017). Additionally, we calculated the magnetic anomaly modulus for the northern South China Sea (see Figure 5(b)), and the results show that its distribution is generally consistent with that of the corrected magnetic anomalies in the high magnetic anomaly belt region, indicating that the remanent magnetism effect in this area is relatively weak. Therefore, conventional magnetic 3D inversion methods, such as the aforementioned mixed L_p norm inversion method, can be used for inversion analysis.

3.2. Physical Property Data

The differences in the physical properties of rocks and ores are the basis of geophysical exploration. A thorough understanding of the physical property parameters of various rocks, ores, and strata in the study area is a prerequisite for inversion interpretation. This study uses rock magnetic measurements and research from the northern South China Sea coastal land and islands (Tanggen Li et al., 1987) [26], as well as the summary report of the comprehensive geophysical survey of the northern South China Sea continental slope area (Shaofei Liu, 2016)[27], to statistically analyze the magnetic reference values for the magmatic rocks, sedimentary strata, and lower crust materials in the area. Geological bodies with similar magnetic parameters or small physical property differences are classified into a single magnetic susceptibility unit. Based on the magnetic susceptibility and corresponding geological types, the study area is divided into six geological units with distinct magnetic susceptibility differences, as follows: (1) intermediate to acidic intrusions (with a large range of magnetic variations and strong magnetism), (2) basic and ultrabasic intrusions (high magnetism), (3) sedimentary strata (low or no magnetism), (4) crust (no magnetism), (5) seawater layer (no magnetism), and (6) mantle (no magnetism). The statistical results are shown in Table 1.

4. Inversion and Results Discussion

4.1. Inversion and Results

Referring to the deep reflection seismic data mentioned earlier (Figure 1), multiple deep reflection seismic profiles reveal depth characteristics resembling ancient volcanic arc structures, providing reliable depth constraints for this study. During the inversion process, the inversion parameters were selected by integrating the depth feature constraints, ensuring that the magnetic source depth of the inversion results in the region was as consistent as possible with the depth of the ancient volcanic arc structure revealed by the deep reflection seismic data. Finally, by selecting the parameters for the mixed L_p norm regularization (p = 0, q = 2, v = 9, αs = 0.3), a 3D magnetic susceptibility model of the crust in the northern South China Sea was obtained. The model has a lateral grid spacing of 0.1° × 0.1°, a depth range of 0-40 km, with a depth step of 2 km. The inversion results are shown in Figure 7.

The model reveals that the planar distribution of the magnetic sources in the northern South China Sea is generally consistent with the high-value regions of the de-clustered aeromagnetic anomalies, presenting an NE-trending high magnetism belt, with the main depth concentration between 20 and 28 km. In addition, the model also shows isolated magnetic sources scattered to the north of the belt, which is consistent with the speculation by Fan et al. (2017) regarding the submarine volcanic intraplate magmatic sources in the northern South China Sea. The inversion residual anomalies are shown in Figure 8, where the relatively small residual values indicate a good fit between the theoretical anomalies of the model and the actual de-clustered aeromagnetic anomaly data.

Only the portion with magnetic susceptibility greater than 0.02 (SI) is shown

4.2. Deep 3D Distribution Characteristics of Volcanic Bodies

To verify the reliability of the model, a multi-channel reflection seismic profile from Gao et al. (2019)[24] (profile line ⑦ in Figure 1) was compared with the model (Figure 9). The results show that the model profile matches well with the deep reflection seismic data, demonstrating the credibility of the model. As mentioned earlier, the Dongsha Rise exhibit significant high magnetic anomalies, which correspond to high-magnetism bodies revealed at depths of 20-30 km in the model. In addition, the deep reflection seismic profile also reveals a Moho discontinuity beneath the Dongsha Rise (at 125 km along the profile), with a high-velocity body in the lower crust at this location. The model inversion results further indicate that the high-magnetic body beneath the Dongsha is located at the bottom of the lower crust, and its planar position and depth are in good agreement with the distribution of the high-velocity body in the lower crust.

The multi-channel deep reflection seismic results, from shallow to deep, include the water layer, sedimentary layer, and the Moho interface.

To more clearly illustrate the spatial distribution of high magnetic susceptibility bodies, four different profiles were selected for visualization, with the magnetic susceptibility distribution and positions shown in Figure 10. The base map represents a horizontal slice of the model at a depth of 25 km, displaying the distribution characteristics of high magnetic susceptibility bodies, primarily composed of intermediate to mafic volcanic rocks. The horizontal slice reveals that the peak magnetic susceptibility of volcanic bodies (approximately 0.04 SI) forms an NE-trending belt along the Dongsha Rise, with large areas of relatively high magnetic material (approximately 0.02 SI) and scattered high-value magnetic material (approximately 0.03 SI) distributed to the north of the belt. The vertical profiles indicate that the volcanic bodies are primarily located in the lower crust (at depths of 20-30 km), suggesting that the high magnetic anomaly zone in this area is caused by deeper and larger-scale volcanic bodies, which is consistent with the findings of Li et al. (2018).

The base map is a horizontal slice at depth=25km

4.3. Discussion on the Origin of the Volcanic Arc

The tectonic evolution of the South China Sea region has long been a focus of research, with many scholars proposing that plate subduction plays a critical role in the geological activities of this area. The OBS2006-3 velocity profile (location shown in Figure 10,Xiaodong Wei et al., 2014[28]) , which crosses the Dongsha Rise, reveals the crustal structure characteristics of the region. The velocity profile indicates that the upper crust is relatively thick, exceeding 10 km, and features an anomalous uplift at the base of the upper crust with a P-wave velocity (Vp) of up to 6.9 km/s. The lower crust is approximately 12 km thick, with a Vp of 7.4 km/s, showing a distinct uplift phenomenon. This uplift corresponds well with the horizontal position and depth of the volcanic bodies revealed by the magnetic susceptibility model, suggesting that crustal thickening and velocity anomalies may be the result of mantle-derived, mafic, high-magnetic, and high-density materials migrating upwards. The intense magmatic activity caused the formation of a high-velocity layer in the lower crust, which extends laterally on both sides.

Qiuju Huang et al. (2005), based on deep reflection seismic data (profile DPSR2002, marked as profile ② in Figure 1), revealed the potential existence of a subduction zone at depth beneath the Baiyun Depression. The boundary between the Dongsha Rise and the Baiyun Depression is considered to be the junction of the subducting plate and the Moho discontinuity. The study found that the volcanic bodies are located above the Moho, at the leading edge of the subducting plate, while the uplift of the high-velocity body in the lower crust suggests upwelling of mantle-derived material. This can be interpreted as the upwelling of mafic material induced by the subducting plate, forming high-temperature magma chambers that gradually evolve into a volcanic arc. Furthermore, to the north of the high-magnetic anomaly belt, the distribution of large areas of higher magnetic material and scattered high-value magnetic anomalies may be the result of subsequent upwelling, diffusion, and localized intrusion of mafic mantle material within the plate.

In addition, drilling data[29] indicate that during the Late Jurassic to Early Cretaceous, the sedimentary environment in the Chaoshan Depression transitioned rapidly from a deep-sea facies to a coastal facies. This suggests that the northern South China Sea experienced significant continental crust uplift during the Late Jurassic, corresponding to subduction of the plate. Overall, the deep magnetic source characteristics closely correspond with the deep reflection seismic and OBS seismic data, indicating that the high magnetic anomalies are consistent with the regional deep structures. The petrophysical data indicate that the high-magnetic materials are primarily igneous rocks, and drilling data reveal Mesozoic subduction in the region. Xu et al. (2016, 2017)[13,14] revealed the presence of Early Jurassic arc-related magmatism in the northeastern South China Sea based on drilling results, primarily manifested by granites and diorites dated to 198-195 Ma. These rocks are characterized by a magnesium-rich, low-potassium tholeiitic basalt series with features including Nb-Ta depletion, low temperatures (678-696 °C), and characteristics of I-type magmatic rocks. Their Sr-Nd isotopic compositions (⁸⁷Sr/⁸⁶Sr_i = 0.705494–0.706623, εNd_t = -0.9 to +2.2) indicate derivation from evolved mantle-derived magmas. Additionally, the enrichment of fluid-mobile elements (e.g., Cs to K and Pb) suggests significant involvement of subduction zone fluids. Similarly, granitic-dioritic rocks found in the southwestern East China Sea, dated to approximately 187 Ma, also exhibit characteristics of arc-related low-temperature magmatism (around 675 °C), with Sr-Nd isotopic signatures (⁸⁷Sr/⁸⁶Sri = 0.705200, εNdt = 1.1) further supporting their origin from evolved mantle-derived magmas. Furthermore, analysis of Paleocene detrital zircons indicates the presence of Early Jurassic arc-related low-temperature magmatism (600-740 °C) in the southwestern East China Sea, suggesting that the provenance may be associated with the Early Jurassic magmatic arc.

According to these findings, it is concluded that the high magnetic anomaly belt in the northern South China Sea is caused by a Mesozoic volcanic arc. This volcanic arc extends in a NE direction, with the main depth between 20 and 28 km, stretching from the Dongsha Rise to Taiwan Island. The volcanic arc near the Dongsha Rise exhibits a larger scale, which is inferred to be the initiation point of the subduction zone. Based on this comprehensive evidence, the study reveals the complexity of the deep structural features in the northern South China Sea and provides new constraints on the geological evolution of the region.