Introduction

Northeast (NE) China contains a large number of intraplate volcanic clusters. Since these volcanoes are far from plate boundaries (Fig. 1a), compared to the existing abundant research and mature theoretical basis of plate boundary volcanic systems, the magmatic plumbing systems and the formation mechanisms of these intraplate volcanoes are still under debate. The Changbaishan volcanic system (CVS) is the largest intraplate volcanic field in NE China, and it is also the most potentially dangerous active volcanic system in all of China, as it has erupted numerous times throughout recorded history1. The broad CVS includes the Changbaishan Tianchi Volcano (CTV) (also known as the Paektu/Baekdu Volcano in North Korea), Wangtiane Volcano (WTEV), and Longgang Volcano (LGV), which contains the largest volume of Quaternary volcanic rocks in China2.

The CTV is a giant composite polygenetic volcano of the CVS, and it produced the Millenium Eruption (AD 946–947), which is one of the largest known post-caldera eruptions on Earth in the past 2000 years1. The CTV is characterized by a large volcanic cone, which is also the highest peak in NE China with an altitude of more than 2500 m (Fig. 1c).

The LGV, in the NW of the CVS, comprises more than 160 volcanoes, including monogenetic trachybasaltic scoria cones, and around a dozen phreatomagmatic centres3. The latest basaltic volcano of LGV, Jinlongdingzi volcano erupted approximately 1700 years ago, indicating the possibility of sustained magma recharge beneath the LGV4. Compared to the solitary giant volcanic cones formed by composite polygenetic volcanoes, monogenetic volcanic activity is characterized by each eruption occurring at a different, unique location through different feeder structures, eventually forming a series of densely distributed volcanic cones. At present, more than 160 maars, volcanic cones, and crater lakes have been documented in the LGV area, which consequently is one of the most eruptively dense volcanic landscapes in China2,4,5, (Fig. 1b).

Previous petrology and geochemistry results6,7,8 suggest that LGV and CTV have a common source at mantle depths. The parental magmas are potassic-trachybasaltic in composition, and the interpretation of a common deep source is supported by previous seismic tomography studies9,10 and by the Magnetotelluric (MT) study11. These studies often detect electrical structure and velocity information below the Moho surface to depths of several hundred kilometers, showing that LGV and CTV share the same low-resistivity and low-speed channel below. However, petrological and geochemical data also suggest that these two volcanoes have undergone different evolutionary processes. CTV is dominated by centralized shield-building potassic-trachybasalt volcanism, trachyte cone building, and late alkaline acid magmatic eruptions (alkali flows and alkaline rhyolite eruptions)5,7,8, which indicates that there are obvious stages of crustal magma chamber evolution in the ejecta. In contrast, LGV contains solely trachybasaltic eruptions, and isotopic ratios and major and trace element data indicate that the erupted magma was derived from a mantle source with minimal crustal contamination5, the ejecta did not go through the stage of crustal magma chamber. Although the CTV and LGV share the same mantle magma source and are within a 100 km distance from each other, they host completely different eruptive styles and products.

Extensive research has been conducted on the lithospheric structure of the two volcanoes, especially the CTV (e.g10,12,13,14,15,16,17). However, due to the differences in resolution, research scope, and scientific issues focused on by various methods, there is currently no consensus for why the two volcanoes have experienced completely different eruption styles. Geochemical data show that intraplate volcanism in NE China is likely affected by processes occurring both in the crust and in the upper mantle2,5. Three-dimensional (3D) MT imaging can effectively reveal the detailed electrical resistivity structure of both the crust and the upper mantle. The MT method is often used to image the subsurface beneath a volcano to identify contemporary melt zones, magma reservoirs and/or crystal mushes18. We therefore collected high-density MT data over these two volcanoes (Fig. 1d). Based on the 3D inversion of these data, we here explore the physical factors that may control the different eruptive styles of these two volcanoes in a unified reference system, constrained by other petrological, geochemical, and geophysical results.

Fig. 1
figure 1

(a) Map showing the location of the study area in a regional setting; the red line indicates the boundary between the Pacific Plate and Eurasian Plate. The base map is quoted from the publicly available Google satellite imagery maps (https://www.google.cn/intl/zh-en/earth/.).(b) Aerial view of LGV (original image modified from previous study3). The red arrows indicate the location of the volcanic cones. (c) Aerial view of CTV (original image modified from previous study3). (d) Locations of the MT stations over LGV and CTV with regional topography around the volcanoes. The red triangles are volcanic cones, the black solid lines are faults, and the blue shaded areas show the locations of Quaternary basalt distribution in the study area6, surface topography from ETOPO1 model (https://www.ngdc.noaa.gov/mgg/global/). The figures are created using Generic Mapping Tools (GMT)19 (v.5.1, https://www.generic-mapping-toolsorg/).

Full size image

Magnetotelluric data

The MT method is one of the most efficient and cost-effective deep-probing geophysical techniques. In recent years, 3D inversion of MT data has been widely applied in volcanic regions around the world (e.g20,21,22,23,24). Figure 1d shows the locations of all MT stations used in this study. The blue triangles represent the latest data collected in 2022, with a total of 68 stations. The black triangles represent data collected in 202025, with a total of 99 stations. The yellow triangles represent data collected in 199712, with a total of 58 stations. The MT stations are distributed within an area of 200 km in the East-West direction and 90 km in the North-South direction, with a total of 225 stations. This study makes use of this most comprehensive and dense MT data obtained in the region to date. The resulting apparent resistivity and phase curves for all stations are included in the Supporting Information, Figures S3–S5, and they show good data quality.

The phase rotational invariant φ2 obtained by phase tensor decomposition26,27,28 can be used to analyze the electrical properties of the subsurface as a function of frequency. Figure 2 shows the phase tensor ellipses at 9 periods from our 225 MT stations. The color of the ellipses represents the geometric mean of the maximum and minimum phase differences between the magnetic and electric field, which is an indicator of the vertical conductivity gradient. As shown in Fig. 2, at shorter periods (< 0.1 s), low φ2 values (< 45) occur in much of the study area, especially in the LGV area, while high φ2 values (> 45) are observed almost exclusively in the CTV area. This pattern indicates that the study area is mainly characterized by relatively uniform high resistivity values at shallow depths, however the CTV requires decrease in resistivity. As the period increases (> 0.1 s), high values of φ2 values are observed in the southwest and central parts of the study area, suggesting decrease in resistivity with depth.

The results of dimensionality analysis (Figs. S6, S7 in the Supporting Information) indicate that the study area is structurally quite complex, with dimensionality ranging from 1D to 3D across the period range of our data. The resistivity structure beneath the survey area is 3D at longer periods, and consequently the observed MT data need to be inverted in 3D.

Fig. 2
figure 2

Phase tensor maps of LGV and CTV at various periods. The figures are created using Generic Mapping Tools (GMT)19 (v.5.1, https://www.generic-mapping-toolsorg/).

Full size image

3D inversion and results

3D inversions were performed using the parallelized nonlinear conjugate gradient algorithm of ModEM29,30. The 3D mesh comprises 66 × 90 × 107 cells in the north-south (X), east-west (Y) and vertical (Z) directions, respectively. The central area of the mesh was 46 × 70 × 107 cells, with a nominal horizontal cell size of 2.5 km (X) × 3.0 km (Y). This portion of the mesh was extended by 10 padding cells in each horizontal direction, with each cell width increasing by a factor of 1.5. In the Z direction, the mesh includes 63 layers that start with a thickness of 20 m at 0 m altitude and increase in thickness by a factor of 1.15 to a total depth of 800 km. To take into account topography, the mesh also includes 44 terrain layers with a constant 50 m thickness above an altitude of 0 m, as detailed in Figures S8–S10 (Supporting Information).

The 3D inversions were run with full impedance tensor elements (Zxx, Zxy, Zyx and Zyy). Error floors were set to an absolute value of 5% and 10% of (|Zxy × Zyx|)½ for the off-diagonal and the diagonal components of the impedance tensor, respectively. The starting model resistivity was set uniformly to 100 Ωm. A ModEM model covariance setting of 0.3 was used, and the ModEM lambda parameter (defining the trade-off between model and data norms) was initially set to 5000 and decreased by a factor of 10 through the NLCG search algorithm. After 86 iterations, a total root mean square (RMS) misfit of 2.11 was achieved. The 3D inversion results are summarized in Figures S11–S15 and Table S1 (Supporting Information).

Our resulting electrical resistivity model is shown in Figs. 3, 4. Figure 3 shows horizontal depth slices of our resistivity model (these depth data values are calculated based on an elevation of 0). Figure 4 shows cross-sections along the profiles L1, L2, L3, L4, L5 and L6 and the profile locations are shown in Fig. 4a. L1 and L2 cross the core of our study area, including both the LGV and CTV; L3 runs through CTV and the northwest of the study area; L4 crosses the crater location of CTV in a southwest-northeast direction; and L5 runs through CTV from south to north. Topography is also plotted above each profile, and the positions of LGV and CTV are marked. These models display the electrical resistivity structure of the crust, lithospheric mantle, and part of the asthenosphere mantle in our study area.

The electrical resistivity information revealed by the MT inversion provides insights into magma plumbing system beneath the LGV and CTV. As shown in Fig. 3a, from the surface to a depth of 16 km, LGV and CTV exhibit significanlty different electrical characteristics. While the upper crust beneath the LGV shows relatively continuous high resistivity values (labeled as R), at 1.5 km beneath CTV, there is a large low-resistivity structure, which we label as C1. Resolution tests prove the robustness of low-resistivity structure C1 beneath CTV (Figs. S24–S26 in the Supporting Information). This low-resistivity structure expands laterally with depth and exhibits a planar distribution, covering an area of up to 20 × 20 km, and it shows a trend of increase with depth. This is similar to the results of recent MT surveys conducted on two volcanoes (Kuju volcanoes, and Unzen volcano) in Japan, where magma can migrate along structural boundaries within a relatively high permeability zone24,31. Structure C1 gradually disappears at a depth of approximately 16 km. From the L4 and L5 profiles, the low-resistivity structure C1 is mainly located beneath the CTV volcanic cone and along its northeast side.

The cross-sectional views in Fig. 4 more clearly display the differences in the upper crust beneath the LGV and CTV. Profiles L1 and L2 show that beneath the LGV, down to the Moho depth, there is generally a continuous high-resistivity zone, with resistivity values generally exceeding 1000 Ωm. There are small, highly localized low-resistivity structures in the middle and lower crust, which we label as C2, C3, and C4. C2 is connected to the low-resistivity structure in the mantle. C3 is in the western portion of our study area, and it may be related to the eruption of the Jinlongdingzi volcano before 1700 years ago; refer to previous study25 for discussion of this possibility. C4 is isolated and exists on the north side of the LGV, without showing a deep connection, and it may represent residual magma from early eruptions. At greater depths, the lower crust in the study area has a heterogeneous low-resistivity structure, and below the Moho depth a quasi-continuous low-resistivity zone gradually appears in the central area of the study area. This low-resistivity structure, which we label as C5, extends downward across the LAB. Resolution tests demonstrate that the MT results constrain structures below 60 km and therefore validate the existence of C5 low-resistivity structure (Figs. S26–S34 in the Supporting Information).

Fig. 3
figure 3

Horizontal slices of the preferred 3D resistivity model at various depths. The white triangles are MT stations, the red triangles are volcanic cones, the black solid lines are faults. The figures are created using Generic Mapping Tools (GMT)19 (v.5.1, https://www.generic-mapping-toolsorg/).

Full size image
Fig. 4
figure 4

Vertical slices of the preferred 3D resistivity model at various cross-sections. (a) Schematic diagram of different cross-sections position, Moho and Lithosphere-Asthenosphere Boundary are from previous study15. (b) Extracted 2D cross-sections from the 3D resistivity model. (c) Geothermal gradient map presented in L2 was created using the data from previous study32. The figures are created using Generic Mapping Tools (GMT)19 (v.5.1, https://www.generic-mapping-toolsorg/).

Full size image

Discussion

The low-resistivity features revealed by MT often correspond to partial melt, which consequently indicate the presence of magma bodies in volcanic areas. The resistivity values within C1 below the CTV are between 1 and 10 Ωm, which would be consistent with a shallow magma system within the upper crust (e.g20,21,22,23,24). Figure 4c presents a cross-sectional view of a geothermal gradient estimate31 (based on surface heat flow data) along the profile L2, which enables a comparison of the geothermal conditions beneath the CTV and LGV. Although this estimated temperature field cannot accurately display the local temperature characteristics below the CTV, it does regionally demonstrate the high temperatures in the crust beneath the CTV, which are very different from the lower crustal temperatures beneath the LGV. Based on borehole observations, the geothermal gradient beneath the CTV exceeds 150° C/km, so the temperature of the C1 exceeds 1000° C at a depth of approximately 7.5 km32. Such a high temperature cannot be explained by fluid advection alone. Instead, we take these high temperatures and our modeled high conductivity values to estimate the melt fraction. Taking the temperature and pressure conditions at a depth of 7.5 km as an example, we calculate the resistivity of pure melt using SIGMELTS33. Geochemical data show that the average SiO2 and Na2O contents in the rocks from the latest eruption of CTV is 65 wt% and 5 wt%, respectively, and the water content is 2 wt%7. The resistivity of the pure melt is therefore ~ 1.2 Ωm. We assume that, in the absence of water and partial melting, the host rock in the upper crust of our study area has a resistivity value of approximately 1000 Ωm. Using the Hashin-Shtrikman upper bound34, we estimate that approximately 12% melt would be required to explain our observed resistivity values beneath the CTV. Given that the magma has a high volatile content7,8, which can further reduce the resistivity value of the pure melt, this melt fraction may be overestimated to some extent. Regardless, we consider this calculated high melt fraction to be reasonable.

Different disciplines have conducted detailed research on the deep magmatic reservoirs beneath the CTV. Although there are still many inconsistencies in these previous results, it has been determined that there are locally high temperatures, low seismic velocities, low resistivities, low densities, and bodies of partial melt within the crust beneath the CTV (e.g10,12,13,14,15,16,17,35). These previously reported anomalous zones correlate well with the location of our structure C1. Based on these correlations and our calculations above, we propose that C1 represents a shallow magma chamber in the upper crust of the CTV. The magma chamber is located directly below and to the northeast of the CTV, with a horizontal distribution of approximately 20 × 20 km. Zhu15 inferred that the magma chamber in the crust of the CTV extends laterally more than 100 km based on a seismic low-velocity zone. Yi8 challenged this conclusion, based on the distribution of volcanic units, and instead proposed that the lateral extent of the CTV magma chamber does not exceed 10 km. Our results support this latter conclusion. The large zone proposed by Zhu15 may be the result of limited seismic stations and insufficient resolution, resulting in an overestimate of the extension of the low-velocity body. In fact, the latest high-resolution seismic results from the CTV constrain the low-velocity region to an area similar to that observed in our MT results36.

The structural characteristics of the nearly continuous high-resistivity body in the upper crust beneath the LGV. The geothermal gradient map32 also shows that the upper crust of the LGV is relatively cold. Geochemical data show that the LGV eruptions have limited crustal contamination5. MT modelling results from Newberry Volcano in USA and Unzen Volcano in Japan identified highly crystalline magma reservoirs with high electrical resistivity values31,37. Our 3D MT data inversion results confirm the absence of an upper-crustal magma chamber (or chambers) beneath the LGV that would facilitate assimilation of crustal material into the eruptible magma, and our resistivity images support an interpretation that Quaternary magmas in the LGV did not experience assimilation and fractional crystallization within shallow crustal magma chambers. There are potentially smaller magma bodies, represented by C2, C3, C4, in the lower crust of the LGV, which may represent mantle-derived melt. However, compared to the LGV, there are much more complex evolution processes operating within the upper crust during the incubation and eruption of magmas in the CTV6. The LGV lacks a magma chamber that would enable such magmatic evolution process in the upper crust, and we consider this to be a key reason for the different eruption style and chemistry between the two volcanoes.

Based on a comprehensive understanding of petrology in our study area7,8, the upper crustal magma chamber of the CTV represents a potential source of volcanic eruptions over the past few thousand years. The CTV experienced a concentrated seismic event between 2002 and 20053,5, with earthquakes mainly concentrated at a depth of 7.5 km below the CTV. We speculate that the shallow seismic event that occurred in the CTV between 2002 and 2005 may be related to the crustal magma chamber C1, which can also be considered as an evidence of continued activity of the CTV.

Petrological and geochemical data suggest that the sources of the LGV and CTV are both potassic-trachybasaltic magmas derived from the mantle5,7,8. Our MT results show a low-resistivity structure C5 extending to the asthenospheric mantle depths (> 60 km) between LGV and CTV. Jones38 conducted a detailed discussion on the reasons for the existence of low resistance bodies in the mantle, proposing hypotheses related to magma or sulfide mineralization, as there were no clear reports on sulfides in the study area. Our results support that magmas of the LGV and CTV may be derived from the asthenospheric mantle, and we consequently propose that this structure C5 represents a mantle-rooted magmatic system that feeds both the LGV and the CTV.

Integrating existing petrological, geochemical, and geophysical results (e.g10,12,13,14,15,16,17,35). , as shown in Fig. 5, we propose a mechanism of formation for the LGV and CTV. These two volcanoes may share a common deep asthenospheric mantle source. However, in the upper crust, there is no shallow magma chamber in the LGV; instead, magma directly rises from the asthenosphere, due to the decompression melting, and passes through the crust to erupt with little to no crustal residence time. Uncontaminated potassic-trachybasalt erupts at the surface to produce a series of densely distributed monogenetic volcanic cones, and the activity of the LGV consequently exhibits uncertain eruption location and small-scale eruption. In contrast, the CTV has an obvious magma chamber in the upper crust. Potassic-trachybasaltic magma derived from the asthenospheric mantle resides in the upper crustal magma chamber for an appreciable amount of time, and the magma in the crust undergoes assimilation and fractional crystallization processes. Magma mixing due to the continuous injection of melt from the mantle into the crustal magma chamber could lead to further evolution of the magma body as well as trigger eruptions due to disturbance of the volatile equilibrium within the chamber. Petrologic and geochemical data support the existence of a shallow crustal magma chamber under the CTV, and the Millenium Eruption demonstrates these sorts of volcanic eruption process7,8. Because there is a shallow crustal magma chamber in the CTV, post-caldera eruption activity is more concentrated. Compared with the densely distributed volcanic cones formed by LGV (Fig. 1c), the magmatic plumbing system revealed in our MT images explains why the CTV has formed a large-scale, single volcanic cone (Fig. 1b).

Fig. 5
figure 5

Proposed conceptual model summarizing eruption mechanisms of LGV and CTV. Generic Mapping Tools (GMT)19 (v.5.1, https://www.generic-mapping-tools.org/) and CorelDraw Graphics Suite 2019 commercial software (https://www.coreldraw.com/cn) were used to create this figure.

Full size image

Conclusions

By performing 3D inversion of an MT dataset that completely covers the LGV and CTV, we have obtained high-resolution electrical resistivity images of the crust and upper mantle in this region. The results reveal that the two volcanoes have distinct magmatic plumbing systems, and this is likely the reason for their different eruptive styles. There is no obvious magma chamber in the shallow crust of LGV. The magma from the mantle rapidly rises through the crust to produce a series of densely distributed volcanic cones. In contrast, there is a magma chamber in the upper crust of the CTV, where the fractional crystallization and mixing of magma has occurred. This magma chamber has facilitated multiple centralized eruptions, and thereby has led to the formation of the large CTV volcanic cone. Our MT results indicate that the differences in magmatic plumbing systems through the crust have controlled the different eruptive styles and geochemistry of the LGV and CTV in NE China.