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).

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.

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.

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| × 10% for the ZXX component, |Zyy| × 10% for the ZYY component, |Zxy| × 5% for the ZXY component, and |Zyx| × 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.

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.

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.

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).

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).

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.

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.

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).

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 = 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).

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 ± 0.067 − 0.102 ± 0.016 w − 4228.5 ± 83.2 − 354.7 ± 18.1 w + 693.6 ± 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.

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.

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 H e / 3 H e 4 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 ( ≥ 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.

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.

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 ± .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.