Driven magmatism and crustal thinning of coastal southern China in response to subduction

The late Mesozoic igneous rocks along the coastal South China Block (SCB) exhibit complex parental sources involving a depleted mantle, subducted sediment-derived melt, and melted crust. This period aligns with the magmatic flare-up and lull in the SCB, debating with the compression or extension in coastal region. Our study employs numerical models to investigate the dynamics of the ascent of underplating magma along the Changle–Nan'ao Belt (CNB), simulating its intrusion and cooling processes while disregarding the formational background. The rheological structure of the lithospheric mantle significantly influences magma pathways, dictating the distribution of magmatism. This work reveals that the ascent of magma in the presence of faults is considerably faster than in the absence of faults, and contemporaneous magmatic melts could produce different cooling and diagenetic processes. Additionally, the influence of pre-existing magma accelerated the emplacement of underplating magma. The magma beneath the fault ascended rapidly, reaching the lower crust within 20 million years, with a cooling rate of approximately ∼ 35 °C Myr⁻¹. Conversely, the thickened magma took 40–50 million years to ascend to the lower crust, cooling at a rate of ∼ 10 °C Myr⁻¹. In contrast, magma without thickening and fault would take a considerably longer time to reach the lower crust. The ascent of magma formed a mush-like head, contributing to magmatic circulation beneath the crust and decreasing crustal thickness. Multiphase magmatism increases the geothermal gradient, reducing lithospheric viscosity and promoting underplating magma ascent, leading to magmatic flare-ups and lulls. Our findings suggest that the Cretaceous magmatism at different times in the coastal SCB may be associated with the effects of lithospheric faults under similar subduction conditions. Boundary compression forces delay magma ascent, while rising magma induces a significant circulation, decreasing the crustal thickness of the coastal SCB. This study provides new insights into the complex interplay of magmatic processes during subduction, emphasizing the role of lithospheric structure in shaping the temporal and spatial evolution of coastal magmatism.

1 Introduction

Magmatism characterized by periodic flare-ups and lulls at convergent plate margins usually manifests a subduction-related origin (Brown, 1994). However, there is no agreement regarding the relationships between magmatism and the roles of subducted slabs and the corresponding subduction styles (Morris et al., 2000; Faccenna et al., 2010). Instead, a contradictory and multifaceted dynamic process emerges due to the complex magmatic composition and considerable time span (Zhou et al., 2006) (Yoo and Lee, 2023). Notably, the absence of high-pressure blueschist, which is typically associated with subduction, has generated intense debate surrounding the accurate timing of initial subduction and variations in slab dip, with proposals ranging from flat-slab subduction to shallow and steep subduction (Suo et al., 2019; Xu, 2023; Su, 2023). Some researchers proposed a model involving Early Permian flat-slab subduction and a Jurassic foundering model in the South China Block (SCB); however, this model lacks substantial evidence of early stage subduction-related magma (Li and Li, 2007). While most research has focused on Jurassic subduction, Early Jurassic intraplate igneous rocks deviate from typical subduction arc-related rocks and display inconsistencies in spatiotemporal distribution during coastward migration (Zhou and Li, 2000; Xu et al., 2017; Li et al., 2019). Researchers contend that the SCB did not immediately experience the influence of the Paleo-Pacific Plate in the Early–Middle Jurassic but experienced intensified activity in the Late Jurassic to Cretaceous (Gan et al., 2021). Gradual steepening of the shallowly subducting slab since the Middle Jurassic is proposed to explain the corresponding flare-up of magmatism in the SCB (Zhou et al., 2006; Mao et al., 2021). In contrast, Xu et al. (2023) interpreted voluminous intraplate silicic magmatism as a response to slab stagnation and coastward migration, overlooking contemporaneous compressional deformation in the Late Jurassic and the end of the Early Cretaceous. In contrast to the coastward migration model, the end of the Early Cretaceous is considered a magmatic lull resulting from crustal shortening due to resubduction of the slab (Wei et al., 2023). The uncertainty lies in whether subduction-induced magma can migrate and intrude concurrently according to these different models. The Mesozoic tectonic magmatism in the SCB was intricately linked not only to lithospheric properties but also to subduction rates and mantle flow (Su, 2023). The transport of intrusive magma spans a significant period, and the pathways of ascending magma play a crucial role in determining the distribution of magmatism during the emplacement process. Unfortunately, lag magma, which is potentially misunderstood as originating from other sources, has received limited attention in current discussions.

The coastal South China Block (SCB) is characterized by the development of a 40–60 km wide, NE–SW-striking ductile shear zone known as the Changle–Nan'ao Belt (CNB) (Cui et al., 2013). This belt comprises gneiss exhibiting evidence of metamorphism at greenschist, amphibolite, and granulite facies (Li et al., 2015). Intrusions of gabbro, diorite, granodiorite, monzogranite, and two-mica granite plutons are also prevalent within the CNB. The U–Pb age analysis of the oldest orthogneiss in the CNB yields a date of 187 ± 1 Ma, with the youngest orthogneiss dating to 130 ± 1 Ma. Additionally, ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}$Ar plateau ages fall within the range of 118–107 Ma (Wang and Lu, 2000). The occurrence of voluminous igneous rocks spans two significant periods, ranging from 143–130 and 110–95 Ma. A debated magmatic lull observed between 130 and 110 Ma is associated with syncollisional orogenesis (Chen et al., 2020; Wei et al., 2023) or post-orogenic extension (Cui et al., 2013; Li et al., 2014; Zhao et al., 2015; Xu et al., 2023) according to the involved deformation and magmatic rocks. Despite these insights, the relationship between intruded magma and orogenesis remains uncertain, mirroring the ambiguous understanding of potential variations in magma migration time in response to tectonic stress. In addressing this uncertainty, our exploration focuses on understanding the emplacement and cooling processes of mantle magmatism and their influence on crustal structure. Additionally, we delve into the question of whether changes in the dynamic background instantaneously produce magmas with different properties that can effectively intrude shallow layers.

2 Emplacement and origin of Cretaceous magma

Cretaceous magmatic rocks cover an area of approximately 117 190 km² on the southeastern coast of China and in the lower Yangtze region of the SCB (Fig. 1) (Liu et al., 2020). The crustal Poisson's ratios range from 0.22–0.26 in the interior to 0.26–0.29 in the eastern coastal region (Guo et al., 2019), implying a high content of felsic minerals and an increasing proportion of mafic minerals from the interior to the coast (Ji et al., 2002). Seismic profiles reveal transparent reflective features of felsic rocks in the upper crust and abundant high-amplitude, short, and isolated reflections of mafic sills in the middle–lower crust (Li et al., 2023). The velocity ratio of P waves to S waves ($V_{\mathrm{p}}/V_{\mathrm{s}}$) on the coast is 1.76, which is slightly greater than that in the interior SCB but lower than the value of 1.79 for mafic underplating of the lower crust (Deng et al., 2019). The coastal Vp values of the lower crust (∼ 6.5 km s⁻¹) are not compatible with the mafic composition (V_p >7.0 km s⁻¹), implying that mafic magma underplating was not common or was removed (Guo et al., 2019). However, high-resistivity anomalies in electrical resistivity profiles indicate local mafic magma underplating (Cheng et al., 2021).

Figure 1 (a) Regional geological map of southeastern China showing distribution of Mesozoic magma. (b) Geological magmatism of the Changle–Nan'ao Belt and corresponding ages (refer to Wei et al., 2023).

The coastal granitic gneisses exhibit varying cooling rates during different periods. These rates were ∼ 35 °C Myr⁻¹ at 130–120 Ma, 13–20 °C Myr⁻¹ at 126–110 Ma, ∼ 10 °C Myr⁻¹ at 110–100 Ma, and ∼ 80 °C Myr⁻¹ at 100–90 Ma (Chen et al., 2002, 2020). The structural pattern of the CNB shows solid-state ductile deformation with temperatures of 300–350 °C in the mylonitic gneiss and deformed volcanic rocks. Syntectonic granitoids exhibit subsolidus magmatic flow (Wang and Lu, 2000; Wei et al., 2015). The temperature and pressure of metamorphic minerals are 540–610 °C and 0.28–0.35 GPa, respectively, on the basis of plagioclase–amphibolite, and 485–640 °C and 0.3 GPa, respectively, based on a mica–quartz schist (Wang and Lu, 2000). The Cretaceous magmatic rocks in the SCB belong to the I-and A-type series, with high-K calc-alkaline to shoshonitic affinities and arc-like features. They are enriched in light rare-earth elements (LREEs) and large ion lithophile elements (LILEs) but depleted in heavy rare-earth elements (HREEs) and high-field-strength elements (HFSEs). They exhibit negative εNd(t) values ranging from −10.1 to −0.3 and variable zircon εHf(t) values ranging from −29.7 to +10.3. The mafic rocks have εNd(t) values ranging from −14.27 to +8.0 and εHf(t) values ranging from −9.5 to +1.9 (Chen et al., 2020). Isotopic data indicate mixed sources, including ancient crust-derived, enriched mantle-derived, and depleted mantle-derived material. Some mafic rocks possibly originated from the melting of the mantle wedge metasomatized by melts from the subducted slab and sediments.

The mixed magmatic source indicates that partial melting was independent of the type of subduction transport across the lithosphere. Subduction-induced melts migrate upwards, and their pathways change depending on the stress conditions, resulting in different time–temperature histories for emplacement and deformation of gneissic magma.

3 Numerical simulation and model setup

The late Mesozoic magmatism in the SCB was triggered by the subduction of the Paleo-Pacific plate (Su, 2023). A substantial volume of magma, originating from the subducted slab, asthenosphere, and the lithospheric base, ascended towards the surface and accumulated at the lower boundary of the lithosphere (Figs. 1 and 2a). To unravel the dynamics of coastal magmatism in response to subduction geometry, two end-member numerical models are considered. Our simulation focuses solely on modeling the intrusion and ascent of magma that has accumulated at the bottom of the lithosphere, without delving into the origins of magma from deeper sources (Fig. 2b). The models are two-dimensional domains (400 km wide and 100 km deep) representing a trench-perpendicular cross-section of a subduction zone, with the trench located to the right of the model (Deng et al., 2019). The underplating magma is represented by a 4–8 km thick body located at 100–80 km depth along the bottom of the model. The 8 km thick magma represents the accumulation of underplating magma due to subduction. A 10 km thick square domain is assigned at 40 km depth, representing previously intruded magma in the lithospheric mantle. The upper 30 km deep rectangle is continental crust. The polygon located between the underplating magma and crust is the mantle lithosphere. A 45 km wide rectangular domain is assigned at the left of the model, representing a lithospheric fault (Cui et al., 2013).

Figure 2 (a) Schematic cross-section illustrating the subduction of the Paleo-Pacific plate. (b) Reference model geometry depicting temperature, density, and viscosity variations with depth (see location in Fig. 1). The boundary conditions are the same between Model 1 and Model 2, except that Model 2 is assigned a horizontal force on the right side.

3.1 Governing equations

The materials of the domain are regarded as incompressible viscous fluids according to the Boussinesq approximation. The models satisfy the following mass, momentum, and energy conservation equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\text{(2)}& {\displaystyle }\mathrm{\nabla }\cdot \left[\mathit{\eta }\cdot \left(\mathrm{\nabla }\mathit{u}+{\left(\mathrm{\nabla }\mathit{u}\right)}^T\right)\right]-\mathrm{\nabla }P+\mathit{\rho }\mathit{g}=\mathrm{0},\text{(3)}& {\displaystyle }\mathit{\rho }{C}_p\cdot ({\displaystyle \frac{\partial T}{\partial t}}+\mathit{u}\mathrm{\nabla }T)=\mathrm{\nabla }\left(k\mathrm{\nabla }T\right)-\mathit{\alpha }\mathit{\rho }{\mathit{v}}_{z}T,\end{array}$$

where u is the velocity field, η is the viscosity, T is the temperature, P is the pressure, ρ is the density, g is gravity, C_p is the specific heat, t is the time, k is the thermal conductivity, α is the thermal expansion coefficient, and v_z is the vertical velocity component (Rodríguez-González et al., 2012). The heat capacity is set to 1000 J (kg K)⁻¹, the ratio of specific heat (gamma) is set to 1, and the thermal conductivity is set at 2.5 W (m K)⁻¹ in all domains of the models (Chapman, 2021). The density of the continental crust varies linearly with depth, increasing from 2600 kg m⁻³ at the surface to 2900 kg m⁻³ at 30 km depth. The mantle lithosphere has a fixed density of 3400 kg m⁻³, whereas the molten magma is modeled with a fixed density of 2800 kg m⁻³ (Chapman, 2021; Su, 2023).

The crustal material has a power law stress–strain rate relationship (Chapman, 2021).

Here, η is the dynamic viscosity and is expressed as follows:

$$\begin{array}{}\text{(4)}& \mathit{\eta }=\left({\displaystyle \frac{d^p}{A\cdot f_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}}\right){}^{\frac{\mathrm{1}}{n}}\cdot {\mathit{\epsilon }}_{\prod }^{\frac{\mathrm{1}-n}{n}}\cdot \exp \left({\displaystyle \frac{E+P_{\mathrm{lit}}\cdot V}{n\cdot R\cdot T}}\right),\end{array}$$

where A is the pre-exponential factor; E and V are the activation energy and volume, respectively; P_lit is the lithostatic pressure; R is the gas constant; and n is the stress exponent. The dynamic viscosity was calculated via the wet quartz flow law of Hirth et al. (2001). The viscosity was calculated at each time step using the temperatures returned from the model, a quartz material parameter (A) of 1.36742 × 10⁻⁵ MPa⁻ⁿ s⁻¹ with a stress exponent (n) of 4, a quartz activation energy (E) of 135 kJ mol⁻¹, a water fugacity ($f_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$) of 1000 MPa, and a strain rate of 10⁻¹⁵ s⁻¹ (Chapman, 2021). The strain applied in the modeling takes a 2D approximation and is based on slab subduction studies (Liu and Currie, 2019). The viscosity was updated after each time step based on the temperature. The mantle lithosphere and the molten magma are modeled with constant dynamic viscosities of 1 × 10²¹ and 1 × 10²⁰ Pa s, respectively. The model fault was assigned a low viscosity of 1 × 10¹⁹ Pa s to represent the CNB (Vissers et al., 1995; Columbu et al., 2015).

The models are run using the three-phase flow, phase field interface option, which accounts for the surface tension between immiscible phases, the contact angles with the walls, and the density and viscosity of each fluid. The three-phase flow model obeys the Cahn–Hilliard equation (Boyer et al., 2010). The fluid motion causes the phase field variable to change from phase to phase, but the sum of all phase field variables ∅_i at each point in the space is 1. Its expression consists of the order parameter of each phase as in the following equation:

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{l}{\mathrm{\varnothing }}_i=\mathrm{phi}i\\ {\mathrm{\varnothing }}_a+{\mathrm{\varnothing }}_b+{\mathrm{\varnothing }}_c=\mathrm{1}\end{array}\right..\end{array}$$

The three phases are continental crust ∅_a, mantle lithosphere ∅_b, and magma ∅_c. Three-phase flow is automatically computed via a phase initialization study step by solving for the geometrical distance to the initial interface. The initialized three-phase flow function is then defined from the analytical steady state solution for a straight fluid–fluid interface.

3.2 Model setup

The resolution of the model is physically controlled but generally represented by a triangular mesh, with side lengths of 2–5 km, a minimum area of 1.8 km² and a maximum area of 3.2 km². The top and sides of the model are no-slip boundaries and are fixed, whereas the bottom of the model is a free-slip boundary in Model 1. The sides of the model are thermally insulated, and the top and bottom are held at constant temperatures of 0 and 950 °C, respectively. The initial temperature of the underplating magma is 1250 °C. The geothermal gradient was assumed to be 3 °C/100 m (Fig. 2b).

In Model 1, magma underplating at varying depths and a lithospheric fault are included to simulate the ascent of magma originating from subduction. In contrast, Model 2 has similar boundary conditions and is assigned a 2 × 10⁹ Pa horizontal force on the right side, representing the compression effect induced by subduction. Model 2 includes a boundary force to simulate magmatism influenced by compression during subduction. The numerical experiments were conducted utilizing the finite-element software COMSOL Multiphysics, which is accessible at https://www.comsol.com (last access: 10 August 2024). The lithospheric thickness was estimated based on the current thickness derived from P-wave velocity measurements (Deng et al., 2019). The density and viscosity of magma were coupled to the thermal model and allowed to vary according to the wet quartz flow law (Chapman, 2021).

4 Underplating magma and circulation

The modeling results elucidated the temporal evolution and migration pathway of magmatism originating from the subsurface. In Model 1 (Fig. 3a–e), underplating magma accumulates at the lithospheric fault and locations of initial thickened magma along the bottom. In Model 1, after 10 Myr the bottom magma was divided into five upwelling magma bodies: 1, 2, 3, 5, and 6. Magma body 1 corresponds to the fault location; bodies 1, 2, and 5 correspond to pre-existing thickened magma bodies; body 4 corresponds to relatively shallow pre-existing magma; and body 3 is located below body 4. At this stage, bodies 1, 2, and 5 are uplifted higher than bodies 3 and 6. As evolution progresses, body 3 reaches the position of body 4 and rapidly achieves heights comparable to those of bodies 2 and 5 (Fig. 3e), whereas body 6 is significantly less uplifted. In Model 2, under compressive stress, four rising magma bodies formed at the bottom. No rising magma bodies formed at the bottom right. The formed magma bodies tilted and ascended to the left (Fig. 3f–j).

Figure 3 Results of Models 1 and 2, illustrating magma upwelling for 10–50 Myr, respectively. (a–e) Underplating magma rising, forming five magma bodies of varying heights in Model 1. (f–j) Underplating magma tilting under right-sided compression in Model 2. The crust, mantle, and magma materials are modeled as phases (fluids) on a dimensionless scale, with values of 0, 0.5, and 1, respectively. The contours denote the flow distribution of mantle fluid.

The magma rises upwards, generating mush-like features in the mantle, and grows laterally beneath the location of initial thickened and shallowly assigned magma (Fig. 3). The ascent of magma follows an up-and-down circulation pattern, driving lower crustal thickening and subsidence in front of the growing magma mush. The underplating magma ascends rapidly from the bottom to the lower crust, crossing a distance of approximately 70–80 km within 20 Myr through the lithospheric fault (Fig. 3), whereas without a fault the magma ascends a distance of only 50–60 km. In Model 2 (Fig. 3f–j), under the influence of compression from the right boundary, the diapiric magma migrates slightly to the left, causing thinning of the crust to the right.

The migration of magma and its temporal pattern are significantly influenced by lithospheric viscosity and temperature (Chapman, 2021). In our models, we specifically examine the effects of pre-existing magma and a lithospheric fault. The models provide cooling histories for five diapiric underplating magmas, which are then compared with observed coastal magmatism (Fig. 4a). The results indicate a more rapid cooling rate for the ascending magma through the fault (∼ 35 °C Myr⁻¹) and a slower cooling rate for underplating magma (∼10 °C Myr⁻¹) (Fig. 4a). The cooling history of other magmas is inconsistent with that of actual igneous rocks, emphasizing the influence of faults in this region. Moreover, the ascent of magma to the location of pre-existing magma is also faster than that without pre-existing magma (e.g., magma 3 and magma 6 in Fig. 3a–e).

Figure 4 (a) Comparison of the observed cooling histories of the CNB magmatic plutons (data from Chen et al., 2020) with time–temperature paths generated by models of rising magma. The green shadow represents the actual cooling age of the CNB magmatic plutons. (b) Zircon Hf isotopes and ages of coastal magmatic rocks in the SCB (data from Li et al., 2023). (c) Sketch illustrating the formation stage of the underplating magma and tectonic background during 80–110, 110–130, and 130–160 Ma. The olive shadow crossing sections (a) and (b) denotes the age range of 110–130 Ma.

5 Discussion

The εHf(t) values of the late Mesozoic igneous rocks along the coastal South China Block (SCB) tend to increase, peaking at positive values during 110–130 Ma (Fig. 4b). These values suggest that these rocks were derived from depleted mantle, subducted sediment-derived melt, and melting crust (Zhao et al., 2015). This period corresponds to the magmatic lull in the SCB, which coincided with a compression phase in the CNB during 130–105 Ma (Wei et al., 2023). Previous researchers attributed this phase to a transitional stage in subduction involving slab foundering, break-off, or steepening (Xu et al., 2023). Intriguingly, they demonstrate the potential for producing underplating magma beneath the lithospheric mantle, resulting in different compositions rising into the crust. Therefore, our models adopt underplating magma to simulate upwards magmatic intrusion and cooling processes, irrespective of the formational background. The rheological structure inevitably influences the pathway of underplating magma as it traverses the thick lithospheric mantle, ultimately dictating the distribution of widespread magmatism (Fig. 3). High viscosities in the lithospheric mantle may limit magma transport, with the effective viscosity of the upper mantle estimated at approximately 1 × 10²⁰–1 × 10²² Pa s in continental China (Shi and Cao, 2008), which decreases with increasing temperature. Considering the widespread magmatism and geothermal activity in the coastal SCB during the Cretaceous, an effective viscosity of 1 × 10²¹ Pa s for the lithospheric mantle is plausible.

It was suggested that once the shear zone network went through, the shear strength of the lithospheric mantle drastically decreased (Vissers et al., 1995). The CNB indicates long-term tectonic shearing activity, which would result in a lower viscosity than that in the interior of the SCB. Therefore, we assigned a low viscosity of 1 × 10¹⁹ Pa s to represent the active fault in the model. If the fault only changes in viscosity and the bottom magma is not thickened, then the magma at the bottom will not upflow. Shear zones control the ascent and emplacement of magmas (Weinberg et al., 2004). This implies that the shear zone should be a pathway for thermal fluid and have higher thermal expansivity (Afonso et al., 2005). Therefore, the weak CNB compared with that of the interior SCB facilitated the emplacement of mantle magma during the magmatic lull.

The rise of underplating magma into the middle crust, with an ascent of 70–80 km, takes 20–25 Myr through a lower-viscosity lithospheric fault, while it takes more than 40–50 Myr without a fault. This time discrepancy aligns with the magmatic lull in the coastal SCB. Thus, magmatism with different ages in the Cretaceous coastal SCB potentially formed through the exploitation of distinct ascent pathways under the same subduction conditions, rather than contemporaneously varying with steepening subduction geometry. Partially molten magma can persist for at least 25 Myr at temperatures exceeding 700 °C (Cavalcante et al., 2018), contributing to the heterogeneous and mixed magmatism observed in the coastal SCB.

The ascent of magma generates a mush-like head, accommodating the rheological structure of the lithospheric mantle and leading to magmatic circulation. These magmas underplate beneath the crust, decreasing the crustal thickness at the head and causing crustal subsidence on both sides. Importantly, pre-existing magma can accelerate the emplacement of underplating magma (Fig. 3). The underplating magma beneath the lithospheric mantle ascends rapidly when pre-existing magma is present. It is possible that multiphase magmatism increases the geothermal gradient in the SCB, reducing lithospheric viscosity and further promoting the ascent of underplating magma and the occurrence of a subsequent magmatic flare-up. In addition, the ascent pathways of magma change under the influence of a boundary force, resulting in increased transport time and delayed magmatic emplacement into the crust (Fig. 4c). Continued compression also contributes to the uplift of the lithospheric mantle, which is associated with the removal of crust, thereby decreasing the crustal thickness (Fig. 4c). This provides a new perspective on the crustal thinning of the coastal regions during subduction.

The deep structure of the late Mesozoic SCB is poorly constrained, resulting in speculative assumptions about key parameters such as fault depth and magma thickness. The model might oversimplify the complex geological features, potentially leading to inaccurate results. Additionally, assuming a uniform crustal thickness may not capture the true variability of the crust. The geometry of the lithospheric faults in the model is simplified, and important details that could affect magmatic processes are neglected. Variations in interpretations of fault characteristics and magma properties contribute to the idealized nature of our model results, which may differ significantly from actual geological conditions. This study aims to use these idealized scenarios to illustrate the complexity and diverse interpretations of magma evolution processes along the southern Chinese coast. Moreover, given the model's two-dimensional nature, it simplifies the intricate three-dimensional processes that likely influence magmatic evolution. Future research will focus on addressing these complexities to provide a more comprehensive understanding.

6 Conclusions

The model results describe the pathways and time spans of underplating magma rising into the crust under the influence of a lithospheric fault, pre-existing magma, and boundary stress. Magmatic flare-ups or lulls are not controlled solely by the slab subduction conditions. The Cretaceous magmatism along the coastal SCB could have occurred under the same subduction conditions, with the CNB facilitating the upwelling and intrusion of underplating magma under various regional stresses. A boundary force delays the ascent of underplating magma, whereas rising magma induces a significant circulation, which decreases the crustal thickness of the coastal SCB.