Tomographic Image of Shear Wave Structure of NE India Based on Analysis of Rayleigh Wave Data

The major scientific purpose of this work is to evaluate the geodynamic processes involved in the development of tectonic features of NE India and its surroundings. In this work, we have obtained tomographic images of the crust and uppermost mantle using inversion of Rayleigh waveform data to augment information about the subsurface gleaned by previous works. The images obtained reveal a very complicated tectonic regime. The Bengal Basin comprises a thick layer of sediments with the thickness increasing from west to east and a sudden steepening of the basement on the eastern side of the Eocene Hinge zone. The nature of the crust below the Bengal Basin varies from oceanic in the south to continental in the north. Indo-Gangetic and Brahmaputra River Valleys comprise ∼5–6-km-thick sediments. Crustal thickness in the higher Himalayas and southern Tibet is ∼70 km but varies between ∼30 and ∼40 km in the remaining part. Several patches of low-velocity medium present in the mid-to-lower crust of southern Tibet along and across the major rifts indicate the presence of either partially molten materials or aqueous fluid. Moho depth decreases drastically from west to east across the Yadong-Gulu rift indicating the complex effect of underthrusting of the Indian plate below the Eurasian plate. Crust and upper mantle below the Shillong Massif and Mikir Hills are at a shallow level. This observation indicates that tectonic forces contribute to the uprising of the Massif.


INTRODUCTION
A long span of convergence between Indian and Eurasian plates has caused subduction, underthrusting, and extrusion resulting in widespread lithospheric deformation (Molnar et al., 1977;Molnar, 1984;Tapponnier et al., 2001;DeCelles et al., 2002;Nábelek et al., 2009). Recently and in the past, the NE Indian region has experienced a frequent occurrence of moderate to large-size earthquakes due to ongoing active tectonic convergence (Ni and Barazangi, 1984;Pandey et al., 1995;Parvez and Ram, 1997;Kayal, 2008;Baruah et al., 2014Baruah et al., , 2016Le Roux-Mallouf et al., 2020). In response to the collision, the lithospheric part of these three plates experiences horizontal shortening and surface uplift. Elevated Tibetan plateau, Indo-Burma Ranges (IBR), and the Himalayas are the result of the collision and crustal thickening, which is associated with mountain building phenomena and crustal shortening (Burchfiel and Royden, 1985;England and Houseman 1989;Murphy and Copeland, 2005;Thiede et al., 2006). The study region includes Brahmaputra River Valley (BRV), Shillong Massif (SM), Mikir Hills (MH), Indo-Burma Ranges (IBR), Bengal Basin (BB), southern Tibet, eastern Himalayas, and part of central Himalayas ( Figure 1A), characterized by the diverse geological setup and a very high variation in topography. The overall study region is demarcated by latitude 22°N to 32°N and longitude 85°E to 98°E. Detailed geological and tectonic studies were carried out by Krishnan (1960), Evans (1964), Desikachar (1974), and Nandy (1980).
In NE India, the Indian plate is supposedly simultaneously underthrusting beneath Eurasia and subducting below Burmese microplate in the north and east, respectively (Kayal, 2008;Sarkar et al., 2013;Kumar et al., 2015;Raoof et al., 2017). Le Dain et al. (1984) proposed that in the early Cenozoic period, the Indian plate subducted below the Burma microplate and presently the hanging lithosphere is being dragged northward. Mukhopadhyay and Dasgupta (1988), Kayal (1989Kayal ( , 1996, and Satyabala (2003) opined that the Indian plate is actively subducting below the Burma plate. GPS measurements in this area and Myanmar suggest EW convergence with a shallow dip of the Indian plate beneath Myanmar, which is locked along a megathrust further down-dip (Gahalaut et al., 2013;Steckler et al., 2016). Steckler et al. (2016) and Gahalaut et al. (2013) suggest that the active subduction has stopped, whereas Steckler et al. (2016) inferred active subduction despite the highly oblique plate motion and presence of thick sediment.
The convergence of Indian and Eurasian plates causes N-S compression mainly within the Himalayas and E-W extension within the Tibetan plateau (Armijo et al., 1986;Angelier and Baruah, 2009;Sarkar et al., 2013). Several workers suggested that low-velocity material is present within southern Tibet as a result of EW extension along several rifts (Kind et al., 1996;Cogan et al., 1998;Jiang et al., 2014) due to clockwise deformation with respect to the northeast corner of NE India. Extension-related normal faults have been identified (Kayal, 2001(Kayal, , 2008Hintersberger et al., 2010) mostly within the Tibetan plateau and also in some cases in the Higher and Lesser Himalayas. Focal mechanism solutions in the eastern Himalayas expose the direction of compression from NNW-SSE to N-S from Arunachal Pradesh, Mishmi hills to Bhutan. However, in the IBR, NE-SW directed average compression is observed (Angelier and Baruah, 2009). IBR is under compression due to oblique subduction of the Indian plate below the Burma plate, where the NE-SW directed stress across the inner and northern arc is switched to the E-W directed stress near the Bengal Basin. Within NE India and SM, Baruah et al. (2016) made an investigation on tectonic stress through inversion of the stress tensor. They observed that compressional stress direction changes from NNW-SSE in the western part of the massif to the NNE-SSW direction in its eastern part. Oblique subduction of the Indian plate underneath the IBR causes compression along the NNE-SSW direction (Baruah et al., 2013(Baruah et al., , 2016. The available data set make it possible to study the crustal and upper mantle structure of Lhasa/southern Tibet and IBR. Eastern Himalayan syntaxis (EHS) is a complex triple junction of Indian and Eurasian plates with the northern end of the Burma plate (Curray, 1989). A study of active tectonics in EHS was done by Holt et al. (1991). They identified an oblique thrust based on focal mechanisms and it shows the movement of the Indian plate directed toward NNE as it underthrusts below this region. Baranowski et al. (1984) found that radial directed slip in the Himalayas causes extension within Tibet; however, it is possibly half of the rate of under-thrusting taking place within the Himalayas. Holt et al. (1991) found that slip vectors are generally not orthogonal to the strike of the Himalayas.
Different workers have proposed different tectonic processes for the formation of SM. Evans (1964) postulated that the Dauki fault is a strike-slip fault along which SM has moved ∼250 km eastward to its present position. Oldham (1899) concluded that movement along a thrust-plane or thrust planes and along secondary thrusts and fault planes may be the cause of the 1897 great Earthquake. Bouguer and isostatic anomalies indicate that high-density crustal/upper mantle strata are present beneath SM (Verma and Mukhopadhyay, 1977). Kailasam (1979) suggested that isostatic adjustment may have caused uplift of SM but positive Bouguer and isostatic anomalies over it imply that it should sink. Dziewonski and Anderson (1983) estimated the P-wave travel time residual beneath the Massif to be −0.57 s, indicating that the subsurface has higher velocity material than the surroundings. Many researchers suggested that the uplift of SM is influenced by active tectonics between NE India and the NE Himalayas (Chen and Molnar, 1990;Khattri et al., 1992;Mukhopadhyay et al., 1993Mukhopadhyay et al., , 1997Najman et al., 2016;Govin et al., 2018). Mukhopadhyay et al. (1997) opined that the Massif has overthrusted southward above the Indian plate along the Dauki fault. Bilham and England (2001) and England and Bilham (2015) proposed that SM uplifted abruptly by ∼11 m during the 1897 Shillong earthquake and suggested that it is a "pop-up" structure, whereas Raoof et al. (2017) suggested that it is the crest of the buckled-up part of the Indian plate. Clark and Bilham (2008) suggested that convergence of India and Eurasia in this part is shared between the Bhutan Himalayas and the Shillong Massif. Strong et al. (2019) stated that the Massif is affected by active uplifting as well as erosion.
From this discussion, it is obvious that there are several postulates about the tectonics and geodynamic characteristics of this area that has led to such a complex crustal/lithospheric structure. Hence, this area needs detailed investigation. The availability of a dense seismic network in NE India makes it possible to carry out such investigations for this region. The study FIGURE 1 | Tectonic map of NE India and the surrounding area (GSI, 2000). Red symbols (shapes same as corresponding curves): location of six points for which Rayleigh wave dispersion curves are shown in region is selected because it shows complex variations in plate motion that may affect the subsurface structure and properties of the lithosphere. In this study, we decipher subsurface characteristics to understand the effects of NS underthrusting and EW subduction in NE India and its surroundings. Most of the geophysical investigations were carried out on collision zones either toward the north or northeast of SM (Singh et al., 2015 and the references therein). However, the southern side of SM and that of BB have not been studied much due to unavailability of sufficient data. Thick sediments within BB make it challenging to delineate the underlying oceanic and/or continental crust (Alam et al., 2013). Sediment thickness and crustal structure of BB have been investigated by Currey (1991), Brune et al. (1992), Mitra et al. (2011), andSingh et al. (2016). The extent of sediment deposition and the nature of the crust in BB are required to be properly investigated. We have estimated the lateral variation in sedimentary layer thickness in BB. Although Bangladesh seismological data were not incorporated in the present work, there are stations in the surrounding part of India that gave good coverage of this zone. We made use of new data sets covering NE India and surrounding areas and tried to solve several issues, such as Moho depth variation, nature of crust and sedimentary layer thickness variation in BB, deformation, and upliftment of SM in response to active NS and EW compression, geodynamic inferences related to oblique subduction of Indian plate below IBR, indentation of EHS and occurrence of low-velocity material in mid-crust of SE Tibetan plateau. To obtain the subsurface crustal and lithosphere S-wave velocity (Vs) image of the study area, surface wave tomography technique developed by Ditmar and Yanovskaya (1987), Yanovskaya and Ditmar (1990), and Yanovskaya et al. (1998) is applied. Tomographic image of the study area is obtained up to ∼90 km depth, which elucidate the geodynamics of NE India through high-resolution data.

DATA
Surface wave data recorded by a network of 20 broadband seismic stations in NE India and three other stations (Kolkata, Bhubaneswar, and Bokaro) operated by the India Meteorological Department (IMD, India) were used in the present study ( Figure 1B; Table 1). The last three stations provide better path coverage for BB. The instruments used in these seismic stations (Table 1) include Trillium-240, STS-2, RT151-120, and CMG40T with natural periods of 240, 120, 120, and 30 s, respectively. The data for 201 earthquakes ( Figure 1B, Supplementary Table S1) used for this study were recorded at local and regional epicentral distances in the range of ∼160-∼1800 km. Surface waves of good signal strength recorded for wave period >1 s generated by earthquakes of focal depth ≤50 km and magnitude ≥5 were used. Shallow focal depth earthquakes generate higher surface wave amplitudes, and therefore even lower magnitude (M∼5) records at regional distances have good surface wave signals. Based on selected stations and event locations, we have obtained maximum 753 source-receiver paths ( Figure 1B; Table 2), which covers a broad range of azimuths and path lengths within the area of investigation. Supplementary Figure S1 shows ray path coverage for periods between 8 and 60 s.
Data processing was performed in four steps: 1) estimation of fundamental mode Rayleigh wave dispersion curves, 2) construction of group velocity maps at different grid points of specified spacing for good resolution, 3) inversion of dispersion curve of each grid for obtaining 1D S-wave (Vs) velocity structure, and 4) construction of 2D and 3D images of Vs.

ESTIMATION OF GROUP VELOCITY
Pre-processing steps were carried out to remove the mean value, linear trend, and avoid spectral leakage (using tapering window) from raw data. The data were also converted into ground motion (displacement) units using the frequency response of the recorders. Fundamental mode dispersion curves of Rayleigh wave were estimated along each path between earthquake epicenter and the seismic station using the multiple filter technique (MFT) developed by Herrmann (2013) and Herrmann and Ammon (2002) based on the methodology of Dziewonski et al. (1969). MFT is a fast and efficient method used for analyzing multiple dispersed signals (Herrmann, 1973;Erduran et al., 2007). In this method, the Fourier transform of the signal selected for a given time window displays amplitude variation with frequency. In the next step, Gaussian filter is applied and followed by inverse Fourier transform to obtain the peak of the spectrum for a given time and frequency. The MFT uses phase-matched filtering to extract the dispersed waveform. It separates the other parts of waveform consisting of body waves, higher modes, and multipath effects (Kolinsky, 2004). Examples of dispersion curves are given in Figure 2.
Selection of range of periods of dispersion curves is based on the natural period of the seismometer, data sampling rate, noise level, and the distance between the epicenter and station along the great circle arc (Yao et al., 2006;Bensen et al., 2007;Yang et al., 2007;Guo et al., 2009;Kumar et al., 2017). The upper limit of the period up to which group velocity can be estimated also depends on the ratio of epicentral distance and the maximum wavelength of the recorded signal. The reliable maximum period of surface wave recorded at a given epicentral distance should have at least three wavelengths within that distance (Yao et al., 2006;Bensen et al., 2007;Yang et al., 2007;Guo et al., 2009). Based on these criteria, we could extract dispersion data for periods between 4 s and 60 s. Wave propagation path coverage for a particular period depends on the availability of good quality data. Path coverage in this study region varies (Supplementary Figure S1), and it has a maximum value at 12 s period ( Figure 1B). In total, maximum 795 Rayleigh wave records were used. Minimum and maximum numbers of paths used for tomographic inversion of group velocity values are 129 and 753 at periods of 4 s and 12 s, respectively ( Table 2).

GROUP VELOCITY TOMOGRAPHY Tomographic Inversion Method
Surface wave dispersion data were used to generate a map of spatial variation of group velocity for different periods within the study region. Surface wave velocity tomography maps are constructed through the inversion technique developed by Ditmar and Yanovskaya (1987) and Yanovskaya and Ditmar (1990). Dispersion curves obtained from the earthquake data are path specific, that is, each one gives average group velocity versus period for that particular path between earthquake epicenter and station location. It gives the group velocity for a particular period at different grid points of the study region. This inversion technique is robust and does not need any initial parameterization or truncation of Taylor series expansion of nonlinear equation because of its straightforwardness (Levshin et al., 1989). The basis functions for the inversion model are formed as a superposition of kernels of the group velocity travel time integral . Measured group velocity is a function of the location described by U(θ,ϕ) on a spherical surface. It is decomposed in different sections across the study region by taking the mean group velocity as a starting value and applying region-specific perturbations. Decomposed reference value U 0 and its location dependence perturbation are expressed by Initially, this method leads to a search for smooth perturbation in group velocity from the reference model U 0 so that U(θ,ϕ) gets fitted by N observed group travel times (t obs i ). Smooth perturbation is based on the weighted least-squares technique. The above approach aims to minimize error function between observed and predicted values at each period and wave type through Eq. 2.
Here, P i is the ith wave path, w i is the associated weight through a map of group velocity U(θ,ϕ), t pred i the predicted group travel time for the same path, and "S" denotes the study region. In Eq. 2, FIGURE 2 | Examples of dispersion curves obtained at four stations. Earthquake source parameters and the recording station names are mentioned on the dispersion curve plot. The range of the data of the dispersion curve used for further analysis is shown with a gray-color curve. The color in a particular contour defines the energy with blue representing minimum energy and its energy increasing trend from blue to cyan to green to yellow to red. Raw waveform is given in the right-hand rectangular block.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 680361 c is the regularization parameter and its selection controls the trade-off between the fit to the data and smoothness of the observed velocity map for each assessed period of the wave Karagianni et al., 2002). The linearized inverse procedure is carried out with very small relative group velocity variations at each iteration for each period and wave type such that δU(θ,ϕ)/U 0 <<1. Tomographic images of Rayleigh wave group velocity were obtained at periods ranging between 4 s and 60 s on 1°× 1°grid spacing. The regularization parameter is chosen such that the velocity range in the group velocity maps is identical or larger than the range of measured group velocities (Yanovskaya et al., 1998). To obtain an appropriate smoothed map, we have used several regularization parameters (c) ranging between 0.05 and 0.4. The lateral variation in the group velocity and associated standard error for 12 s period for the different six regularization parameters are shown in Supplementary Figures S2, S3, respectively. Decrease in c gives a sharp solution with large standard (solution) error and vice-versa. An increase in c enhances the smoothness of the solution region (Supplementary Figure S2). Finally, the value of c 0.3 gives a relatively smoothed map with lesser solution error.
For best results in a tomographic study, the area of investigation should be covered uniformly by ray paths from all directions. Unfortunately, that is very rarely possible. Uneven distribution of ray paths may introduce some anomalous values to the solution. The solution is obtained after rejecting the data of a path with anomalous values categorized based on regularization parameters. The anomalous value represents a large deviation from the average velocity of a particular localized region. This error may be due to unresolvable lateral velocity variations, measurement errors, multipath effects, etc. (Yanovskaya and Kozhevnikov, 2003). The quality of the solution for the given regularization parameter is measured by comparing the initial and the remaining (unaccounted) mean square travel time residuals (Yanovskaya et al., 1998). Here, the initial means the travel time residual before the inversion based on the mean velocity of the data. The value of standard error (σ) was used for data selection. If an individual value of the travel time residual shows the value laying more than 3σ, the corresponding data were rejected from the data set and further recalculation of the solution is done. The observed values are almost comparable to the a priori errors resulting in fine velocity maps ( Table 2). The standard errors associated with the evaluated group velocities are ≤0.07 km/s. Ditmar and Yanovskaya (1987) and Yanovskaya et al. (1998) used two parameters, viz., mean size and stretching (ε) of the averaging area to prioritize lateral resolution. Yanovskaya (1997) adopted a technique similar to the one formulated by Backus and Gilbert (1968) for the 1D problem. This is formulated for 2D inversion through the calculation of the size of the averaging area S(x,y) for each point in different directions. Obtained resolution approximated by an averaging area over an ellipse centered at the point with S max (x,y) and S min (x,y) as the largest and smallest axes, respectively, over the area S (x,y) (Karagianni et al., 2002;Yanovskaya and Kozhevnikov, 2003). The mean size of the averaging area, L (S max + S min )/2, is the measure of the resolution for a particular point and these values for all grid points was used to calculate the lateral resolution of the area under investigation. The resolution depends on the path coverage between the epicenter and station as well as its azimuthal distribution (Gonzalez et al., 2007). The lateral resolution of the mean of the average area is less than 200 km in most of the regions under investigation ( Figure 3A and Supplementary Figure S4A). The lateral resolution for the present data set was evaluated for different periods between 4 and 60 s. Its variation for a set of four periods (8, 12, 30, and 60 s) is shown in Figure 3A. Plots for the remaining periods are shown in Supplementary Figure S4A. It is noticed that the resolution is good between the 8 s and 30 s period. The evaluated 2D resolution dimension mainly portrays the propagation path geometry and its density as shown in Figure 1B; Supplementary Figure S1.
The stretching parameter (ε) indicates how uniform the spatial distribution of the path is (Gonzalez et al., 2007). The stretching (ε) of the averaging area is obtained by 2(Smax (x,y)-Smin(x,y))/ (Smax (x,y)+Smin(x,y)), which describes the quality of the source-receiver path coverage (Kumar et al., 2019). Small values of the stretching imply that the paths are uniformly distributed in all directions, whereas large values indicate a preferred orientation of paths (Yanovskaya, 1997). Areas with stretching values ≤ 1.0 are considered to have a uniform azimuthal distribution of path and resolution is nearly the same along all directions. Figure 3B and Supplementary Figure S4B show stretching for the same periods for which resolution has been plotted. Tomographic images for group velocity and Vs are shown later for only the part where the resolution, as indicated in Figure 3A; Supplementary Figure  S4A, is good.

Tomographic Images
Rayleigh wave group velocity tomograms for 16 different periods were obtained at 4, 6,8,10,12,16,21,27,30,34,38,42,46,50,55, and 60 s to measure the spatial variation. Tomographic images are shown for the area where the spatial resolution is good (i.e., the mean size of the averaging area L ≤ 200 km). This encloses the area covered by the dense ray paths ( Figures 1B, 3A, and Supplementary Figure S1). Local dispersion curves obtained through tomography for six grid point locations ( Figure 1A) are shown in Figure 4 that indicates a high variation of group velocities at different periods from one place to another. Eight group velocity tomograms at 6, 10, 16, 24, 30, 38, 50, and 60 s are shown in Figure 5, and the rest are plotted in Supplementary Figure S5.
Rayleigh wave group velocity varies between ∼1.80 km/s and ∼4.15 km/s for periods ranging between 4 s and 60 s ( Figure 5; Supplementary Figure S5). The data set is enough to investigate the crustal and uppermost mantle structures of the study region. Lateral variations of group velocities represent different geological and tectonic features. A shorter period of the Rayleigh wave group velocity map is very sensitive to Vs in the uppermost crust, whereas a longer period has deeper depth sensitivity depending on its wavelength. Lateral variation of group velocity can be discussed based on various geotectonic features. According to Mitra et al. (2006), group velocities at <22 s periods represent upper crust, between 25 s and 35 s are mainly influenced by lower crust and up to 60 s represents upper mantle structure. Figure 4 shows that dispersion curves obtained for grid nodes ( Figure 1A) located in the BB, IBR, and BRV have lower group velocity values at lower periods compared to those located at SM, STD, and Lhasa block. This indicates the effect of sediments at shallower depths in BB, IBR, and BRV regions. The rate of increase of group velocity value at periods greater than 20 s is very high except for the Lhasa block, indicating the presence of a low-velocity medium in lower crustal depth in the Lhasa block.
At periods ≤10 s, the entire BB, Indo-Gangetic Plain (IGP), eastern Nepal Himalayas, IBR, and parts of southern Tibet  Figure S5) indicate the presence of sedimentary rocks. As group velocity at a higher period is affected by the Vs of the medium at greater depth, this shows that sedimentary rocks are present in BB and the southern part of the IBR up to greater depth. This observation for BB is supported by earlier investigations (Verma and Mukhopadhyay, 1977;Mitra et al., 2006;Acton et al., 2010), but for IBR, this is new information. Shifting of low group velocity from BB to IBR is supportive evidence of oblique subduction (Mitchell and McKerrow, 1975;Curray, 1989;Ni et al., 1989) of the Indian plate below the Burma plate. A higher velocity gradient in the western part of the BB relative to the eastern part of it shows that the basement depth increases from west to east. Low group velocity (from 4 s to 16 s) in BRV and IGP is also evidence of the presence of sediment at shallow depth (Kumar et al., 2018a). Kayal (2008) reported that the average sediment thickness in BRV is ∼4 km. No anomalous change is found in this zone and a continuous increase in group velocity is observed, indicating a gradual increase in medium velocity with depth. The Singhbhum craton and Rajmahal trap in the Indian shields in the southwestern part of the study region are represented by relatively higher group velocity up to 30 s. A prominent increase in group velocity from 21 s along the NE direction of the Indian continent to SMM and EHS is observed.
It may reflect upward crustal buckling at the middle of the study region as suggested by Raoof et al. (2017).
At all periods, SMM shows up as a zone of high group velocity indicating the presence of high-velocity crystalline rock bodies at all depth levels. A sudden decrease in group velocity is observed across the Dauki fault from north to south up to 21 s period. At periods >30 s, southern Tibet has lower group velocity compared to eastern Himalayas, NE India, and BB, indicating that at a deeper depth, the former region comprises lower velocity material. This corroborates the presence of thick crust in Tibet. This is also in accordance with the suggestion of the presence of partial melts and/or aqueous fluids in the midcrust in Tibet Cogan et al., 1998;Hetényi et al., 2011;Jiang et al., 2014). Such low-velocity zones in southern Tibet are observed in between 27 s and 50 s. Acton et al. (2010) also observed this type of change of velocity pattern between 50 s and 70 s and inferred that the crust below the Tibetan plateau is thick. Variation of group velocity all along the Himalayas, especially at lower periods, indicates the possible variation of crustal thickness along this tectonic trend as suggested by Raoof et al. (2017). In Nepal Himalayas, the low-velocity at short period has also been observed by Guo et al. (2009). The Gangetic plains adjacent to this region also have relatively lower group velocity with high gradients up to 16 s, which is similar to observations of Mitra et al. (2006).

INVERSION FOR CRUSTAL STRUCTURE
Tomographic maps provide averaged and smoothened group velocities estimated at discrete points at 1°× 1°spacing within the area under investigation. These are described as the local dispersion curves of different grid points (ex. Figure 4). The smooth group velocity versus period data were inverted to obtain an average 1D Vs model at each grid node using CPS3.30 (Herrmann and Ammon, 2002;Herrmann, 2013). The lateral variation at each depth level from 4 to 90 km was obtained by contouring these values as per the procedure described by Kumar et al. (2019). In this process, the dispersion curve extracted from the group velocity maps at each grid point was inverted iteratively using the damped least-squares method. Selected grid points to the part of the study area where the resolution is at least 200 km (see Figure 3, Supplementary Figure S4) were inverted. The available path densities (Table 2, Supplementary Figure S1) show that the resolution for 8 s-30 s periods covers a very wide area (Figure 3, Supplementary Figure S4). For each data set, a fixed number of iterations were carried out. To avoid the trade-off between model variance and data variance, the selection of the damping factor is an important parameter to obtain the unbiased shear wave velocity structure. The damping factor is evaluated from the L curve of the trade-off between data variance and model variance for different damping values between 0.01 and 10 for single iteration inversion (Kumar et al., 2017) Figure S6). As an example, the optimum damping value obtained for the dispersion curve of one grid point is 1.1. The selected damping factor gives the smallest value which provides a balanced fit between the model and data from the L-curve FIGURE 4 | Local dispersion curves for six different grid points (locations shown in Figure 1A).

(Supplementary
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 680361 (Gubbins, 2004). This procedure is adopted for obtaining the damping factor of each grid point. To start inversion, an initial assumed velocity model is required (Yanovskaya et al., 1998;Bhattacharya et al., 2013;Motazedian and Ma, 2014). The inversion is mainly performed in two steps, in the first step, 30 different initial models are selected from the study region and adjoining parts of the Himalayas and IBR (Kumar et al., 2019 and references within). These initial models are used to invert each local dispersion curve for each model to obtain a priori Vs model. To overcome the dependency of the inverted model on the initial model, we constructed the average 1D Vs model at each grid node for each tomographic dispersion curve at the respective grid in this manner. In the next step, 250 randomly positive and negative anomalies were inserted in constructed average 1D Vs model for a wide range of perturbation at different depths and velocities. Thereafter, the inversion is performed for all these generated models at each grid nodes to obtain the respective output Vs models.
This brings the inverted model close to the local structure where perturbations minimize the dependency of inversion on the initial model. To follow the dispersion curve, a gradual increase in layer thickness and velocity is considered with depth. This process was carried out for each grid node within the study area separately by using the respective average dispersion curve of that grid. An example of initial and inverted models for the grid node at 24°N, 90°E are shown in Figures 6A,B, respectively. The best velocity model obtained from the weighted average of 250 output models is represented by a solid blue line ( Figure 6B), which is based on the weighted average values for different depths from 250 output models. Our result indicates that the final model converges very well from a wide range of initial models and, therefore, the output is independent of initial values.

SHEAR WAVE VELOCITY STRUCTURE
Inverted 1D Vs models at all nodes are used to obtain lateral 2D variation (Figure 7, Supplementary Figure S7), and combining these 2D models, 3D velocity structure is obtained. 2D Vs map indicates ∼2.8-3.0 km/s average velocity at 4 km depth in the southern Tibet and Lhasa block. It increases gradually to ∼3.3 km/s at 6 km depth, to ∼3.4-3.6 km/s at 10 km, and ∼3.6-3.8 km/s at 24 km depth. The average velocity values at different depths match well with the 1D velocity results of Monsalve et al. (2008). A low-velocity zone is observed at the mid-crustal level (∼27-50 km) for southern Tibet, which is similar to the observations of Jiang et al. (2011), Sun et al. (2010, and Kumar et al. (2018b;. The ultra-low velocity (∼3.3 km/s) zones in the mid-crust of southern Tibet have been suggested to be due to partial melting or aqueous fluid (Cogan et al., 1998;Makovsky & Klemperer, 1999;Unsworth et al., 2005;Jiang et al., 2014). Low-velocity derived in the present study is supported by previous geophysical studies, for example, observation of low-velocity anomaly (Kind et al., 1996;Hung et al., 2010;Hetényi et al., 2011;Kumar et al., 2019), bright spot Makovsky et al., 1996), high electric conductivity (Chen et al., 1996;Wei et al., 2001), magnetotelluric study (Unsworth et al., 2004;Arora et al., 2007), and high heat flow (Francheteau et al., 1984). Part of the Nepal Himalayas (between 85°E-∼88°E) is located within the study area. Top ∼8 km here contains low-velocity material (<2.8 km/s) that has a velocity lower than that in NE Himalayas at the same depth level (Figure 7, Supplementary  Figure S7). Down to ∼10 km depth, there is a variation in Vs along the Himalayas from Nepal to NE Himalayas. A positive velocity gradient is observed from 4 to 18 km depth (ranging from ∼2.4 to 2.6 km/s to ∼3.6-3.8 km/s) in the Nepal Himalayas. In NE Himalayas, this value varies from 3.0 to 3.6 km/s in the same depth range. Velocity gradient up to 18 km is relatively higher in the Nepal Himalayas as compared to that in the NE Himalayas and south Tibet (see also Raoof et al., 2017).
In the upper part of the crust down to ∼21 km depth, Vs in BB is low compared to other parts of the study area (Figure 7, Supplementary Figure S7). This contrast is very high for the uppermost crust and distinctly visible between ∼6 and 8 km depth indicating thick sediment depositions in the basin. The margin of low-velocity sediment deposition is demarcated by Dauki fault in the north. Although Vs increases gradually with depth, roots of low-velocity are visible down to ∼21 km depth. Curray et al. (1982) and Brune and Singh (1986) also suggested that sediments/meta-sediments are present down to ∼22 km depth in this region. Variation of sedimentary layer thickness in the BB observed by us is also seen by Singh et al. (2016). We find Vs representative of the crustal basement at depth below ∼24 km. At depth ≥24 km, the Vs below BB is similar to that of oceanic crust (Sowers and Boyd, 2019) The basement also dips in the SE direction.
There are numerous active thrust faults in IBR (Le Dain et al., 1984). Low-velocity up to a depth of ∼10 km is observed below the entire IBR. This could be because of several reasons, viz., presence of cracks, fluid, partial melting, or underthrusting of sediments. However, as Vs in this zone is similar to that observed in the adjacent BB, underthrusting seems most plausible. The low-velocity zone shifts southward along IBR. It is observed down to ∼15-∼18 km depth within the southern part of IBR. Below this depth, the velocity values are representative of crustal material. Therefore, the thickness of sediments below IBR increases from north to south. The high-velocity medium below 50 km depth may represent the subducting lithosphere as proposed by Raoof et al. (2017).
The BB, SMM, and BRV are trapped between two tectonic arcs. There is a very sharp change in velocity across the Dauki fault that separates BB from SMM with a very clear signature from ∼4 to ∼18 km depth (Figure 7, Supplementary Figure S7). This shows that the Dauki fault separates the deep-seated sediments of BB from the crystalline rocks of SMM. Low Vs in the upper ∼4 km in IGP and BRV can be associated with thick sediment deposits (Figure 7, Supplementary Figure S7) (see also Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 680361 Mitra et al., 2006;Suresh et al., 2008;Kumar et al., 2017). It is observed that sedimentary layer thickness varies in the IGP and especially in BRV. Within IGP and NE corner of BRV sedimentary rocks are observed up to 6 km depth, whereas within BRV north of SMM, its thickness is less than 4 km (Figure 7, Supplementary Figure S7).
To get a better perspective about the medium characteristics, S-wave velocity variation along eight profiles ( Figure 1A) are presented (Figure 8). Profile AA ' (Figures 1A, 8A) is aligned SSW-NNE which crosses the IGP, Nepal Himalayas, and Tibetan plateau. Prominent tectonic features along this section are HFT, MBT, MCT, STD, and ITSZ. The surface wave is not sensitive to sharp discontinuities like Moho but gives absolute velocity for a depth section across the discontinuity. Kumar et al. (2019) correlated the Rayleigh wave tomography data with the receiver function of the Himalaya-Karakoram-Tibet region to evaluate 4 km/s S-wave contours of tomography results as the approximation of Moho discontinuity. Acton et al. (2010) obtained a 4.1 km/s contour of Moho based on the Rayleigh wave tomography data of the Indian continent and Tibet. Based on the receiver function analysis, Bora et al. (2014) obtained 4 km/s S-wave velocity for Moho in a large part of NE India. In the present study, Moho is marked at 4 km/s (Figure 8). Figure 8A shows that crustal thickness increases gradually from south to north which clearly shows the effect of underthrusting of Indian plate below the Himalayas. In the south beneath IGP, the Moho depth is ∼30 km, which increases to ∼80 km beneath Tibet. In the Nepal Himalayas, interfingering of high-and low-velocity material indicates repeated underthrusting and overthrusting of crustal and mantle material, respectively. Based on the receiver function analysis, Acton et al. (2011) observed imbrication of the crust FIGURE 7 | Lateral variation of Vs at different depths. Abbreviations of tectonic features are as in Figure 1A.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 680361 in the same region. This further supports our observation. A very low-velocity medium north of ITSZ near the Kung-Co rift at a depth range of 30-50 km probably represents a zone of partial melting and/or the presence of aqueous fluid (Makovsky and Klemperer, 1999;Unsworth et al., 2005;Klemperer, 2006). Our observation along this profile is showing good agreement with the finding of the presence of ramp-flat-ramp geometry of the Main Himalayan Thrust (MHT) in the Sikkim Himalaya using joint inversion of receiver functions and surface wave dispersions by Paul and Mitra (2017). Profile BB ' (Figures 1A, 8B) starts from the BB; passes through the SM, NEH, and Lhasa block; and reaches the Bangong-Nujiang Suture (BNS) of south Tibet. The topography is highly variable (minimum and close to mean sea level at BB and BRV). The Moho depth is also highly variable from ∼33 km beneath BB, gradually increasing toward the north up to the frontal Himalayas. Around SM, the crustal thickness is ∼40 km. In the Himalayas from MBT to MCT, the Moho depth suddenly increases and mirrors the increase in topography height. In the Bhutan Himalayas, the Moho depth obtained in this study is similar to that obtained by Singer et al. (2017). Further north, it gradually increases to ∼75 km depth. A very striking result is the low-velocity in the uppermost crust down to ∼20-22 km in the BB for more than 200 km in the southern part. This suggests the presence of thick sediment deposits. The uppermost crustal velocity of SM is higher than that of BB and BRV. The boundary between the sedimentary rock and upper crust (dashed line) clearly shows that the Dauki fault  Figure 1A.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 680361 juxtaposes sediments of BB and crystalline rocks of SM. This also suggests that the Massif is either a pop-up structure (Bilham and England, 2001) or the crest of buckled-up Indian plate (Raoof et al., 2017). Profile CC' (Figures 1A, 8C) passes through the IBR, SM, MH, NEH, and southern Tibet and shows abnormal Moho depth in the central part of the profile. In the southernmost part, this discontinuity is marked at ∼40 km depth with a gradual increase toward the north for about 300 km passing through IBR. South of Dauki fault, a sudden bend is detected in the Moho discontinuity with a sudden uplift in a localized zone of 200 km beneath MH. This observation is similar to Profile BB', where the sudden uplift is observed in the boundary between the sedimentary rock layer and the crystalline crust and supports the hypothesis that MH, like SM, is either a pop-up structure or it sits atop the crest of the buckled-up Indian plate. Moho depth is observed between ∼30 and 32 km beneath the MH. Borah et al. (2016) obtained crustal thickness of 32-36 km north of SM based on the S-wave velocity model. The occurrence of thinned crust (30-32 km) beneath SM is also reported by Mitra et al. (2018) and Priestley et al. (2019). The sedimentary layer thickness progressively decreases from IBR to NEH. The low-velocity zone is seen in the mid-to-lower crust beneath Tibet and close to ITSZ that is also visible in Profile HH' (Figure 8H).
Profile DD' (Figures 1A, 8D) starts from the west of BB and continues up to the Sagaing fault. An ∼10-km thick layer of sedimentary rocks observed in the western part of BB up to EHZ increases down to ∼20-21 km in the east beneath BB and IBR. Underthrusting of sedimentary wedge below the IBR is observed in both east and west. Moho depth varies between ∼30 and 40 km.
Profile EE' (Figures 1A, 8E) is along WNW-ESE from the IGP up to IBR. Within BB, low-velocity in the uppermost part is due to the thick sediments. In the west of the IBR, an ∼20-21 km thick pile of sediment is visible. Moho is undulating and its depth is ∼30 km at the WNW end and ∼50 km at the ESE end.
Profile FF' (Figures 1A, 8F), from the southwest of BB, crosses the BB, passes through SM, a long stretch of BRV, and then perpendicularly traverses the syntaxial bend of EHS. It is the longest section that stretches over 1,450 km and crosses many geological features and tectonic discontinuities. In the BB, it aligns nearly along the strike of EHZ. The profile crosses the Dauki fault to the south of SM and DKF to its northeast. Along the EHS, the section passes through the Mishmi Hills and several NW-SE trending faults. It shows low-velocity in the uppermost crust at both the sides of the SM related to thick sediment deposits of the BB and BRV. In the SM, the sedimentary rock layer thickness is very low. Crustal thickness varies between 35 and 70 km with the thinnest crust beneath the Massif and below the BRV close to the Massif and thickest in the northeast end beneath the EHS. Borah et al. (2016) obtained a crustal thickness of 36-40 km in the Assam valley using the S-wave velocity model. In BB, its value is between 40 and 45 km and decreases toward the northeast. Its thickest geometry in the EHS indicates subduction of the Indian plate. The increase in crustal thickness on either side of the Massif supports the observations of Raoof et al. (2017) suggesting upward crustal buckling beneath the Massif.
Profile GG' (Figures 1A, 8G) is oblique to the Himalayan arc starting from the central Himalayas, crossing a narrow zone of the BRV, SMM, and further east perpendicular to IBR, then it crosses into the Burma plate. In the uppermost crust within the Himalayas and BRV, low-velocity is related to metasedimentary and sedimentary deposits of the Tethys and BRV, respectively. Borah et al. (2016) observed a 44 km thick crust beneath the lesser Himalayan region. Our result shows that the crustal thickness varies between 45 and 60 km with maximum thickness below the Himalayas.
Profile HH' (Figures 1A, 8H) covers the Lhasa block and EHS. It starts from south Tibet near the Kung-Co rift, passes through the Yadung-Gulu rift, ITSZ, Mishmi Hills, and ends within the EHS. The topography in southern Tibet is high and the Moho depth is ∼70-75 km for the initial 500 km of the profile. Crustal thickness suddenly decreases near the Yadong-Gulu rift and after that, it gradually decreases to ∼40 km near the EHS. Beneath the Tibetan plateau/Lhasa block, two broad low-velocity zones are detected within ∼30-40 km for a distance of about 300 km. Other workers Kind et al., 1996;Nelson et al., 1996;Cogan et al., 1998;Hetényi et al., 2011;Jiang et al., 2014) also observed low-velocity medium in this region. Priestley et al. (2008) and Mitra et al. (2005) using the receiver function estimated the Moho depth variation along an NS profile starting from 22.60°N and 91.25°E which is close to the profile BB' of this study ( Figure 8B). In Figure 8I, the Moho depth reported by them is superposed on profile BB'. Our estimated Moho depth is shallower in most parts of the profile than their estimates. Although the Moho position estimated by them does not exactly match that of ours, the general trend of Moho depth variation is similar. Beneath the Higher Himalayas, the Moho depth suddenly dips toward north; similar observations are reported by Singer et al. (2017) and Mitra et al. (2005).

Bay of Bengal and Indo-Burma Range
The Bay of Bengal is characterized by the thickest sediment deposition in the whole world. Huge sediment deposition occurred in BB during the early-mid Miocene following collision between India and Eurasia (Alam et al., 2003). The large and quick influx of clastic sediments within the basin from the Himalayas in the north and IBR in the east led to sudden subsidence of the basin. At this stage, deep-marine sedimentation at the deeper part of the basin and deep to shallow marine sedimentation in the eastern part of the basin occurred (Alam et al., 2003). Mitra et al. (2018) observed the presence of 18-20 km thick sediment underlain by crystalline crust beneath the central BB. Our results show that sedimentary layer thickness is ∼10 km in the western part of BB and it starts increasing rapidly east of EHZ with a maximum value of ∼21 km at the eastern end of BB (Figures 7, 8, Supplementary Figure  S7). The margin of low-velocity sediment deposition is demarcated by the Dauki fault in the north and IBR in the east. Although Vs increases gradually with depth, low-velocity strata are visible down to ∼21 km depth. Thick sediment deposition in the BB camouflages basement configuration and makes it difficult to find the exact location of the boundary between oceanic and continental crusts (Alam et al., 2003). The oceanic crust contains only the mafic rock layer, whereas the continental crust has both the felsic layer overlying the mafic layer. Mafic rocks usually have higher P/S wave velocity than felsic rocks. Also, oceanic crust is thinner than the continental crust. It is observed that the crust underlying the thick sedimentary cover of BB has a thickness varying between ∼10 and 18 km and Vs ∼3.8 km/s (Figures 7, 8B,C, Supplementary Figure S7), which is similar to that of oceanic crust (Gonzalez et al., 2007;Pasyanos and Walter, 2002). North of the Dauki fault, a crustal layer having an S wave velocity similar to that of continental crust (Vs ∼3.4 km/s) overlying this layer is also observed and the total thickness of these two layers increases northward (Data and Estimation of Group Velocity in Figure 7). Hence, north of the Dauki fault crustal structure is similar to the continental crust. Based on these observations, we propose that BB is underlain by oceanic crust, and north of the Dauki fault, the crust is continental (Figures 7, 8B,C, Supplementary Figure S7). Hence, we conclude that the Dauki fault represents an area where the transition of oceanic crust toward its south to continental crust toward its north occurs. Thus, we find clear demarcation between the oceanic and continental crusts. Based on paleomagnetic, gravity, and deep seismic sounding data, Talwani et al. (2016) too found the zone of transition from oceanic to continental crust near this area. Brune and Singh (1986) and Kumar et al. (2018a) found that average Vs in BB changes from oceanic type in its southern part to more continental type toward the north. Based on previous works using Rayleigh wave close to this study region (Acton et al., 2010;Bora et al., 2014;Kumar et al., 2019), 4.0 km/s contour is taken as the Moho discontinuity.
Numerous active thrust faults have been reported in Burma (Le Dain et al., 1984). Low-velocity up to a depth of ∼10 km is observed below the entire IBR. With increasing depth, this lowvelocity zone shifts southward along IBR with a gradual increase in velocity. It is observed up to ∼15-∼18 km depth within a very small portion representing the thickness of sediments as they match with velocity observed in the adjacent part of BB ( Figures  7, 8, Supplementary Figure S7). Below this depth, the velocity values are representative of crustal material. This shows that the depth up to which sedimentary rocks are present below the IBR increases from north to south. High-velocity medium below 50 km depth in this region may represent the subducting lithosphere (Figures 7, 8, Supplementary Figure S7) as reported by Raoof et al. (2017). Based on local earthquake tomography over a part of the IBR and Burma plate, Zhang et al. (2021) too find evidence of eastward subduction of the India plate below the Burma plate.

Shillong-Mikir Massif, Indo-Gangetic Plain, and Brahmaputra River Valley
Low Vs in the upper ∼4 km in the IGP and BRV is associated with thick sediment deposits. Evaluated Moho depth, based on 4 km/s velocity contour, varies between ∼30 and 40 km (Figures 7, 8).
Thick sediments observed by us and highly oblique plate motion suggested by Kayal (2008) between India and IBR suggest the presence of active shallow dipping and locked megathrust faults . We find a high Vs structure between ∼50 and ∼75 km (Supplementary Figure S7) within the BB formed due to crustal-scale buckling. At deeper depth, a high-velocity gradient shifted eastward and toward the IBR (Figure 7).
The SMM and BRV are sandwiched between two tectonic arcs where large earthquakes occur. Vs in SMM is representative of high-velocity crystalline rock bodies near the surface, whereas an ∼6 km thick sediment is present in the BRV (Figure 7). A sharp change in velocity across the Dauki fault separating the BB from the SMM is observed, especially between ∼6 and ∼18 km depth. This shows that the Dauki fault separates the deep-seated sediments of BB from the crystalline rocks of SMM, where sediment thickness is negligible. High Vs beneath SMM reveals uplifted crust/uppermost mantle caused by it being in a vice-like grip between the Himalayas and IBR ( Figures 1A,  8B,C,F).
It has been suggested that upliftment of SM is caused by N-S compression due to Himalayan collision and E-W directed compression due to the Indo-Burman subduction (Chen and Molnar, 1990;Mukhopadhyay et al., 1997;Rao and Kumar, 1997;Najman et al., 2016;Govin et al., 2018). Upliftment of the SM is linked with kinematic changes in response to a collision between the Himalayas and the Indian plate and subduction of the Indian plate below the Burma plate cumulatively (Clark and Bilham, 2008). The crustal structure of NE India using the P-wave receiver function (Mitra et al., 2018) suggested that the SM uplifted due to thrust faulting on the Dauki fault, continent margin paleo rift fault, and back-thrusting Oldham fault. Strong et al. (2019) stated that the Massif consists of an actively uplifting block with associated continuous erosion processes.

Nepal Himalayas, Northeastern Himalayas, and South Tibet
A variation in Vs along the Himalayas from Nepal to NEH is noted down to a depth of ∼10 km (Figures 7, 8). Velocity values at each of 12-24 km, 33-45 km, and 65-90 km depth range all along this zone are uniform. At other depth ranges, variation in velocity is observed along the trend of the Himalayas. This may be the effect of complex structural geometry in this region. Velocity gradient down to 18 km is relatively higher in the Nepal Himalayas than that in the NEH and south Tibet. Koulakov et al. (2015) estimated that the Moho depth varies both along and across the tectonic trend of the Nepal Himalayas and suggested that the first is due to the crumpling of the crust caused by the stress along the tectonic trend due to counterclockwise rotation of the Indian plate after the first contact with the Asian plate in the NW Himalayas. Observed lateral velocity variation along the tectonic trend in the NE and NW Himalayas (Koulakov et al., 2015;Raoof et al., 2017;Raoof et al., 2018) and our results in this study also support this contention.
Variation in Vs (Figures 7, 8) in southern Tibet and the Lhasa block matches well with the 1D velocity model estimated by Monsalve et al. (2008) for this region. In the mid-crust (∼27-50 km), ultralow-velocity (∼3.3 km/s) zone is observed which is also supported by other researchers using different geophysical studies (Francheteau et al., 1984;Brown et al., 1996;Chen et al., 1996;Kind et al., 1996;Makovsky et al., 1996;Nelson et al., 1996;Wei et al., 2001;Unsworth et al., 2004;Arora et al., 2007;Hung et al., 2010;Jiang et al., 2011). This may indicate partial melting or the presence of aqueous fluid. It is observed that the zones having partial melt/aqueous fluids do not form a continuous layer in southern Tibet but appear rather localized (Figures 7, 8). Hetényi et al. (2011) too observed that zones having partial melts/aqueous fluids are discontinuous in nature.

Eastern Himalayan Syntaxis
EHS is the frontal part of the sandwiched Indian plate between the Eurasian and Burma plate. Sediment thickness of ∼6 km is observed toward the northeast of BRV. An abrupt change in crustal thickness is observed in EHS between the F′ and H′ edges ( Figures 1A, 8F,H) of profiles FF′ and HH′, respectively. These two points are on the two sides of EHS, but they are quite close to each other. However, our study shows that the crust is 70 and 45 km thick near the F' and H' edges, respectively. Significant variation in crustal thickness within the northeast part of the EHS suggests that the nature of crustal deformation in the EHS is complex. Crustal thickness across the Tidding Suture is 55 km (Hazarika et al., 2012). An increase in crustal thickness within the EHS can be interpreted as an effect of indentation geometry of the Indian plate from the southwest direction. Convergence of Indian and Eurasian plates generates N-S and E-W compression within the NE India, leading to crustal buckling across the NE axis as shown in this work. CONCLUSION a) Bengal Basin is underlain by oceanic crust. Further north below the Shillong Massif/Brahmaputra River Valley, the crust changes to continental type. b) Thickness of BB sediments increases from west to east. This inference is similar to previous estimates of sediment thickness for this area. The dip of the basement of this sedimentary basin increases suddenly to the east of the Eastern Hinge Zone. c) Sediments are also getting underthrusted at the southern part of the Indo-Burma Ranges that lie within the study area. d) In the Indo-Gangetic Plains and Brahmaputra River Valley, 5-6 km thick sedimentary layers are present. This observation is also reported by previous workers. However, north of SMM sedimentary layer is very thin. e) Across the Dauki fault, low-velocity sediments of Bengal Basin are juxtaposed with high-velocity granitic/metamorphic rocks of the Shillong Massif and Mikir Hills. In the Shillong Massif and Mikir Hills, sedimentary layer cover is very thin or absent and crystalline rocks having higher velocity are present near the surface. f) Shillong Massif and Mikir Hills lie atop the Indian plate that has buckled up because in this area, the Indian plate is in a vice-like grip between tectonically active Himalayan ranges toward the north, below which it is underthrusting, and the Indo-Burma Ranges toward the east, below which it is subducting (Grujic et al., 2018). g) The Moho depth varies not only from south to north but also along the tectonic trends of the Himalayas. h) Average Moho depth increases from ∼40 km in the south to ∼70 km in the northern part of the study area smoothly in most places but abruptly at other places. i) In the eastern part of the Nepal Himalayas, the signature of repeated underthrusting and overthrusting is observed, whereas such a pattern is not observed further east along the trend of the Himalayas. j) For southern Tibet/Lhasa block, two important findings are as follows 1) the Moho depth decreases from ∼70 km west of the Yadong-Gulu rift to ∼60 km toward the east which may be the effect of collision geometry to control the underplating of Indian plate in this section of southern Tibet/Lhasa block (Shi et al., 2015); 2) few low-velocity zones are observed in velocity sections at mid-crustal levels near or below the rift zones in the southern Tibet/Lhasa block. This could be due to the presence of zones having partially molten material and/or aqueous fluids. This observation is similar to many seismic and magnetotelluric studies in this area.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.