Unveiling the Evolution Journey from Pangea to Present Himalayan Orogeny with Relation to Seismic Hazard Assessment

Abstract

The objective is to understand incessant seismic activities in Northwest and Central Himalayan regions. GPS data acquired (2017–2020, Nepal; 2015–2019, Uttarakhand) from 65 GNSS stations are used to generate velocity solutions with respect to International Terrestrial Reference Frame 2014 & Indian fixed reference frame to determine the site’s precise position. These velocities are further used to calculate the strain rate and prevailing convergence rate by the respective Triangulation method and Okada’s formulation. The estimated mean maximum and minimum principal strain rate are 12.19 nano strain/yr. and − 102.94 nano strain/yr. respectively. And the respective mean shear strain and dilatation are 115.13 nano strain/yr. −90.75 nano strain, which implies that Higher Himalaya observes high compression rate compared to Outer and Lesser Himalayan region. Estimations have also elucidated presence of extensional deformation in the Northwestern part of the Himalayan arc. Accordingly, in Central Himalaya, paleoliquefaction investigations have deciphered turbidites, confirming that the seismic ruptures did not reach the surface during the 2015 Gorkha earthquake. The best-fit locking depth of 14 km and convergence rate of 21 mm/yr. (Nepal) & 18 mm/yr. (Uttarakhand) are obtained. The strain budget analysis indicates that Northwest and Central Himalaya can beckon a megathrust earthquake in the future.

1. Introduction

The collision between the Indian and Eurasian Plate gave birth to the world’s youngest orogeny or fold mountain known as Himalaya. The Himalayan orogeny is very young compared to Aravallis and the Appalachian Mountains. During the mid-Permian, a shallow sea known as Tethys was present in the latitudinal area now covered by the Himalaya. In this era, the Pangea had also started splitting up into small landmasses. In this course, the Northern Eurasian landmass or Angara and the Southern Indian landmass or Gondwana deposited a reasonable amount of sediments into the Tethys, reducing the gap between the Eurasian and Indian landmass. The mountain building process took ∼55 Ma, which started in the Upper Cretaceous period and continued until the Middle Miocene. In the last stage, the formation of Siwaliks also took place that continued from Early Miocene to Neogene age [1, 2].

The convergence activity in the Himalayan arc has made it very vulnerable or important to high seismic activity [3], leading to megathrust earthquakes in the brittle part of the crust that has accumulated elastic strain for a prolonged period. Also, the complex rheology of Himalayan wedges and crustal structure is another reason for the seismic hazard potential. Past studies have shown that Uttarakhand and Western Nepal had been the most potential zone for megathrust earthquakes for 200–500 years [4, 5, 6, 7]. Therefore, in the past half-century, the Northwest Himalaya was ruptured by various small and large earthquakes. Similarly, in the Central Himalaya, the recent Gorkha earthquake of ∼Mw >7 that took place in Kathmandu, Nepal on 25 April 2015 resonated with the 1833 megathrust earthquake, thereby indicating an ongoing convergence rate of ∼14–20 mm/yr. (also consistent with the estimated value) as shown in Section no. 6 & 7 in the Himalayan region [8].

Figure 1.

A schematic diagram illustrating the study area along the northwest-central Himalayan stretch and the locations of the 65 GNSS stations. The square on the left figure is enlarged on the right highlighting the northwest Himalaya or Uttarakhand region along with the GNSS stations that are located in this realm. Note the white vectors indicate the velocity with respective to ITRF2014 and the white circles on it defines the error ellipse. While the red vectors represent the velocities derived with respect to India fixed reference frame (Jade et al. [19]) and the red circles on it define the seismicity of varying magnitude. The yellow lines feature the three Main thrusts of Himalaya; Himalayan frontal thrust (HFT), Main boundary thrust (MBT) and Main central thrust (MCT).

Figure 2.

The strain rate analysis using the triangulation approach. The red and blue vectors indicate maximum and minimum principal strain rate respectively. Note that the blue vectors are highly dominant towards the south eastern part of Nepal region indicating a high rate of compression. While in the northwest part of Uttarakhand region the red vectors are much dominant thus indicating an extension. The maximum and minimum values estimated are tabulated in Table A1.

A dense network of GNSS stations with ∼ > 5 years long period also accomplices in crustal deformation study and help in detecting the slip rate asperity on decollement to unravel the zones/areas with a high risk of seismicity. A detailed study is needed to understand the factors and precursors that trigger megathrust earthquakes. This research work estimated the surface deformation to assess the future seismic hazard potential in Northwest-Central Himalayan region using a close-knit network of 65 GNSS stations (both permanent and campaign).

2. Tectonic setting

Himalaya is located along the southern fringe of Tibetan plateau with a stretch of ∼2500 km length and 250–300 km width, is bounded by Nanga Parbhat (Indus gorge) in northwest and Namcha Barwa (Tsangpo gorge) in the Northeast. The morphologic and structural framework of Himalaya is classified as Kashmir-Punjab of length ∼ 550 km, Kumaun-Garhwal of length ∼ 320 km, Assam segments of length ∼ 400 km and Nepal of length ∼ 800 km [9]. At the east and west syntaxes, the Tsangpo River separates Himalaya from the Indo-Burmese range in the east, and the Indus River separates Hindukush in the west, respectively. The Himalaya with peaks high resides on compressed, thrusted, folded and 70 km thick (twice the thickness from the normal continental crust) continental crust, witnessed through geophysical surveys (seismic reflection and gravity) [10, 11, 12]. MFT lies adjacent to Indo-Gangetic Plain a foreland formed due to convergence between the Indian and Eurasian Plate. The Main Himalayan Thrust acts as a décollement or detachment that dips at a shallow depth. The South Tibetan Detachment, the Main Central Thrust, the Main Boundary Thrust, the Himalayan Frontal Thrust coalesce and meet [13, 14]. The Formation of the Himalaya took place in five distinct phases, which are described in Table 1 relative to their geology and age modified after [5, 15, 16].

Evolution Phase	Age of formation	Tectonic/geological process	Geology (metamorphism grade, Fossils)
Trans-Himalayan Uplift	55–35 Ma	Upliftment due to thrusting on South-directed thrust such as Gangdese Fault, Transportation of sediments by Indus and Tsangpo river, rapid cooled and eroded igneous rocks, formation of major basins such as Kargil, Kailash and Lhasa	Granitic and volcanic rock intruded by metamorphic and sedimentary rock of the South Tibetan block
Tethyan Himalaya Uplift	45–35 Ma	Extensive folding, uplift of granite gneiss domes, formation of North Himalayan thrust, ‘antecedent’ pattern of Himalayan rivers such as Indus, Satluj, Kali, Ganga, Gandaki, and Brahmaputra	Cambrian-Eocene sandstone, shale and limestone (contains fossils which includes cretaceous age ammonites), Miocene Leucogranite, Post early middle Eocene-Biotite grade and lower; Early Miocene and younger- Kyanite-sillimanite grade
Higher Himalaya Uplift	24–17 Ma	Uplift of the metamorphic rocks buried at 20–25 km depth along the MCT, separation of Tethyan Himalaya from Higher Himalaya along STD that is characterized by tectonic extension, gravitational gliding and back-folding of the sediments of Tethyan Himalaya and upliftment and exhumation of the Higher Himalaya, Klippe formation, deposition of sediments in the foothills of Himalaya	10–20 km thick schist and Gneiss, Miocene leucogranite, Central Crystalline Zone, Post-Early middle Eocene-garnet grade
Lesser Himalaya uplift	11–7 Ma	Thrusting across MBT, began of Monsoon, rapid erosion of Himalaya, excessive rate of sedimentation in the Siwalik basin,	Quartzite, marble, slate, schist and gneiss, Eocene-Chlorite grade and lower
Neotectonics	2.6–0 Ma	Uplift of Siwalik range along HFT, rapid erosion, deep and narrow gorges, abundant of coarse-grained fluvial sediments	Sandstone, mudstone (abundant of mammal fossils), Unmetamorphosed

Table 1.

Evolution phase of the Himalaya with associated tectonic processes and geology.

3. GPS velocity measurements within close-knitted GNSS stations

Himalaya has a history of seismicity or earthquakes of very high magnitude that lead to devastation and casualties. To assess the seismic hazard potential, the interseismic strain that is getting accumulated for quite a long period is evaluated with the aid of a close-knitted network of continuous and campaign GNSS stations [5, 6, 17, 18]. These stations are located exquisitely along the Northwest and Central Himalayan stretch spanning the three major thrusts (HFT, MBT and MCT, Figure 1). The Continuous GNSS stations or CORS (Continuously Operating Reference Stations) are well equipped with robust & firmly constructed receivers from Trimble and Leica, Choke-ring geodetic antenna, <10° elevation angle and the time taken for data recording is at 15 s sampling interval. Consequently, the campaign or episodic GNSS stations are equipped with a Zephyr antenna. The data is recorded for 3–4 days repeatedly at regular intervals so that the hydrological mass distribution and seasonal variation effect in this region diminish accordingly. The data is acquired from 65 stations for 2015–2019 (Uttarakhand region) and 2017–2020 (Nepal region). GPS data acquired is processed in GAMIT/GLOBK (vs. 10.71) software along with nearby IGS (International GNSS Services) stations to calculate the precise position of these 65 stations [20]. It takes two steps to complete the GPS data processing. In the first step, the daily relative position of an individual station/site is estimated with respect to clock errors and orbital parameters of the satellite, along with error predicted due to ocean-tidal effect, ionospheric electron content and atmospheric water vapor. The Finite Element Solution (FES) 2012 and International Earth Rotation system (IERS) 2010 models are employed to reduce the respective ocean tidal and Earth tidal effects. While the Global Pressure (GPT) model and Global Mapping Function (GMF) correct the tropospheric delay, which results due to wet and dry water mass [21]. The final or the second step of this process uses GLOBK modules to analyze the loosely constrained solution for obtaining the site velocity with respect to International Terrestrial Reference Frame 2008 (ITRF08) and ITRF14 Uttarakhand and Nepal region, respectively. The velocities obtained from ITRF08 is converted to ITRF14 using the HTDP software (geodesy.noaa.gov/TOOLS/Htdp/Htdp.shtml). Also, the velocities for each station is calculated with respect to fixed Indian reference Frame/Euler Pole as suggested by [19] (Table 2, Figure 1). The velocities estimated with respect to ITRF14 vary from 34.70 ± 1.5 (BADR) to 52.42 ± 0.07 (BRN2), thereby indicating crustal shortening rate in the Higher Himalaya is very high compared to the stations that are located in the Outer Himalayan region. Accordingly, the velocities obtained from India Plate fixed reference frame, vary from 23.3 (CHLM) to 0.1 (DEHR). Such variation signifies that the surface deformation in the Himalayan stretch occurs due to the thrusting of the Indian Plate below the Eurasian Plate.

Figure 3.

The colored patches show the dilatation at the Centre of the respective triangles as described in Table A1. The prominent red patch in northwest of the Uttarakhand region defines extension (with reference to Figures 2 and 4, here the rate of extension is high). And the prominent blue patch in the Southeastern part of Nepal region indicates high compression rate.

SITE	LONG	LAT	East_vel	Error(east)	North_vel	Error(north)	Fixed Indian (North)	Fixed Indian (East)
BADR	79.49	30.74	27.4	1.2	21.3	1	−10.7	−3.6
BUGG	77.89	30.08	35.9	1.1	30.6	0.9	−2.5	0.3
DEHR×	78.04	30.34	34	0.3	35	0.4	0.1	0.1
GNFC×	78.05	30.45	33.3	0.4	33	0.5	−1.9	−0.5
HHIM	78.16	29.95	31.1	1	34.5	0.8	0.1	−1.1
KARN	79.21	30.26	28.9	1.2	30.1	0.9	−5.1	−5.4
MAND×	79.27	30.45	32.9	0.7	26.1	0.6	−9.1	−1.2
MHIM	79.89	30.68	27.7	1.2	23.7	1.3	−11.6	−6.5
MMHI	77.88	30.18	31.3	0.8	32	0.8	−1.9	0.1
MUSS	78.09	30.45	31.1	0.9	30.4	0.8	−2.8	−1
PAND	79.55	30.63	33.5	1.2	24.1	1	−11.2	−0.6
PANT×	79.48	29.01	36.2	0.5	35	0.4	−0.3	1.1
PIPA	79.43	30.44	29.4	0.9	28.6	1.1	−6.6	−4.8
RAJP	78.09	30.38	27.4	1.2	32.1	0.8	−0.6	−6.2
ROOR	77.88	29.85	34.8	1.1	31	1	−2.6	1
SAHA	77.56	29.93	34.9	1	33.8	1.2	−0.9	0.9
TAPO	79.64	30.49	27.8	1.1	25.6	1.2	−9.7	−6.5
THIM	79.20	30.49	29.5	0.9	29.6	1	−5.6	−4.6
BRN2	87.27	26.51	39.06	0.05	34.97	0.06	−1.94	0.16
CHLM	85.31	28.20	31.97	0.04	13.89	0.05	−22.67	−5.42
CHWN	84.38	27.66	36.86	0.14	32.6	0.14	−3.78	−0.56
DLPA	82.81	28.98	35.42	0.1	24.7	0.12	−11.35	−0.77
DNC4	85.24	28.07	28.18	0.21	17.67	0.24	−18.88	−9.27
DNSG	83.76	28.34	35.85	0.04	27.46	0.05	−8.80	−1.00
GRHI	82.49	27.95	38.31	0.37	29.88	0.4	−6.10	1.60
HETA	85.05	27.41	37.39	0.06	33.66	0.07	−2.85	−0.37
JIR2	86.18	27.65	35.5	0.04	24.13	0.05	−12.59	−2.46
JMSM	83.74	28.80	35.25	0.04	24.33	0.06	−11.92	−1.32
KKN4	85.27	27.80	34.09	4.05	29.53	4.83	−7.03	−3.52
KUGE	85.53	27.61	35.95	0.13	33.57	0.13	−3.04	−1.84
LHAZ	91.10	29.65	46.38	0.04	15.63	0.05	−21.84	7.96
LMJG	84.57	28.17	34.49	0.04	25.86	0.05	−10.56	−2.70
NAST	85.32	27.65	32.44	0.07	34.25	0.08	−2.32	−5.27
NPGJ	81.59	28.11	38.34	0.05	34.49	0.07	−1.29	1.99
RMJT	86.55	27.30	38.21	0.04	29.38	0.05	−7.41	−0.05
SYBC	86.71	27.81	38.35	0.04	21.72	0.04	−15.10	0.31
BDRI	79.49	30.74	27.56	0.3	23.67	0.1	−5.3	12.6
MALA	79.88	30.68	30.1	0.1	22.9	0.2	−11.6	−6.5
DHAR	78.78	31.04	29.4	0.1	24.2	0.1	−13	−1.4
MUNS	80.25	30.06	30.7	0.9	25.9	1.8	−11.6	−4.9
DRCL	80.5	29.73	30.3	0.3	26.5	0.5	6.3	9.7
PANG	80.69	29.97	34.1	0.8	21.5	0.9	−10.1	−5
PHOL	78.42	30.95	29.9	0.4	27.7	0.8	−8.7	−4.6
HRMN	79.33	30.38	30.6	0.2	27.3	0.2	−8	−5.3
BAIJ	79.60	29.93	31	1.3	27.5	1.7	−7.5	−3.6
MUN2	80.24	30.06	31.13	0.6	27.96	0.3	−4.3	−7.1
AULI	79.56	30.53	33.1	0.4	27.7	0.5	−7	−4.2
GHUT	78.74	30.53	30.96	0.2	30.88	0.1	−6.6	−4.5
GUTU	78.74	30.53	32.1	0.1	31.2	0.2	−6.7	−1.8
BDKD	78.66	30.55	32.2	0.2	31.1	0.2	−5.5	−3.2
PITH	80.28	29.57	32.1	0.5	31.4	0.5	−4.2	−3
SUNR	78.13	30.79	31.4	0.1	32.3	0.3	−2.7	−3.5
PITH2	80.28	29.57	31.55	0.5	32.77	0.3	−3.3	−2.6
AGAR	78.33	30.21	34.2	0.1	32.8	0.2	−3.2	−2.5
WIHD	78.01	30.32	32.67	0.1	34.59	0.2	−2	−2.9
KUNR	78.40	30.46	33.5	0.5	33.8	0.4	−2.6	0.5
KHIR	78.88	30.16	33.8	0.1	33.6	0.1	0.7	−2.5
KORA	79.59	29.64	35.8	2.1	31.8	0.9	−0.3	−1.7
NAIN	79.45	29.38	33.7	0.2	34.2	0.3	−0.2	−1.6
NAJI	78.33	29.59	34.8	0.7	33.1	0.7	−1.5	−0.5
LOHA	80.08	29.42	34	0.3	34.3	0.2	−0.8	−1.1
NNKM	79.80	28.94	33.2	0.7	35.2	0.7	−1.3	−0.2
KOTA	79.09	29.64	34.6	0.1	34.3	0.2	−0.3	−1.2
RATH	78.60	30.81	47.1	0.7	29.1	0.2	−0.2	−0.6
GUPT	79.07	30.53	44.4	0.6	40.7	0.5	−0.2	−0.1

Table 2.

Estimated velocities from both continuous and campaign stations with respect to ITRF 2014 and Indian fixed reference frame [19].

The starred (×) stations indicate the permanent ones.

4. Strain rate estimation on a localized scale using the triangulation method

The estimated site velocities depend on the deformation rate and the translation motion [22]. The translation motion does not depend on the position in a localized zone and has a similar direction and magnitude for a particular reference frame. Whereas the deformation rate indicates a change in the displacement ($∆u_i$) or position of the site secularly, which is explained in 2D by the following equation.

$$∆u_i=L_{\mathit{ij}}∆x_j\text{for}i,j=x,y$$E1

where $∆x_j$ is the change in the site position and $L_{\mathit{ij}}$ is the displacement gradient tensor. The $L_{\mathit{ij}}$ or displacement gradient tensor is decomposed into asymmetric and symmetric parts (Eq. 2), thereby representing the rotation rate and the infinitesimal strain rate.

$$L_{\mathit{ij}}={\epsilon }_{\mathit{ij}}+{\omega }_{\mathit{ij}}=\frac{1}{2}\left(\frac{{\partial u}_i}{{\partial x}_j}+\frac{{\partial u}_j}{{\partial x}_i}\right)+\frac{1}{2}\left(\frac{{\partial u}_i}{{\partial x}_j}-\frac{{\partial u}_j}{{\partial x}_i}\right)$$E2

Within a close clustered network of GNSS stations that are distributed along the Himalayan stretch, the estimation of strain rate (on a local scale) and the rate of rotation can be carried with the aid of the Triangulation method [23]. In this method, the horizontal field velocity $({\dot{u}}_x$, ${\dot{u}}_y$) determines the deformation rate at the centroid of the triangle at three nonlinear stations. Each of the horizontal velocity component is expressed in terms of deformation rate (${\dot{\epsilon }}_{\mathit{xx}},{\dot{\epsilon }}_{\mathit{xy}},{\dot{\epsilon }}_{\mathit{yx}},{\dot{\epsilon }}_{\mathit{yy}}$) and translation motion ($t_x,t_y$) with the recognized initial site location ($x_o,y_o)$. So, the velocity field at three different stations gives six equations (i.e. 2 for each site) alongwith six unknown, and it is expressed by the following equation.

$$\left(\begin{array}{c}{}_{.}{}^1{\dot{u}}_x\\ {}_{.}{}^1{\dot{u}}_y\\ \begin{array}{c}{}_{.}{}^2{\dot{u}}_x\\ {}_{.}{}^2{\dot{u}}_y\\ \begin{array}{c}{}_{.}{}^3{\dot{u}}_x\\ {}_{.}{}^3{\dot{u}}_y\end{array}\end{array}\end{array}\right)=\left(\begin{array}{ccc}1& 0{}_{.}{}^1{x}_o{}_{.}{}^1{y}_o& \begin{array}{cc}0& 0\end{array}\\ 0& 100& \begin{array}{cc}{}_{.}{}^1{x}_o& {}_{.}{}^1{y}_o\end{array}\\ \begin{array}{c}1\\ 0\\ \begin{array}{c}1\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}0\\ 1\\ \begin{array}{c}0\\ 1\end{array}\end{array}& \begin{array}{c}{}_{.}{}^2{x}_o\\ 0\\ \begin{array}{c}{}_{.}{}^3{x}_o\\ 0\end{array}\end{array}& \begin{array}{c}{}_{.}{}^2{y}_o\\ 0\\ \begin{array}{c}{}_{.}{}^3{y}_o\\ 0\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}0\\ {}_{.}{}^2{x}_o\\ \begin{array}{c}0\\ {}_{.}{}^3{x}_o\end{array}\end{array}& \begin{array}{c}0\\ {}_{.}{}^2{y}_o\\ \begin{array}{c}0\\ {}_{.}{}^3{y}_o\end{array}\end{array}\end{array}\end{array}\right)\left(\begin{array}{c}{t}_x\\ t_y\\ \begin{array}{c}{\dot{\epsilon }}_{\mathit{xx}}\\ {\dot{\epsilon }}_{\mathit{xy}}\\ \begin{array}{c}{\dot{\epsilon }}_{\mathit{yx}}\\ {\dot{\epsilon }}_{\mathit{yy}}\end{array}\end{array}\end{array}\right)$$E3

The above equation and its solution for three non-collinear GPS stations can be written in the form of

$$\dot{u}=\mathit{Gm}$$E4

$$m=G^{-1}\dot{u}$$E5

where $\dot{u}$ is a 1 × 6 column matrix of known instantaneous displacement/velocity vectors at three stations, m is a 1 × 6 column matrix of unknown translation vector and the deformation gradient tensor, and G is a square matrix of 6 × 6 order, and G⁻¹ is its inverse matrix which includes zeros, ones and the location vector coordinates of the three GNSS stations.

In a triangulation network, the translation gradient, rotational gradient, maximum shear strain rate and dilatation is estimated using the derived model parameters [24].

The translation gradient is calculated as:

$$\mathit{\text{Speed}}=\sqrt{t_x^2+t_y^2}$$E6

The rotation rate is expressed as:

$$\dot{\omega }=\left(\begin{array}{cc}0& \frac{\left({\dot{\mathrm{\epsilon }}}_{\mathrm{xy}}-{\dot{\epsilon }}_{\mathrm{yx}}\right)}{2}\\ \frac{\left({\dot{\mathrm{\epsilon }}}_{\mathrm{yx}}-{\dot{\epsilon }}_{\mathrm{xy}}\right)}{2}& 0\end{array}\right)$$E7

2D Lagrangian strain rate tensor described as;

$${\dot{\epsilon }}_{\mathit{ij}}=\left(\begin{array}{cc}{\dot{\epsilon }}_{\mathit{xx}}& \frac{{\dot{\epsilon }}_{\mathit{xy}}+{\dot{\epsilon }}_{\mathit{yx}}}{2}\\ \frac{{\dot{\epsilon }}_{\mathit{yx}}+{\dot{\epsilon }}_{\mathit{xy}}}{2}& {\dot{\epsilon }}_{\mathit{yy}}\end{array}\right)$$E8

The magnitude and orientation of maximum $(e_1$) and minimum ($e_2$) principal strain rate tensor. Also, compression and extension in a regime are determined by the respective negative and positive values of principal strain rates. And the maximum infinitesimal shear strain (${\gamma }_{\max }$) is expressed by the following equation.

$${\gamma }_{\max }=e_1-e_2=2\sqrt{{\left(\frac{{\epsilon }_{\mathit{xx}}-{\epsilon }_{\mathit{yy}}}{2}\right)}^2+{\left({\epsilon }_{\mathit{xy}}\right)}^2}$$E9

The first invariant of the 2D strain tensor is the areal strain or dilatation ((ellipse area-circle area)/circle area) and is explained by the following equation.

$$\text{Area strain rate}/\text{Dilatation}=e_1+e_2$$E10

5. Inferences from localized strain rate estimation

The strain rate and rotation rate are estimated at the centroid of 89 triangular zones of the 65 sites forming an angle of 30° or more (Table A1). The maximum strain rate varies from −47.90 to 195.22 nano strain/yr. with an azimuth variation from 23.21°N to 161.02° N. Similarly, the minimum strain rate varies from −333.606 to −8.9 nano strain/yr. with an azimuth variation from 178.55 °N to 0.07 °N. The value of maximum shear strain lies between 13.47 to 301.51 nano strain/yr., and the angle of rotation range from −9.9 E-0.7 to 9.12E-0.7 (°/yr). Along with this, the change in dilatation or areal strain differs between 136.07 to −381.405 nano strain. From these observations, it can be stated that the stations (mainly station number 71–86) located in the Higher Himalaya (South East) of the Nepal region experiences a high rate of compression (Figure 2, Table A1), clockwise motion with a mean shear value of 115.13 nano strain/yr. (Figure 4) and negative dilatation (highlighted by blue to yellow color, Figure 3). On proceeding towards the Northwest Himalaya, the stations located there observe a low compression rate (especially along the MBT boundary). At the same time, the cluster of stations (mainly station number 36–50) located in the Higher Himalaya of Uttarakhand region (Northwest) reflects positive dilatation (highlighted by dark red color, Figure 3), anticlockwise movement (Figure 4) and positive strain rate values (Figure 2) or extension.

S.NO.	Rotation ± uncertainty (degrees/yr)	Max horizontal extension (e1H) (nano-strain)	Azimuth of S1H (degrees)	Min horizontal extension (e2H) (nano-strain)	Azimuth of S2H (degrees)	Max shear strain (nano-strain)	Area strain (nano-strain)	Long (°)	Lat (°)
1	5.36E-08	0.755184	106.9466	-78.8757	16.94656	79.63087	-78.1205	82.01544	28.97287
2	-4.04E-07	2.578359	103.6807	-99.4337	13.68069	102.012	-96.8553	81.31119	29.27868
3	-4.47E-07	3.897905	104.6112	-97.7991	14.6112	101.6971	-93.9012	81.00014	29.33331
4	-8.06E-07	12.2127	110.9687	-87.5555	20.96874	99.76821	-75.3428	80.89607	28.98761
5	-6.17E-07	11.00359	108.7238	-91.4848	18.72384	102.4884	-80.4813	80.59397	29.04022
6	-1.21E-06	23.8691	108.2314	-99.8852	18.23138	123.7543	-76.0161	80.58626	29.16166
7	-1.21E-06	17.19301	111.0656	-105.442	21.06559	122.6353	-88.2493	80.46837	29.25081
8	-4.80E-07	9.787359	107.9763	-88.2259	17.97631	98.01327	-78.4386	80.21492	29.23759
9	5.15E-07	17.46843	114.1938	-75.7052	24.1938	93.17367	-58.2368	81.94994	28.89188
10	3.83E-07	17.42987	117.4941	-77.7848	27.49407	95.21467	-60.3549	80.85491	29.36348
11	2.50E-07	13.2956	114.4983	-70.6135	24.49831	83.90915	-57.3179	81.87103	29.00414
12	5.06E-07	6.954867	118.3468	-72.1839	28.34675	79.13882	-65.2291	80.06962	29.26691
13	8.85E-08	15.74241	111.9557	-76.4613	21.9557	92.20373	-60.7189	81.34453	29.01988
14	-2.13E-07	5.76327	115.0188	-61.9755	25.01882	67.73874	-56.2122	81.6495	28.86617
15	-3.89E-07	6.863233	114.4877	-70.1774	24.4877	77.04068	-63.3142	81.33987	28.92157
16	2.63E-07	22.81866	112.4271	-77.6835	22.42707	100.5022	-54.8649	81.64041	28.94797
17	-1.50E-07	23.51346	116.8895	-70.4583	26.88947	93.97175	-46.9448	81.56862	28.89528
18	1.99E-07	19.41619	118.9625	-63.8279	28.96255	83.24404	-44.4117	81.87805	28.83936
19	-3.99E-07	22.79689	117.718	-45.9614	27.71804	68.75828	-23.1645	80.9715	29.17439
20	-4.13E-07	20.65258	119.0297	-67.1152	29.02969	87.76776	-46.4626	81.50303	28.84527
21	-6.32E-08	18.17555	121.5551	-57.198	31.5551	75.3736	-39.0225	81.81234	28.7895
22	8.85E-08	15.74241	111.9557	-76.4613	21.9557	92.20373	-60.7189	81.34453	29.01988
23	-3.27E-07	12.22406	115.5537	-69.9911	25.5537	82.21512	-57.767	81.17633	28.92427
24	1.84E-06	27.0679	122.0977	-38.0549	32.09772	65.12281	-10.987	79.67635	29.27331
25	-4.26E-07	-42.8061	161.0264	-56.2841	71.02641	13.47797	-99.0902	79.78683	29.4118
26	5.10E-07	0.804883	119.2044	-91.7348	29.20436	92.53966	-90.9299	79.8613	29.56377
27	2.93E-07	9.469466	117.993	-77.5553	27.99301	87.02474	-68.0858	80.44904	29.06996
28	-4.64E-07	10.72079	108.4077	-67.2098	18.40767	77.93056	-56.489	80.2173	29.25883
29	7.92E-07	14.78235	120.7276	-68.7109	30.7276	83.49326	-53.9286	79.65318	29.73983
30	4.85E-07	68.3212	122.4197	-116.211	32.41966	184.5321	-47.8897	79.83462	29.96318
31	-1.16E-06	27.92151	109.0378	-119.019	19.03778	146.9404	-91.0974	79.91971	29.96864
32	7.26E-07	-11.1356	139.2819	-198.881	49.28187	187.7449	-210.016	78.94783	30.32082
33	6.86E-07	-22.1811	142.0982	-57.3622	52.09815	35.18111	-79.5433	78.6589	30.27227
34	9.12E-07	-40.4595	88.55891	-60.9931	178.5589	20.53359	-101.453	78.62523	30.38797
35	-1.32E-07	44.33593	112.9551	-91.363	22.95505	135.6989	-47.0271	79.89045	29.91392
36	-1.16E-06	19.0765	135.5199	-81.4348	45.51989	100.5113	-62.3583	78.75035	30.43268
37	1.42E-06	7.120198	84.47429	-134.444	174.4743	141.5642	-127.324	79.67702	30.6391
38	-7.74E-07	9.237582	122.8693	-159.482	32.86931	168.7197	-150.244	78.93326	30.62309
39	-1.59E-06	-2.95515	127.2037	-106.236	37.20367	103.2811	-109.191	78.77293	30.7093
40	-6.95E-06	42.29392	23.21754	-174.236	113.2175	216.53	-131.942	78.7894	30.4542
41	-6.08E-06	195.2231	16.63236	-59.1451	106.6324	254.3682	136.078	78.85412	30.1137
42	-1.98E-06	119.6129	31.89535	-25.1974	121.8954	144.8103	94.41548	79.01866	29.94196
43	-5.04E-07	56.26886	150.589	-17.0123	60.58895	73.28111	39.25661	78.497	30.28397
44	3.95E-06	113.3468	126.2494	-103.518	36.24939	216.865	9.828563	78.44253	30.27166
45	1.83E-07	13.50632	159.8563	-8.96657	69.85626	22.47289	4.539758	78.26003	29.90215
46	2.17E-07	12.03511	154.3127	-12.1777	64.3127	24.21281	-0.1426	78.1004	29.99925
47	2.62E-06	47.60433	171.3537	-21.4498	81.35375	69.05415	26.1545	78.09562	30.2528
48	5.16E-07	-14.7983	116.051	-97.2158	26.051	82.41755	-112.014	78.40788	30.64831
49	-8.64E-07	5.991775	110.0822	-55.0814	20.0822	61.07316	-49.0896	79.02436	30.1249
50	-1.19E-07	29.3643	120.2667	-105.102	30.26673	134.4663	-75.7377	78.70181	30.65096
51	-9.85E-07	25.4069	126.6595	-60.1971	36.65949	85.60399	-34.7902	82.30015	28.3516
52	-1.36E-07	-2.76345	106.9935	-78.4504	16.99351	75.68695	-81.2138	82.72469	28.48492
53	-6.24E-08	-3.59438	104.1042	-81.4359	14.10417	77.84151	-85.0303	83.2421	28.20207
54	-6.48E-07	-7.67581	108.3235	-61.3681	18.32355	53.69227	-69.0439	83.30321	28.37135
55	-3.49E-07	-5.66473	108.4554	-74.8367	18.45539	69.17195	-80.5014	82.99606	28.43045
56	1.61E-09	-3.39801	104.8787	-82.0956	14.87865	78.6976	-85.4936	82.93539	28.261
57	5.42E-07	-7.34975	108.8526	-48.5599	18.8526	41.21012	-55.9096	83.02396	28.42762
58	3.15E-07	0.965328	112.7242	-49.3059	22.7242	50.27122	-48.3406	83.01617	28.58104
59	2.06E-07	-6.9486	116.768	-50.2089	26.76797	43.26029	-57.1575	83.295	28.37257
60	-9.76E-07	11.55383	130.9644	-62.7837	40.96443	74.33758	-51.2299	83.54244	28.38635
61	-7.73E-07	-5.95005	136.1987	-64.1231	46.19874	58.17304	-70.0731	83.52032	28.34257
62	5.14E-07	-11.7672	110.0825	-48.6664	20.08247	36.89918	-60.4335	83.53281	28.25053
63	8.67E-07	-22.0907	98.77466	-48.6519	8.774663	26.56119	-70.7426	83.54998	28.20266
64	2.47E-07	-14.3185	100.9803	-54.7164	10.98025	40.39787	-69.0349	83.83905	28.19027
65	4.03E-07	-21.2864	87.3828	-68.9126	177.3828	47.62626	-90.199	83.85607	28.14237
66	6.75E-07	-6.42711	100.1659	-61.7546	10.1659	55.32747	-68.1817	83.76443	28.0608
67	6.77E-07	-9.64322	97.05148	-64.5422	7.051483	54.89894	-74.1854	83.54023	28.14369
68	-1.16E-07	-12.0751	113.4563	-44.9703	23.4563	32.89521	-57.0454	83.60176	28.31282
69	8.06E-07	-9.65619	99.36524	-67.3229	9.365238	57.66676	-76.9791	83.54729	27.99028
70	-1.41E-06	1.666171	111.5565	-227.039	21.55649	228.7048	-225.372	84.71106	27.99063
71	-1.77E-06	1.997162	113.5982	-185.794	23.59824	187.7912	-183.797	84.4879	28.075
72	1.82E-07	2.16203	107.8961	-121.065	17.89611	123.2265	-118.902	84.67058	27.75272
73	1.18E-06	0.762812	115.1201	-242.47	25.12014	243.2332	-241.708	84.89394	27.72048
74	4.45E-07	-19.6028	104.9755	-121.595	14.97547	101.9927	-141.198	84.74588	27.88156
75	1.12E-06	-47.9061	94.27342	-126.353	4.273422	78.44712	-174.259	84.76245	27.83363
76	-3.05E-06	22.5541	117.7994	-232.354	27.79936	254.9081	-209.8	84.97985	27.93243
77	-1.75E-06	-14.9844	127.9858	-273.508	37.98576	258.5236	-288.492	84.95728	27.88894
78	-5.29E-07	9.224338	109.923	-139.637	19.923	148.8609	-130.412	85.27155	27.75039
79	-4.47E-07	21.88203	104.3644	-123.731	14.36442	145.6126	-101.849	85.44619	27.80415
80	1.53E-06	-7.35197	97.29582	-308.868	7.295815	301.5161	-316.22	85.07896	27.83218
81	4.42E-07	-47.7993	90.07481	-333.606	0.074808	285.8065	-381.405	85.00853	27.84473
82	4.87E-07	6.189758	104.0459	-253.488	14.04587	259.6779	-247.298	85.29538	27.84613
83	3.13E-07	10.3839	105.5967	-244.154	15.59665	254.5374	-233.77	85.41778	27.72976
84	1.55E-06	30.97821	99.49511	-228.562	9.495111	259.5402	-197.584	85.69234	27.81402
85	7.23E-07	7.642276	102.2797	-225.582	12.27966	233.2245	-217.94	85.51755	27.76057
86	2.64E-06	40.61847	108.5187	-236.438	18.51867	277.0563	-195.819	85.49479	27.71716
87	-2.53E-06	48.25904	100.9263	-128.556	10.92632	176.8148	-80.2968	86.44012	27.33266
88	3.76E-07	60.1606	92.38606	-135.83	2.386064	195.9902	-75.669	86.48312	27.59231
89	-1.19E-06	10.99028	100.9509	-113.997	10.95087	124.9877	-103.007	86.34737	27.25283

Table A1.

The estimated maximum & minimum principal strain rate, maximum shear strain rate, rotation rate, dilatation with their corresponding azimuth and coordinated by Triangulation approach.

Figure 4.

The maximum shear strain rate is estimated by the triangulation approach. Note the green and blue lines represent right lateral (RL) and left lateral (LL) respectively. The respective estimates and direction of the angular rotation is described in Table A1.

Figure 5.

A 2D second invariant map estimated by the triangulation approach showing the spatial distribution of seismicity. Note the hollow red circles are the seismic events ranging in magnitude from 4 to >5.

Figure 6.

The rate of ongoing convergence in the Nepal region is estimated along the arc parallel or fault parallel direction. Note the green arrows are the estimated convergence rate while the red arrows define modeled convergence rate. The hollow blue circles indicate seismic events (2017–2021) of magnitude varying from 4 to >5.

The second invariant of the strain rate field, where the second invariant is defined as;

$$\dot{}{{\epsilon }^2}_{\mathit{\phi \phi }}+\dot{}{{\epsilon }^2}_{\mathit{\theta \theta }}+2\dot{}{{\epsilon }^2}_{\mathit{\phi \theta }}$$E11

where ˙ε_φφ, ˙ε_θθ and ˙ε_φθ are the horizontal components of the strain rate tensor. When the second invariant is estimated, the whole tectonic regime experiences moderate to high areal dilatation (color varies from yellow to orange to red, Figure 5), leaving out certain portions of the Eastern or NNE region (blue to yellow color, Figure 5). The 2D second invariant to the strain rate tensor has been estimated independently from the principal components of strain (Figure 5) to establish a relation between the rate of deformation, seismicity and the interseismic strain accumulating for a prolonged period [25, 26]. The map of the second invariant also focuses on the spatial distribution of earthquakes (Figure 5).

Figure 7.

The rate of ongoing convergence in the Uttarakhand region is estimated along the arc parallel or fault parallel direction. Note the green arrows are the estimated convergence rate while the red arrows define modeled convergence rate. The hollow blue circles indicate seismic events (1960–2021) of magnitude varying from 4 to >6.

In an orogen of convergence, observation of such paradox situation becomes quite inquisitive [27]. Past studies have shown that prominent extensional deformation exists in an orogen of convergence like Andean cordillera, Scandinavian Caledonides etc. [28, 29]. In the Himalayan realm, intra-continental tectonic features, prominent folds, and several other deformational geological structures indicate a compressional tectonic area. But in recent studies, it has been observed that normal faults and many other extensional geological structures are also distributed in different Himalayan sectors [30], but these are not consistent with ongoing seismicity and thrusting in the southward direction. Also, the presence of large scale lineament, shifting and tilting of river, subsidence of older rocks, upliftment of river terraces & active thrusts and normal faults overriding the young Holocene sediments demonstrates that extensional and compressional features co-exist with each other in Himalaya [31, 32].

6. Estimation of ongoing convergence rate and its future implications

The rate of convergence along the active plate boundaries determines the seismic potential and the rate of seismicity occurring today. Previous studies have stated that currently, MFT is the most active fault because here, the crustal shortening is taking place due to the plate movement and evidence of GPS vectors [33]. The site velocities (obtained with respect to ITRF) are distributed along the strike direction (fault parallel or strike slip fault) (Figures 6 and 7) and across the strike direction (fault normal, or oblique fault) (Figures 8 and 9), and the structural trend of MFT is taken as 303° and 288° in the Nepal & Uttarakhand region. It is assumed that the fault normal velocity as dip slip and fault parallel velocity as strike slip on the plate interface, ‘i.e.’ MHT (Main Himalayan Thrust). In fault normal velocity, the uncertainty is calculated by applying the following formula

Figure 8.

The rate of ongoing convergence in the Nepal region is estimated across the arc normal or fault normal direction. Note the green arrows are the estimated convergence rate while the red arrows define modeled convergence rate. The hollow blue circles indicate seismic events (2017–2021) of magnitude varying from 4 to >5.

Figure 9.

The rate of ongoing convergence in the Uttarakhand region is estimated across the arc normal or fault normal direction. Note the green arrows are the estimated convergence rate while the red arrows define modeled convergence rate. The hollow blue circles indicate seismic events (1960–2021) of magnitude varying from 4 to >6.

$${\sigma }_{\mathit{\text{horz}}}=\sqrt{\frac{{\sigma }_E^2{\sigma }_N^2}{{\sigma }_N^2{\mathit{cos}}^2\theta +{\sigma }_E^2{\mathit{sin}}^2\theta }}$$E12

where ${\sigma }_E$ is the standard deviation of the east velocity, ${\sigma }_N$ is the standard deviation of the north velocity, and $\theta$ is the bearing of the convergence vector, counterclockwise from east. At the site (GNSS stations) location, the subsurface deformation generally depends on several fault parameters: rake, dip, strike, fault length, width, etc. This is due to finite subsurface dislocation, and it can be explained by the back slip approach as expressed by [34] in the following expression

$$\mathrm{u}=\mathrm{Gs}$$E13

where u is the surface deformation vectors, s is the vector of parameters, and G is the Green’s functions matrix computed using the semi-analytical formulation published by [35].

The Okada dip slip dislocation approach [35] estimates the surface deformation at each GNSS site due to a finite rectangular dip slip dislocation on MHT for an elastic half-space. In this approach, it is presumed that the frontal portion of MHT is completely locked to some extent and aseismic creeping occurs with a uniform slip rate at downdip zone of MHT. After that, reduced chi-square uncertainty is calculated because a grid search analysis was required to estimate the best fit value of slip rate and locking depth. The expression used for calculating the chi-square uncertainty is as follows;

$${\chi }_v^2=\frac{1}{{\left(n-P\right)}^2}\sum _{i=1}^n{\left(\frac{r_i}{{\mathit{f\sigma }}_i}\right)}^2$$E14

where n is the number of observations, P is the number of free parameters, $r_i$ is the residual between observed and calculated velocity, ${\sigma }_i$ is the formal data uncertainty, and f is the data scaling factor. The scaling factors helps to avoid the inclusion of any additional uncertainty and prevent influencing various other types of data. Also, the formal GPS velocity uncertainty estimates are generally underestimated by a factor of 2 to 5 [36].

7. Inferences from ongoing plate convergence and seismicity assesment

The plate convergence rate (fault normal) or the elastic strain accumulated in the brittle frontal part of MHT is 21 mm/yr. and 18 mm/yr. for Nepal and the Uttarakhand region, respectively, with 14 km locked depth and dip angle 10° (Nepal) and 7° (Uttarakhand) (Figures 10 and 11). But the fault parallel motion is relatively small (3–4 mm/yr) and varies inconsistently (Figures 6 and 7), so the fault normal motion is only considered. The slip deficit in the Uttarakhand region is calculated as 3.9 m, indicating that a significant amount of strain had accumulated since the last magnitude earthquake (1803 earthquake) [37]. And the estimated slip deficit corresponds to seismic moment ∼4.9 × 10²¹ Nm and M_w ∼ 8.3 by taking fault length 412 km and assuming its width as 111 Km. Therefore, the strain budget analysis in the Uttarakhand region signifies that this area bears the potential to produce a megathrust earthquake, and the results are consistent with [6, 16, 38, 39]. While the Nepal region was struck by a recent high magnitude (Mw >7) Gorkha earthquake that immensely disturbed the tectonic activity of this region. It has been observed that the 2015 seismicity took place along the planar MHT and also gives good information about ground motion frequencies, the direction of rupture etc. It is also presumed that earthquakes with M_w > 8 can rupture the surface devastatingly. But the earthquakes with M_w between 7 and 8 rupture mainly within depth and rarely reach the surface. Also, Paleoliquefaction investigations have deciphered turbidites that formed due to 2015 seismic activity from the Rara Lake, which presents alternative evidence for ruptures that did not reach the surface [40, 41, 42, 43, 44]. Therefore, it can be concluded that interseismic strain is building up within the Nepal Himalayan region to beckon the next calamitous earthquake in the near future.

3D reference model (for both Nepal and the Uttarakhand region) illustrating the predominant occurrence of earthquakes ranging from low to high magnitude is developed using the parameters that have been applied in Okada’s formulation. This model discusses the evolution of slip rate spatiotemporally and the related seismic events that happened in and around the Main Central Thrust fault, thereby indicating that the locations of hypocentre fall in the Higher Himalayan region (Figures 10 and 11) [41, 45, 46]. But previous studies have shown that many historic megathrust earthquakes have ruptured partially and did not reach the surface [7, 47]. Therefore, through magnitude moment calculations and 3D reference model, it can be presumed that a megathrust event similar to the Gorkha earthquake or even more catastrophic than that will occur instantaneously.

Figure 10.

A 3D reference model showing the cross section of the Himalayan geometry in Nepal region. Note the locked width is 96 km, locked depth is 14 km, dip angle 10° and slip rate estimated is 21 mm/yr. the hollow blue circles represent seismic events from 2017 to 2020 ranging between 3 to >5. Also, the concentration of seismicity is high near the Main central thrust, thereby indicating hypocentre of these earthquakes falls in and around the higher Himalaya.

Figure 11.

A 3D reference model showing the cross section of the Himalayan geometry in Uttarakhand region. Note the locked width is 111 km, locked depth is 14 km, dip angle 7° and slip rate estimated is 18 mm/yr. the hollow blue circles represent seismic events from 2015 to 2021 ranging between 3 to >5. Also, the concentration of seismicity is high near the Main central thrust, thereby indicating hypocentre of these earthquakes falls in and around the higher Himalaya.

8. Conclusion

The surface deformation rate is measured along the Northwest -Central Himalayan belt using 65 GNSS stations located close to the wide spatial distribution. GPS data acquired (2015–2019, Uttarakhand and 2017–2020, Nepal) from these GNSS stations are used to estimate the velocity solution and determine the precise position of each site with respect to ITRF2014 and Indian fixed reference frame as suggested by Jade et al. [19]. These estimated velocities are then used to calculate the strain rate and prevailing convergence rate of the Himalayan stretch using the Triangulation method and Okada formulation. From the triangulation approach, the mean observed values of maximum principal strain are 12.19 nano strain/yr., minimum principal strain is −102.94 nano strain/yr., maximum shear strain is 115.13 nano strain/yr. and dilatation of −90.75 nano strain. Although the velocity solution from ITRF2014 implies an increase in convergence rate towards the Higher Himalaya, but results from triangulation have also featured comparatively high extensional deformation in the Northwest (Uttarakhand region) than in the Nepal region. Accordingly, using the grid search analysis, the best fit convergence calculated is 21 mm/yr. and 18 mm/yr. in Nepal and Uttarakhand regions. The locking depth is 14 km in both the regions with dip angle 10° (Nepal) and 7° (Uttarakhand). Considering the last megathrust earthquake in the Uttarakhand region (1803 earthquake), the slip deficit estimated is ∼3.9 m corresponding to ∼4.9 × 10²¹ Nm seismic moment and M_w > 8. Whereas in Nepal, after the 2015 Gorkha earthquake (M_w 7.3), the tectonic regime is disturbed greatly and has significantly impacted the infrastructure and the GNSS stations present here. Secondary evidences like turbidites from palaeoseismic investigation show that 2015 seismic activity has not ruptured the surface to a megathrust front. This implies building up interseismic stress waiting to release into a calamitous earthquake in the near future.

Acknowledgments

The authors are grateful to the Director of the Indian Institute of Remote Sensing for his consistent encouragement in carrying out this study. We acknowledge those scholars who have assisted us in the field surveys and those who have helped directly and indirectly in this work. We are thankful to Dr. Robert W. King (MIT) for providing GAMIT/GLOBK 10.70 software for processing GPS data. The authors heartily acknowledges the Editor-in-chief and Reviewer for helping in exemplifying their work exquisitely.