Our study area includes most parts of the western Qiangtang terrane, the west end of the BNS, and the northern part of the western Lhasa terrane. We used 97 MT stations (Figure 1), including eight long-period MT data and 89 broadband MT data, collected by the SinoProbe (Dong et al., 2013) and INDEPTH-MT projects (Jin et al., 2007). The SinoProbe broadband MT (BBMT) data were collected using the Phoenix MTU-5 system, while the INDEPTH-MT BBMT data were collected using both the Phoenix MTU-5 system and the EMI MT24 system for more than 24 h. The long-period MT (LMT) data were collected using the LVIV Lemi-417 system for more than 7 days per site. All data were collected using at least two electric channels and three magnetic channels.

Time series were analyzed and processed using robust statistic methods (BBMT: (Egbert and Booker, 1986; Egbert, 1997); LMT (Varentsov et al., 2003)) to estimate the MT transfer functions of BBMT data with a period range of 0.01–3,000 s and LMT data with a period range of 0.01–10,000 s. Because of the lack of human activity around the stations, the data quality of most stations was pretty good.

To date, only one other 2D MT study involved in the western Qiangtang terrane has been published (Jin et al., 2007), and this prior study focused on the electrical structure of the BNS. The dataset presented in this paper will provide the first 3D electrical structure focused on the western Qiangtang terrane and its surrounding areas.

As the inversion of MT data is greatly simplified when applying a 2D assumption, it is necessary to get an overview of the MT data with its dimensionality. The phase tensors are calculated by using the algorithm given by Caldwell et al. (2004) in different frequencies (in Figure 2). In 1D configurations, the phase tensor ellipse is a circle. In 2D, the major axis direction of the phase tensor ellipse may be parallel or perpendicular to the electrical major axis direction. The filled color in the ellipse shows the 2D deviation angle β. If |β| > ±3°, it indicates significant 3D tectonic features (Booker, 2014).

In the period of 0.1–1 s, most of the data have a 2D deviation angle |β| = 0 and the ellipses have unequal lengths of the long and short axes, indicating a 2D tectonic feature. The data with a 2D deviation angle |β| > 3 are mainly distributed near the suture zone, indicating inhomogeneous electrical characteristics within the suture zone. In the period of 10–100 s, there is an increase in the amount of data with a 2D deviation angle of |β| > 9, and the length of the long and short axes of the ellipse is unequal. The direction of the long axis of the polarized ellipse of the data located at the tectonic boundary is consistent with the area’s direction. The closer to the boundary, the larger |β|, the larger it is. 3D tectonic features are evident in the moderately deep region of the study area. In the period of 1,000 s, most of the data have a 2D deviation angle |β| > 9, indicating strong 3D features in the deeper part of the study area. Therefore, the data in the study area require 3D inversion to establish the electrical structure of the area.

The 3D inversion code package named ModEM (Kelbert et al., 2014) with the L-BFGS algorithm (Dong and Egbert, 2019) was carried out with the 117 full-impedance tensor MT data. Error floors of them were set to 5% of |Z_xy×Z_yx|^0.5. To obtain further details of the deeper structure, we chose four frequency points between 0.1 and 1 s and 25 frequency points between 1 and 10,000 s evenly. The “core” area was discretized into a 10 km mesh in the north–south (x) and east–west (y) directions and a 50 m-thick uppermost layer with a step factor of 1.1. Therefore, a 60 × 63 × 66 mesh was finally formed. The start and initial models are uniform half-spaces of 100 Ω•m.

The preferred model was obtained after 209 iterations with the overall normalized-root-mean-squared (nRMS) misfit of 2.12, which has an initial nRMS misfit of 23.94. The distribution of the nRMS misfits of each station is illustrated in Figure 3. Overall, the data are fitted well, indicating that our preferred model is well-constrained and described satisfactorily by the observed data. For the impedance tensor, all the main features of the observed data are well reproduced by the calculated data by using our preferred model. A detailed comparison between the observed and calculated data can be found in Supplemental Material.

As the previous 2D result has shown (Jin et al., 2007), our preferred model also finds resistivity anomalies in the upper crust and conductivity anomalies in the mid-lower crust. The horizontal resistivity maps reveal a more intricate and complex geometry. To avoid potential false displays caused by site location and regularization mechanisms, we determined the confidence region at every depth to be circular, centered on each station’s location with a radius equal to the depth (opacity). The remaining regions are included for reference only (sub-transparent). According to the sensitivity test, the model is considered reliable for depths above 85 km (see Supplemental Material).

The main features of the resistivity model are the two conductivity anomalies in the mid-north of the Qiangtang terrane (C1) and the Lhasa terrane (C3), which are separated by resistivity anomaly under the BNS and south of the Qiangtang terrane (R) as shown in the longitude (N-S) resistivity section (Figures 4E–G). C1 and C3 extend from 20 km to the Moho (Zhang et al., 2014; Gilligan et al., 2015), with C1 extending deeper and C3 reaching some areas beneath the BNS. C1 is continuous only in the Qiangtang terrane’s northern region, while C3 is discontinuous in an east–west direction under the BNS. West of C1 and C3, there is a conductivity anomaly C2 in the junction of the KF, BNS, and LGCF (Figures 4A–C). The figure illustrates that the high conductivity anomaly C2 is closer to the normal faults rather than the strike-slip faults on the surface.

Our results reveal, for the first time, the electrical structure beneath the LGCF, which is believed to be the western extension of the Altyn Tagh fault and the western boundary of the eastward material movement. The sensitivity test identified a discontinuity between C1 and C2 beneath the LGCF (see Supplemental Material). Previous research using aeromagnetic (Xiong et al., 2016) and gravity data (Yin and Xu, 2017) found no clear boundary between the Tianshuihai subterrane and the Qiangtang terrane. The seismic result also suggests no significant discontinuity along the Moho surface beneath the LGCF, and the Vp/Vs result exhibits no significant variation between either side of the fault (Priestley et al., 2008; Wang et al., 2015). These findings suggest that the LGCF may not cut the entire crust.

The presence of high resistivity anomalies in the deep is indicative of the crystalline basement. The occurrence of a low-velocity or high-conductivity zone in the Tibetan Plateau crust is generally attributed to the presence of melt or salt-bearing fluids, which can significantly weaken the crust. In fact, most earthquakes occur within highly resistant bodies or at the boundary between highly resistant and low-resistant bodies (Figures 4E–G). The highly conductive Qiangtang mid-lower crust exhibits almost no seismic activity, suggesting the rocks in this region may be too weak for brittle damage. Conversely, earthquakes occurring near the BNS exhibit deeper focal depths, implying that the rocks are characterized by a higher strength profile in comparison to the Qiangtang terrane. This distinctive contrast in mechanical properties may potentially be linked to the distribution of ophiolites within the BNS (Wang et al., 2016). Seismicity is notably scarce below 40 km, suggesting that the rocks at this depth lack the ability to undergo brittle deformation.

The conductivity of the deep crust has been used to estimate fluid fraction and viscosity based on laboratory experiments in calibrating internal relationships. Our discussions below will focus on using the high conductivity zone (especially C1 in the Qiangtang terrane) to constrain the rheology of the crust, which can further help to understand surface deformation in the western Qiangtang terrane and its surrounding areas.

Similar to other regularized geophysical inversions, the ModEM MT inversion method aims to find a model that fits the data, while also maintaining smooth spatial variations between different layers, which is accomplished through the use of a smooth prior model and regularization to that prior (Kelbert et al., 2014; Dong and Egbert, 2019). Therefore, it is better to use the vertical conductance describing a conductor than the individual thickness and conductivities (Bai et al., 2010; Unsworth, 2010). Here, we integrate the vertical conductance values from depths of 20–75 km, which show that the conductance of C1 ranges from 5,500 to 16,500 S with 55 km thickness (Figure 5). In other words, the effective conductivity of C1 ranges from 0.1 to 0.3 S/m.

The high-conductivity zones observed in the mid-lower crust of the Tibetan Plateau, particularly when coinciding with low-velocity zones, are often interpreted as indicative of partial melt (Zeng et al., 2015; Li et al., 2020). Although the relationship between the high-conductivity zone and melt existence is still debated, it is worth considering for our study area, as volcanic activity has been found in the northern Qiangtang terrane from 5 Ma to near 1 Ma (Williams et al., 2004). Thus, we assume that C1 in our study area is indicative of the presence of melt.

To determine conductivity in partial melting zones, the region can be treated as a two-phase medium with solid and liquid phases, and conductivity can be estimated using the Modified Archie’s law empirical formula (Glover et al., 2000): σ e f f = σ r 1 − Φ log ⁡ ⁡ 1 − Φ m log ⁡ ⁡ 1 − Φ + σ f Φ m .

σ e f f is the effective bulk conductivity, σ r is the rock matrix’s conductivity, σ f is the fluid or melt conductivity, and m is the cementation exponent, where m=1 means that melts are highly interconnected and m=2 means that melts are isolated. Experiments have shown that the primary factors influencing the conductivity of the rock matrix and pure melt are temperature and pressure. In the absence of density and geothermal heat flow data, the depth–temperature relationship is established based on the average temperature model of the Tibetan Plateau (Shen, 1991), while the depth–pressure relationship is established using the PREM crustal density model (Dziewonski and Anderson, 1981) (see in Figure 6).

At this temperature and pressure condition, rock matrix conductivity can be estimated (Hashim et al., 2013). Pure melt conductivity (Guo et al., 2017) is determined from the latest volcanic composition, which is medium-acidic trachyandesite in this area (Wei et al., 2017). Water content can affect melt fraction estimates (Guo et al., 2017). The highest recorded water content in the western Qiangtang terrane is ∼6% (Williams et al., 2004). Due to the recent volcanic activities (<5 Ma) in the study area (Williams et al., 2004), we assume the melts are highly interconnected in the mid-lower crust. Thus, the relationship between melt fraction and conductivity at different depths was established using this water content (Figure 7A).

Estimations (Figure 7A) show that if the depth is larger, the melt fraction is lower for achieving a specific effective conductivity. At depths below 40 km, the estimated melt fraction is approximately 1%–7% for an effective conductivity of 0.1 S/m and approximately 4%–20% for an effective conductivity of 0.3 S/m. At 30 km depth, 0.1 S/m conductivity yields >20% melt, while higher conductivity requires more melt. Shallower depths need more melt for 0.1 S/m conductivity. Combined with the vertical profile, this estimate may indicate a low melt fraction (<7%) in the lower crust and a region of very high melt fraction (>20%) in the middle crust.

Our result suggests a higher estimate of melt fraction in the mid-lower crust compared to the previous seismic result (<3%) (Li et al., 2020). It should be noted that the previous study estimated the melt fraction of the entire crust, resulting in an average value. In contrast, our study provides a more detailed description of the melt distribution. The previous seismic result (Li et al., 2020) also indicated the presence of lower Vs in the mid-crust, implying a higher melt fraction in that region, which is consistent with our result.

Moreover, the extrusive igneous rocks in the northern part of the Qiangtang terrane are relatively young, with ages less than 5 Ma and some even younger than 0.3 Ma (Williams et al., 2004), indicating that C1 in the Qiangtang mid-crust could potentially be a high melt fraction zone, as previous studies proposed chambers at 10–25 km depth in this region (Yu et al., 2014). While it is possible that other volatile fractions may exist that could produce the same electrical conductivity anomalies with lower melt fractions (Hacker et al., 2014), this aspect is beyond the scope of this paper and will be investigated in future research. To estimate the rheology of the region, the current estimates are adequate.

The estimated melt fraction provides a qualitative understanding of the rheological properties of the region. Experimental studies have indicated that a melt fraction of more than 10% will significantly reduce the crustal strength (Rosenberg and Handy, 2005; Costa, 2021). The mid-lower crust’s shear viscosity can now be estimated based on the melt fraction using the pore elasticity equation (Schmeling, 1986; 1985; Schmeling et al., 2012) and an average viscosity of 2.0×10²⁰ Pa•s (Shinevar et al., 2015), assuming that the melt occurs in spherical or tubular inclusions (Figure 7B). Depending on the assumed melt geometry, both bulk and shear viscosity may be lower than 3×10¹⁸ Pa•s at different melt fractions: ∼50% for spherical molten pockets, ∼20% for tapered molten tubes, or 7% for ellipsoidal molten films. In other words, areas where the melt fraction is less than 7% are unlikely to exhibit viscous flow behavior. Therefore, C1 in the lower crust mainly deforms ductility and retains strength, with some patchy areas beneath the BNS and north of the Qiangtang terrane exhibiting viscous deformation. The mid-crust’s high-conductivity zone weakens and deforms viscously, while the rest of the crust deforms through ductility.

Based on the analysis from the previous section, the viscosity surpasses 3×10¹⁸ Pa s, in most of the western Qiangtang terrane. Therefore we believe that the primary deformation type for the regional crust should be ductile, except for the scattered high-conductivity anomalies representing viscosities below 3 × 10¹⁸ Pa•s. Moreover, the electrical model shows that LGCF does not fully cut through the crust, which indicates that the western Qiangtang and the Tianshuihai subterrane may remain connected in the deep crust (Chevalier et al., 2012). This also correlates with the GPS results showing consistent rotation directions for the two blocks (Wang and Shen, 2020). However, it is not easy to relate the rheology structure of the western Qiangtang terrane to the complex surface major tectonic distribution.

Unlike the situation in the central plateau, the craton-like Tarim Basin in the study area serves as an oblique boundary condition, owing to its angled orientation with respect to the northward subduction of the Indian plate (Wang and Shen, 2020). On the surface, the distributions of the strike-slip and north–south normal faults in the western Qiangtang terrane are much closer to each other compared to their distance in the central plateau (Taylor et al., 2003). The V-shaped conjugate strike-slip faults, such as Bue Co conjugate faults and Aishui Co conjugate faults, as well as KF and LGCF (Taylor et al., 2003) all exhibit an east-facing “<” orientation, indicating the occurrence of eastward motion in the region. However, the strike-slip motions of the KF and LGCF are less than 3 mm/yr (Wright et al., 2004; Raterman et al., 2007; Chevalier et al., 2017; 2016), while Aishui Co and Bue Co conjugate faults exhibit strike-slip rates ranging from 2.1–4.1 mm/yr (Taylor and Peltzer, 2006). These velocities are much lower than the strike-slip rates in the central Tibet Plateau (e.g., 6–11 mm/yr for Dong Co faults) (Taylor and Peltzer, 2006) and lower than the northward motion rate of the Indian plate (40 mm/yr) (Wang and Shen, 2020). This may imply that the strike-slip faults have limited contribution to the regional deformation in the study region, which further suggests that the surface deformation may be associated with the deep structures in the region.

Comparing surface deformation and deep rheological structures derived from electrical structures, it is observed that the weak, viscous zones in the mid-crust scattered inside the Qiangtang terrane correlate with the location of the normal faults on the surface (Figure 8). Previous petrographic analyses indicated that the melt responsible for the crustal weakening of the area may be linked to mantle upwelling (Guo et al., 2014; 2006; Wei et al., 2017). Therefore, the mid-lower crust of the western Qiangtang terrane may be divided into a series of small, dispersed ductile blocks separated by N–S directed viscous deformation areas. As the Indian Plate continues its northward subduction, the oblique compressional stress between the Tarim and Lhasa terranes would certainly lead to the eastward extrusion of the triangle-shaped west Qiangtang. Apart from the E–W expansion of the weak/viscous mid-lower crust material that contributes to the formation of the NW–SE-directed normal faults, the viscous region may also act as weak boundaries that allow the individual clockwise rotation of the ductile blocks in between. Indeed, these clockwise rotations may help facilitate the sinistral strike-slip faults in the north and south boundaries of the west Qiangtang terrane (Li et al., 2021). The combination of the extensional and rotational deformation may also explain the enigmatic combination of sinistral strike-slip and extensional earthquake focal mechanisms in the western Qiangtang terrane.