Journal of Geography  & Natural Disasters
Open Access

Seismological Evidence for Intra-Crustal Low Velocity and Thick Mantle Transition Zones in the North-West Himalaya

Chinmay Haldar^*, Narendra Kumar and Kalachand Sain ^*Correspondence: Chinmay Haldar, Department of Geography, Wadia Institute of Himalayan Geology, Dehradun, India, Email:

Author info »

Introduction

Since 55 Ma years, a significant topographic variation has occurred due to the continent-continent collision and convergence between the Indian and Eurasian plates in the Himalayan mountain belt. The main topographic variations in the Himalayan region occurred due to the crustal shortening, development of ample folds and thrusts, crustal deformation, magmatic activity, and uplift of the Tibetan Plateau. Due to the presence of highly dense material in the Indian plate, it under thrust below the Eurasian plate along a northerly dipping detachment surface called the Main Himalayan Thrust (MHT), and this detachment surface separates the down going Indian plate from the overriding Himalayan wedge. The collision of Indian and Eurasian plates accumulates a lot of strain energy along the MHT zone. The release of this strain energy causes earthquakes of varying magnitudes in the Himalayan region, such as the 1905 Kangra Earthquake (Mw 7.8) the 1934 Bihar- Nepal earthquake (Mw 8.4), the 1991 Uttarkashi Earthquake (Mw 6.8), the 1999 Chamoli Earthquake (Mw 6.3), the 2001 Kunlun Earthquake (Mw 8.1), the 2005 Kashmir Earthquake (Mw 7.6), the 2008 Wenchuan Earthquake (Mw 8.0) and recent the 2015 Gorkha-Nepal earthquake (Mw 7.9). Many geophysical studies reveal the presence of a mid-crustal ramp in the NW Himalaya. Most of the earthquakes beneath the Kumau-Garhwal Himalaya occurred due to this ramp geometry of the MHT beneath the Munsyari. Several geophysical and geological studies proposed that the locked zone along the southern part of the MHT can generate great earthquakes (>Mw 8). Bilham suggested that few segments of the Himalayas are mature enough to produce a great earthquake (Mw>8.0). Recently, Pappachen proposed that the Kumaun-Garhwal, situated in the northwest Himalayan region, falls under a high seismic risk zone and can produce a great earthquake[1].

The intra-crustal Low Velocity Layer (LVL) is crucial to understanding any region's seismic activities and geodynamic evolution. Previous researchers based on the LVL study came up with a vital indication to understand lithospheric deformation, nature of MHT, and geodynamic development of the himalaya. But the hypothesis related to the genesis of LVL is under debate among geoscientists [2]. Some studies have been carried out to understand the potential causes for the existence of the LVL in a few parts of himalaya the eastern and central parts of the himalaya found a fluid-filled LVL beneath the Southern Tibet detachment using the converted wave technique. Hazarika, et al. and Kanna and Gupta also showed evidence of LVL beneath the sub and higher himalaya. Magneto-telluric studies support the presence of LVL beneath the Himalaya, accounting for the low resistivity zone at the mid-crustal layer. Magneto-telluric studies support the presence of LVL beneath the Himalayas, accounting for the low resistivity zone at the mid-crustal layer [3]. This low resistivity zone in the upper crust may be due to the partial melt, aqueous fluid, or graphite in this region. But, the structure of the intra-crustal LVL, crust-mantle boundary (Moho), and upper mantle discontinuity remain unexplored [4].

Using the converted wave technique, the current research will provide new inferences of intra-crustal LVL, Moho geometry, and upper mantle discontinuities beneath the kumaun-garhwal region. Due to the complex structure of the Himalayas and to avoid the overlapping effects of conversion points, a Common Conversion Point (CCP) stacking technique is used to understand the upper mantle geometry better. The current information from past studies is incorporated to shed light on this region's geological and geophysical behavior. Based on the inversions results, the present research shows a good correlation between the known geological and tectonic features with the upper crustal LVL below the North-West (NW) Himalaya [5].

Materials and Methods

Receiver function

Data recorded at seven broadband seismological stations beneath the Kumaun-Garhwal, NW Himalaya from 2016 to December 2019, operated by the Wadia Institute of Himalayan Geology (WIHG), Dehradun, has been used here. The seismological units consist of a nano metrics trillium 240 broadband sensor with a flat velocity response between 0.004 Hz and 35 Hz and a sensitivity of 200 V/m/s connected with a centaur digitizer with a sample rate of 100 samples/second. The Global Positioning System (GPS) receivers synchronize the data [6]. The Receiver Function (RF) method is a converted wave technique widely used to investigate the sub-surface structure below the receiver. Langston proposed the RF technique based on the spectral ratio method to understand the crustal structure by isolating conversions from different discontinuities [7]. On the other hand, Vinnik proposed a similar approach to explore the upper mantle discontinuities using delay and sum techniques. We used teleseismic data to understand the study region's crust and upper mantle structure. We took data recorded from all the back azimuths and sorted these events with their magnitude greater than ≥ 5.5 Mb and epicentral distance in 300-900 [8]. We have visually selected the waveform data for the quality purpose, which has a perfect P phase in the vertical (Z) component. The recorded three-component waveforms, i.e., the Z, N (North- South), and E (East-West), have been rotated using the backazimuth to get a new coordinate system, i.e., into Z (Vertical), R (Radial), and T (Transverse). We have rotated these Z-R-T components into the L (P), Q (SV), and T (SH) by using the theoretical incidence angle estimated from the maximum polarization direction of the P-wave. We have chosen only those events whose differences in theoretical and waveform incidence angle are less than 50 and 300 in the case of back azimuth. We did not consider the effects from the free surface to estimate the angle of incidence [9]. We have applied a low pass filter with a corner frequency of ̴1 Hz to avoid the scattering from shallow depths below the stations. After rotation, we have estimated the time domain P-Receiver Functions (PRFs) using the deconvolution of the radial components with their respective vertical ones [10]. This deconvolution technique has eliminated the source and propagation path effects. Sometimes, the crosscorrelation method can be used in place of the deconvolution technique. However, Kumar et al., retrieved the structure using a summation technique without the deconvolution techniques. We have used the spiking deconvolution method, a least-squares inverse filter and estimated this filter by choosing a time window of 10 s before the theoretical onset of P and 90 s after. If we choose a more extended time window, the waveform may contain multiple reflections instead of the source function [11]. We applied the moveout correction to the PRFs using an epicentral distance of 67˚ corresponding to a ray parameter of 6.4 s/deg or 0.057 s/km to eliminate the source-receiver distance effect. In the actual field seismological data, the Ps converted phases generally have a low amplitude (compared to the synthetic ones). We have stacked all the derived individual moveout corrected PRFs to increase the signal-to-noise ratio or amplitude (Figure 1) [12].

Figure 1: a: Location of broadband seismic stations (demarcated by red triangles) superimposed over the topographic map of NW Himalaya. Inset shows the study region within the Indian subcontinent. b: The distribution of earthquake epicenters (red circle) of all events used in the present study.

Neighbourhood Algorithm (NA) inversion method

The Neighborhood Algorithm (NA) inversion scheme is a stochastic inversion method. NA method mainly depends on the combination of parameter space [13]. Inversion of PRFs using the NA scheme is a complex nonlinear waveform matching method that overcomes the non-linearity and non-uniqueness problems where the linearized inversion method either fails or over-reliance on the starting model. Twenty-four variables parametrize the shear wave velocity structure of the crust [14]. This method provides the velocity structure of six horizontal layers, i.e., sedimentary layer, basement, upper-crust, middlecrust, lower-crust, and upper mantle beneath each station. Each layer is characterized by four parameters in each layer, i.e., the thickness of the layer, shear wave velocity at the top and bottom of the layer, and the ratio of P to S wave velocity at each layer [15]. The shear wave velocity profile is computed using the direct search method from a multi-dimensional range of acceptable wide velocity gradients between the two velocities at each layer's top and bottom parts. NA inversion scheme depends on the Voronoi tessellation concept, where the division of the initial model or resample has taken place in terms of Voronoi cells, i.e., N_r (cells resampled number for every iteration) and Ns (models generation number).

The synthetic RF is computed using the Thomson Haskel matrix scheme where the observed misfit is constant within each Voronoi cell. In our study total of 200 iterations have been carried out where the model generation number is 100 (N_s) and the number of cells resampled is 2 (N_r) [16]. This method produces more than 20,000 shear wave velocity models that match the observed and synthetic receiver functions with every inversion. From these 20,000 models, we got 1000 best-fitting models, with the final velocity model having the least misfit. The stack of PRFs from all the stations are inverted using the NA inversion scheme to estimate the crustal shear wave velocity structure of the Kumaun Garhwal Himalaya. Since the study region lies in the Himalayas, the initial model for crustal velocity has been taken from the previously published literature as a constraint to the inversion scheme [17]. This type of practice saves the analytical time to get the best model by fitting the synthetic receiver functions with the observations. The misfit error between the observed and a set of acceptable synthetic receiver functions is relatively low for each inversion. The results estimated from the stacked PRFs and NA inversion is summarized in Table 1 [18].

Table 1: Name of all seven (7) seismological stations along with depth of Low Velocity Layer (LVL), Crust-Mantle Boundary (M) and V_p/V_s. The conversion times from P-to-s from the upper mantle discontinuities have shown. LH-Lesser Himalaya, SH-Sub Himalaya and HH-Higher Himalaya, geo-geological.

Results and Discussion

Intra-crustal low-velocity layer and crust-mantle boundary

In Figure 2 we have shown the individual and stack PRFs at the PITG, MUNS, THRI, GRGA, BHTW, RWSR, and HDWR stations, where we have arranged the stations with increasing order of distance from Himalayan Frontal Thrust (HFT). The number of receiver functions that have been stacked is 47, 13, 25, 14, 41, 24, and 31 for station HDWR, THRI, RWSR, PITG, BHTW, GRGA, and MUNS. From the stack PFRs, we can see two clear phases, i.e., one negative (grey color) and one positive (black color). In the present study, the first positive (black) polarity observed just after 0 sec below the four stations, namely HDWR, BHTW RWSR, and GRGA, indicates the presence of a sedimentary layer (Figure 2). The sedimentary effect has also been included in the NA inversion modeling for mapping the Moho [19]. Previous studies also demonstrated the effect of sedimentary deposition in the Himalayan region. However, the objective of the present research mainly aims to estimate the structure of the crust and upper mantle in the study region. The positive conversions occurred from the crust-mantle boundary, and the negative conversion happened at the intra-crustal LVL zone. Both these conversions are observed from the individual and stacked receiver function traces. Clear, and coherent negative P to S conversion from LVL is observed, whose conversion time varies from 1.4 s to 3.0 s. In contrast, the positive conversion time from the crust-mantle boundary varies from 4.8 s to 6.8 s. To quantify the layered structure beneath a station in terms of the P and S velocity variations, we have used the neighborhood algorithm inversion technique (Figure 3) (detailed in the data and methodology section). The misfit values between the observed and synthetic PRFs for all the stations (Figure 4) are smaller than the misfit value of 1.42, as Sambridge regarded the misfit values for stations HDWR, THRI, RWSR, PITG, BHTW, GRGA, and MUNS are 1.02, 1.13, 1.07, 1.07, 0.69, 0.95, and 1.12. In this iterative process, by matching the observed and synthetic PRFs, we get ~1000 best-fitting models and a final average model. Inversion results show that the crustal thickness beneath station HDWR, THRI, RWSR, PITG, BHTW, GRGA, and MUNS are 54.1 km ± 2 km, 50.4 km ± 2 km, 47.0 km ± 2 km, 44.2 km ± 2 km, 48.8 km ± 2 km, 52.5 km ± 2 km and 49.0 km ± 2 km respectively. The inverted velocity models obtained from the inversion of the receiver functions are shown in Figure 3. Poisson’s ratio (σ) is related with V_p/V_s (γ) ratio by a relation σ=1/2 × (1-1/(γ²-1)). We observe high V_p/V_s corresponding to a high Poisson’s ratio. The LVL zone is located at depths 20.6 km, 10.9 km, 9.0 km, 24.5 km, 10.6 km, and 11.0 km beneath stations THRI, RWSR, PITG, BHTW, GRGA, and MUNS, respectively, with a high value of V_p/V_s ratio. Different studies have also reported LVL in a few places in the Himalayas. They interpreted that this LVL is probably due to the shear heating in the ductile regime with the elevated temperature that creates a conducive environment for fractional melt. The high Poisson's ratio (or V_p/ V_s) is responsive to the mineralogical composition, the existence of fluid content, and fractional melt. The higher V_p/V_s ratio values in the LVL around the Kumaun-Garhwal Himalaya indicate the presence of fluid or partial melt in the upper crust. The LVL obtained in the study area are well consistent with the previous studies in different parts of the Himalaya i.e., NW Himalaya- Ladakh Garhwal Himalaya have also suggested the LVL beneath the high Himalayan orogenic fold and thrust belt and gangetic plains by using the inversion of the rayleigh wave group velocity.

Similar low velocities in the top crust have also been reported earlier in the NW Himalaya-Ladakh the Doda-Kistwar region and the Garhwal Himalaya. The magnetotelluric study has provided evidence of the high conductivity zone with resistivity lower than 3 Ω m between 15 km and 20 km depth due to the availability of free aqueous fluids. While other studies based on tomographic inversion also indicate low V_P and high V_p/V_s around the Kumaun-Garhwal region. Another reason for this LVL zone could be the shear heating caused by the exhumation of the greater Himalayan leucogranites and decompression melting. The results from receiver function inversion suggest a thick crust below stations HDWR (in the sub-Himalaya), GRGA, and THRI (in the lower Himalaya) due to the presence of sedimentary cover, producing a delay in the arrival of Ps phases (Figure 5). Moho variation for stations BHTW, THRI, and HDWR are 448.8 km ± 2 km, 50.4 km ± 2 km, and 54.1 km ± 2 km between the sub- Himalaya and lesser Himalaya (Figure 6). Our results are consistent with the moho depth of 48 km beneath the lesser Himalaya, as reported [20].

Figure 2: (a) Move out corrected receiver functions with one second low pass filtered. Stations are arranged with increasing order of distance from HFT. Black and gray colors denote positive and negative peaks, respectively. Red lines show P to S conversions times (Pms) from Moho. (b) Stacked receiver functions of all seismic stations. The receiver functions stacked are 47, 13, 25, 14, 41, 24, and 31 for station HDWR, THRI, RWSR, PITG, BHTW, GRGA, and MUNS. Two bounding grey lines on both sides of the PRF represent PRF amplitude added and subtracted with 2 sigma standard errors.

Figure 3: Shear wave velocity models at different stations using Neighborhood Algorithm (NA) method located in the study region in the Northwest Himalaya. The solid red line on the left side represents the V_p/V_s ratio. Stations' names are shown in bold at the top right side of each plot. Best fit 1000 models (green area) having the least misfit between the observed and synthetic receiver functions are shown within the total 20000 models (grey) generated by this study.

Figure 4: Match between observed and synthetic waveform for all stations (the station names shown in the right-hand corner). The blue-colored dot line represents the synthetic waveform obtained from the output velocity model, and the red-colored solid line represents the observed waveform.

Figure 5: Red and black crosses are the P to S conversion points corresponding to the depth at 410 km, and 660 km discontinues. Red triangles show the location of broadband seismic stations.

Figure 6: Observed stacked receiver functions (move out corrected) after filtered with a corner frequency of 3 s reveal the upper 410 km and 660 km discontinuity for five stations. Two bounding grey lines on both sides of the receiver functions represent receiver functions amplitude added and subtracted with 2 sigma standard errors. Red and blue colors represent the positive and negative conversion. The Black and green dash line shows the P_s conversion time from 410 km and 660 km discontinuity. Stations are arranged with increasing order of distance from HFT. b. 2D migrated image of the receiver functions showing the disposition of upper mantle 410 km and 660 km discontinuities. Red and blue colors represent the positive and negative conversion [21].

Mantle transition zone

The upper mantle discontinuities are related to the significant mineral phase transition. The zone between 410 km to 660 km is called the Mantle Transition Zone (MTZ), which provides information about the mantle's rheology, thermal state, and deep dynamics, which control the plate tectonics. The phase transition occurs at 410 km due to α-olivine to β-wadsleyite (Mg2SiO4) with a positive clapeyron slope. In contrast, the 660 km discontinuity occurs due to the phase transition from spinel to perovskite and magnesiowustite, with a negative clapeyron slope. Thermal equilibrium studies based on the clapeyron slope suggest that change in pressure and temperature can occur a difference of 5 km-25 km in the MTZ. In the case of colder MTZ, uplift of 410 km discontinuity and deepening of 660 km discontinuity occurs, which causes the thicker MTZ. For hot MTZ, deepening of 410 discontinuity and uplift of 660 discontinuity occurred, which caused thinner MTZ. The thickness between these two discontinuities decreased by 20 km-30 km, increasing temperature by 200°C. The P to S arrival times for the 410 km discontinuity and 660 km discontinuity are typically detected at 44.1 s and 68.1 s. To avoid overlapping conversion points at 410 km and 660 km discontinuities, we have used the Common Conversion Points (CCP) stacking technique. Stacking by single stations is meaningful if piercing points from different stations are not overlapping in the depth range studied. If several stations are closely spaced, and piercing points overlap at the study depth, summation by common regions of piercing points may be more meaningful. In this case, piercing point locations for conversions at a defined depth are computed, and all traces with piercing points inside a defined area are summed. CCP stacking is a converted phase back-propagation stacking method. In the CCP stacking process, all the traces with a common conversion point are stacked. As we know that the conversion time from 410 km to 660 km can vary significantly with the velocity structure, summation using the Common Conversion Points (CCP) will be more meaningful. Here we migrate the data from time to depth domain by using backprojection of converted phases along the incoming ray path with a 1-D velocity model. Figure 5 shows the spatial distribution of P to S conversion points at 410 and 660 km depths beneath all the stations. The stack PRFs and PRFs amplitude with depth (image map) have shown in Figure 6. We have taken only those phases for Ps conversion from 410 km and 660 km, which are beyond the ± 2 standard error limit. The Ps time conversion from 410 km and 660 km discontinuity beneath stations GRGA, HDWR, BHTW, PITG, and RWSR are beyond the error limit, while stations THRI and MUNS are within the error limit. Conversion times at 410 km vary from 42 s to 44 s, while those for the 660 km vary from 66 s to 68 s. P to S conversion times from 410 km discontinuity beneath all the stations show early arrival w.r.t the global average values estimated from the IASP91 velocity model. The thickness of the transition zone in the study area is higher than usual, indicating the significant water contamination which caused the upliftment of the 410 km discontinuity.

This upliftment of the 410 km discontinuity is feasible with a wet transition zone due to the relatively low temperature and water contamination and melting on the top of the 410 km discontinuity. The wet transition zone may be due to water content on the order of 1 wt. %. But, in a specific hydrous condition, Smyth and Frost suggest that the width of the olivine to wadsleyite phase transition will be a few tens of kilometers in place of the typical few kilometers. The thick mantle transition zone in the study region is due to the subduction of the Indian plate beneath the Eurasian plate, which causes the low temperature in the Mantel Transition Zone (MTZ). Earlier studies suggested that the cause of the shallow depth of 410 km discontinuity may be due to mantle upwelling. Observed 2 s early conversion from the 410 km discontinuity beneath central India, western Tibet, and a few parts of the Himalaya. Singh and Kumar observed a comparatively large transition zone thickness in the eastern Himalayas based on the early conversion time. Recently, proposed a thick MTZ beneath North-East Himalaya due to the subduction of the Indian plate beneath the Eurasian plates [22].

Conclusion

The study presents the seismic images of the crust, the intracrustal LVL, and the upper mantle discontinuities beneath the Kumaun-Garhwal Himalaya. The modeling result suggests that the crust-mantle boundary varies from 42 km to 54 km. The LVL varies from 11 km to 24 km, with a high Poisson's ratio beneath all stations indicating the presence of fluid in the upper crust. The LVL associated with a weak zone may reactivate the pre-existing faults to generate seismic activity in the upper to the mid crust. The early P to S conversion time at 410 discontinuity estimated from the stacked PRFs and 2D migrated image of PRFs indicates the thick transition zone in the study region, which may be due to the low temperature or the presence of water content in the MTZ. This cold and thick MTZ indicates that the Indian plate would be cooler by a few hundred degrees Celsius than the ambient mantle at this depth. Our observation suggests that the consumed part of the Indian shield has extended up to the mantle transition zone.

Acknowledgements

The authors are thankful to the director of the Wadia institute of Himalayan geology, Dehradun, India, for providing support and all necessary facilities and permission to publish this work. We also thank the SERB, Department of Science and Technology, government of India for funding the research project. Seismic Handler software is used for seismological data analysis. Generic mapping tool is used for plotting.

References

1. Acton CE, Priestley K, Mitra S, Gaur VK. Crustal structure of the Darjeeling-Sikkim Himalaya and southern Tibet. Geophy J Intern. 2011;184(2):829-852.

   [Crossref] [Google Scholar]

2. Ammon CJ, Randall GE, Zandt G. On the nonuniqueness of receiver function inversions. J Geophy Res Solid Earth. 1990;95(B10):15303-15318.

   [Crossref] [Google Scholar]

3. Ammon CJ, Zandt G. Receiver structure beneath the Southern Mojave block, California. Bull Seismol Soc Ame. 1993;83(3):737-755.

   [Google Scholar]

4. Arora BR, Unsworth MJ, Rawat G. Deep resistivity structure of the Northwest Indian Himalaya and its tectonic implications. Geophys Res Lett. 2007;34(4).

   [Crossref] [Google Scholar]

5. Avouac JP. Mountain building: From earthquakes to geological deformation. Dynamic processes in extensional and compressional settings. Treat Geophys. 2007;6:377-439.

   [Crossref] [Google Scholar]

6. Avouac JP, Meng L, Wei S, Wang T, Ampuero JP. Lower edge of locked main Himalayan thrust unzipped by the 2015 Gorkha earthquake. Nature Geosci. 2015;8(9):708-711.

   [Crossref] [Google Scholar]

7. Avouac JP, Ayoub F, Leprince S, Konca O, Helmberger DV. The 2005, Mw 7.6 Kashmir earthquake: Sub-pixel correlation of ASTER images and seismic waveforms analysis. Earth Planet Sci Lett. 2006;249(3-4):514-528.

   [Crossref] [Google Scholar]

8. Bai L, Liu H, Ritsema J, Mori J, Zhang T, Ishikawa Y, et al. Faulting structure above the main Himalayan thrust as shown by relocated aftershocks of the 2015 Mw7.8 Gorkha, Nepal, earthquake. Geophys Res Lett. 2016;43(2):637-642.

   [Crossref] [Google Scholar]

9. Benz HM, Vidale JE. Sharpness of upper mantle discontinuities determined from high frequency reflections. Nature. 1993;365(6442):147-150.

   [Crossref] [Google Scholar]

10. Berkhout AJ. Least squares inverse filtering and wavelet deconvolution. Geophysics. 1977;42(7):1369-1383.

    [Crossref] [Google Scholar]

11. Bina CR, Helffrich G. Phase transition clapeyron slopes and transition zone seismic discontinuity topography. J Geophy Res Solid Earth. 1994;99(B8):15853-15860.

    [Crossref] [Google Scholar]

12. Bollinger L, Henry P, Avouac JP. Mountain building in the Nepal Himalaya: Thermal and kinematic model. Earth Planet Sci Lett. 2006;244(1-2):58-71.

    [Crossref] [Google Scholar]

13. Bostock M, Hyndman R, Rondenay S, Peacock S. An inverted continental Moho and the serpentinization of the Cascadia forearc mantle. Nature. 2003;417:536-538.

    [Google Scholar]

14. Bostock MG. Green's functions, source signatures, and the normalization of teleseismic wave fields. J Geophy Res Solid Earth. 2004;109(B3).

    [Crossref] [Google Scholar]

15. Bouchon M, Vallée M. Observation of long supershear rupture during the magnitude 8.1 Kunlunshan earthquake. Science. 2003;301(5634):824-826.

    [Crossref] [Google Scholar]

16. Borah K, Kanna N, Rai SS, Prakasam KS. Sediment thickness beneath the Indo-Gangetic plain and Siwalik Himalaya inferred from receiver function modelling. J Asian Earth Sci. 2015;99:41-56.

    [Crossref] [Google Scholar]

17. Burdick LJ, Langston CA. Modeling crustal structure through the use of converted phases in teleseismic body wave forms. Bull Seismol Soc Amer. 1977;67(3):677-691.

    [Crossref] [Google Scholar]

18. Caldwell WB, Klemperer SL, Rai SS, Lawrence JF. Partial melt in the upper middle crust of the northwest Himalaya revealed by Rayleigh wave dispersion. Tectonophy. 2009;477(1-2):58-65.

    [Crossref] [Google Scholar]

19. Caldwell WB, Klemperer SL, Lawrence JF, Rai SS. Characterizing the main Himalayan thrust in the Garhwal Himalaya, India with receiver function CCP stacking. Earth Planet Sci Lett. 2013;367:15-27.

    [Crossref] [Google Scholar]

20. Chen J, Gaillard F, Villaros A, Yang X, Laumonier M, Jolivet L, et al. Melting conditions in the modern Tibetan crust since the Miocene. Nature Comm. 2018;9(1):3515.

    [Crossref] [Google Scholar]

21. Chen WP, Molnar P. Seismic moments of major earthquakes and the average rate of slip in central Asia. J Geophys Res. 1977;82(20):2945-2969.

    [Crossref] [Google Scholar]

22. Chen WP, Molnar P. Focal depths of intracontinental and intraplate earthquakes and their implications for the thermal and mechanical properties of the lithosphere. J Geophys Res Solid Earth. 1983;88(B5):4183-4214.

    [Crossref] [Google Scholar]