Seismic Imaging of a Magma Chamber and Melt Recharge of the Dormant Datong Volcanoes

Abstract

The Datong volcanoes erupted episodically from ca. 0.74 to ca 0.10 Ma with the last eruption at ca. 100 ka, and ever since has remained dormant. It remains elusive whether the Datong volcanoes are being recharged and would become active again in the near future. Here, we image the crustal structures beneath the Datong volcanoes using local earthquake data recorded at a newly deployed dense Datong Nodal Array. The high‐resolution 3‐D model reveals, for the first time, the geometry and melt fraction of a shallow crustal magma chamber at the depths of 10–20 km beneath the Datong volcanoes. Here, we suggest a magmatic system beneath the Datong volcanoes that could transport basaltic melt from the uppermost mantle to the shallow crustal magma chamber beneath the Datong volcanoes. Further, the melt‐poor crustal magma chamber (melt fraction of 5%–10%) can explain why the Datong volcanoes are inactive at present, and future eruptions can be expected if melt recharge from a mantle source is sustained.

Introduction

Volcanic activities often take place over several life-cycles. The repose stage of a volcano can vary from days to tens of thousands of years (Passarelli & Brodsky, 2012; Schmandt et al., 2019). However, what controls the episodicity of volcanic activities and how a dormant volcano is re-activated still remain enigmatic. Knowledge of magma chamber geometry and melt fraction is important for understanding such magmatic activities and monitoring the unrest of volcanoes. The Quaternary Datong volcanoes locates at the eastern margin of Ordos block, which preserve an Archean craton keel, and northern end of the Shanxi rift in the Northern Trans-North China Orogen (TNCO).

The Datong volcanoes, consisting of about 30 calderas, are one of the largest Quaternary intraplate volcanic groups in China (Chen et al., 1992) (Figure 1a), which can be divided into two different parts (Xu et al., 2005). The eastern Datong volcanoes, distributed along the Liulingshanqian fault (LLSQF, Figure 1b), are widely covered by basaltic lava flows with a thickness of ∼3–25 m. The volcanic rocks in the eastern group are mainly composed of tholeiitic basalts (Xu et al., 2005). The western Datong volcanic group, including at least 13 volcanic cones, are composed of volcanic clastic debris with lavas occurring at the base. Potassium-Argon dating suggests that the volcanism in the western Datong volcanoes began in the Late Pleistocene (∼0.4 Ma), which commenced later than the eastern group (Chen et al., 1992). Most of volcanic rocks in the western group are tholeiitic basalts. Xu et al. (2005) suggested that these alkali and tholeiitic basalts in the Datong volcanoes are the consequence of low-degree partial melting of upwelling asthenosphere and the reaction of asthenosphere-derived magma with the heterogeneous lithospheric mantle.

Figure 1. Geological map of the Datong Volcanic filed. (a) IMSZ, Inner Mongolia suture zone; WNCC, the western NCC; TNCO, the trans-North China Orogen; ENCC, the eastern NCC. Red combination logo: Datong volcanoes (DTV). The red rectangle shows the study area. Black lines represent tectonic lines. Red box in the inset shows the location of the study area. (b) Major geological features in our study area. YG-TZF, Yanggao-Tianzheng fault; KQF, Kouquan fault; CLSF, Cailiangshan fault; LLSQF, Liulingshanqian fault; SXR, Shanxi Rift; LLSM, Liulingshan Mountain; CLM, Cailiang Mountain; YMM, Yunmen Mountain. Purple dotted line separates the western and eastern DTV. Blue triangles: portable stations and yellow triangles: permanent stations. Black dots: earthquakes between 2010 and 2019; White star: 1989 Ms 6.1 earthquake. The red volcano logos indicate the Datong volcanoes; Black lines represent faults.

The Datong volcanic eruption took place within several life-cycles from 0.74 to 0.10 Ma with a relatively consistent repose interval of ∼0.2–0.3 Ma (Zhai et al., 2011), which is almost a longest end-member in terms of repose time. Although the Datong volcanoes are commonly known as dormant without any apparent volcanic activities at present, seismic images from surface wave (Li et al., 2018; Xu et al., 2020) and teleseismic body wave (Lei, 2012) tomography revealed a low velocity anomaly in the upper mantle beneath the Datong volcanoes. A recent magnetotelluric study revealed significant middle crustal conductors, which may indicate the presence of crustal magma bodies (Zhang et al., 2016). However, it is still not yet known what causes the regular eruption cycles and such a long stage of repose time in Datong volcanoes, and whether there still exists a crustal-uppermost mantle magmatic system that would transport melt from uppermost mantle to crustal magma chambers and fuel the next eruption in the future.

In this study, to address these questions, we deployed a dense nodal seismic array (DNArray) in the Datong volcanic field. Using local earthquake data recorded by this array, we construct a high-resolution 3-D crustal velocity model of Datong volcanoes. The unprecedented high-resolution 3-D model provides seismic evidence for the presence of partial melt and outlines detailed distributions of magma chambers beneath the currently dormant Datong volcanoes. These observations are crucial for deciphering the crustal magmatic system and understanding the mechanisms that control the life-cycles of volcanoes with long repose times, like Datong volcanoes.

Data and Methods

We collect seismic data recorded by a total of 165 seismic stations, including 13 permanent stations from China National Seismic Network and 152 newly deployed short period portable seismic stations from the DNArray (Figure 1b) that are deployed at two stages between September and December in 2018 and between June and July in 2019. The mean interstation spacing is about 5 km.

Detecting Earthquakes Using Match and Locate Method

We utilize the M&L (Match and Locate) method developed by Zhang and Wen (2015) to detect local earthquakes in the study area. The M&L involves computing the running cross-correlograms between the template waveforms and continuous waveform at each station to detect events and scanning potential locations of the detected events around the template events. Vertical components of continuous seismic data recorded by DNArray are first resampled to 100 Hz and band-pass filtered between 2 and 12 Hz. We then use 4 s waveform segments (1 s before and 3 s after the arrival time of P waves) of template events to cross-correlate with the continuous data. Finally, an event is detected when mean correlation coefficient (CC) and median absolute deviation (MAD) values of the stacked correlogram exceed the defined thresholds (CC ≥0.2 and MAD ≥10). A total of 28 local earthquakes (blue circles in Figure 2) from the China Earthquake Administration (CEA) catalog are used as templates in the M&L search. Several examples of seismic waveforms at different stations from a detected local earthquake are presented in Figure S1 in Supporting Information S1. Finally, we detect 1,075 local earthquakes (red circles in Figure 2) using the M&L method, in which a total of 702 detected events are recorded by at least five stations and are used in seismic tomography.

Figure 2. Map showing location of catalog earthquakes (blue circle) and detected earthquakes (red circle). Blue triangles represent the portable stations and yellow triangles represent the permanent stations. The detected earthquakes are about 30 times the catalog earthquakes.

Local Earthquake Tomography

We manually pick 20,358 P wave arrivals from all these 730 local events (702 new detected events and 28 catalog earthquakes). In addition, we also include local earthquakes between 2010 and 2019 reported by the CEA catalog in the tomography. The catalog included 14,323 P-wave arrival times from a total of 2,632 events that are recorded by 13 permanent CEA stations. We then plot all P-wave travel times according to epicentral distances and fit the data with a least squares linear trend. The outliers that fall outside 2SDs of the linear trend are discarded. Finally, we obtain 31,130 P-wave arrival times (Figure 3) for the local earthquake tomography.

Figure 3. P wave travel times as a function of distance. (a) All travel times from earthquakes including detected earthquakes and catalog earthquakes during 2010 and 2019. Blue lines represent the two standard deviations of travel times. (b) Travel times used in this study after outliers discarded.

The double-difference earthquake location algorithm (HypoDD) (Waldhauser & Ellsworth, 2000) is used to relocate these local earthquakes before tomographic inversion, which is an efficient method to determine high-resolution hypocenter locations. The relocated local earthquakes are presented in Figure 4, most of which are distributed along the KQF (Kouquan fault) to the west of Datong volcanoes (Figure 4).

Figure 4. Three-dimensional visualization of the locations of earthquakes. Red dots denote earthquakes. (a) Before relocated. (b) After relocated.

We then invert all the selected P wave arrival times to obtain a 3-D P-wave velocity model of the Datong volcanoes. We adopt the Fast-Marching tomographic technique (de Kool et al., 2006; Rawlinson et al., 2006) to iteratively invert for the 3-D P-wave velocity model. The method simultaneously solves for source parameters (location and origin time) using the updated 3-D velocity model at each iteration. Fast marching tomography uses a grid-based Eikonal solver to trace wavefronts in a 3-D domain.

The study area is 180 × 200 km wide and vertically extends from the surface to 30 km depth. We parameterize our study area into discrete grid nodes with 2 km spacing in latitude and longitude and 1 km in depth, respectively. We also set the upper and lower bounds at −1.8 km depth (above the surface) and 30 km depth, respectively, to account for the boundary effect. We run six iterations to obtain the final velocity model. The ray-paths for all data used in the tomography is presented in Figure S2 in Supporting Information S1. After six iterations, the Root Mean Square (The RMS is mainly affected by the recovery of velocity model in this work) of travel-time residual reduces from 0.878 to 0.230 s, an 89% variance reduction (Figures 5, S3 and S4 in Supporting Information S1).

Figure 5. Travel time residual distribution before (a) and after (b) inversion.

Results Resolution Tests

We conduct checkerboard tests to evaluate the resolution of our 3-D P-wave velocity model. In the checkerboard resolution test, we use the same 1-D reference model as the earthquake relocation. We set the magnitude of the input velocity perturbations to ±0.5 km/s relative to the reference 1-D model. The size of checkerboard anomalies is 10 × 10 km laterally and 5 km in depth (Figure 6). In the tests, we use the same source-receiver geometry and regularization parameters as those used in the field data inversion. We add 0.25 s error of travel times to permanent stations and 0.1 s error of travel times to portable stations in the resolution tests. Figure 6 show the results of checkerboard resolution tests. The tests demonstrate that our tomographic model has reasonably good resolution that can resolve the observed low velocity body at depths of 0–20 km, although some smearing appears at some margins of the study area. More importantly, the amplitude of the input anomaly can also be well recovered at depths of 0–20 km beneath the Datong volcanoes, which is important for the melt fraction estimation discussed in the following section.

Figure 6. (a–f) Map view of checkerboard test at the depths of 0, 3, 8, 14, 20, and 26 km. Red boxes outlines the area for which we present our 3D model as shown in Figures 8–10. (g–j) Cross-section images of resolution test results. The sizes of checkerboard anomalies are 10 × 10 km laterally and are 5 km in depth. Velocity anomalies are recovered very well at depths above 20 km depth.

We also compare the results of checkerboard tests using ray paths solely recorded by permanent stations (Figure 7a) or solely by the portable array (Figure 7b). When solely using permanent stations, the upper crustal structure cannot be recovered and the velocity anomalies are significantly underestimated from inversion. The recovered models suggest that the inclusion of DNArray in the local earthquake tomography significantly improves the resolution in the upper crust (0–10 km depths).

Figure 7. Checkerboard tests based on data recorded solely by permanent CEA stations (a) and solely by the portable array (b). The recovered models suggest that the inclusion of DNArray in the tomography significantly improves the resolution at the upper crust (0–10 km depths). The meaning of the red boxes is the same as those in Figure 6.

Three-Dimensional Velocity Model

Our 3-D velocity model is concentrated in the high-resolution area from our checkerboard test (red boxes in Figures 6 and 7). Figure 8 shows horizontal slices of our 3-D P-wave velocity model at various depths (Figure S5 in Supporting Information S1 shows horizontal slices with consistent color bars). At 2 km depth, low velocities are widely distributed along the Shanxi rift, which is covered by relatively thick Quaternary deposits (Figure 8b). The mountainous areas, including the Cailiang Moutain (CLM) and Yunmen Mountain (YMM), are characterized by broad high velocities. A pipe-like, prominent high velocity anomaly (>6 km/s) is observed right at the cones of Datong volcanoes (Figure 8b).

Figure 8. (a) Blue lines show profiles of cross-sections shown in Figure 9. (b–f) Horizontal slices of the 3-D P wave velocity model at different depths.

At 6 km depths, broad high velocities are found beneath the Datong volcanoes and the mountainous areas, whereas low velocities are found in the southwestern Shanxi rift. However, at 10 km depth, low velocities emerge beneath the Datong volcanoes, and YMM (Figure 8d). In the upper-to-middle crust (Figure 8e), our tomographic model exhibits fork-shaped low velocities, with one branch extending in the northeastern direction and the other in the northwestern direction from the Datong volcanoes. High velocities are found beneath the southwestern Shanxi rift. At greater depth (18 km) (Figure 8f), broad low P-wave velocities are observed beneath the Datong volcanic cones, the Cailiang (CLM), and Yunmen Mountains (YMM).

Figure 9 presents three vertical cross-sections of our 3-D P-wave velocity model with their surficial profiles delineated in Figure 8a. All these vertical transects traverse the cones of Datong volcanoes. The local earthquakes within 5 km of each profile are also presented along profiles.

Figure 9. (a–c) P wave velocities along three cross-sections A–A′, B–B′, and C–C′. F1: LLSQF; Black dots: earthquakes within 5 km of each transect. (d) Theoretical compressional velocities of granite as a function of porosity filled with basaltic melt and fluid contents. Blue stars denote velocities of 5.0 and 5.5 km corresponding to melt fraction of 5%–10%.

Beneath the mountains and volcanic cones, there is almost no sedimentary cover; while the Shanxi rift, home to the volcanoes, is covered by thick sedimentary deposits. The sedimentary deposit is much thicker beneath the southern than that beneath the northeastern end of the Shanxi rift (Figures 9a–9c). A high velocity body is imaged beneath the volcanic cones at 0–5 km depths. Furthermore, a small low velocity body is shown at 5–10 km depths below this high velocity body.

The most prominent feature of the tomographic model is a large low velocity body at the depth range of 10–20 km beneath the Datong volcanoes. The low velocity body is outlined by a 5.5 km/s contour, and the center of the low velocity body has a velocity as low as ∼5 km/s. Furthermore, a high velocity layer with a thickness of 3–5 km is imaged in the upper-to-middle crust. The SE-NW trending A–A′ profile (Figure 9a) shows that the high velocity layer is connected with the high velocity at the LLSQF to the southeast. Furthermore, a strong low velocity body at depths of 10–20 km is also found ∼20–40 km away from the recently active Datong vents (0.1 Ma) (Figure 9a), extending across the LLSQF and is located roughly beneath the prior Datong vents (the eastern Datong volcanic group). The lateral distribution and depth extent of the low velocities beneath the Datong volcanoes generally agree with conductive bodies observed by MT survey (Zhang et al., 2016, Figure S6 in Supporting Information S1).

From the NE–SW B–B′ profile along the rift (Figure 9b), we observe that the spire-shaped low velocity is symmetrical with the summit right beneath the Datong volcanoes in the upper-to-middle crust (5–20 km depths). The geometry of the low velocity anomaly is more clearly depicted in the 3-D perspective plotted in Figure 10a. The NNW–SSE cross-section C–C′ (Figure 9c) reveals a high velocity anomaly beneath the LLSQF extending from ∼5 to ∼16 km depths.

Figure 10. (a) 3D visualization of the upper-to-middle crustal magma chamber beneath the Datong volcanoes. (b) S wave velocity perturbation down to a depth of 100 km from ambient noise and teleseismic surface wave tomography (data from Li et al., 2018), which shows the mid-to-lower crustal and uppermost mantle low velocity anomalies beneath the Datong volcano. Black line indicates Moho depth. (c) The same absolute versus model from surface wave tomography (Li et al., 2018) but down to 250 km depth, which highlights the mantle flow and mantle upwelling that finally feeds the Datong volcanoes. Black line indicates the depth of Moho and dash line indicates the depth of LAB.

To further illustrate that the upper-to-middle crustal low velocity body is a robust feature, we conduct a custom resolution test. We design a model with a conical frustum-shaped low velocity (5 km/s) at depths between 10 and 20 km covered by a high velocity (6 km/s) layer with 5 km thick. The recovered velocity model shows that main feature of the input velocity anomaly can be well recovered (Figure 11). The absolute value of low velocity can be well recovered in the western area although it is slightly underestimated in the west. To test whether clusters of events (Figures 2 and 4) bias the inversion results, we randomly select 10% events of this cluster (black box in Figure 2) to carry out the inversion. The resulting model velocity shown in Figure S7 in Supporting Information S1 is in good agreement with the model using all events, suggesting that the inversion is not biased by the uneven source distribution.

Figure 11. (a) A synthetic 3D model with a cone-shaped low velocity (5 km/s) at depths of 10–20 km covered by a high velocity (6 km/s) layer of 5 km thick. (b) Recovered model from our resolution test, showing that the main geometry as well as the amplitude of the input velocity model can be well recovered.

Discussions Shallow Crustal Magma Chamber and Melt Fraction

The resulting tomographic model reveals prominent low velocity anomalies (5–5.5 km/s) beneath the Datong volcanoes at depths of 10 and 20 km (Figure 9). The most common granitic-granodioritic rocks in the upper crust have a typical P-wave velocity of 5.83–5.92 km/s at a temperature range of 250–500°C and a pressure range of 0.2–0.48 GPa, which are the temperate and pressure values at ∼10–20 km depths (Christensen & Mooney, 1995). If we assume that velocity reduction beneath the Datong volcanoes is solely due to the elevated temperature, then the observed velocity reduction (0.33–0.42 km/s) in the middle crust beneath the Datong volcanoes would require a high temperature of over 1100°C, according to the granite temperature coefficient ((∂Vp)/(∂T)) of −0.39×10−3 km/s/°C (Christensen and Mooney (1995). Such an estimated temperature is significantly higher than the subsolidus temperature of the upper crust (∼800°C) (Sawyer et al., 2011). Thus, we suggest that the elevated temperature alone cannot explain the observed low P-wave velocity and therefore, partial melt is required to account for the upper-to-middle crustal low velocity beneath the Datong volcanoes. To estimate melt fraction beneath the Datong volcanoes (see details in the supplemental material and Figure 9d), we adopt the relationship between P-wave velocity and melt fraction proposed by Chu et al. (2010), which assumes an equilibrium state melt-saturated porous media. It should be noted that the melt-fraction relationship developed by Chu et al. (2010) was specific to the composition estimated for Yellowstone rhyolite. We thus modify the relationship by replacing rhyolite with basalt for Datong volcanic system and consider the effects of temperature beneath the Datong volcanoes on P-wave velocity (see the supplemental material). Based on this calculation, a low velocity anomaly of 5.5 km/s at depths of 10–20 km corresponds to a melt fraction of ∼5% (Figure 9d). Our estimate of melt fraction from seismic velocity is generally consistent with the recent estimate based on a MT study (Zhang et al., 2016), in which 10 Ωm corresponds to over 6% melt fraction. We suggest that the low velocity zone encompassed by 5.5 km/s contour in the upper-to-middle crust represents the shallow magma chamber beneath the Datong volcanoes which is now filled with at least over 5% melt fraction. In addition, the fact that the magma chamber becomes wider with depth (Figure 10a) supports the notion that magma commonly accumulates in the middle crust, where melt differentiation further takes place (Cashman et al., 2017; Dufek & Bachmann, 2010; Sparks et al., 1977). At the center of the magma chamber, P-wave velocity reaches as low as 5 km/s, implying at least 10% melt fraction in the center of the magma chamber (Figure 9d). Note that the resolution tests show that the absolute value of low velocity at 10–20 km depths is only slightly underestimated (∼5%), particularly in the eastern study area. This underestimation can only lead to a slight underestimate of melt fraction ∼2% in the magma chamber.

The major part of the magma chamber with a ∼5% melt fraction could be in a super-solidus state that lacks melt interconnectivity (Sparks et al., 2019), while only the inner part of the magma chamber with melt fraction >10% exhibits a melt-poor, mushy state with melt pockets interconnected. Thus, the observations of the melt-poor magma chamber (Figures 9 and 10a) could explain why the Datong volcanoes are currently lacking features of volcanic activities, such as high surface heat flow and active sulfide and carbon dioxide emissions.

Here, the estimated 5%–10% melt fraction of the shallow magma chamber beneath the Datong volcanoes is far below the melt fraction at the critical point of an eruption (∼40%) (Gelman et al., 2014; Parmigiani et al., 2014). In comparison with those of active volcanoes, like the Weishan volcano (15.3%, Gao et al., 2020) and Yellowstone (∼32%, Chu et al., 2010), the melt fraction of the Datong volcanoes is much smaller. Furthermore, our tomographic model shows that the magma chamber beneath the Datong volcanoes is covered by a high velocity layer (>6 km/s) with a thickness of 3–5 km. The high velocity layer (Figure 9) could be the crystallized igneous rock, which outlines the top boundary of the crustal magma chamber. In particular, the pipe-like high velocity right beneath the volcanic cones could be the frozen conduit that magma ascended through during past eruptions.

It has been suggested that earthquakes in the Datong volcanoes are mainly controlled by faults (e.g., Wang et al., 2002). The distribution of earthquakes, including Ms 6.1 main shock, mainly occurred at the high velocity anomalies at the depths of 5–20 km beneath the fault zone, and only a few earthquakes occurred at the low velocity anomalies in the magma chamber (Figure 9). We suggest that most of earthquake in this area, including the Ms 6.1 mainshock, are mainly controlled by faults (e.g., LLSQF), but a few earthquakes occurred at the low velocity anomalies may be related to the shallow magma chamber beneath the Datong volcanoes. However, an in-depth study of the controlling factors of these earthquakes occurring in the magma chamber, which is beyond the scope of this study, will be of great significance to the study of the magma system beneath the Datong volcanoes.

Upper Mantle Dynamics and Uppermost Mantle-Crustal Magmatic System

It is generally accepted that a shallow magma chamber is short-lived without melt recharge from a deeper source (Solano et al., 2012). Small magma bodies in the upper crust achieve thermal equilibrium within 1 ka and only those large igneous systems could have lifetimes up to 100 ka (Cashman et al., 2017). Given that the last eruption of the Datong volcanoes was ca. 100 ka (Zhai et al., 2011), it would be expected that the whole magma chamber is almost completely crystallized and exhibits as normal or even high velocities (Daniels et al., 2014; Flinders et al., 2018), like the high velocity layer that covers it as shown in Figure 9. Thus, it is unlikely that the melt in the upper-to-middle crustal magma chamber is the remnants of melt intrusions at ca. 100 ka. We thus suggest that the small melt fraction in the upper-to-middle crustal magma chamber is the consequence of a recent melt recharge from deep sources in the uppermost mantle. The episodic eruption of the Datong volcanoes with a relatively consistent interval (0.2–0.3 Ma) from the Oligocene to Pleistocene further implies episodic influx of upper mantle-derived melt. Basaltic melt supply from the upper mantle is mainly controlled by mantle dynamics as well as melt segregation and transportation in a shallow magmatic system (Cashman et al., 2017).

Recent seismic imaging using surface wave reveals an eastward dipping asthenospheric upwelling, which originates from a depth of ∼200 km and finally reaches to the uppermost mantle beneath the Datong volcanoes (Figure 10) (Li et al., 2018). Such localized asthenospheric mantle upwelling is likely to be induced by the shear-driven convection of an eastward mantle flow at the craton edge in response to the far-field effects of continuous Indian-Eurasian continents converge and retreat of Pacific slab (Guo et al., 2016; Li et al., 2018; Li et al., 2020). The asthenospheric upwelling may cause thermal erosion to the overlying continental lithosphere, resulting in a thinned lithosphere beneath the Datong volcano (Figure 10). Significant low S-wave velocity (>−6%) in the uppermost mantle (∼60 km) may imply the presence of partial melt. The melt supply rate from mantle is likely to be slow and episodic. Our study further demonstrates that basaltic melt from the uppermost mantle finally ponds in a shallow magma chamber between 10 and 20 km depths beneath the Datong volcano.

Thus, on the basis of these observations, we suggest a magmatic system beneath the Datong volcanoes that may transport basaltic melt from the uppermost mantle to the shallow crustal magma chamber as shown in Figure 10. The complex and deep-seated faults systems, like the LLSQF and CLXF, in the Shanxi rift may control the melt migration in the crust and the surface distribution of volcanic activities, like the migration of eruptions between eastern and western volcanic groups. Thus, it is possible that the Datong volcanoes are now at a stage between two eruption cycles, and future eruptions could be expected if melt recharge from mantle source can be sustained. On the other hand, the most significant low velocity bodies in the upper crust locate roughly beneath the eastern Datong volcanic vents, suggesting that the next stage of volcanic activities is likely to occur in the eastern volcanic group.

In summary, the dense seismic array covering the entire Datong volcanic region provides us with an unprecedented opportunity to decipher the geometry and melt fraction of a shallow magmatic system beneath the Datong volcanoes (Figure 10). Here, we have observed a magmatic system beneath the Datong volcanoes that relate upper mantle partial melting zone to the shallow magma chamber. The melt-poor crustal magma chamber and small volume melt recharge by the deep source in the uppermost mantle beneath the Datong volcanoes since its last eruption should provide important understanding of the life-cycles of intraplate volcanoes such as the Datong volcanoes.

Conclusions

In this study, we have obtained an unprecedented high resolution 3-D velocity model of the upper and middle crust beneath the Datong volcanoes using seismic data from a dense seismic array. Our new velocity model helps better illuminate the crustal magmatic system beneath the Datong volcanoes. We observed a shallow crustal magma chamber at the depths of 10–20 km beneath the Datong volcanoes with a melt fraction of 5%–10%, and the small melt fraction of the magma chamber could explain why the Datong volcanoes are inactive at present. This magma chamber is covered by a high velocity layer with a thickness of about 5 km, which could be the frozen conduit that magma ascended through during previous eruptions. Given that the last eruption of the Datong volcanoes was at ca. 100 ka, we suggest that the small melt fraction in the upper-to-middle crustal magma chamber is the consequence of a recent melt recharge from deep sources in the uppermost mantle. In summary, we suggest a magmatic system within beneath the Datong volcanoes that could transport basaltic melt from the uppermost mantle to the shallow crustal magma chamber, and the complex and deep-seated faults systems in the Shanxi rift may control the melt migration in the crust and the surface distribution of volcanic activities.

Acknowledgments

The authors thank all the people who participated in the field deployment of seismic arrays for the DNArray. This work is supported by the NSFC grants (41974052, 42122027 and 41774052). The authors thank Editor Dr. Paolo Diviacco and two reviewers Dr. YoungHee Kim and Dr. Pieter-Ewald Share for their constructive comments and suggestions that greatly improved the manuscript.

Data Availability Statement

Data are available from https://zenodo.org/record/4680179.