Deep Electrical Structure of the Qilian Orogenic Belt with Dynamic Implications for the Northeastern Margin of the Tibetan Plateau: Revealed by 3D Magnetotelluric Inversion Using Unstructured Tetrahedral Elements

The northeastern margin of the Tibetan Plateau is distant from the Indian-Asian collision front. This margin consists primarily of the Qaidam, the Qilian orogenic belt, and Alxa block (a part of the North China Craton) from south to north. The area extending from the south of Tianshan to the north of Kunlun Mountain (e.g., the Tarim and Qaidam basins, Qilian Shan, and Alxa block) constituted a plate known as “the western plate” in the Late Proterozoic [1, 2]. NE-SW extension led to the opening of the North Qilian Ocean during the Neoproterozoic to Mid-Late Cambrian. The North Qilian Ocean underwent northward subduction and closed along the North Qilian suture during the Late Ordovician. The collision of the Qilian-Qaidam and Alxa blocks occurred due to the compressive stress regime associated with the closure of the North Qilian Ocean, and the North Qilian subsequently entered an orogenic stage. The surface expression of the Qilian orogenic belt was strongly denuded after the Late Devonian owing to continuous extension following collisional orogenesis. The Paleozoic orogeny in the Qilian area ceased completely in the Carboniferous. The Indian and Asian plates began to collide during the Cenozoic, with the strong compression of the Indian plate activating extensive orogenesis [3, 4].

The Tibetan Plateau has undergone an enormous amount of uplift since the Cenozoic [5, 6]. The outward growth and expansion of the plateau have been modeled using several tectonic mechanisms, such as continuous shortening [7–9], crustal channel flow [10, 11], and rigid plate compression [5]. The Qilian orogenic belt, which is located in the front-most zone of outward expansion of the Tibetan Plateau, can be divided into three parts: the complex North Qilian, a transitional Middle Qilian, and the relatively simple, high-resistivity South Qilian [12–14]. Several irregularly shaped, unconnected “fragments” of extremely high-resistivity blocks in the mid-upper crust and a stratified low-resistivity layer in the mid-lower crust comprise the electrical structure of the crust in this area [15, 16]. Gao and Cheng [17] defined Kuantan Shan-Hei Shan fault (KHF) as the northern boundary thrust (NBT), which is similar to the thrust structure of the main Himalayan thrust (MHT) belt along the southern margin of the Tibetan Plateau, indicating northward thrusting of the Tibetan Plateau along this boundary. Furthermore, the subduction position of the Asian plate relative to the northeastern margin Tibetan Plateau is still disputed, particularly in the northern margin of the Qaidam basin and the lower part of the Qilian orogenic belt [18–20].

The MT inversion results of Liang et al. through the Qilian orogenic belt have found large-scale high-resistivity bodies dipping south under the Alxa block and isolated conductive bodies in the middle-lower crust of the orogenic belt [21]. Since the degree of freedom of 3D inversion of single profile is significantly increased, there is a risk of introducing false anomalies [22]. To reduce this risk and further study the mechanism governing the outward expansion and growth of the Tibetan Plateau, we added 35 additional stations to build a 3D dataset (red dots in Figure 1(b)), providing more constraints for the model. In addition, we used unstructured tetrahedral elements for 3D inversion, which can fit more complex geological structures [23].

The MT array mainly transected the Qilian orogenic belt, extending from the northern margin of the Qaidam basin to the southern margin of the Alxa block and passing through Muli, Tianjun, and Gaotai cities. A total of 128 MT stations (blue and red dots in Figure 1(b)) were recorded to form an approximate 3D array, of which the new stations were recorded in 2019. The total length of profile L1 is 550 km (white solid line in Figure 1(b)), with 93 stations, and the average station spacing was 5 km. The electrical structure of this profile has been published by Liang et al. [21]. The MT data were collected using the Phoenix V5-2000 system (Phoenix Geophysics, Canada), with an average recording time of >20 h. The three mutually orthogonal magnetic field components (Hx, Hy, and Hz) and two mutually orthogonal horizontal electric field components (Ex and Ey) were recorded for each station. The $x$⁠, $y$⁠, and $z$ subscripts indicate the north-south, east-west, and vertical directions, respectively. The original MT time series data were transformed into the frequency domain after data collection, with the frequency-dependent transfer function computed using a statistically robust algorithm [26]. Finally, the MT data up to 2000 s was obtained.

Figure 2 shows the apparent resistivity and phase curves for typical measuring stations that sample the different geological units along the transect (yellow triangles in Figure 1(b)). It can be seen that the sounding data are generally with relatively high quality.

The MT data dimensionality should be ascertained prior to the MT inversion. We used phase tensors to assess the dimensionality of the resistivity structures. A phase tensor can be drawn as an ellipse defined by the minimum phase (⁠$φmin$⁠), maximum phase (⁠$φmax$⁠), azimuth of the reference axis (⁠$α$⁠), and skew angle (⁠$β$⁠). When the $β$ is greater than 3°, it indicates that the underground structure is three-dimensional features [27]. Figure 3 shows the phase tensor ellipses and corresponding absolute skew angles (⁠$β$⁠) at the MT stations for four observational periods (0.1, 1, 100, and 1000 s). The 0.1 and 1 s periods generally possess smaller absolute $β$ values (<3°), which imply that the shallow regional electrical structure primarily consists of 2D features. Larger absolute $β$ values (>3°) are observed at the 100 and 1000 s periods for some of the MT stations in the Qaidam basin and Qilian orogenic belt, which are indicative of 3D features (Figure 3). We therefore conclude that a 3D inversion approach is applicable for our MT data.

The MT data were inverted using the new 3D inversion program-GEM3D [23, 28]. In this new inversion code, the model is discretized by unstructured grids to better fit the complex underground structures. In addition, this program uses the limited-memory quasi-Newton (L-BFGS) inversion algorithm, which avoids the calculation of Hessian matrices and sensitivity matrices, greatly reduces the memory requirements, and improves the efficiency of 3D inversion [23, 29].

The full impedance tensors (⁠$Z$⁠) of 128 MT stations in the frequency range of 10–0.0005 Hz with 20 frequencies in total constituted the input data. For the off-diagonal elements, the error floor was set at 5% of $Zxy×Zyx1/2$⁠. Since the diagonal elements were more susceptible to noise [30], we set the error floor at 20% of $Zxy×Zyx1/2$ for the diagonal elements. Unstructured tetrahedral elements were used for mesh generation (Figure 4). The overall grid size was $8000km×8000km×4000km$ in $X×Y×Z$ direction. It contains 2000 km of air layer in $Z$ direction. The total number of elements was 1203521. The mesh in the central portion of the domain, which includes MT stations, was finer. We called it the core region. The grid size of this region was $500km×400km×100km$ in $X×Y×Z$ direction. The area near the MT stations was meshed with a tetrahedral element size of 5-20 m. The part of the core region not near the MT stations was meshed with a tetrahedral element size of 0.1-10 km. The final three-dimensional mesh was generated using the TetGen mesh generator [31]. The initial model was set to a uniform half space of 100 Ω·m. And the initial $λ$ was set to 100. When the relative change in overall root mean square (RMS) between two consecutive iterations was below 2%, it was then divided by 10. After 110 iterations, the RMS reduced from 18.11 to 2.97. Figure 5 shows the distribution of RMS. The RMS of most stations is less than 3, indicating that our inversion model is well constrained by the observed data. The apparent resistivity phase fitting curve of each station on profile L1 is shown in Figure S1 of the Supplementary Materials. The Niblett-Bostick method [32] was used to estimate the penetration depth of each station in L1 profile (Figure S2 in Supplementary Materials).

The electrical structural model of profile L1 obtained by the inversion is shown in Figure 6; the red areas (C1–C3) represent relatively low-resistivity areas, and the blue areas (R1–R3) represent relatively high-resistivity areas. We also marked the Moho depth in this region based on the deep geophysical results from a previous study [33]. Figure 7 shows 4 map view slices at depths of 5, 20, 40, and 60 km, respectively. Because there are additional stations in the Qaidam basin, compared with the previous electrical structure of Liang et al. [21], our results show that the high-resistivity body R1 in the Qaidam basin is larger in scale and may extend below the Moho. In addition, in the northeast of the study area, we have arranged new stations along both sides of the profile L1, which provides more constraints for the model. The new inversion result shows that there is a medium- to high-resistivity body R3 under the Moho of North Qilian and Hexi Corridor, and the shallow part is connected with the high-resistivity body R2 under the Alxa block (Figure 6(b)). The inversion results at 5 km depth show that the shallow part of the Qilian orogenic belt is characterized by high resistivity (Figure 7(a)), which may be related to the surface exposed Ordovician-Silurian rock bodies and arc associations, ophiolitic melanges, and low- to high-grade metamorphic rocks [34]. The most obvious characteristic of the electrical structure is that the middle-lower crust of the Qilian orogenic belt possesses low-resistivity bodies C1, C2, and C3. The Alxa block is a relatively stable block with a high-resistivity crust.

In addition, we also conducted a sensitivity test on the high-resistivity anomaly R3, by replacing the resistivity value in the rectangular area at the depths 60-100 km with 10 Ω·m (Figure 9(a)). The RMS of the modified model was 3.19, and the RMS_Update of the stations above R3 has changed significantly (Figure 9(b)).

The subduction position of the Asian plate relative to the Tibetan Plateau is still disputed. In our result, we observe that high-resistivity body R3 is located primarily under the North Qilian and Hexi Corridor (Figure 6(b)). The deeper section of R3 wedges southward into the lower part of the mantle of the North Qilian, forming a high angle. The shallow section is associated with high-resistivity basement R2 of the Alxa block, which has manifested into a wedging feature that begins at the hard North China Plate and extends to the lower part of the relatively weak North Qilian orogenic belt. The electrical structure of adjacent MT profiles also found high-resistivity bodies wedged southward under the Hexi Corridor and Alxa block [13, 16]. The results of seismic tomography show that there is high-velocity anomalies inserted southward by the Alxa block under the Qilian orogenic belt [20, 37]. The anisotropy of the orogenic belt reflects that it is affected by the southward subduction of the lithosphere of the northern block [38]. Besides, the results of receiver function show that there is a south-dipping LAB under the Alxa block and the North Qilian [19]. A south-dipping reflection interface has been identified in a deep seismic reflection profile along the North Qilian-Hexi Corridor [39]. The reactivation of the Haiyuan fault during the Late Miocene may be affected by the southward subduction of the North China Plate [19, 40, 41]. All these results indicate that the contribution of the North China Plate to the northeastern lateral growth of the Tibetan Plateau is not due to passive blockage but rather active wedging to the side of the Tibetan Plateau. R3 may represent subducted Asian lithosphere. The results of Liang et al. also found traces of southward subduction of the Alxa block [21]. On this basis, we believe that the Asian lithospheric mantle has been subducted under the North Qilian orogenic belt.

The electrical structure of the Qilian orogenic belt is obviously different to that of the Qaidam basin. Low-resistivity bodies C1–C3 are present in the middle and lower crust. Horizontal slices at different depths also show that there are low-resistivity bodies in the depth range of 20-60 km below the Qilian orogenic belt (Figure 7). Recent seismic studies have indicated that there is a low-velocity zone that corresponds to the low-resistivity zone in the middle and lower crust of the Middle and North Qilian [19, 42]. The genesis of low resistivities is usually related to graphite, sulfide minerals, partial melting, and fluids [43, 44]. Graphite films do not remain intact at high temperatures [45]. The temperature of the lower crust below the North Qilian can reach 600–840°C [46], which would result in unstable graphite films. The distribution range of low-resistivity anomalies that are produced by sulfide minerals is limited [43], whereas low-resistivity bodies C1–C3 cover a large area in our imaged electrical structure. Multiple lines of evidence have indicated that the low resistivities in the middle and lower crust of the Tibetan Plateau are closely related to partial melting and/or aqueous fluids [10, 47, 48]. In the front of the subduction of the Indian plate, there are also low-resistivity bodies in the middle and lower crust of the South Lhasa terrane, which may be a partial melting caused by tectonic heat [49]. The far-field effect of the convergence of Indo-Asian since the Cenozoic can affect the northeast of the Asian plate [50, 51]. Affected by the southward subduction of the Asian plate and the far-field effect of the northward advancement of the Indian plate, the Qilian orogenic belt is in a compressional setting and the crust undergoes napping deformation [5, 19, 21]. In addition, the middle-lower crust of the Qilian orogenic belt is rather felsic. Under the same temperature and pressure conditions, felsic rocks are generally weaker than intermediate/mafic rocks [52], which will weaken the middle-lower crust [53]. Zhao et al. believe that shear heating is the origin of low-velocity zone in the middle and lower crust of the Qilian orogenic belt [53]. Besides, the crustal nappe deformation generates frictional heating [49]. Therefore, the tectonic heat may be enough to partially melt the middle and lower crust, leading to lower resistivities.

The underlying causes of the Tibetan Plateau’s uplift are somewhat controversial [5, 7, 11]. The potential existence of crustal flow in the northeastern Tibetan Plateau also remains a controversial topic [54]. Our inversion results indicate that the crust of the Qaidam basin shows a high-resistivity R1. The Qaidam basin was a part of the Precambrian craton, which was called the Qilian-Qaidam Craton [3]. R1 may represent the rigid crust of the craton. Due to the blocking effect of the rigid crust in the Qaidam basin, crustal flow does not penetrate beneath the Qilian orogenic belt. Our results show that there are only local and small-scale high conductors under the Qilian orogenic belt, which is significantly different from the large-scale interconnected high conductors under the Qiangtang and Songpan-Ganzi blocks, where crustal flow is widely developed [10, 47, 48].

According to the Archie’s formula, the melt fraction of C2 is estimated in Figure 10, which is in the range of 2–5%, and is far less than the melt fractions of high conductors below the Qiangtang and Songpan-Ganzi blocks [57]. Besides, laboratory studies have shown that crustal flow is possible only when the melt fraction is greater than 5%, because in this case, the strength and viscosity of rock will be reduced by an order of magnitude [44, 55, 58]. It therefore seems highly unlikely that channel flow is responsible for the enhanced thickness of the lower crust.

Earlier studies have revealed the existence of a low-velocity layer in the Qilian orogenic belt [19, 42] that may be a weak layer in the crust [59]. Recent studies have also revealed that seismic anisotropy exists in both the mid-upper crust and upper mantle of this region [60, 61], such that the presence of low-resistivity bodies may indicate that the deformation in the upper and middle-lower crust is decoupled in this region. The lithospheric mantle of the North China Craton was continuously squeezed into the lower part of the Qilian orogenic belt owing to continuous compression between the North China Craton and the northeastern margin of the Tibetan Plateau. This caused the decoupled middle-lower crust driving the upper crust northward, as well as strong deformation in the upper crust and an outward thrust, which resulted in the uplift of the Qilian orogenic belt and generated a series of thrust fault [34, 62, 63].

We acquired an MT array across the Qilian orogenic belt in the northeastern margin of the Tibetan Plateau, and the unstructured tetrahedral element is used for 3D inversion of MT data to constrain the deep electrical structure of the region. Similar to the previous result of Liang et al., we found traces of southward subduction of the Asian plate, but the new electrical structure shows that there are medium- to high-resistivity bodies below the Moho of the North Qilian and Hexi Corridor, which indicates that the Asian lithospheric mantle has been subducted southward beneath the North Qilian orogenic belt. We also found that there are low-resistivity bodies beneath the Qilian orogenic belt, but this may be partial melting caused by tectonic heat. The melt fraction of the low-resistivity body calculated by Archie’s formula is 2-5%. Finally, we believe that there is no crustal flow beneath the Qilian orogenic belt. The lithospheric mantle of the North China Craton has been continuously squeezed into the lower part of the Qilian orogenic belt owing to continuous compression and convergence between the northeastern margin of the Tibetan Plateau and the North China Craton. The decoupled middle-lower crust has subsequently driven the northward movement of the upper crust, such that the Paleozoic sedimentary strata in the Qilian orogenic belt have undergone strong deformation and outward thrusting, resulting in its uplift.

MT data used to support the 3D electrical structure is included in this article. Please ask for permission for further use of the data.

The authors declare there are no competing interests.

This study was sponsored by the National Natural Science Foundation of China (Grants 41590863, 41430213, and 41504076), National Key R&D Program of China (Grant 2017YFC0601305), Central University Basic Scientific Research Business Expenses Special Funds, and Jilin Province Science and Technology Development Plan Project (Grant 20180101093JC). Some figures were created using the GMT software package [64].