3-D electrical structure and tectonic dynamics in the Yangbajing area based on the array magnetotelluric data

Lei, Lulu; Jin, Sheng; Dong, Hao; Wei, Wenbo; Ye, Gaofeng; Zhang, Letian

doi:10.3389/feart.2022.1089675

ORIGINAL RESEARCH article

Front. Earth Sci., 05 January 2023

Sec. Solid Earth Geophysics

Volume 10 - 2022 | https://doi.org/10.3389/feart.2022.1089675

This article is part of the Research Topic Episodic Plate Destruction and Construction in Southeast Asia: Observations, Modeling, and Case Studies View all 14 articles

3-D electrical structure and tectonic dynamics in the Yangbajing area based on the array magnetotelluric data

Lulu Lei¹^†

Sheng Jin^1,2*^†

Hao Dong^1,2

Wenbo Wei^1,2

Gaofeng Ye^1,2

Letian Zhang^1,2

¹School of Geophysics and Information Technology, China University of Geosciences, Beijing, China
²Key Laboratory of Geo-detection of Ministry of Education, Beijing, China

The well-known N-S-trending fault in the Yangbajing area plays a crucial role in the tectonic evolution of the Tibetan Plateau. Previous researches on a few E-W geophysical profiles suggested that the eastern shear at the base of the upper crust and/or lithosphere deformation brought on by asthenosphere upwelling are the major causes of the Yadong-Gulu rift’s creation. Here we propose a 3-D electrical resistivity model derived from the magnetotelluric (MT) array data spanning the Yadong-Gulu rift (YGR), and the distribution of temperature and melt fraction is estimated by the experimental calibrated relationships bridging electrical conductivity and temperature/melt fraction. The result reveals that the Indian slab subducted steeply in the east of the Yadong-Gulu rift, while Indian slab may have delaminated with a flat subduction angle in the west. The temperature distribution shows that the upper mantle of the northern Lhasa terrane is hotter than that of the southern Lhasa terrane. This is likely the result of mantle upwelling caused by either the subduction of the Indian slab or thickened Tibetan lithosphere delamination. Moreover, the strength of the mid-lower crust is so low that it may meet the conditions of the local crust flow in the west-east direction. The local crustal flow and the pulling force from the upwelling asthenosphere jointly contributed to the formation of the Yadong-Gulu rift. These main factors exist in different stages of the evolution of the Yadong-Gulu rift.

1 Introduction

Since the collision of Indian and Eurasia about 50 Ma ago, at least 2,000 km of convergence has been accommodated by thickening the crust and elevating the Himalayan-Tibetan Plateau (Yin and Harrison, 2000). How the high topography and the thickened crust of the Tibetan Plateau formed has been explained by several geodynamic models. They usually include the lateral eastern extrusion of continental lithosphere along several plateau bounding strike-slip faults (Tapponnier et al., 2001), a thin viscous sheet where the Tibetan Plateau undergoes distributed shortening and crustal thickening (England and McKenzie, 1982) and ductile flow in the middle-lower crust, that decouples deformation between the upper crust and mantle lithosphere (Clark and Royden, 2000).

As the continent-continent collision, a variety of tectonically and seismically active faulting zones are formed in different regions of Tibetan Plateau, with NW-SE and NE-SW trending conjugate strike-slip faults in central Tibet and N-S normal faulting systems in southern Tibet (Pang et al., 2018), and these rifts always suggest generally east-west extension of the Tibet (Molnar and Tapponnier, 1978). However, the mechanism of these extensional structures is still unclear. According to previous studies, two main genetic models have been proposed. These include as follows: 1) model related to the gravity collapse after a maximum uplift and crustal thickening of the plateau or delamination of the lithospheric mantle (England and Houseman, 1989). 2) model related to the regional stress filed and boundary conditions, which consists of lateral extrusion model (Tapponnier et al., 1982), arc bending model (Klootwijk et al., 1985), radioactive spreading model (Murphy and Copeland, 2005), oblique convergence model (McCaffrey and Nabelek, 1998) and Pacific plate rotation model (Yin, 2000). The broad similarities in the history of volcanism, the age of rift initiation, and the trend and direction of extension of rifts in Tibet, Lake Baikal, and Shanxi were emphasized in the model of the Pacific plate rotation model (Yin, 2000). Furthermore, Yin (2000) proposed that as the Indian northward subducted, all of eastern Asian experienced mantle upwelling beginning at ∼40–35 Ma leading to thermal weakening of the lithosphere and eventual rift development at ∼8–4 Ma. The Ultrapotassic N-S trending dikes and adakitic intrusives usually accompany the W-E extension in southern Tibet implying a magma source in the lithosphere mantle that also triggered melting of an eclogitic lower crust (Hou et al., 2004). Despite the insights provided by previous studies, the mechanism of the formation and evolution of the N-S trending faults in the southern Tibet is complex and controversial, and it may be not controlled by merely one factor.

As one of the most famous rifts in the southern Tibet, the formation of the Yadong-Gulu rift (YGR) has been studied by seismic studies (e.g., Zhang et al., 2013; Tian et al., 2015; Liu et al., 2019) and magnetotelluric study (e.g., Wang et al., 2017). Evidenced by the observation of shear wave splitting, Zhang et al. (2013) suggests that the onset of the N-S normal faulting does not indicate the gravitational collapse of the Tibetan lithosphere, and the shear rifting in the YGR is attributed to an eastern shear at the base of the upper crust. However, Tian et al. (2015) proposed that the continued extension of the YGR is mainly due to lithospheric deformation in response to asthenospheric upwelling and far-field W-E extensional stress. Wang et al. (2017) and Liu et al. (2019) indicate that the underthrusting of the Indian plate may cause slab tearing of the Indian slab and the upwelling of the asthenosphere, which could also drive E-W extension of the crust and contribute to the formation of the YGR in southern Tibet. However, these results are mainly reported by some E-W profiles in some specific regions. To better understand the evolution mechanism of the YGR, we did a 3-D magnetotelluric study by using MT array data from Sinoprobe, deployed in the Yangbajing region. Compared with the previous work of 3-D electrical profile across the YGR, this work covers a wider range near the Yadong-Gulu rift, which may provide us with a more comprehensive view of the electrical structure of this area. Also, more widely distributed electrical structure indicates the essential role of the mid-lower crustal conductors in the N-S rifting. In addition, the adding of the long-period MT data provides more reliable constraints for deep electrical structure.

The thermal state of the upper mantle of the lithosphere can help us understand 1) the origin and evolution of the lithosphere. 2) Lithospheric structure changes such as tectonic shortening, rifting, and mantle convection. 3) Relationship between shallow characteristics and deep dynamic processes. 4) Location of important mineral deposits and so on. The lithospheric thermal state is mainly acquired by four geophysical ways. We can use gravity data, surface heat flow data (SHF) and elevation to obtain a temperature model that can fit the observed data within the error range (e.g., Kaban et al., 2010). However, this method ignores the influence of factors such as the compressibility coefficient, phase transformation of rocks, and heterogeneity of composition in the upper mantle, which may lead to some differences between the thermal structures obtained by other methods. Thermal state of the lithosphere can also be established based on seismic data (usually shear wave). This method is usually dependent on the empirical relationship between seismic wave velocity (usually shear wave) and petrological parameters (Priestley and McKenzie, 2006). Although this is a direct method to detect the mantle structure, the mantle structure information obtained from seismic data cannot clearly distinguish the influence of thermal structure and composition structure in the mantle and thus this method still has a certain degree of ambiguity (e.g. Schutt and Lesher, 2006). Since the electrical conductivity of solid aggregates is exponentially sensitive to temperature through an Arrhenius relationship, Magnetotellurics (MT) has the potential to provide constraints on the lithospheric thermal state. Recent studies have shown encouraging results towards linking conductivity, composition, water content, temperature and melt fraction (e.g. Sheng et al., 2020).

In this work, from the derived 3-D electrical structure, we estimated the lithospheric thermal structure. Thermal state in this region can help us have a knowledge of the strength of the lithosphere, so as to further understand the formation of YGR and the tectonic and magmatic activities.

2 Magnetotelluric data and analysis

2.1 Magnetotelluric data

The magnetotelluric (MT) data in this study are located in the longitude of 90°E to 93°E and the latitude of 29°N to 32°N (Figure 1), which are from the China Magnetotelluric Standard-Grid Network (SINOPROBE, Dong and Li, 2009) acquired by China University of Geosciences, Beijing (CUGB) from 2011 to 2013. There is a total of 119 MT data in the study area, including 107 broad-band MT data (BBMT) and 12 long-period MT data (LMT).

FIGURE 1

FIGURE 1. (A) Topographic map of study area showing location of broad-band MT stations (blue triangles) and long-period MT stations (white triangles). The green lines mark the locations of the electrical resistivity profiles. QT, Qiangtang terrane; NL, northern Lhasa; CL, central Lhasa; SL, southern Lhasa; TH, Tethys-Himalaya terrane; BNS, Banggong-Nujiang suture; YGR, Yadong-Gulu rift; RCR, Riduo-Cuona Rift; SNMZ, Shiquanhe-Yongzhu-Nam Tso Mélange Zone; LMF, Luobadui-Milashan fault; IYS, Indus-Yarlung Zangbo suture (B) The extent of our study area (white rectangle). (C) Location of the thermal heat flow points (purple circles), hot springs (red, yellow, blue, and green circles), porphyry Cu-Mo deposits.

tHigh-quality broadband MT (BBMT) data were collected by Phoenix MTU-5 instruments over a period range of approximately .003–3,000 s and High-quality Long-period MT (LMT) data were collected by LVIV Lemi-417 long-period MT (LMT) instruments over an optimized period range of approximately 10–15,000 s. The incoherent noise was removed by synchronously recording the orthogonal electric field (E) and magnetic field (H) of different sites with GPS based on the geomagnetic coordinates (so called remote reference), and the data quality was improved. Robust statistic methods (LMT: Egbert and Booker, 1986; BBMT; Egbert, 1997) were used to analyze and process the MT time series files and estimate the frequency-dependent transfer function. The BBMT and LMT data were merged into single responses for each station. The apparent resistivity and phase curves of representative MT sites in the study area are shown in Figure 2. Furthermore, because of the extremely low environmental and human noise in Tibet, almost all of the MT data are of good quality.

FIGURE 2

FIGURE 2. Examples of apparent resistivity and phase curves for the sites 9029s50, 9030s30, 9030s40, 9031n50, 9031w30, 9232e10 (marked by the white circles in Figure 1A). The red triangles and blue circles represent the off-diagonal elements of the impedances. The solid lines (red, blue) show the predicted responses for the preferred model. The error bars are corresponding to the measurement errors of the observed data.

2.2 Dimensional analysis

Before the MT data are used to employ inversion, it is needed to analyze whether the 2-D inversion or 3-D inversion is preferred for the MT data in the study area. Given that the phase tensor (PT) parameters can generally avoid the frequency—independent electric field galvanic distortion from heterogeneities (Caldwell et al., 2004) and need no prior dimensionality assumption, we use the phase tensor decomposition method to evaluate the dimensionality of the MT data. The equivalent geoelectrical strike direction of the MT data is indicated by the ellipse long axis, and the skew angle β denoted by the filled color of the ellipse can indicate the dimensionality at different periods (Caldwell et al., 2004). It is generally indicated that the MT data show complex 3-D structures in the isotropric case when the absolute value of β (|β|-values) is more than 3° (Booker, 2014). The more the absolute value β changes, the stronger three-dimensionality exhibits.

As illustrated in Figure 3, in the period of .01–0.1 s, the |β|-values are less than 3° and the ellipses are nearly circular, which indicates that the shallow structure exhibits the quasi-2D feature. In the period of 0.1 s, the |β|-values are less than 3° in the northern Lhasa terrane and near the IYS, which indicates generally 1-D or 2-D feature, but in the other areas, the |β|-values are more than 3° which shows the 3-D feature. The results of the phase tensor in the high-frequency period show that the shallow structure is pretty complex which may due to the fault and rift zone. The flatted ellipses with light color show that the consistent dimensionality pattern in deep structure and 3-D dimensionality become increasingly stronger with the increase of depth. Furthermore, the lithospheric-scale suture zones probably dominate the dimensionality in the deep structure. As suggested by the dimensionality above, the MT data in the survey necessarily call for a 3-D inversion and interpretation.

FIGURE 3

FIGURE 3. The phase tensor ellipses for periods. (A)0.01 s, (B)0.1 s, (C)1 s, (D)10 s, (E)100 s, (F)1,000 s. The colors of the ellipses indicate the skew angle $β$ .

3 The 3-D inversion of the MT data

3.1 Details of the inversion and preferred model

Here, the ModEM inversion algorithm (Egbert and Kelbert, 2012; Kelbert et al., 2014), which utilizes the staggered-grid finite-difference approach to solve Maxwell’s equations and the non-linear conjugate gradient method to invert the MT data, is employed to invert the full impedance tensors (Z) in the study. A total of 32 frequency points in the range of .1–1,0000 s are used in the 3-D inversion. The initial model was a uniform half-space with resistivity of 100 $Ω ∙ m$ . In the initial model, there are 107 cells in the N-S direction, 85 cells in the E-W direction, and 71 cells in the vertical direction (including seven air layers). In the vertical direction, the thickness of the first layer was 100 m and the space of mesh was increased by a factor 1.1. The horizontal grid had a pacing of 4 km (in both directions). In the 3-D inversion, we set error floors as |Zxx| $\times$ 10% for the ZXX component, |Zyy| $\times$ 10% for the ZYY component, |Zxy| $\times$ 5% for the ZXY component, and |Zyx| $\times$ 5% for the ZYX component. Note that 5% error in impedance is equivalent to approximately 10% in apparent resistivity and 2.86° in phase. The model covariance smoothing parameters of the X and Y directions are chosen to be .3 and the initial damping parameter lambda is 1,000. After 170 iterations, the normalized-root-mean-square (nR.M.S.) misfit was reduced from 16.8 to 1.77. The misfit was less than 3 for almost all of the sites (Figure 4), which indicates that the model responses fit well with the observed data.

FIGURE 4

FIGURE 4. Site-by-site RMS of the 3D inversion. The misfit is <3 at almost all sites.

Supplementary Figures S1, S2 are the comparations of the horizontal pseudo-slices of the observed data and 3-D inversion response at different periods. The apparent resistivity shows pretty good agreement between the observed data and 3-D inversion response at all periods while the phases of both are also consistent in the periods of less than 1,000 s. Meanwhile, several sites have minor difference in the impedance phase in the period of 1,000 s and 1,0000 s. It is further proved that the 3-D electrical resistivity model is reliable.

3.2 Inversion results

Our 3-D electrical resistivity model (see Figure 5) shows that a large-scale resistive layer with the uneven bottom interface is at depths of less than 20 km in the upper crust, which is separated by several small-scale, isolated conductors. Five conductors marked C1, C2, C3, C4, and C5 are distributed beneath the resistive layer in the Lhasa terrane, respectively. Furthermore, a resistor (R1) is located in the mid-lower crust of the TH terrane in the south of the IYS, the northern boundary of which is also uneven.

FIGURE 5

FIGURE 5. Horizontal slices of the 3-D electrical resistivity model for depth of 10, 20, 30, 40, 50, 60,70, 80, and 100 km. See Figure 1 for labels.

FIGURE 6

FIGURE 6. Cross sections of 3-D electrical resistivity model, approximately oriented East to West (line 29, line30, line 31, line 32 marked in Figure 5). Inverted triangles: magnetotelluric measurement stations; Red lines: Yadong-Gulu rift; Cyan circles: major Miocene porphyry Cu-Mo deposits (PCDs); Yellow dashed open circles: possible supply magma chambers. The Moho depths are from Li et al. (2014)

FIGURE 7

FIGURE 7. Cross sections of the 3-D electrical resistivity model, approximately oriented South to North (line 90, line 91, line 92 marked in Figure 5). Inverted triangles: magnetotelluric measurement stations; Bold dashed black lines: possible locations of the main Himalayan thrust fault (MHT, Wei et al., 2001)

A series of sensitivity tests were carried out to verify whether the MT data has good constraints on these electrical anomalies. As a result, the 3-D electrical model is frozen with a block of 100 $Ω ∙ m$ beneath a series of specified depths (e.g. 40, 60, 80 and 100 km; Supplementary Figures S3–S7). Compared with the responses of the original inversion model, conductors C1, C2, C3, C4, and C5 are well constrained by the MT data.

Our electrical structure also shows three channel-like conductive zones that link the main conductors in the middle-lower crust. The specific locations of these conductive zones are illustrated in the Supplementary Figures S8–S10 and they were marked as A, B, and C. In order to verify the existence of the three conductors, we conducted sensitivity tests of conductors A and B in the depth range of ca. 50–80 km, by replacing the resistivity of the corresponding space with 30 $Ω ∙ m$ . And we also conducted a sensitivity test of conductor C in the depth range of ca. 60–80 km, by replacing the resistivity of the corresponding space with 20 $Ω ∙ m$ . The comparations between the calculated responses of the sensitivity tests and the original responses indicate that conductors A, B and C may exist in their corresponding depths with a little lower conductivity (20–30 $Ω ∙ m$ ) compared with the surrounding conductors.

In addition, Supplementary Figures S11A, B, C show the comparisons between the cross-section electrical structure of line 30 (this study) and the previous result of Wang et al. (2017). An obvious difference is whether a conductor with its bottom at a depth of approximately 77 km is in the east of the YGR. Note that the distribution of the MT data covering the conductive zones in the work of Wang et al. (2017) is denser than that in this study. In order to discuss what causes the difference between these two models, some sensitivity tests were carried out. Firstly, based on the scale of conductor C6 of the result of Wang et al. (2017), three groups of the modified models were established. In different depth ranges, we replaced the same area as C6 in our electrical resistivity model with the blocks with the resistivity of 3, 10, and 30 Ω·m respectively (Supplementary Figures S11–S13). Moreover, we also conducted another group of 3D inversion only with BBMT data (Supplementary Figures S14, S15). Generally, LMT data always has a relatively better constrains on the deep electrical structure. Compared with these three 3D inversion results and sensitivity tests, we found that involving of LMT data provides a better constrains on the bottom of the conductors and thus, there may not exist conductor C6 beneath the YGR in the depth range of 60–77 km or the C6 exists with a resistivity greater than 30 Ω·m. The difference in electrical structures between this work and Wang et al. (2017) may be due to the lack of LMT data which may provide more information about the deep earth. Furthermore, compared with the profile data, array data provide more information about the electrical structure of the study area and avoid the effect on the electrical resistivity profile from the conductors around.

Magnetotelluric studies often reveal anomalously high electrical conductivity in the mantle wedge above subducting plates (e.g., Worzewski et al., 2010; Pommier and Garnero, 2014; Ichiki et al., 2015), in the deep crust beneath young fold belts, such as the Himalaya (e.g., Gaillard, 2004; Le Pape et al., 2012). These anomalies are usually attributed to the presence of partial melts, brine-bearing aqueous fluids, or other high conductive phases, such as graphite (Wei et al., 2001; Unsworth et al., 2005; Jin, 2009), granulite (Fuji-ta et al., 2004), gabbro (Kariya and Shankland 1983), quartzite (Shimojuku et al., 2012) and granite (Olhoeft, 1981). Furthermore, recent petrological study of Guo et al. (2022) proposed that dehydration melting of garnet amphibolite contributes the conductivity comparable to the high conductivity anomalies. The possible origin of the high-conductivity layers (.1–1 S/m) widely distributed in the mid-lower crust of the Tibet plateau, is a long-standing and controversial issue (Nelson et al., 1996; Wei et al., 2001; Unsworth et al., 2005; Le Pape et al., 2012).

Among the factors that may lead to the high conductivity anomalies in the mid-lower crust of the Tibet plateau, metallic minerals, such as iron and copper sulphide, and graphite films can be easily ruled out due to the limited spatial distribution and low stability at the geological time-scale (Yoshino and Noritake, 2011). Furthermore, the conductivities contributed by typical rocks, such as graphite, granulite, gabbro, quartzite and granite are much lower than the high-conductivity anomalies (Yang et al., 2012). Therefore, Aqueous fluids or crustal-derived and/or mantle derived melts, or combination of both are always used to explain the high-conductivity anomalies in the Tibet plateau, because the conductivity of the aqueous fluids and melts are high enough for the crustal high-conductivity anomalies. Even a small amount of fluid can cause an order of magnitude reduction in the bulk resistivity, and in the viscosity (e.g., Rosenberg and Handy, 2005), but corresponding changes in seismic velocity are small-and may go undetected (Watanabe, 2013).

Aqueous fluids in the crust and upper mantle of the subduction zones usually contain some dissolved salts, especially the Nacl (Kawamoto et al., 2013; Frezzotti and Ferrando, 2015), which greatly enhance conductivity (Sinmyo and Keppler 2017). However, a brine-bearing, crustal-rock model is only applicable in regions with low geotherm (Guo et al., 2022).

The melts especially the partial dehydration melting of amphibolite and garnet amphibolite under lower crust conditions are also a potential candidate model for explaining seismic anisotropy and is also the origin of the low-velocity-high-conductivity zones (LV-HCZs) in the mid-lower crust of Tibet, especially for the following three reasons: The first is that amphibolite dominates the rock type in the region where LV-HCZs exist (Hou et al., 2017). The second is that the melts of amphibolite are of relatively high conductivity and low velocity; The third is the temperature is high enough to initial the dehydration melting of the amphibolite. In addition, phase equilibrium experiments have revealed that the dehydration melting of amphibolite stars at nearly 1123 K at 1.5 Gpa (Wang H. et al., 2019). Therefore, we reasonably consider that dehydration melting of amphibolite is likely to occur in the geotherm of the Tibet (Guo et al., 2022).

In our electrical resistivity model, conductors are located in the three different subterranes of the Lhasa terrane. Widely distributed shallow conductive layers (<5 km) in the Lhasa terrane may be considered to be a reflection of the geothermal activities. We notice that the relatively high-conductive structure in the western segment (west of 90°E) of the northern Lhasa subterrane is at depths of approximately 30–80 km with resistivity of 1–20 $Ω ∙ m$ . The conductor marked C1 corresponds to a low Vs. velocity anomaly in the Vs. velocity model (Huang et al., 2020), which was considered to be the weakened structure in central Tibet. In addition, the shallow part of the C1 may be a combination of aqueous fluids or saline fluids overlying a zone of partial melting, which greatly increase the conductivity of C1. A conductive upper mantle beneath C1 was thought to be consistent with the volcanic rocks in this area, which may be generated by the delamination of the Indian lithospheric mantle (Wei et al., 2001). Conductor C5 is corresponded with the LVEL (low-velocity zone in the eastern Lhasa terrane) which is at depths of 20–40 km, in the 3-D shear wave velocity model (Huang et al., 2020). Previous MT profile interpreted the conductive layer in the mid-upper crust is partially molten (Wei et al., 2001; Unsworth et al., 2005; Jin, 2009; Bai et al., 2010; Xie et al., 2016; Wang et al., 2017) and overlain by a thin fluid layer (Li et al., 2003). Moreover, the southern boundary of conductor C5 is approximately in the south of YZS and the northernmost part ends in ∼30.5 $° N$ . Conductor C4 is mainly located in the southern Lhasa, to the west of the YGR, and it is coincident with the conductor cross the Yarlung-Zangbo suture in the work of Wei et al. (2001). The geothermal activities due to aqueous fluids (Hochestein and Regernauer-Lieb, 1998) in this area may be associated with the high-conductivity feature of C4 in the shallow part (<15 km). While the deeper part of the conductor is thought to be the partial melt generated by the asthenosphere from deep mantle.

The widely distributed resistive cover identified in the shallow crust (∼10 km) might be the granite and multiple volcanic materials in the Lhasa terrane (Xu et al., 2006). The resistivity body in the Tethys-Himalaya terrane, such as R1 and R2 may represent the northern subducted Indian plate, which has been verified in the previous studies (Wei et al., 2001; Gao et al., 2014; Shi et al., 2015; Xie et al., 2017; Sheng et al., 2020; Sheng et al., 2021), while the resistive layer above the MHT presumably represents the high Himalaya crystalline rocks (Wei et al., 2001). In addition, the thick resistive layers (NR, MR1, MR2, and SR) in the Lhasa terrane are explained as the upward extrusion and cooling process of deep magma, which reflects the extensive and strong magmatic events under the Tibet plateau (Mo et al., 2009).

4 Thermal characteristic of lithospheric upper mantle and mid-lower crust

4.1 Thermal state of the lithospheric upper mantle

It is generally assumed minerals’ electrical conductivity is dependent on pressure, temperature, water content and composition, which can be described by the Arrhenius equation (Logan, 1982). Considering the different conduction modes in the minerals, the petrologically calibrated relationships between electrical conductivity and temperature of the main upper mantle minerals (55% Olivine, 28% orthopyroxene, 11% clinopyroxene and 6% garnet) (Peslier, 2010) (Dai and Karato, 2009; Yoshino et al., 2009; Yang et al., 2011; Zhang et al., 2012) and the basalt melt (Tyburczy and Waff, 1983) were utilized to calculated the lithospheric temperature. Moreover, Berryman (1995) and Jones et al. (2009) concluded that the Hashin-Shtrikman (HS) extreme bounds (Hashin and Shtrikman, 1962) are the most reliable approach for estimating the electrical conductivity of multiphase materials. Here, under the condition of the lithospheric mantle, we use the upper and lower bounds of the HS boundary to calculate the solid and solid-liquid mixed lithospheric mantle respectively. Besides, melt fraction of the lithospheric mantle was calculated based on thermo-dynamic modelling as well as experimental investigations of peridotite melting, and hydrous equilibria conducted by Katz et al. (2003).

It is worth noting that since MT data are better at constraining the conductance than in resolving the electrical conductivity (EC) with depth, we prefer to use bulk conductivity estimated by conductance and corresponding thickness rather than EC to estimate the thermal state of the upper mantle. Moreover, Water content, carbon dioxide content and pressure are very important for calculating the temperature and melt fraction of the lithospheric mantle (Peslier, 2010; Zhang et al., 2021). The maximum water content in the Tibetan plateau, according to Vozar et al. (2014) and Zhang (2017), was 200 ppm, but Le Pape et al. (2015) found that the water content in northern Tibet is roughly 100 ppm. Thus, we calculated the upper mantle thermal state with water content of 0 ppm (dry mantle) and 150 ppm (wet mantle with average water content) respectively, given that the ratio of water to carbon dioxide mass of the oceanic island basalt derived mantle (OIB) and basaltic mantle from mid ocean ridge was .82 ( $C_{w} C O_{2}$ =.82 $C_{w} H_{2} O$ ) (Sifre et al., 2014). In our study, based on the density model in the Tibet of Bai et al. (2013), the pressure is estimated by the formulation of $P = ρ_{i} ∙ g ∙ h_{i}$ (Sun et al., 2013) (Figure 8).

FIGURE 8

FIGURE 8. (A,B) The density model (modified from Bai et al., 2013) and pressure model calculated based on the $P = ρ_{i} ∙ g ∙ h_{i}$ (Sun et al., 2013), of Line 90 approximately oriented South to North. (C,D) The density model and pressure model calculated based on the $P = ρ_{i} ∙ g ∙ h_{i}$ , of Line 91 and Line 92 approximately oriented South to North.

FIGURE 9

FIGURE 9. (A–H) Distribution of the temperature and melt fraction at depths of 70, 80, 90, and 100 km [ $C_{w} (H_{2} O)$ =0 ppm]. (A–H) Distribution of the temperature and melt fraction at depths of 70, 80, 90, and 100 km [ $C_{w} (H_{2} O)$ = 150 ppm].

The temperature model of the upper mantle shows that temperature increase slightly with the increasing depths and the maximum temperature of the lithospheric upper mantle can reach to ca. 1,142°C and 1,087°C with a condition of dry mantle and wet mantle respectively. However, the melt fraction of the upper mantle appears to decrease with the increasing depth and increasing temperature. The maximum melt fraction of the upper mantle can reach to 12.39% at 70 km when the mantle is dry, while the maximum melt fraction is only 1.90% at 70 km when the mantle water content is 150 ppm.

The thermal state difference between the west and east of the YGR is so obvious that we speculate that the differences in underground metamorphism, material migration, volcanic earthquake activity between the two regions are partially caused by the asymmetry thermal state and eventually form the unique surface features.

4.2 Temperature of the lithosphere

Thermal state of the lithosphere controls its mechanical and may affect the lithospheric dynamics. Based on the steady heat conduction equation, the data of heat flow and element abundance of U, Th, K in the Tibetan plateau are used to calculate the crustal temperature in/beneath the Tibet plateau (Wang and Huang, 1990; Jiang et al., 2019; Čermák and Laštovičková, 1987 and reference therein). According to the study of Zhang and Zhao (2003), the depth of the Curie isotherm in the Yangbajing area is approximately 15–18.6 km, with the temperature of 580°C. The least square method was adapted to combine the calculated temperature of lithospheric upper mantle with the Curie temperature and crustal temperature to obtain the lithospheric temperature (Figure 10). The results show that crustal temperature in the east of YGR is relatively higher than that in the west of YGR, while the lithospreric mantle temperature in the two regions is similar.

FIGURE 10

FIGURE 10. (A) Lithospheric temperature in the north west of the YGR estimated by using average terrestrial heat flow of 95.5 $m W / m^{2}$ . (B) Lithospheric temperature in the south east of the YGR estimated by using average terrestrial heat flow of 106 $m W / m^{2}$ .

4.3 Moho temperature

The Moho discontinuity marks the boundary between the crust and the mantle and it plays a key role in mass transport. Constraining the depth and shape of the Moho beneath the Lhasa terrane may help elucidate the geometry of the subducting Indian plate beneath Lhasa Terrane (Shi D. et al., 2020). Moreover, the degree of coupling between upper and mid-lower crust and/or between the crust and lithospheric mantle is mainly controlled by the Moho temperature acting as a critical factor in influencing lithospheric rheology (Liao and Gerya, 2017). Here, by using the temperature calculation method same as the method mentioned in Section 4.1, the Moho temperature is estimated (Figure 11C). In this study, Moho depths refer to a compilation of receiver function (Li et al., 2014; Figure 11B) and water content in Moho depth is also similar to the “wet lithospheric mantle.” The average Moho temperature in our study area is greater than 1,200°C with relatively higher Moho temperature greater than 1,300°C beneath the Yadong-Gulu rift and Beng Co fault (BCF). This primarily results from the greater crustal thickness there. In contrast, lower Moho temperature west of the NL might at first seem counterintuitive and this primarily reflects thinner crustal thickness there. Moreover, previous study shows that when a crustal thickness of 70 km and linear temperature distribution with a constant strain rate of $10^{- 15} S^{- 1}$ , the threshold value Moho temperature for the occurrence of decoupling between the brittle upper crust and ductile mid-lower crust is estimated to be less than 500°C (Li et al., 2021). The higher Moho temperature, the larger degree of decoupling occurring in the crust (Liao and Gerya, 2017). Thus, the high Moho temperature beneath the YGR and BCF may indicate that the decoupling in these areas is strong, which may be related to the basal shear of the conjugate strike-slip zone (BCF) (Yin and Taylor 2011).

FIGURE 11

FIGURE 11. (A) A curved section of the 3-D electrical resistivity model along the Moho (B) Moho depth from a compilation of receiver function (Li et al., 2014); (C) Estimated Moho temperature [ $C_{w} (H_{2} O)$ =150 ppm; $C_{w} (C O_{2})$ =.82 $C_{w} H_{2} O$ , density-P], Purple dashed rectangles mark the locations of Moho offset from Zhang and Klemperer (2005), Tian et al. (2015) and Shi Z. et al. (2020). The brown filled circles mark high ${}^{3}{H e /} {}^{4}{H e}$ ratio anomaly sample points (Klemperer et al., 2022). The white circles indicate intermediate-depth earthquakes ( $\geq 70 k m$ ) from the high-quality (HQ) determinations (Chu et al., 2009)

4.4 The strength of the mid-lower crust

Archie’s law is considered to be a feasible method to estimate the melt fraction of the crust when the melt in the crust is interconnected in the previous studies (e.g., Bai et al., 2010; Rippe and Unsworth, 2010; Le Pape et al., 2015; Wang et al., 2017; Comeau et al., 2018). However, conventional Archie’s law is only valid for one conducting phase distributed within a non-conducting phase, which make it inapplicable in instance where the rock matrix has a significance conductivity such as for clay-rich rocks and in calculation involving partial melting (Glover et al., 2000). The modified Archie’s law (see Eq. 1; Glover et al., 2000) was put forward and it describes the electrical behavior of the tow-phase porous medium extremely well, improving greatly on the conventional Archie’s law.

σ_{e f f} = σ_{1} {χ_{1}}^{p} + σ_{2} χ_{2}^{m} (1)

p = \frac{\log (1 - χ_{2}^{m})}{\log (1 - χ_{2})} (2)

χ_{1} = 1 - χ_{2} (3)

Where $σ_{e f f}$ is the effective conductivity of the rock consisting of conductivities of pure melt $σ_{1}$ and conductivities of the rock matrix $σ_{2}$ . $χ_{1}$ is the volume fraction of melt and $χ_{2}$ is the volume fraction of rock matrix. The exponents (m and p) are constants describing the connectivity of the two-phase porous medium/fluid system. In this study we assume that the two-phase system consists of solid phase and melts in the crust, and m which is equivalent to n in the conventional Archie’s law is estimated to be 1.30 (Li et al., 2003; Unsworth et al., 2005). Since conductance can better reflect the distribution of conductive zones, conductance can be preferred to calculate the bulk conductivity ( $σ_{e f f}$ ) of mid-lower crust, which allow us to effectively evaluate whether locally deformation has occurred in the crust. The conductance is calculated for upper crust (0–30 km), middle crust (30–45 km) and lower crust (45–70 km) respectively (See Figure 12).

FIGURE 12

FIGURE 12. A map showing the distribution of conductance (A) Horizontal distribution of conductance at depths of 0–30 km. (B) Horizontal distribution of conductance at depths of 30–45 km. (C) Horizontal distribution of conductance at depths of 45–70 km.

FIGURE 13

FIGURE 13. The melt fraction in the mid-lower crust. (A) Horizontal distribution of melt fraction in the middle crust (∼30–45 km). (B) Horizontal distribution of melt fraction in the lower crust (∼45–70 km).

Previous studies have shown that melt conductivity can vary in a certain range. Some scholars argued that when the melt was assumed to be well interconnected, the pure melt resistivity was assumed to be .1–.3 $Ω ∙ m$ in the southern Tibet (Li et al., 2003; Unsworth et al., 2005), while another study in the southeastern Tibet shows that the pure melt resistivity is .05–.2 $Ω ∙ m$ (Dong et al., 2016). In this study, we use the general empirical formulation of peraluminous granite conductivity proposed by Guo et al. (2018) to calculate the pure melt conductivity under specific crustal temperature, pressure and water content conditions.

\log σ = 3.205 (\pm 0.067) - 0.102 (\pm 0.016) w - \frac{4228.5 (\pm 83.2) - 354.7 (\pm 18.1) w + 693.6 (\pm 37.8) P}{T} (4)

Where $σ$ is the electrical conductivity in S/m, $T$ is the temperature in K, $P$ is pressure in Gpa, and $w$ is the water content in wt%. Moreover, water content plays a key role in the process of melting (Unsworth et al., 2005), and previous studies have estimated water content of melt from different tectonic background (e.g., Patiño Douce and Harris, 1998; Hashim et al., 2013; Wang et al., 2013; Comeau et al., 2016). For example, more than 4wt% water content was consumed to exist in the lower crustal melt and water content of 5-6wt% (Patiño Douce and Harris, 1998) or probably >6wt% (Hashim et al., 2013) in the Himalaya granite magmas. Moreover, the maximum water content of crust in southern Tibet was deemed to be 8.4wt% (Guo et al., 2018). Based on previous studies mentioned above, the water content of the crust is assumed to be 7 $w t %$ in this study and thus the pure melt conductivity is estimated as 1.76–2.19 S/m under middle crustal condition (1173 K, 1.2 Gpa), and 1.55–2.02 S/m under lower crust condition (1273 K, 1.6 Gpa).

In the middle crust, the partial melting mainly occurred in the northwest of the northern Lhasa terrane and southeast of the central with more than 6% melt fraction and southern Lhasa terrane with minimum melt fraction of 10%, which shows a pretty consistent with the hot springs in the study area. In contrast, the melt fraction in most areas of the lower crust exceeds 10%. The effective viscosity would be reduced by an order of magnitude, when the melt fraction was larger than 5% (Rosenberg and Handy, 2005; Unsworth et al., 2005). The distribution of the melt fraction in the mid-lower crust shows that the mid-lower crust in our study area is relatively weak with a maximum melt fraction of approximately 20 vol%. It is well known that melt may introduce a drop in wave velocity. Moreover, Liu et al. (2019) indicates that there is a subtle drop in S-wave velocity (>12%) in the mid-lower crust, which may be due to the weak mid-lower crust in this area. Our result sheds lights on that possible local crustal flow may exist in our study area, even though more evidences are needed to support our deductions.

5 Interpretation and discussion

5.1 Subduction of the indian plate

The northward dipping resistors R1 and R2 located in the mid-lower crust and upper mantle beneath the Tethys-Himalaya terrane, which is also shown as the resistive zones (Wei et al., 2001; Wei et al., 2009; Xie et al., 2016; 2017; Sheng et al., 2020) and high-velocity zones (Liang et al., 2016; Huang et al., 2020), are indicated as the Indian lithospheric plate. Previous studies (Liang et al., 2016; Wang et al., 2017; Li and Song, 2018) indicated that the subduction angles of the Indian slab moving northwards are various in this region. Taking the YGR as a boundary, the subduction angle of the western segment is as different as than that of the eastern segment.

In our electrical resistivity model, the north-dipping resistor R1, with its northern boundary at approximately 29.5°N at depths of 30–80 km, is beyond IYS and does not extend to LMF in the west of the YGR. The resistor R2 in the mid-lower crust does not extend to IYS, which is consistent with the work of Sheng et al. (2020). No obvious resistive zones are between resistors R1 and R2. The resistive features mentioned above may indicate the different locations of the Indian lower crust (ILC) in different regions, which may result from the different onset locations of the Indian subduction because of the arched nature of the Himalayas and the IYS (Rosenbaum et al., 2008).

It is reported that the Indian lithospheric mantle may detach from the Indian crust and continue to move northwards (Kosarev et al., 1999; Shi et al., 2015). Previous studies have revealed that central Lhasa subterrane and the eastern segment of the northern Lhasa subterrane (east of 90^°E) were considered to be subterranes with ancient basement rock, while the southern Lhasa and western segment of the northern Lhasa subterranes were juvenile blocks (Zhu et al., 2011; Hou et al., 2015). Our works reveal that in our study area, conductive layers widely distribute in the Lhasa terrane while the eastern portion of the northern subterrane (ENL) and the west portion of the central subterrane (WCL) exhibit moderate resistive. In addition, the ENL and WCL correspond well to the high-velocity, ancient, rigid and stable areas (Huang et al., 2020). Comparing with the WCL, the ECL appears to be relatively high-conductivity and highly molten from middle to lower crust. The maximum melt fraction of the ECL can reach to ∼20% and the lithospheric upper mantle temperature in this region is ∼1,400°C, indicating that the ancient subterrane has undergone active magmatism.

Based on the discussion on the ENL, it is possible that ENL represents the remnant of ancient lithosphere of the Lhasa terrane (Zhu et al., 2011), which may play a key role in blocking the subduction of Indian slab, resulting in steep-angle subduction and short-distance advancement in the east of the YGR.

The continuous subduction of the Indian mantle beneath the Lhasa terrane caused the disturbance of the asthenosphere and upwelling of the hot materials to heat the mid-lower crust, resulting in the partial melting of the mid-lower crust (Li et al., 2008; Kind and Yuan, 2010). The possible subduction Frontier of the Indian lithosphere could be sketched according to the northern edge of the conductive layers in the mid-lower crust. Thus, combined with the electrical resistivity model and previous studies (Li and Song, 2018; Liu et al., 2019), we sketch a modified subduction model (Figure 14). Taking the YGR as the boundary, the subduction angle of the Indian slab in the east of the YGR may be steeper than that in the west of the YGR, while the Indian slab in the west of the YGR advances more northward than that in the east of the YGR. The sufficient strain localization triggered by the distinct variations in the subduction geometry of the Indian lithosphere may be one of the factors that caused the slab tearing. In addition, Liu et al. (2019) reported that the thickness of the Tibetan lithosphere is ∼30 km greater in the west of the YGR than in the east. The Indian lithosphere has an offset of ∼70 km across the YGR, which results in the deeper Moho depth in the west of the YGR than in the east (Liu et al., 2019). The unique structure of the Indian lithosphere led to the upwelling of the asthenosphere at different depths in west-east direction and formed the conductors C1 in the mid-lower crust and C5 in the mid-upper crust. Thus, the inherited structures of the Indian lithosphere may become another factor resulting in the tearing of the Indian lithosphere.

FIGURE 14

FIGURE 14. 3-D interpretative diagram of the current crustal and lithospheric mantle structure beneath YGR, showing the subduction mode of the Indian slab. MFT, Main Frontal Thrust; MBT, Main Boundary Thrust; MCT, Main Central Thrust; STDS, Southern Tibetan detachment system; JF, Jiali fault.

The uneven Moho surface (Figure 11B) corresponds well with the Hf mapping (Hou and Wang, 2018), which may imply that the Lhasa terrane has undergone the crustal thickening, mainly resulting from the wedge of the Indian slab and emplacement of the large-scale mantle sourced magma during the co-collision (Mo et al., 2007; Zhu et al., 2017). The thicker crust, with high Hf value beneath YGR suggests the double thickening effects and it may represent the juvenile crust caused by the subduction of the Indian slab. The small-scale conductor C3 to the north of the BNS may indicate the juvenile thickened crust caused by the bi-subduction of the Bangong Co-Nujiang (Hou et al., 2015), while the large scaled thickening beneath the YGR may due to the active tectonic activities caused by the N-S trending fault disturbing the asthenosphere and forming the extensive emplacement of the magmas in the mid-lower crust. However, the thickening of the crust in the area of the conductive layer C1 is not as prominent as the surrounding region. It is generally known that rolling back and delamination may occur as the subduction of the Indian slab and thereby triggering upwelling of the asthenosphere, forming large-scale magmatic activities (Chung et al., 2005; Ji et al., 2016; Niu et al., 2017). Thus, we indicate that the Indian mantle delamination may occur during the underthrust of the Indian slab along 90°E, which has been also verified by the teleseismic waves recorded by a temporary broadband network in approximately the same region (Kosarev et al., 1999), resulting in the not thickened crust in this area. The delaminated Indian slab resulted in the slab rebounding and forming the lower subduction angle of the Indian slab (Hou et al., 2020).

Localized Indian lithospheric mantle delamination into upwelling asthenoshperic mantle along 90°E can also well interpreted the distribution of the ${}^{3}{H e /} {}^{4}{H e}$ ratio anomalies (Figure 11B,; Klemperer et al., 2022) and ultrapotassic adakites (Guo et al., 2015), as well as the negative aeromagnetic anomaly (He et al., 2007) along 90°E, presumably a result of demagnetization from mantle heat flow. Furthermore, the distribution of the Linzizong volcanic rocks (Zhu et al., 2015), mainly formed by the asthenosphere mantle metasomatized by subduction plate fluid, may approximately outline the locations of the Frontier of the Indian lithospheric mantle. The northern limitation of the Indian lithosphere west of the YGR is at ∼30°N while the Frontier of the Indian lithosphere east of the YGR is at ∼30.5, which is consistent with the present-day Frontier of the Indian slab identified by shear-wave anisotropy (Chen et al., 2015). The lack of the Linzizong volcanic rocks in the northern and central Lhasa terrane along 90°E may mainly due to the delamination of the Indian lithosphere.

In addition, phase transitions may occur during the injection of the Indian lithosphere. The base of the crust may be eclogitized with the low geotherm at the early stage of mountain building (Nabelek et al., 2009; Zhao et al., 2010), and the followed major extensional event with magma intrusion from asthenosphere may cause the increasing temperature, rocks re-equilibrating of the rocks within the granulite facies in the lower crust (Richardson and England, 1979). Jackson et al. (2004) indicated that earthquakes near the Moho may relate to the infiltration of hydrous fluids and brittle fracture that accompanies a transformation of metastable granulite to eclogite. In the southern Lhasa terrane, the obviously thick crust suggests the granulite lower crust beneath the YGR. Moreover, dehydration of the subducting Indian lithosphere may generate fluids in/above the Indian lithospheric mantle (Zhu et al., 2015), thereby make the granulite-eclogite transformation process (Jackson et al., 2004). Combined with all these factors above, the intermediate-depth earthquakes ( $\geq 70 k m$ ), which are approximately at the present-day northern edge of the Indian slab, in the west of the YGR may also result from the increased strain-rate (Wang X. et al., 2019) associated with the deforming of the Indian slab into the ramp-and-flat geometry.

5.2 Formation of the Yadong-Gulu rift

Yadong-Gulu rift (YGR) is one of the most important structures in our study area. The thermal state, geophysical features and even structures on both sides of the YGR show obvious asymmetry. The different depths of C1 (20–80 km) and C5 (20–60 km) are mainly due to the inherited structure of the Indian lithosphere, which leads to the upwelling of the asthenosphere at different depths in west-east direction. Meanwhile, tearing windows were formed as the tearing happened and asthenosphere upwelled through the tearing windows, heating and weakening the overlying crust and resulting in the mantle convection. It is worth noting that the conductive layers do not spatially match the surface rifts. We speculated that in the process of asthenosphere upwelling, mantle convection may provide east-west driving force on the Tibetan crust, which may cause the extension of the Tibetan and eventually formed the YGR. The formation mechanism is similar to the pattern in the Tian et al. (2015).

In addition, previous works have proposed that middle crust flow or lower crustal flow was pervasive in the Tibetan plateau (Nelson et al., 1996; Cogan et al., 1998; Clark and Royden, 2000; Huang et al., 2020). In our study, large-scale conductive layers are located in the mid-lower crust and we do a set of sensitivity tests to verify the connections between the isolated conductors. It seems that the isolated conductors (C1, C2, C3, C4, and C5) were interconnected in the middle and lower crusts, but whether the local crust flow exists in these areas still need to be verified from other physical characteristics. Based on the parameters of the Archie’s law in this region, the melt fraction of 5%, acted as the minimum threshold of the reduction of the effective viscosity by an order of magnitude, is estimated by the bulk conductivity of at least .038 S/m (26.5 Ω m), which may contribute to the local “crust flow.” When the resistivity is less than 26.5 $Ω ∙ m$ , the melt fraction will be greater than 5%, resulting in the effective viscosity reduced by an order of magnitude (Unsworth et al., 2005) and thus, middle and lower crustal flow is possible to exist locally. The 3-D crustal resistivity model (Figure 15) is established with the maximum resistivity value of 26.5 $Ω ∙ m$ between depths of 20–70 km. Areas between conductor C1 and conductor C5 show strong connectivity, which indicates that local crustal flow only occurs cross the YGR and the connection C provides a channel for the materials migration between the east and the west. In contrast, the resistivity sensitivity tests revealed that connections A and B may also exist in the mid-lower crust, but the melt fraction in these connections is not enough to meet the rheological conditions.

FIGURE 15

FIGURE 15. 3-D view of the crust electrical structure of the study area. The maximum resistivity of this model is 26.5 $Ω ∙ m$ between depths of 20 and 70 km. Tectonic processes are marked schematically. IL, Indian Lithosphere; YGR, Yadong-Gulu rift.

In addition, the temperatures on the both sides of the YGR are calculated in the Section 4.2, which shows that the lithospheric temperature is relatively higher in the west of the YGR than that in the east. As a result, connection C is not only a channel for material migration, but also a pathway for heat convection. Under the condition of bidirectional (N-S compression and E-W extension) lithospheric deformation, the widely distributed mid-lower crustal weak zones resulting from dehydrate melting or wet melting play a crucial role in the formation of the N-S trending rifts (Pang et al., 2018). Eastward crustal flow occurred due to the pressure difference in the west-east direction (Clark and Royden, 2000), which could create a west-east basal shear on the Tibetan crust that facilitates surface rifting. The high Moho temperature (See Figure 11) also indicates that strong decoupling exists between the brittle upper crust and ductile mid-lower crust. And thus, the west-east direction crustal flow also makes contributions to the dextral slip of the Beng Co fault to some extent, resulting in the faster ∼EW extension rate of up to 3–6 mm/yr near the northern end of the rift (Gulu) compared with only 1.3 $\pm$ .3 mm/yr in the south (Yangbajing) (Chevalier et al., 2020).

Both mantle convection and local crust flow lead to the west-east slip of the YGR supporting the focal mechanisms along YGR (Zhu and Helmberger, 1996; Chen and Yang, 2004; Baur, 2007; Ekstrom et al., 2012). The genesis of the YGR has the characteristics of multi-stage and multi-factor.

Comparing with the previous viewpoint that the rift zones were just restricted to the upper crust (Masek et al., 1994; Nelson et al., 1996; Cogan et al., 1998), the estimated melt fraction and the electrical resistivity structure of the lower crust and upper mantle suggest that northern and central YGR may cuts throughout the entire crust and mantle lithosphere for the local migration of the lower crust and the hot materials upwelling from the upper mantle through the rift zones. Furthermore, instability analysis of rift spacing (Yin, 2000), the ultrapotassic and adakitic magmatism along rift zones (Williams et al., 2001; Chung et al., 2005), body wave tomography studies (Ren and Shen, 2008; Liang et al., 2011; 2012; Huang et al., 2020), surface wave tomography (Jiang et al., 2014) provided the direct evidences to support that the YGR originate from deep lithosphere.

6 Conclusion

In this study, we have presented the results from a 3-D magnetotelluric study covering the Yadong-Gulu rift. Based on the electrical structure, the experimentally calibrated relationships bridging electrical conductivity and temperature/melt fraction were applied to construct the thermal and rheological structures of lithosphere.

On basis of the electrical structure, thermal structure as well as rheological structure, our study proposed the subduction model of the Indian slab which shows a steeper subduction angle in the east of the YGR and an Indian slab delamination in the west of the YGR. Moreover, subduction of the Indian slab may disturb the asthenosphere and trigger the mantle flow beneath the Lhasa terrane and the upwelling asthenosphere from the deep may widely melt the mid-lower crust of the corresponding area.

The estimated Moho temperature suggests that the decoupling between the brittle upper crust and ductile mid-lower crust is so strong that they may contribute to the formation of the YGR and BCF. Furthermore, the derived rheological feature shows that the mid-lower crust is so weak that it meets the rheological conditions to generate local crust flow from west to east under a pressure difference beneath the YGR. It is the local crustal flow that play an essential role in pulling the Indian and Tibetan crust and forming the YGR.

What’s more, the combination of the granulite lower crust beneath the Lhasa terrane, the fluids generated by dehydration of the Indian slab and the extremely high strain-rate for the special locations contributed to the occurrence of the intermediate-depth earthquakes along the Frontier of the delaminated Indian slab.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions

LL: data curation, methodology, visualization, software, writing—original draft. SJ: writing—review and editing, project administration, funding acquisition. HD: writing—review and editing. WW: writing—review and editing. GY: writing—review and editing. LZ: writing—review and editing.

Funding

This study is funded by the Second Tibetan Plateau Scientific Expedition and Research Program (2019QZKK0701), China Scholarship Council (202006400054) and the National Key R&D Program of China (2016YFC0600301).

Acknowledgments

We thank Gary Egbert and Anna Kelbert for providing their 3-D MT inversion code Modem. Special thanks must go to our field crews and students; this study would not have been possible without their efforts. We used the GMT software package (Wessel and Smith, 1998) for production of some figures.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.1089675/full#supplementary-material

References

Bai, D., Unsworth, M. J., Meju, M. A., Ma, X., Teng, J., Kong, X., et al. (2010). Crustal deformation of the eastern Tibetan plateau revealed by Magnetotelluric imaging. Nat. Geosci. 3 (5), 358–362. doi:10.1038/ngeo830