This study investigates the upper mantle deformation pattern beneath the Indo-Eurasia collision zone utilizing the core-refracted (S(K)KS) phases from 167 earthquakes recorded by 20 broadband seismic stations deployed in the Western Himalaya. The 76 new shear wave splitting measurements reveal that the fast polarization azimuths (FPAs) are mainly oriented in the ENE-WSW direction, with the delay times varying between 0.2 and 1.7 s. The FPAs at most of the stations tend to be orthogonal to the major geological boundaries in the Western Himalaya. The average trend of the FPAs at each station indicates that the seismic anisotropy is primarily caused due to strain-induced deformation in the top ~200 km of the upper mantle as a result of the ongoing Indo-Eurasian collision. A contribution from the mantle flow in the direction of the Indian plate motion is possible. The mantle strain revealed in the present study may be due to a combination of basal shear resulting from plate motion and ductile flow along the collision front due to compression.

The Himalaya-Tibet orogen is a result of a collision between the Indian and Eurasian plates that has been ongoing since 50 Ma. An effective way of deciphering the deformation associated with such geodynamic processes is through measurement of seismic anisotropy using the core-refracted (SK(K)S) phases. These phases carry information about the nature of azimuthal anisotropy along the ray path from the core-mantle boundary (CMB) to the surface beneath the station (e.g., [1]). However, the lower mantle above the D layer is considered to be isotropic because the deformation is associated with diffusion creep (e.g., [1]). This implies that a major contribution to anisotropy comes from the first 400 km of the upper mantle [2].

Seismic anisotropy is an intrinsic property associated with the earth’s medium that characterizes a change in the velocity of a seismic wave with direction [3, 4]. This causes the shear wave to split into two components while passing through an anisotropic medium. These two components are orthogonally polarized and travel with different velocities, with the component that travels parallel to the fast direction moving faster than the other. Shear wave splitting is one of the widely used techniques to understand the dynamic processes in the earth’s mantle due to past and present deformations [5]. To quantify anisotropy, the splitting parameters phi (Φ) and delay time (δt) are calculated. The parameter δt provides information about the thickness of the anisotropic layer and its strength, and Φ denotes the polarization of the fast wave. Seismic anisotropy in the earth’s upper mantle results due to the preferred crystallographic orientation (“lattice preferred orientation”) adopted by its constituent minerals (e.g., olivine) in response to deformation [6]. For large strains, the LPO, which is represented by the direction of the fast axis, lies parallel to the flow direction, and for small strains, the LPO is rotated with respect to the flow direction, indicating that seismic anisotropy is extremely sensitive to the sense of shear in the upper mantle [7].

The NW Himalaya is rigidly bounded at its southern end by the east-west trending Himalayan Frontal Thrust (HFT) and to its north by the Indus-Tsangpo suture zone (ITSZ). The Himalayan orogeny in the study area can be broadly divided into four major lithotectonic units (Figure 1): Outer (Sub) Himalaya (SH), Lesser Himalaya (LH), Higher Himalayan Crystalline (HHC), and Tethys Himalaya Sequence (TH), from south to north [8]. The major tectonic divisions are segregated from each other by NW-SE trending fault zones [9] (Figures 1 and 2), which from south to north include the Himalayan Frontal Thrust, the Main Boundary Thrust, the Main Central Thrust, and the South Tibetan Detachment. The study area covers the major seismotectonic units that can be termed as the Garhwal-Himalayan region, the Kangra-Chamba sector, and the Kinnaur segment of the NW Himalayan region (Figure 2). These lithospheric units are reported to be extremely heterogeneous and structurally complex due to the ongoing Indo-Eurasian collision.

Numerous shear wave splitting studies have been performed for the Indo-Eurasian collision zone, which infer a complex deformation pattern beneath the region (Figure 1) [1019]. The existing data sets show that seismic anisotropy is different in the diverse tectonic units like the Indian subcontinent, Northeast Himalaya and its foothills, Western Himalaya, and Southern Tibet [14, 1826]. Shear wave splitting studies for the Indian subcontinent are interpreted in terms of the absolute plate motion- (APM-) related strain in the mantle due to the simple progressive shear at the base of the lithosphere [14, 23, 24, 27].

In the NE Himalaya and its foothills, the FPAs are mostly oriented along the NNE-WSW direction and/or are parallel to the strike of the orogens ([18, 25, 26]). In the Northern and Central India, the FPAs are mostly oriented in accordance with the plate motion direction of the Indian plate. In addition, the effect of the Indo-Eurasian collision is also observed [25] (Figure 1). In southern India, the effect of both plate motion and frozen anisotropy is observed [14, 2729]. The anisotropic orientation changes from ENE-WSW direction in the NW and Western Himalaya to ~E-W direction to the north of the Indus-Tsangpo suture zone [13, 16, 17, 22, 25, 30, 31]. These authors suggest that the major contribution of anisotropy is mainly due to strain induced by the Indo-Eurasia collision and the asthenospheric flow. In Southern Tibet, the lateral variation of FPAs is interpreted in terms of underthrusting of the Indian plate in the Himalaya-Tibetan zone [20]. The region is also associated with weak anisotropy due to the existence of LPO fabric, with a subvertical axis of symmetry [17, 21].

The NW region of India, an area covering Garhwal and Himachal Pradesh, has been hit by four destructive moderate to great earthquakes since the beginning of the 20th century. These are the Kangra earthquake of 1905 (Mw 7.8), the Kinnaur earthquake of 1975 (Ms 6.8), the Uttarkashi earthquake of 1991 (Mb 6.5), and the Chamoli earthquake of 1999 (Mb 6.4). These seismic activities manifest the large-scale subsurface deformation and weak zones, which can produce significantly larger events in the Himalayas [32], which could as well have happened in the past [23, 24, 3336]. Therefore, to have deeper insights into the ongoing deformation beneath these tectonically unstable zones, a 20-station broadband seismic array was deployed in the NW Himalayan region. In the present work, we determine the anisotropic nature of the NW Himalayan (India) region using the core-refracted SK(S)S phases. The cause, source, and the associated geological structures responsible for anisotropy of the NW Himalayan region are discussed in the context of this study.

In this study, data acquired from 20 broadband seismic stations operated by the Wadia Institute of Himalayan Geology (WIHG), Dehradun, during the period 2007–2013, are used (Figure 2; Table S1). The units include Trillium-240 broadband sensors having a flat velocity response between 0.004 and 35 Hz. The sensors are connected to 24-bit Taurus digitizers set to record the data at 100 samples/s. To perform shear wave splitting analysis, the SK(K)S phases are extracted from high quality (S/Nratio5.0) waveforms of 167 earthquakes of magnitude5.0 in the epicentral distance range of 85° to 145°. The study resulted in 76 shear wave splitting and 28 null measurements (Figure 3).

The splitting parameters phi (Φ) and delay time (δt) are measured assuming that the orientation of the fast symmetry axis is predominantly horizontal. These parameters are estimated through a grid search method [5], by correcting the observed components for the anisotropy effect, achieved through minimization of energy on the corrected transverse component. To enhance the dominant period of SK(K)S phases, the waveforms are band-pass filtered between 0.05 and 0.25 Hz using a Butterworth filter. The splitting parameters are estimated using the SplitLab software, which uses a manual selection of window that further improves the quality controls on the seismogram [37]. SplitLab handles large datasets efficiently, with easy quality checks for the seismograms. The software enables estimation of parameters using three different methods, namely, Rotation correlation [38], Minimum energy [5], and eigenvalue [5].

The splitting parameters presented in this study are obtained using the Minimum energy method, since it provides a more stable and robust solution [39]. The splitting parameters are selected based on these criteria: (1) the desired phases (SKS and SKKS) should be well identified on the seismic waveforms, (2) the energy corresponding to the transverse component of the waveform should be minimum after applying the anisotropic correction, and (3) the initial particle motion should be elliptical in nature and should get linearized after the correction (Figure 4). Null measurements obtained from this study can be attributed either to the presence of multiple anisotropic layers beneath the study area [40] or to the presence of an isotropic layer [4] (Figure S1). Null measurements may also indicate that the incoming polarization of the SK(K)S wave is aligned with the fast axis (or perpendicular to the fast axis).

Each splitting measurement is carefully checked and sorted into three categories, namely, “good,” “fair,” and “poor,” based on the values and errors of phi (Φ) and delay time (δt). For a measurement to be categorized as “good,” the 2σ errors should be ≤±15° for phi (Φ) and ≤±0.3 s for δt, and the deviation observed for the estimated splitting parameters calculated by ΦRCΦME should be within ±7°, and δtRC/δtME should be between 0.9 and 1.1. For the event to be categorized as “fair,” the 2σ errors should be ±30° for phi (Φ) and ±0.6 s for δt, and the deviation observed for the estimated splitting parameters calculated by ΦRCΦME should be within ±15°, and δtRC/δtME should be between 0.7 and 1.2 [25].

In this study, a total of 76 well-constrained shear wave splitting measurements and 28 null measurements are obtained (Table S2). The delay times vary from 0.2 s to 1.7 s, with a majority of them between 0.5 s and 0.7 s (Figure 5(a)). The orientations of FPAs show some variation but are mainly in the ENE-WSW direction (Figure 5(b)). The individual shear wave splitting parameters and null measurements are plotted at the respective stations (Figures 6(a) and 6(b)). At stations DBN, NHN, KHI, and PULG located in the western part of the foothills of NW Himalaya, the FPAs are oriented in the NE direction (Figure 6(a)), in accordance with the APM of the Indian plate in a no-net rotation frame [41]. In the central and northern parts of the Main Central Thrust, the FPAs are mainly oriented in the ENE-WSW directions, while to the south of it, the FPAs are oriented in the NE direction. Station KAZA yields the maximum number of splitting and null measurements, and the FPAs are mainly oriented in the ENE-WSW direction. Overall, the FPAs seem to be consistent across the different tectonic regions in the NW Himalayan segment. A good coherence is observed in both the fast axis orientations and the delay times, from west to east in the NW Himalaya, with small variation on a local scale. To examine any back azimuthal variations, the splitting and null measurements are plotted as a function of back azimuth and incidence angle (Figure 7). It can be seen that the splitting measurements are mostly from the south-eastern back azimuth and are predominantly oriented in the ENE-WSW direction. We obtained the splitting measurements mostly from a back azimuthal range of 90° to 120°, the maximum from the back azimuth of ~100°.

In the NW Himalaya, the FPAs do not reveal much variation among the stations. The average delay time is 0.8 s with a cluster around 0.5 s (Figure 5(a)). The FPAs in the NW Himalaya are mostly oriented in the NE and ENE-WSW direction, the latter being the predominant orientation. However, the major structural trends are predominantly in the NW-SE direction, orthogonal to the absolute plate motion direction of the Indian plate in a no-net rotation frame ([41] (Figures 1 and 2). From south to north in the Western Himalaya, the majority of FPAs are orthogonal to the major structural trends and subparallel to the absolute plate motion of the subducting Indian plate. Shear at the base of the Indian lithosphere might result in the orientation of mantle fabrics along the plate motion direction. The regional stress axis computed by inverting the focal mechanisms shows a NE-SW orientation ([42]; Figure S2), in concurrence with the ongoing lithospheric deformation in the NW Himalaya. On the other hand, the ENE-WSW orientation reflects lithospheric strain due to a convergence of ~14-21 mm/yr inferred from GPS measurements in this part of Himalaya [43, 44]. Previous studies in the NW Himalaya also observed ENE-WSW orientated FPAs [22, 25]. The FPA at station KHI in the western part is along the plate motion direction, while in the eastern part, station GTU shows an orientation along the ENE-WSW direction. Similar observations are found at stations in the Lesser Himalaya. In the Higher and Tethys Himalaya, majority of the stations have ENE-WSW direction as the FPA, except for station PULG where it is oriented in the absolute plate motion direction. The ongoing India-Eurasia plate collision has reoriented the a-axis of the olivine in the underlying Indian plate lithosphere, which has caused the FPAs to align ENE-WSW for the NW Himalaya, which was initially along the absolute plate motion direction of the Indian plate [40].

In the Ladakh Himalaya that lies to the north of the study region (Figure 1), the crustal deformation pattern inferred through local shear wave and the Ps splitting is found to be complex [45]. The FPAs estimated from shear wave splitting are oriented along the strike of the trend (NW-SE) for the Karakoram fault and the Indus suture zone, whereas in the Ladakh Himalaya (Figure 1), they are oriented along the plate motion direction of the Indian plate. The FPAs derived from Ps splitting are found to be consistently oriented along the strike of the faults with a slightly smaller average delay time of ~0.15 s [45]. These observations are explained in terms of anisotropy in the middle and lower crusts. This suggests that the crust and mantle are decoupled, and the degree of deformation due to the collision is found to be larger in the crust than in the upper mantle. It is further suggested that the source of SKS anisotropy beneath the NW Himalaya might lie in the lithospheric mantle (excluding the crust) and asthenosphere. It seems that the observed anisotropy from the lithospheric mantle mainly reflects the anisotropy of the underthrusting Indian plate.

There is a significant variation in the nature of shear wave splitting from east to west in Southern Tibet ([21, 4650]). Further, the southern part of Tibet might be associated with slab tear and lateral variation of subduction of the Indian lithosphere [20]. The central and western parts of the torn slab propagate parallel to the convergence direction [20, 51, 52]. They further suggest that the central and western Indian lithosphere subducts at a shallow dip angle which favors the Indian plate to slide to a greater distance. The modeling of SKS waveforms also suggests that the observed anisotropy in central Tibet is due to underthrusting of the Indian lithosphere beneath Eurasia and is associated with a flatter axis of symmetry rather than a steeply dipping symmetry axis [52]. As the western Indian lithosphere slides easily, the effect of collision should be comparatively small. This implies that the degree of deformation should be small in the western Indian lithosphere.

A recent study in the Eastern Himalaya also found that the FPAs are mostly oriented along the arc curvature, lending credence to the hypothesis that lithospheric deformation is the primary source of anisotropy [26]. Further west in the Sikkim Himalaya, strike parallel FPA orientation follows [19]. Similar orientations are observed to the further east in the Assam foredeep and the Himalayan front [18, 25]. Such a correlation with the arc curvature is not observed for NW Himalaya. Further, due to the higher convergence rate, the amount of shortening is found to be more for the eastern Himalaya [49]. This might cause the FPAs in Arunachal Himalaya to orient along the strike of the orogen, suggesting that the coherently deformed mantle is under compression [26]. The amount of shortening decreases towards west. It seems that in the NW Himalaya, the degree of deformation due to collision seems to be stronger for the crust and its impact decreases gradually deep down in the upper mantle.

This study presents results of shear wave splitting analysis of core-refracted phases (SK(K)S) at 20 broadband seismological stations across the Western Himalaya. The delay times range between 0.2 and 1.7 s, with a cluster around 0.5 s, and the fast polarization azimuths (FPAs) are predominantly oriented in the ENE-WSW direction. The FPAs are orthogonal to the strike of the major thrust faults at the leading edge of the India-Eurasia plate collision. The observed FPAs reveal deformational strain in the lithosphere ensuing from the ongoing Indo-Eurasian collision. The compression in the Indo-Eurasian collision zone of NW Himalaya might be causing the a-axis of olivine in the underthrusting Indian lithosphere to get aligned in the ENE-WSW direction to minimize the strain. There is not much variation in the orientation of FPAs from south to north. However, there is also a possible contribution from the asthenospheric mantle flow due to the motion of the Indian plate.

The seismic data supporting the results are available with the Wadia Institute of Himalayan Geology (WIHG), Dehradun, India, internal data server and is restricted for access and can only be accessed with prior permission from the Institute competent authority.

The authors declare that they have no conflicts of interest.

The authors are extremely grateful to the Director, Wadia Institute of Himalayan Geology, Dehradun, for supporting this work in all respects. This work is carried out under a collaborative project named SIVAH of WIHG, Dehradun, and CSIR-NGRI, Hyderabad. The data used in this paper are collected by the Geophysics Group of WIHG. All supporting staff of WIHG, Dehradun, are accredited for their necessary help and support during the work. The first author acknowledges the Department of Science and Technology, Government of India, for providing Inspire fellowship for conducting her doctoral research.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).

Supplementary data