Velocity Structure and Deep Earthquakes beneath the Kinnaur, NW Himalaya: Constraints from Relocated Seismicity and Moment Tensor Analysis

The optimum 1D velocity model is calculated for the Kinnaur sector of the NW Himalaya utilizing the arrival time information of the local earthquakes (137 no.) recorded with 12 broadband seismic network within the azimuthal gap of ≤ 180 ° . This optimum 1D velocity model is a ﬁ ve-layer model and ranges from the surface to 90 km in the shallow mantle. P velocity varies from 5.5km/s to 8.6 km/s in the crust and upper mantle, and S-wave velocity varies between 3.2km/s and 4.9km/s for the same range. When we relocated the earthquakes with the Joint Hypocenter Determination program incorporating the optimum 1D velocity model, it resulted in a lower RMS residual error of 0.23 s for the hypocenter locations compared to initial hypo71 locations. A total of 1274 P and 1272 S arrival times were utilized to compute station delays. We observed positive variations in P-station delays from -0.19s below the PULG station to 0.11s below the SRHN station. Similarly, for S-station delays, we observed negative delays at each individual site from -0.65s at LOSR station to -0.16s at the SRHN station. This large variation in P- and S-station delays corresponds to the 3D nature of the subsurface below the Kinnaur Himalaya. The relocated seismicity is clustered along the STD fault at sub-Moho and Moho depths ranging between 40 km and 80 km. The seismicity distribution aligned across the strike of the STD and along the strike of the Kaurik-Chango fault (KCF) can be attributed to the cross-fault interactions of the KCF and the STD fault in the area. We also observed bimodal depth distribution of seismicity in the Higher and Tethys Himalayas. The occurrence of earthquakes down to a depth of ~ 0-40km and 50-80 km in the study area can be interpreted in terms of stress contribution from interseismic stress loading associated with the India-Eurasia collision tectonics. The presence of hypocentres in the shallow mantle at 120km depth highlights the strength of the mantle, which seems to be deforming in a brittle manner below the region. The computed focal mechanisms exhibit generally the ﬂ exing of the Indian plate below the Lesser Himalaya with shear parallel to the strike of the MCT and extension orthogonal to it. This study shows deformation over the entire crust and shallow upper mantle levels, with di ﬀ erential stress conditions. Thus, we can consider the crust and the shallow upper mantle down to depths of 120km to be seismogenic in nature and is capable of producing the microseismicity beneath the Kinnaur Himalaya.


Introduction
The collision between the India and Eurasia continental plates has created the~2500 km long, acute shaped Himalayan Range [1,2]. The ongoing underthrusting of the Indian Plate below the Tibetan Plateau starting~55 Ma ago has resulted in the high elevation, the formation of fold-thrust, and crustal shortening of more than 1500 km [3]. Seeber and Armbruster [4] proposed that the convergence in the Himalaya is accommodated by the Main Himalayan Thrust (MHT), a decollement acting as a root to all three splay faults: the Main Central Thrust (MCT), the Main Boundary Thrust (MBT), and the Himalayan Frontal Thrust (HFT). The MHT is divided into two major regions: one seismogenic region that accumulates strain during the interseismic period, and one aseismic, creeping region that does not contribute to interseismic strain accumulation. The seismic portion of the MHT lies beneath the Siwalik and Lesser Himalayas while the aseismic part of the MHT lies beneath the Higher Himalayas. A ramp structure, referred to as the mid-crustal ramp (MCR), joins the two sections of the MHT. The MCR has a dip angle of~16° [5,6]. A locking depth along the Himalaya is denoted by a line of microseismicity that closely follows the 3.5 km elevation contour marking the transition from the locked to the creeping upper surface of the Indian Plate. The geometry of the interseismic locking zone likely controls the microseismicity distribution and coseismic slip on the MHT during great earthquakes. The Himalayan collision zone hosts major active normal fault systems that strike approximately perpendicular to the trend of the orogeny [7][8][9][10][11][12][13][14][15][16][17]. Earthquake moment tensor analysis indicates the existence of normal faults and other extensional features in different parts of the Himalaya [6,[18][19][20][21][22], and the existence of these extensional structures within the greater compressional regime has led to the complex distribution of stress and strain which controls the resulting active deformations of the MHT [23]. This variation in the mechanism of deformation over the MHT at different depth levels resulted in bimodal focal depth distribution of microseismicity along the Himalayan arc [24,25].
An accurate seismic velocity model is essential for understanding the active tectonic deformations and evaluating earthquake hazard. Although the seismic structure of the Kangra and Garhwal sector of the Himalaya have been extensively studied [26][27][28], little is known of the Kinnaur sector. In this study, we utilize recordings from a 12 broadband seismographs located in the Kinnaur region to determine a one-dimensional (1D) compressional and shear wave velocity model and use this model to evaluate the dynamics of the deformation within the source zone of the 1975 Ms 6.8 Kinnaur earthquake. The epicentral region of the Kinnaur earthquake is still experiencing intense seismicity indicated by the high concentration of small to intermediate-magnitude earthquakes ( Figure 1).

Tectonic Setting
From south to north, the Himalaya is divided into four major geologic units: the Siwalik Himalaya (SH), the Lesser Himalaya (LH), the Higher Himalaya (HH) and the Tethys Himalaya (TH). Each of these units is bounded to the south by a north dipping thrust fault: the Himalayan Frontal Thrust (HFT), the Main Boundary Thrust (MBT), the Main Central Thrust (MCT), the South Tibetan Detachment (STD) and the Indus Tsangpo Suture Zone (ITSZ). The Kinnaur region lies between the MCT and the STD faults in the Higher Himalaya. All the major thrusts (HFT, MBT, MCT, STD, and ITSZ) roots into a common detachment or decollement surface, referred to as the Main Himalayan Thrust (MHT) [29]. The MHT detachment accommodates 30%-50% of the~30 mm/yr India-Eurasia convergence [30,31]. A number of smaller scale, tectonic features also exist in the Kinnaur Himalaya including the Kaurik-Chango Fault (KCF), (Figure 1) [17,32]. The MCT for Kinnaur and Himachal Himalaya forms a 2 km thick, continuous, folded, southwest-directed thrust shear zone, which remained active during the Early and Middle Miocene age [33].

Data and Methodology
3.1. Seismic Network. Between 2008 and 2013, a 12-station broadband (BB) seismic network was operated by the Wadia Institute of Himalayan Geology (WIHG) in the Sutlej Valley ( Figure 1). This seismic network was established to study the ongoing microseismic activity occurring in the source zone of the 1975 Kinnaur earthquake [34]. Each site was equipped with a Trillium-240 broadband sensor having a velocity response between 0.004 and 35 Hz and a 24-bit Taurus digitizer (100 samples/s). Table 1 summarises the station information used in the network.

Preliminary
Locations from Hypo71. Earthquakes with a minimum of five P-and five S-phases were filtered with a 1-5 Hz bandpass Butterworth filter and the arrival times picked utilizing the Seisan software [35]. We initially located these events using Hypo71 [36] and the velocity model of Parija et al. [37]. This resulted in a total of 266 earthquakes with an average root-mean-square (RMS) residual error up to 0.8 s. To improve the earthquake locations, we chose the 137 best constrained events for carrying out a travel-time inversion for determining an improved 1D velocity model for the region.
3.3. Minimum 1D Velocity Model (VELEST). An improved 1D velocity model was determined using the VELEST routine [38][39][40]. This method jointly inverts for a 1D velocity model with improved epicentral locations and station corrections. This method solves the forward problem by ray tracing and the inverse problem employing damped leastsquares. The VELEST method comprises two steps: the first step involves computing the complete forward problem and the second solves the inverse problem iteratively to obtain the new 1D velocity model along with the relocated hypocentres and station corrections. The selected 137 earthquake epicentres are bounded between 31.5°and 32.5°in latitude and 77.5°and 78.8°in longitude covering an area of approximately 120 km 2 . The earthquake events selected for the inversion have an azimuthal distribution ≤ 180°and root mean square (rms) residuals of travel time of 0.8 s. The initial hypocentres were mostly shallower than 80 km with a few deeper hypocentres. The initial hypocentres were estimated utilizing a priori velocity model [37] with a velocity constraint maximum up to the depth of 30 km. A set of five different 1D velocity models available for the Himalaya are used as a priori velocity model to compute the final minimum 1D velocity model. The reasons for selecting the five different velocity models for computing the best 1D model for the study area is that these are the only 1D velocity 2 Lithosphere models that are available for the NW Himalaya and adjoining regions. These models are obtained utilizing the local earthquakes that are mostly concentrated in the crust and also range down to the upper mantle levels. Details of the initial velocity models are provided under Table S1 to Table S5 in the supplementary materials. The initial velocity models utilized for the travel time inversion were (1) 1D velocity of Parija et al. [37] for the northwest Himalaya, (2) 1D velocity of Kanaujia et al. [26], (3) 1D velocity of Mahesh et al. [41], (4) 1D velocity model of Monsalve et al. [24] for Nepal, and (5) 1D velocity model of Monsalve et al. [24] for southern Tibet. For each initial velocity model, the process of inversion was terminated when average rms did not change significantly in successive iterations. First, each of the initial velocity models were jointly inverted with a hypocenter and the output 1D velocity models were obtained with a different number of iterations that range between 5 and 10. The initial rms residual error associated with the hypocenter distribution for each of the individual initial and final velocity models was recorded from the first and the last iterations, respectively. We noted that for the 1D model of Parija et al. [37] yield an average rms residual error of 0.27 s in the final iteration in hypocenter locations, which was initially up to 1.04 s during the first iteration. Again similarly, the other two initial models for the Garhwal Himalaya by Kanaujia et al. [26] and the Mahesh et al. [41] resulted in the average rms residual error in the hypocenter locations of 0.42 s and 0.29 s, respectively, in their respective final iterations. Then lastly, the other two initial velocity models of Monsalve et al. [24] for the Nepal  Himalaya and the South Tibet were considered for carrying out the inversion using the same procedure in VELEST and the inversion resulted in 10 iterations for each input velocity model and ended up with a minimum average residual error of 0.66 s and 0.24 s, respectively, for the hypocenter locations. After considering all the five 1D velocity models and carrying out the joint inversion with the VELEST algorithm, we could finally succeed with getting an initial velocity model that can serve as an input velocity model for computing the minimum 1D velocity structure beneath the Kinnaur Himalaya. This final input model is the 1D model of Monsalve et al. [24], which was obtained for the neighbouring Southern Tibet region. Table 2 represents the quality check parameter (average root mean square error) residual for the iteration velocity models.
3.4. Joint Hypocenter Determination (JHD) Method. The joint hypocenter determination (JHD) technique [42] is an efficient location technique that quantifies unmodeled lateral variations in the velocity parameters that are unaccountable in a 1D velocity models used to locate seismic events. This method was first utilized and proposed for earthquake location by Douglas [42] and later improved by Dewey [43], Pavlis and Booker [44] and Pujol [45]. The 1D velocity models often assumed in hypocentre determination are an oversimplified version of the real earth structure. In reality, a velocity structure is quite complex and this complexity reduces the accuracy of the hypocentres which can lead to incorrect conclusions of tectonic structures and their deforming patterns. Thus, it is necessary that we incorporate an appropriate velocity model that has a minimum variation from the actual velocity model for the region. The propagation path can be subdivided into three segments: (1) near the earthquake epicentre, (2) near the recording stations, and (3) along the intermediate ray path between the source and receiver. The variations in the first of these are accounted for by considering groups of earthquakes and performing simultaneous location. The other two errors are accounted for by the term "station corrections" after accounting for irregularities near the earthquake hypocentres [45]. Hypocentre relocation using JHD is an iterative process involving two steps. First, station corrections are calculated with initial hypocentres; second, hypocentres and origin times are determined utilizing the station corrections using a least square regression approach. This process is repeated until a stable solution is achieved. Station corrections observed for both P-and S-wave velocities can be both positive and negative depending upon the lithological variations beneath respective stations. Positive station corrections correspond to lower velocity anomaly and negative station corrections higher velocity anomaly along the ray path, respectively. Here, we have utilized the JHD method to relocate the earthquake hypocentres with the 1D VELEST velocity model obtained in Section 4.2. The final JHD relocations of the 137 earthquakes resulted in an average rms of 0.23 s which is a marginal improvement over the hypocentre location with the VELEST algorithm. This suggests that the 1D VELEST Pand S-wave velocity model is suitable for locating the hypocentres in the Kinnaur region.

Source Moment Tensors (MT).
To better understand the seismotectonics of the Kinnaur region, we calculated the source mechanism for the moderate earthquakes (M ≥ 3:0). The event locations used in computing the moment tensor solutions are from the JHD relocations. The ISOLA code (https://www.mathworks.com/products/ MATLAB) is used for computing the source mechanism using a full waveform inversion approach [46]. Retrieving the waveform inversion approach using the entire waveform is more reliable and strongly helps in overcoming the difficulties in source characteristics and understanding the complex tectonics of the study region [6]. The method utilizes the iterative deconvolution method by Kikuchi and Kanamori [47] modified for regional distances and newly encoded by Zahradnik and Plesinger [48]. This method utilizes the discrete wavenumber method [49] to calculate the Green functions in a given 1D velocity model. This study is based on the single-point source and deviatoric inversion. The moment tensor inversion follows a grid search over a set of trial source positions and time shifts for obtaining the optimal centroid position, time, and corresponding moment tensor with minimum misfit. This misfit results from the correlation between the observed and the computed synthetic seismogram. The results obtained are represented by the percentage of double-couple (DC) and the percentage of non-double-couple (non-DC) components, which is considered a compensated linear vector dipole (CLVD). The results obtained are stated in terms of the double-couple component of the deviatoric solution, represented by the scalar moment, strike, dip, and rake.
In this study, we selected earthquakes at local to regional distances from the recording stations. The waveforms are first preprocessed and only the recording stations with high signal-to-noise ratio and greater variance reduction in terms of the synthetic-observed waveforms fit are considered for inversion. Table 3 summarizes our MT results.

Location Results with the Five Different Starting Models.
The minimum residual velocity model was determined using 137 events located within the seismic network with azimuthal gap ≤ 180°. Here, we have first computed the travel times (calculated) of the hypocentres utilizing the initial velocity models and the error in hypocenter locations and origin times are reflected in terms of misfit (or residuals). This misfit is generally calculated by the RMS (root-meansquared) norm. The RMS value corresponds to any possible combination of hypocentres, velocity model, and the station corrections. First, we started with locating the earthquake hypocentres with the 7-layer 1D velocity model by Parija et al. [37], which resulted in an initial average rms residual error of 1.04 s in travel times in the first iteration. We continued with the joint inversion of hypocentres and velocity model until five iterations under simultaneous mode, and with each iteration, this rms misfit significantly reduced, and finally, it got to 0.27 s for the last iteration. From this inversion, we obtained a 1D velocity model in the final iteration that comprised of seven-layer velocity parameter for both P-and S-wave down to 30 km depth from the surface. The velocity model has a velocity layering thickness of 5 km from the surface down to 30 km depth. The earthquake hypocentres located with the computed final output velocity model calculated from the initial velocity model of Parija et al. [37] were placed in a depth range extending from 1.5 km down to 120 km depth. This study assigns a deep Moho beneath the Kinnaur sector, which can be approxi-mately~60-70 km.
After this, we tested with another available velocity model for the nearby Tehri region of the Garhwal section of the northwest Himalaya proposed by Kannuajia et al. [26] utilizing the local seismicity. The model determined in the VELEST analysis consists of six layers between the surface and 24 km depth. The crust is 46 km thick. The model also observed low velocity at 12 and 14 km depth, respectively, which it ascribed to the fractured basement thrust representing the upper surface of the Indian plate [26].
We started our inversion with this velocity model, which resulted in an initial average rms misfit of 0.92 s in the hypocenter location for the first iteration. This rms misfit in the hypocenter location reduced to 0.42 s for the last final iteration no. 5. This final rms error in hypocenter parameters is comparatively higher by about 0.15 s than the final output model obtained from inverting the seven-layer velocity model of Parija et al. [37]. This final model obtained placed the hypocentres between 1.5 km and 120 km, respectively. We observed that there is not much variation in depths for the final hypocentres with both the 1D velocity models of Parija et al. [37] and Kannuajia et al. [26]. However, this study reported a Moho depth at near about 50 km below the Garhwal Himalaya region as observed from the sharp velocity contrast at this depth. The insignificant changes observed in the rms residual error associated with hypocenter parameters compelled for a consideration of another local P-and S-wave 1D velocity model in the area. Therefore, we considered the 1D velocity model calculated by Mahesh et al. [41] for the joint inversion of hypocenter parameters and velocity model under simultaneous mode in the VELEST program. The inversion process was terminated when the final solution was reached after five iterations. The rms residual error associated with the hypocenter parameters reduced from 1.13 s in the first iteration to 0.29 s for the final fifth iteration. This model does not report any low velocity within 20 km of the crust during the inversion. This model also calculated a depth in the wide range extending from very shallow depth of 0.1 km to 130 km below the area. After completing the joint inversion of hypocentres with the three different initial P-and S-wave velocity structures, we obtained three different output velocity models from each of the initial velocity models and coupled hypocenter inversion.
In the next step, we fused all the three resultant velocity models to form a single input velocity model that has nonuniform layer spacing of 2 km, 4 km, 6 km, and 8 km, respectively. We have also tried to look into the dependency of the misfits on the introduction of low velocity at a depth of 10 km (case 1) and in absence of a low-velocity layer at the same depth of 10 km (case 2). For the first trial, we introduced a uniform layer spacing at 2 km down to 10 km depth in the upper crust. We observed that the rms error in the hypocentres and origin times reduced from 0.42 s in the first iteration to 0.24 s for the third and final iteration. In the next step, we tried with case 2, which is not having a low velocity at 10 km depth interface. This iterative inversion resulted in an average rms residual error of 0.43 s and 0.28 s for the first and the last iterations, respectively. The layer spacing remained the same as it was for case 1. We also tried with different combinations of layer spacing at upper and midcrustal levels but ended up with increasing rms residual error in hypocenter locations and origin times.
After utilizing the three above preliminary velocity models and their resultant fusion model for inversion, we observed that the hypocentres are present down to deeper levels maybe in the shallow mantle level for the Kinnaur Himalaya as reported earlier for the Nepal and southern Tibet region [24,52].
In addition to that, we know that all the three initial velocity models and their resultant fusion model are extremely shallow in nature and constrain velocity up to 46 km, which is well above the Moho depth in the area. As this analysis and the preliminary locations of hypocentres in the area signifies a number of sub-Moho and shallow mantle events, we decided to proceed with the joint inversion of hypocenter locations and the 1D velocity model by Monsalve et al. [24] for the Nepal and the Southern Tibet area, which had velocity constraint down to shallow mantle depths of 90 km.

Development of VELEST Model.
First, we considered the 1D velocity model by Monsalve et al. [24] for the Nepal Himalaya as the initial input velocity model for inversion. This model resulted in an average rms misfit of 0.25 s in the hypocenter locations for the tenth and final iteration, which was earlier up to 1.20s for the first iteration. Then, we tried utilizing the second velocity model proposed for 6 Lithosphere the southern Tibet region. This model registered the minimum rms residual error in the hypocenter locations in its final tenth iteration, which is 0.24 s (Figure 2). We decided to take this model as our final input velocity model (  (Table 5). This optimum 1D velocity model is a five-layer model and ranges over a depth of 0 km down to 90 km in the shallow mantle. P velocity varies from 5.5 km/s to 8.6 km/s in crust and upper mantle, and and S-wave velocity varies between 3.2 km/s and 4.9 km/s for the same region ( Figure 3, Table 5).

Variations in Lateral Velocity and Station Corrections.
Another important parameter that gets constrained along with the 1D velocity structure is the station corrections, which corresponds to the lateral variations in the velocity structure beneath respective stations. This station correction provides information about the nature of the subsurface below individual seismic stations. Station corrections can be either positive or negative depending upon the properties of the materials constituting the subsurface and are represented in the form of station delays. Positive variations are observed where the actual velocity is less than the predicted one and vice versa. We computed station delays at 10 BBS stations spread over an area of 200 × 200 km 2 for the Kinnaur Himalaya. A total of 2546 P and S arrival times were utilized to compute the station delays below 10 seismic stations in the study area. We observed positive variations in P-station delays from -0.19 s below the PULG station to 0.11 s below the SRHN station. Similarly, for S-station delays we observed negative delays at each individual site from -0.65 s at LOSR station to -0.16 s at SRHN station ( Table 6). This large variation in P-and S-station delays corresponds to the 3D nature of the subsurface below the Kinnaur Himalaya. Figures 4(a) and 4(b) completely depicts the contour map of P and S delay at various stations lying within the array.

Seismicity from JHD Relocations.
Studying the spatial and temporal distributions of seismicity is an essential step towards understanding the complex geodynamics of the Himalayan subduction zone. Seismicity is the reflection of the geometry of deformation going on along the fault structure in the subsurface. Thus, the accurate hypocenter locations would help us to image the depth at which the fault is deforming along its plane. The hypocenter parameters of all the 137 earthquake events were relocated using the JHD, which incorporated the estimated optimum 1D velocity model and station corrections. The relocated seismicity   shows a tight clustering of earthquake events in the Higher Himalayas for the study region ( Figure 5). The relocated seismicity has a magnitude range (Ml) between 2.0 and 5.6. We observe a significant decrease in uncertainty of hypocenter parameters as compared with the initial locations obtained with Hypo71. The uncertainty in longitude, latitude and depth was limited to 5 km maximum, which was initially up to 15 km for Hypo71 locations (Figure 6). The depth distribution of the relocated hypocentres shows a very dense concentration throughout the entire crust down to Moho depth, which from this study is estimated to lie at 70 km approximately in the region. The hypocentres are clustered along the STD fault at sub-Moho and Moho depths range between 40 km to 80 km. The seismicity distribution aligned across strike of the STD and along strike of the Kaurik-Chango fault (KCF) can be attributed to the cross-fault interactions of the KCF and the STD fault in the area. The seismicity in the area is mainly influenced by the deformation of the Kaurik-Chango fault, which produced the last significant earthquake of Ms 6.8 in the region. The clustered seismicity is showing an E-W extension at the mid to lower crustal level from 40 to 60 km centroid depth. The dense clustering of earthquake hypocentres at lower crustal depths signifies that the lower crust has enough strength to host the microseismicity. This observation is significant as for other adjacent regions of Garhwal Himalaya or Nepal Himalaya, it is proposed that most of the seismicity is hosted by the upper crustal stress and the lower crust is ductile in nature. We also report five earthquakes at shallow mantle levels lying between 80 km to 120 km depth below the region. The location of earthquake hypocentres at this shallow mantle level advocates for the brittle nature of the    [6]. Earthquake residing within shallow focal depth beneath the upper part of the Higher Himalaya corresponds to the ongoing active deformation over the midcrustal ramp [53], which is dipping at anglẽ 20°-30°and is present depth range of ∼15-20 km [5,54]. Arora et al. [34] marked seismicity as a possible proxy to identify the geometry of the under-thrusting Indian plate beneath the Himalayan wedge. Absence of this seismicity to south of the MCT illuminates the fact that the detachment is almost flat in this locality with a dip angle ranging from 2.5°to 4° [55,56]. The presence of lower crustal earthquakes in the study region is due to dry granulite rocks, a dominant constituent of the subducted Indian crust. It became brittle when deformed under conditions corresponding to the Eclogite stability field [57]. This method is based on the iterative deconvolution of Kikuchi and Kanamori [47], modified for regional distances and newly encoded by Zahradnik and Plesinger [48]. We computed the moment tensor inversion of eight moderate earthquakes utilizing the ISOLA code, which is a full waveform inversion approach [46]. In the current study, we carried out the moment tensor analysis of eight earthquakes with a magnitude (Ml) ranging between 3.0 and 4.9 in the Kinnaur Himalaya region. The moment tensor calculations were performed utilizing the high signalnoise ratio and frequency range of 0.01-5.0 Hz as this is the lowest frequency with highest signal-noise ratio. The maximum fit between the observed and the synthetic seismogram is measured in terms of its correlation. We have considered modeling in low frequency, because at low frequency the results are independent of the accurate constraints of the crustal structure like if the velocity is varying 1D or 3D, it does not affect the modelling. We consider that if the total variance reduction is positive for all seismic stations used in the inversion then the solution is accepted along with the calculated DC%. Here, some events have a higher CLVD%, which is also earlier observed by Parija et al. [6], the reason they mentioned for such observations is the presence of fluids or geothermal structures, and we would add up to that it may be because of mantle materials getting trapped within the crust.
In this study, we consider a frequency range between 0.01 Hz and 3.0 Hz for inversion. The broadband seismic data of minimum three to maximum of six stations were utilized for computed the focal mechanism solutions. The computation of moment tensors follows certain steps. The first step involves the data preparation, where the data is converted to adequate format like the SAC or PITSA before being considered for further analysis. Here, we have converted our SEED data to SAC by utilizing the wave-tool program encoded in the SEISAN software package. Then, we have provided a 1D velocity model, which we have estimated from the present study using the crust-mod program of the ISOLA code. This model is utilized further to generate the Green functions utilizing the green pretool that is utilized for generating synthetics for inversion. The inversion was carried out in the frequency domain having a frequency range of 0.01-3.0 Hz with cosine tapering applied to both ends. Additionally, for tectonic interpretation, we have also included the five earthquakes recorded by the USGS and reported in Parija et al. [6] in the Kinnaur Himalaya. These focal mechanisms are mainly dominated by thrust and normal fault plane solutions and only a few strike-slip mechanisms. In all the cases, we obtained stable moment tensor estimates ( Table 3) that reproduce well near-regional waveforms with broad azimuthal distribution. The depth distribution for all the estimated and reported earthquakes is between 5 km and 60 km. This shows that this region is characterized by shallow and deep deformation mechanisms in contrast to other regions in the Himalaya, which are only characterized by shallow deformations. The MT solutions of the eight earthquakes determined from the present study plus those from the previously reported by [6], USGS GCMT catalogue, and [50,51] are described in Table 3. An example of the observed and synthetic match for the regional waveforms of earthquake having a magnitude of M 4.9 in the NW Himalaya along with its correlation with the Double couple (DC %) has been shown in Figure 7.

Discussion
Our study provides new and interesting insights into the seismotectonic of the Kinnaur Himalaya region. The improved 1D velocity model has resolved the layers down to 90 km in the upper mantle and assists in constraining the focal depths of events occurring in this collision region. Additionally, it provides a distinct discernment about the seismicity pattern, which is associated with the MCT and STD fault zone along with the KCF rift. We observed bimodal depth distribution of seismicity in the Higher Himalayas: one cluster down to 50 km depth and another from 50 km to 100 km depth. The occurrence of earthquakes down to 50 km in the study area can be interpreted in terms of stress contribution from interseismic stress loading 10 Lithosphere associated with the India-Eurasia collision tectonics.
Another reason for the occurrence of seismic activity is the region down to 50 km is presence of the Kaurik-Chango rift structure, which extends further south into the Higher Himalaya and plays a significant role in stress distribution in the Himalayan seismic belt (HSB) [34,58]. The earthquake hypocentres present below the Higher Himalaya down to lower crustal depths of 50 km suggests that the entire crust is seismogenic in nature [6,59]. This lower crustal distribution of seismicity indicates the interaction between the secondary fault structures with the Main Himalayan Thrust (MHT) system. The moment tensor focal mechanisms computed along with the reported earlier shows E-W extension mechanism for deep lower crustal earthquakes up to 60 km. Our computed mechanisms for the area are in agreement with what was earlier reported by Jackson    Figure 8: A summary of the tectonic processes across the Kinnaur Himalayas determined in this study. The seismicity (black dots) is from the JHD analysis. The red focal mechanisms are the eight moment tensor solutions determined in this study, the blue focal mechanisms are from Parija et al. [6], the brown focal mechanisms from the USGS GCMT, and the purple focal mechanisms are from Dziewonski [50,51]. The coloured sections of the focal mechanisms denote the compressional quadrant. The NE-SW transect AB is marked in the figure. 11 Lithosphere [60] for the southernmost Himalaya and foreland. The earthquake hypocentres lying at a depth from 50 km to 100 km suggests that the shallow upper mantle beneath the Higher Himalayas is subjected to high differential stresses between the lower crust and the uppermost mantle as observed for Ganga basins and southernmost Tibetan plateau [52]. The deepest earthquake of the region is at a depth of 120 km, located in the upper mantle, and it deforms by brittle manner this may be due to the presence of hydrous materials at temperatures less than 600°C as observed for lower crust [25,61,62]. This characteristic of the deformational process is marked opposite to the rest of the Himalaya towards the south i.e. the Indo-gangetic Plain, Siwalik Himalaya and Lesser Himalaya. Towards the south of the STD fault, we observe other deep crustal earthquakes as compared to shallower ones along the MCT with extensional E-W trend due to normal faulting beneath Higher Himalayas. This variation in earthquake depth in the study region can be due to the spatial variation of stress among the crust and the upper mantle between the two adjacent tectonic units.
We have obtained 14 moment tensor focal mechanism solutions for the study area ( Figure 8). Out of these 14 focal mechanisms, one for the Lesser Himalayas, five for the Higher Himalayas, and eight are reported for the Tethys Himalayas. The centroid depth reported for the earthquakes are between 5 km and 60 km, respectively, for the three different tectonic units. The focal mechanisms along with their source parameter details are described in Table 3. The focal mechanism is mostly dominated by normal faulting with N-S compression and E-W extension for the study region.
The focal mechanism with corresponding ID no. 5 obtained for the Lesser Himalaya has a shallow focal depth of 5 km and the P-axis is oriented N-S with shallow plunge and T-axis is subhorizontal with deep plunge. This characterizes ongoing shear parallel to the strike of the Main Central Thrust (MCT) and E-W extension. This also accounts for a flexed Indian crust below the Lesser Himalayas in the region. This trend is observed for four other earthquakes with ID. No. 4, 6, 8, and 12 having focal depths between 5 km and 25 km, respectively (Figure 9). But an exception to this is an earthquake at a deeper depth of 60 km with ID  Figure 9: Seismicity and focal mechanisms along the cross-section AB indicated in Figure 8. The profile crosses the major tectonic segments starting from the Main Central Thrust (MCT) in the south to South Tibetan Detachment (STD) in the north. This cross-section characterizes the style of deformation beneath the Kinnaur section of the NW Himalaya. The seismicity is shown as red solid circles and its uncertainty in depth as green solid bars. The maximum depth of the investigation below the cross-section is 120 km. The colour of the compression quadrants for the focal mechanisms is the same for both the map view as well as for the cross-sectional view. The focal mechanisms are represented as side-on views.
12 Lithosphere no. 7 showed the opposite trend with the pension axis being oriented N-S and the pressure axis being oriented E-W. This shows that the Indian plate flexes at lower crustal level down to 60 km below the Higher Himalayas. Further moving to the north of the South Tibetan Detachment (STD) fault, we have about eight focal mechanisms represented by ID no. 1, 2, 3, 9, 10, 11, 13, and 14, respectively. They have a centroid depths range between 5 km and 60 km, respectively. All the focal mechanisms show similar orientations of the Pand T-axis irrespective of their focal depths. They show compressional P-axis-oriented N-S and tensional t-axis E-W oriented for the entire Tethys Himalaya region. This observation accounts for a shear and E-W extension at the base of the Indian lithosphere as reported from upper mantle anisotropy [63]. These observations highlight the fact that the subduction process is related to the earthquakes at Moho depths for the Kinnaur Himalaya region.

Conclusion
In the present study, we have utilized the seismic events recorded through the local broadband seismic network operated in the Kinnaur Himalaya to derive the velocity structure of the region by travel time inversion of P and S waves, respectively. The obtained minimum 1D velocity model can be used as an input model to derive the 3D velocity model below the region. This analysis utilizing the newly computed minimum 1D model resulted in relocated hypocentres throughout the entire crust down to Moho depth, which from this study is estimated to lie at 70 km approximately in the region. The hypocentres are clustered in the Higher and Tethys Himalayas along the STD fault at sub-Moho and Moho depths range between 40 km and 80 km. The seismicity distribution aligned across strike of the STD and along strike of the Kaurik-Chango fault (KCF) can be attributed to the cross-fault interactions of the KCF and the STD fault in the area. The seismicity in the area is mainly influenced by the deformation of the Kaurik-Chango fault, which produced the last significant earthquake of Ms 6.8 in the region. The presence of hypocentres in the shallow mantle highlights the strength of the mantle, which seems to be deforming in brittle manner below the region. The computed focal mechanisms exhibit generally the flexing of the Indian plate below the Lesser Himalaya with shear parallel to the strike of the MCT and extension orthogonal to it. This study shows deformation over the entire crust and shallow upper mantle levels, with differential stress conditions. Thus, we can consider the crust and the shallow upper mantle down to depths of 120 km to be seismogenic in nature and is capable of producing the microseismicity beneath the Kinnaur Himalaya.

Data Availability
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.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.