Estimation of active surface deformation using PSInSAR technique of the Central Himalayan region

Abstract

The Kumaun Himalaya is considered as the most active part of the Central Seismic gap in the Indian Sub-continent. In this paper, we presented active surface deformation rates of the Kumaun region from February 2017 to February 2021 using the Persistent Scatterer Interferometric Synthetic Aperture Radar (PSI) technique. The cumulative displacement that occurred during the span of 4 years is ±55 mm, whereas the Line of Sight (LOS) deformation velocity rate ranges is ±7 mm/yr. Apart from PSI, we also estimated the b-value of the Kumaun region from 1803 to 2021 (399 events) and its value is 0.54 ± 0.03. A distinct NE-SW trend of b-value is observed where earthquakes with M > 6 occurred towards NE of this trend. The PSI-derived deformation reveals that the central part of the Inner Lesser Himalaya along with the Main Central Thrust (MCT) zone is dominated by uplift. The zone between the Munsiari Thrust (MT) and MCT in the central region shows the maximum uplift ranging 5–7 mm/yr which exactly lies above the mid-crustal ramp of the Main Himalayan Thrust (MHT). Our results are well corroborated with available observations of geodetic strain and peak ground acceleration values. However, the deformation patterns and high-velocity rates in the central part of the study area between MT and MCT indicate the accumulation of high stress.

Keywords:

• Deformation
• PSInSAR
• Kumaun Himalaya
• Uttarakhand
• India

1. Introduction

Late Paleocene-Eocene inter-continental docking of the Indian and Eurasian plates resulted the occurrence of a large-magnitude earthquakes along 2500 km long strike distance of the lofty Himalayan terrain (Molnar and Tapponnier 1975; Patriat and Achache 1984; Thakur 1992; Rowley 1996; Hodges 2000). This vital process which shortens the lateral spreading of the lithosphere has been ongoing since this collision (Bilham et al. 1998; Avouac 2003). It is also conspicuous from the convergence along the Himalayan arc, Tibetan plateau, and Tien Shan orogen due to the eastward rise of the earth crust and southward transposition at the eastern syntaxes (Searle et al. 1987; Sinha 1989; Molnar and Lyon-Caent 1989; Wang et al. 2001; Zhang et al. 2004; Wang and Barbot 2023). These processes are still active and at present, perpetually exposing the Himalayan arc to natural hazards, particularly earthquakes, and landslides. The inter-continental convergence between India and Eurasia has led to the generation of several rapidly developing strain zones, thrusts, faults, and highly fractured and jointed rock formations in the Himalayan terrain which caused instability in the region due to seismic activities (Ambraseys and Bilham 2000; Banerjee and Bürgmann 2002). Recent studies on the Himalayan deformation suggest that southern Tibet has advanced towards India with a speed of ∼16–18 mm/yr by sliding over the top of the under thrusting Indian plate (Bilham et al. 1997; Jouanne et al. 1999; Jade et al. 2007, 2014, 2017; Banerjee et al. 2008; Ader et al. 2012; Dumka et al. 2014, 2018; Stevens and Avouac 2015; Yadav et al. 2019; Dal Zilio et al. 2020; Ghavri and Jade 2021). About 10–20 mm/yr of differential shortening rate is suggested for the entire length of the Himalayas (Jade 2004; Jade et al. 2014, 2017). The Himalayan orogenic belt has been stuck by several devastating earthquakes in the past viz., The 1897 Shillong earthquake M_w > 8; the 1905 Kangra earthquake M_w 7.8, the 1934 Bihar-Nepal earthquake M_w > 8, 1950 Tibet-Assam earthquake M_w 8.6, the 2005 Kashmir earthquake M_w 7.6, and the recent most 2015 Gorkha earthquake M_w 7.8 (Seeber and Armbruster 1981; Pandey and Molnar 1988; Ambraseys and Bilham 2000; Ambraseys and Douglas 2004; Rajendran and Rajendran 2005; Negi and Paul 2015; Bilham 2019; Ghavri and Jade 2021). Furthermore, the studies suggest that after the 1950, Tibet-Assam earthquake M_w 8.6, there is a quiescence for the occurrence of the great earthquake in the Himalayan arc. Various researchers predicted the recurrence of similar magnitude earthquakes in the central Himalayan region, which could be more devastating than the previous earthquakes because of the speedy rise in inhabitants and tremendous increase in infrastructures (Bilham et al. 2001; Bansal and Verma 2013; Ghavri and Jade 2021; Srivastava et al. 2013; Stevens and Avouac 2015, 2016).

The active tectonic belt of Kumaun Himalaya is experienced several large magnitude earthquakes of M ≥ 6 (Khattri et al. 1994; Kayal 1996; Mandal et al. 2000; Rajendran et al. 2000; Rastogi 2000; Kayal et al. 2003) (Figure 1a and b). Previous paleoseismic studies along with historical data indicate the occurrence of late Holocene earthquakes in the Central Himalaya, which had occurred around 1222–1422 AD (M_w ∼8.5), 1st September 1803 (M_w ∼7.7), 1842 (M_w ∼7.2), 1852 (M_w ∼7) implying active thrusting in central Himalayan region (Wesnousky et al. 1999; Kumar et al. 2006; 2015; Mahesh et al. 2015; Negi et al. 2017; Prasath et al. 2017; Rajendran et al. 2018; 2020). Furthermore, valley-wide tectonic-geomorphic field investigations along the major structures indicate the active nature of thrusts and faults during the Late Quaternary including the Holocene period (Valdiya 1988, 2001; Pant et al. 1992, 2007; Valdiya et al. 1992, 1999; Pandey and Pandey 2004; Luirei et al. 2006; Patidar et al. 2007; Kothyari and Pant 2008; Kothyari 2010, 2014; Kothyari and Juyal 2013; Kothyari and Luirei 2016; Kothyari et al. 2017a, 2017b, 2018, 2020a, 2020b; Tiwari et al. 2021). The elucidation of the deformation mechanism from the northernmost part of the Kumaun Himalaya to the foothills thus constitutes a serious problem, both to address the as yet inscrutable earth processes of continental collision operating at the highest rate, as well as to provide numerical estimates of seismic hazards to form a rational basis for the design of land use plans as well as of safe habitats and critical facilities in the region.

Figure 1. (a) Geological map of the Kumaun Himalaya (Modified after Valdiya (1980)). (b) The seismotectonic map of Kumaun Himalaya shows the distribution of earthquakes from 1800 to 2021 (seismic data source: National Centre for Seismology (NCS), New Delhi; https://seismo.gov.in/MIS/riseq/earthquake/archive).

The Persistent Scatterer Interferometric Synthetic Aperture Radar (PSI) is an advanced InSAR technique and is an effective tool that is capable of monitoring and estimating the earth’s surface deformation up to millimeter scale precision (Crosetto et al. 2016). The most important aspect of the InSAR is to estimate the ground deformation over a larger areal extent without physical access in compression to other geodetic measurement techniques (Kandregula et al. 2021; Kothyari et al. 2021; Dumka et al. 2022; Joshi et al. 2022; Suribabu et al. 2022; and references therein). The InSAR technique is successfully implemented in various fields of geosciences viz., subsidence monitoring (Ferretti et al. 2000), landslide monitoring (Ferretti et al. 2001; Hilley et al. 2004), seismology (Massonnet et al. 1993; Massonnet et al. 1994; Dixon 1995; Peltzer et al. 1996, 1999; Massonnet et al. 1996; Carnec and Delacourt 2000), glaciology (Mohr and Madsen 2000), hydrology (Borgeaud and Wegmüller 1997; Rodriguez and Martin 1992; Lakhote et al. 2023), ecology (Dixon 1995; Ludwig et al. 2000) and volcanology (Massonnet et al. 1995; Salvi et al. 2004). The PSI method was proposed by Ferretti et al. (2000, 2001), wherein it successfully eliminates the residual topographic error and atmospheric artifacts overcoming the limitations of its antecedent technique (i.e. Differential InSAR). However, the integrated analysis of the PSI and geodetic measurement can improve the accuracy of ground deformation patterns and come up with reasonable estimates.

In this study, we estimated the surface deformation along the active thrusts of the Kumaun Himalaya region using PSI method and conducted tectonic-geomorphic field investigations for neotectonic evidence of the major thrusts. Further, we have estimated the b-value of the earthquakes that occurred in the region from 1803 to 2021 to identify the stress accumulation zones in the region. Further, we correlated our results with GPS results, Geodetic strain, and Peak Ground Acceleration (PGA) values for better evaluation/assessment. Eventually, we provided neotectonic evidence for each thrust to understand the geodynamics of the region and the results will potentially aid in a better assessment of seismic hazard.

2. Geology of the study area

The fold-thrust belt of the Central Himalayan arc of the Kumaun region covers approximately 21,100 km² area. The major thrusts in the Kumaun Himalaya such as Main Frontal Thrust (MFT), Main Boundary Thrust (MBT), Ramgarh Thrust (RT), South Almora Thrust (SAT), North Almora Thrust (NAT), Main Central Thrust (MCT), Pindari Thrust (PT) and Trans-Himadri Fault (THF) respectively are well exposed in the region (Figure 1a). The MCT, MBT, and MFT merge at 15 to 20 km depth level with the top of the under-thrusting Indian plate beneath the Himalaya (basal decollement zone) (Auden 1935; Gansser 1964; Misra and Sharma 1967; Jain 1971; Valdiya 1976; Srivastava and Mitra 1994; Richards et al. 2005). The region is categorized into four main lithotectonic units beginning from its southern emergence above the plains to the Indus-Tsangpo Suture Zone (ITSZ) i.e. Sub-Himalaya, Lesser Himalaya, Higher Himalaya, and Tethys Himalaya (Gansser 1964; Valdiya 1979, 1980). North of MFT, lies the ruggedly youthful Siwaliks with an average height of 1500–2000 m which in turn over-thrusts the Lesser Himalaya (LH) along the MBT. The rocks of the LH zone are comprised of Cambrian and Proterozoic meta-sedimentary and are over-trusted by crystalline of the Higher Himalaya (HH) (Figure 1a). However, the MCT is categorized by Munsiari Thrust (MT) and Vaikrita Thrust (VT) which delimits the litho-units of LH and HH (Figure 1a). The HH represents several peaks between the MCT and Trans Himadri Fault (THF) constituting the crystalline rocks. The THF is also called as Malari-Martoli Fault in the NW Himalaya and acts as a boundary between HH from the Tethys Himalayan sedimentary sequence (Valdiya 1980, 2005; Yin 2006). The Tethys Himalayan boundary is demarcated amid the THF and ITSZ in southern Tibet (Valdiya 1980). The ITSZ comprises of belt of ophiolites which indicates the northern limit of the Himalayas.

The lithological distribution of Kumaun Himalaya is characterized by the presence of thick crystalline sheets along their nappe such as Almora, Baijnath, Askot, and Chhiplakot (Valdiya 1980; Paul 1998) (Figure 1a). The Lesser Himalaya that was tectonically detached from the Siwalik Himalaya by MBT encompasses four lithotectonic units e.g. autochthonous unit, Krol, Ramgarh, and Almora nappes (Figure 1a). The rocks of Tejam Group are characterized by white, brownish, purple, siliceous, and coarse-grained crystalline marble with intercalation of white talcose schist. Furthermore, the NAT and SAT demarcate the outer and inner boundaries of the crystalline zone of Almora Nappe (Figure 1b). The Munsiari Formation consists of mylonitised augen gneiss, amphibolite, sericite schist, and garnet-bearing schist. The Munsiari Formation is separated from the Tejam sedimentaries and Berinag Quartzite by the MT in the south (Valdiya 1979, 1980, 1981) (Figure 1a). The lithotectonic setup of the region is structurally represented by the presence of a tectonic duplex (Saklani and Bahuguna 1983; Johnson 1986). The thick slab of 10–15 km of basement complex is thrust, giving rise to the ‘Central Crystalline Zone (CCZ)’ which forms the root of the Tethyan sequence of Late Precambrian to Upper Cretaceous age (Valdiya 1980; Talukdar et al. 2020) (Figure 1a).

3. Seismotectonics

Occurrences of low to moderate-intensity earthquakes are partially releasing accumulated strain beneath the central Himalayan region (Bilham et al. 1998; Mahesh et al. 2015; Kothyari et al. 2021). About 30% of the convergence rates are consumed by central Himalaya (Molnar and Tapponnier 1975) and ∼ 17.7 ± 2 mm/yr is being accommodated by slip rates of large magnitude earthquakes (Molnar 1990). The seismo-tectonics of the Himalayas is mainly controlled by slip along the Main Himalayan Thrust (MHT) that appears as MFT (Berthet et al. 2013; Hammer et al. 2013). The basal detachment surface (MHT) beneath the Himalayas is believed to be sealed which is accommodating a slip deficit at a rate of 14 ± 1 mm/yr (Seeber and Armbruster 1981; Srivastava and Mitra 1994; Banerjee and Bürgmann 2002). In the Nepal Himalaya, the slip rate observed is 5 mm/yr, however, the central sector of the Himalayan arc seems to be locked and accumulating a third of the convergence slip (Bilham and Gaur 2000). GPS study carried out in western Nepal from 1995 to 1999 indicates a vertical velocity of 12–15 ± 3 mm/yr (Jouanne et al. 1999). The steady aseismic southward approach of southern Tibet and Tethys Himalaya (∼50–100 km wide) continues for a few hundred years at a rate of ∼18 mm/yr replenishing the drained strain budget (Bilham 2005). When the slip deficit (horizontal compression) behind the locked part of the Himalayas exceeds the frictional strength of the locked surfaces, the later is driven forward by the amount of accumulated slip (Bilham and Wallace 2005). Numerous moderate earthquakes are detected in the region, of which the more reliably documented 19th-century events are those of 1803 (M∼7.5), 1833 (∼M7.7), 1991 (∼M 6.6), and 1999 (∼M 6.8). Many smaller earthquakes, of course, continue to occur as the strained budget in the region approaches failure (Rajendran and Rajendran 2005) (Figure 2). However, the low accuracy of their locations has precluded their identification with geological faults, and a credible tectonic model explaining their mechanism and predicting future events is yet to be developed. Small and moderate earthquakes release only a small percentage of the accumulated slip leaving the main budget to be released by a great earthquake, as testified by the ubiquitously observed frequency-magnitude relationship (b-value) of earthquakes (Aki 1965; Utsu 1965; Bender 1983; Bilham and Gaur 2000; Wiemer 2001; Bilham 2005; Woessner and Wiemer 2005).

Figure 2. (a) b value curve of earthquake events occurred in Kumaun Himalaya show Mc = 3.2, histograms of no of earthquake vs. depth and earthquake vs. magnitude. (b) Spatial distribution of b value in Kumaun Himalaya. The earthquake epicentral distribution is highlighted by white circles.

The ESE-WNW concentration of seismicity in the central sector of Kumaun Himalaya suggests the ongoing strain-building process beneath the region between the MCT and MBT (Seeber and Armbruster 1981). The faults/thrusts crossing from the northeastern border of the Kumaun Himalaya is discharging accumulated stress in the form of earthquakes (Pant and Paul 2007) (Figure 2a). Most of the seismic events with low magnitude were situated in the area encompassed by the MCT and NAT (Pant and Paul 2007; Kothyari and Pant 2008; Joshi and Kothyari 2010; Kothyari and Luirei 2016; Talukdar et al. 2020 and reference therein; Pudi et al. 2021). This indicates that the accumulation of strain is influenced by moderate to low earthquake (M ≥ 3 and M ≤ 1.5) activity at a depth of 10–30 km (Paul and Pant 2007; Pudi et al. 2021). Another NE-SW trend of low-magnitude earthquakes (M 1.5 to 3) is also observed (Figure 2). However, near the Thal-Bageshwar area strike parallel trend is visible. Specifically, across the NAT zone, more than 100 earthquake events of magnitude M ≥ 3 and ∼20 earthquake events of magnitude M ≤ 1.5 are observed (Pant and Paul 2007). The strain accumulation is going below the Lesser Himalayan domain which is indicative of a possible large-magnitude earthquake in the future.

4. Methodology

To establish the satellite line-of-sight (LOS) time series PSI deformation, the data was collected from February 2, 2017 to February 10, 2021 (4 years) with a 12 day repeat cycle. We processed 121 images of Sentinel-1A, SLC (Single Look Complex), C-band (5.6 cm wavelength), Interferometric Wide (IW) swath of dual polarization (VV + VH), path-56 and frame-93 acquired in TOPS (Terrain Observation with Progressive Scan) mode in ascending direction. The data is downloaded from the Alaska Satellite Facility (ASF) (https://search.asf.alaska.edu) website. The TOPS mode acquires Sentinel data in three sub-swaths encompassing a total coverage area of 250 km with a ground cell range of 5*20 m. The variation of perpendicular distance of each image acquisition from the single reference scene has been provided in Figure 3a. The perpendicular baseline distance for used datasets were kept within the optimum baseline range of < i.e. 300 m. The Kumaun region was covered in sub-swath IW1, therefore, it is processed only, for all 108 images. The United States Geological Survey (USGS) Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM) of 30 m resolution is used for the elimination of the topographic phase, flattening of the interferogram, and ellipsoidal corrections and time series analysis have been performed using the SARPROZ software (Perissin 2015).

Figure 3. (a) Normal baseline graph of the Sentinel data sets used in present study. (b) PSI methodology flow chart applied for estimation of LOS displacement.

For the estimation of b-value, we have acquired the earthquake data from National Centre for Seismology (NCS), New Delhi (https://seismo.gov.in/MIS/riseq/earthquake/archive) since January, 1800 to July, 2021. We used MATLAB base Z-map software developed by Weimer (2001) for b-value estimation.

4.1. PSI methodology

In interferometry, the differences in phase changes among two radar (Sentinel-1A) images of a target area are obtained at different times with similar polarizations (Kandregula et al. 2021). The phase changes (ɸ) in the interferogram are commensurate with the wavelength (λ) and the distance between two sensors (dR). dR relies on the target al.titude analogous to a reference surface area (dR_topo), target deformation (dR_disp) occurs due to tectonic, climatic, or anthropogenic activity, atmospheric phase disturbances (dR_atmos) caused by the tropospheric signal delays and noise (dR_noise) produced by geometrical and temporal decorrelation as given in Equation (1)(1) $$\mathrm{ɸ}=2\mathrm{\pi }\frac{\mathrm{\text{dR}}}{\left(\frac{\mathrm{\lambda }}{2}\right)}=\frac{4\mathrm{\pi }}{\mathrm{\lambda }}({\mathrm{\text{dR}}}_{\mathrm{\text{topo}}}+{\mathrm{\text{dR}}}_{\mathrm{\text{disp}}}+{\mathrm{\text{dR}}}_{\mathrm{\text{atmos}}}+{\mathrm{\text{dR}}}_{\mathrm{\text{noise}}})$$(1) (Wasowski and Bovenga 2014; Yazici and Gormus 2020) (1) $$\mathrm{ɸ}=2\mathrm{\pi }\frac{\mathrm{\text{dR}}}{\left(\frac{\mathrm{\lambda }}{2}\right)}=\frac{4\mathrm{\pi }}{\mathrm{\lambda }}({\mathrm{\text{dR}}}_{\mathrm{\text{topo}}}+{\mathrm{\text{dR}}}_{\mathrm{\text{disp}}}+{\mathrm{\text{dR}}}_{\mathrm{\text{atmos}}}+{\mathrm{\text{dR}}}_{\mathrm{\text{noise}}})$$(1)

The two distinguished methods of interferometry are the Differential InSAR and PSI. In the DInSAR technique, deformation can be obtained by eliminating topographic noise by employing an effective DEM accompanied by precise orbital information. However, the DInSAR technique certainly has some limitations because of the noise produced by incoherent deviations in the target backscattering (dR_noise) and atmospheric noise (dR_atmos). To overcome these limitations, Ferretti et al. (2001) introduced the PSI technique. The PSI technique processes a large number (more than 25) of SAR scenes acquired over a certain period that identifies radar targets wherein provide a backscattered phase signal that can be quantifiable through time so that velocities of the surface deformation can be estimated in the sensors LOS (Yazici and Gormus 2020).

In this study, the PSI method was implemented using the SARPROZ software. The spatial correlations in interferogram phases are utilized to detect mild-phase variations in terrain by employing multi-temporal SAR data in SARPROZ (Kandregula et al. 2021). The flow chart depicting the methodology of the PSI method is provided in Figure 3b. The single reference scene from the data sets is automatically chosen by the SARPROZ software, after scrutinizing for any substantial weather (atmospheric) fluctuations. The user can choose to process a subset of the image or the whole image can be analyzed. After importing all the scenes (reading the metadata), based on the spatial and temporal baseline the SARPROZ automatically chooses the single reference image for coregistration purpose and downloads the orbital files of each scene and the DEM (Hooper et al. 2004; Hooper 2008) (Figure 3b). Then the area of interest (AoI)/subset of AoI can be selected. The co-registration of all the secondary images with the single reference image is accomplished by adopting an optimum tie point interval technique (Zou and Chen 2019). The accuracy of co-registration greatly depends on the tie point interval which consists of tie point density. The autocorrelation coefficient is extracted in range and azimuth is disintegrated into frequencies ranging from low to high using a wavelet transform. Eventually, the optimum tie point interval can be decided by an amalgamated analysis of fluctuations in wave amplitudes along with distance (Zou and Chen 2019). In this study, the 1900 tie-point interval is adopted as the optimum interval which is estimated based on trial and error. The obtained results from the co-registration process are cross-checked and ensured of any positional errors. The next step after the co-registration process is the permanent scatterers (PS) generation. The SARPROZ software automatically selects the PS points based on the amplitude stability index by giving an appropriate threshold value. In our analysis, we have assigned a threshold value of 0.8, as it yields the abundant number of PS points. To improve the coherence and eliminate any atmospheric artifacts, a feature called Atmospheric Phase Screening (APS) is available to estimate accurate surface displacement rates (Cheng et al. 2012; Kandregula et al. 2021). After APS, the obtained results are geocoded, and the results can be viewed in Microsoft Excel (spreadsheet format) and also in Google Earth. The uncertainties and standard deviation occurred during the processing of data set are provided in supplementary Figure 1.

4.2. b-value estimation

The earthquake probability genesis and its frequency and magnitude can be mathematically established by the b-value (Gutenberg and Richter 1944; Godano et al. 2014). The b-value can be estimated by the following relation: (2) $$\mathrm{\text{Log}}N=a-\mathit{\text{bM}}$$(2) where N is the cumulative number of shocks, and a and b are constants. ‘a’ depends upon the level of seismicity (log number of events with M = 0), and slope ‘b’ reflects upon the heterogeneity of the medium and level of stress (Mogi 1962; Scholz 1968; Wyss 1973; Urbancic et al. 1992; Wiemer and Benoit 1996; Wiemer and Wyss 1997; Kayal et al. 2012). Further, an explanation regarding the b-value is provided in Supplementary Document 1.

5. Results

5.1. b-value

We analyzed ∼399 earthquake events from 1800 to 2021 to estimate the b-value for Kumaun Himalayan. The estimated b-value is 0.54 ± 0.03, the a-value is 4.042 and the Magnitude of completeness (M_c) is 3.2 (Figure 2). A lower b-value (< 1) prevails in the Kumaun Himalayan region. Further, the comparatively high a-value indicates enhanced seismicity (Figure 2).

To understand the spatial variation of the b-value in the Kumaun Himalaya region, a color contour map has been prepared (Figure 2). The Kumaun Himalaya shows a very low b-value ranging between (0.45–0.6). The low b-value range of 0.45–0.48 is visible in the Lesser Himalayan domain which indicates a build-up of stress in this region (Figure 2). A distinct NE-SW trend of b-value ranging from 0.46 to 0.52 is observed (Figure 2). Two major earthquakes 1945 (M 6.5) and 1958 (M 6.3) occurred along this trend. Further, magnitudes greater than 6 and many moderate to low magnitude earthquake events are observed to the right side of this trend.

5.2. PSI

The cumulative time series displacement has been estimated from February 7, 2017, to February 7, 2021 (Figure 4). We have presented the results of time series cumulative displacement for every four months covering the 4 years (Figure 4). A total of 119,152 PS points were extracted from the PSI analysis. The positive values indicate that the distance between the target region and the SAR sensor in its LOS is reducing, therefore implying an uplift, while the negative values represent that the distance between the target and radar sensor in its LOS is increasing, hence denoting subsidence. The total cumulative displacement that occurred during this period is ± 55 mm (Figure 4). The average cumulative displacement observed is ± 40 mm with a standard deviation ranging from 2.7 to 7.95 mm. From Figure 2, we observed three main zones of active uplift i.e. (a) Uplift between PT, MCT, and MT; (b) Uplift towards NE of THF (c) patches of uplift between Pithoragarh and Champawat. In these three zones, the uplift gradually increases through time which indicates stress buildup (Figure 4. Prominent subsidence is observed at AT near Gwaldam and to the south of Munsiari at the lower reaches of the Goriganga river and in the Nepal region.

Figure 4. Sentinel base time series of cumulative deformation pattern estimated from 7^th Feb 2017 to 10^th Feb 2021 evaluated for every four months covering the 4 years span.

Figure 5a shows PS points extracted from the PSI analysis and their velocities. During the 4 years, the LOS mean deformation velocity is estimated to be ±10mm/yr with a 1σ deviation ranging from 0.55 to 5.63 mm (Figure 5a and b). Only those PS points were considered whose 1σ deviations ranging from 0.55 to 0.8 mm for better precision. Except for a few regions, the deformation in the study area lies in the range of ±7 mm/yr (Figure 5a). The central part of the study area around Bageshwar locality is relatively stable or shows very low deformation as compared to other regions. The Goriganga and Kali Rivers show subsidence majorly along their flow paths (Figure 5a and b). The rivers Saryu, Ramganga and Goriganga exhibit uplift in between VT/MCT and MT. Panar River shows uplift along its path whereas, the Ladiya and Lohawati Rivers are dominated by subsidence trend along their paths (Figure 5a and b).

Figure 5. (a) PSI scatterer map shows LOS velocity (mm/yr) of Kumaun Himalaya, (b) especial distribution of LOS velocities (mm/yr). The positive and negative velocities are represented by different colour codes. The permanent GPS stations occupied by Dumka et al. (2018) are highlighted by white triangles. The elevation swath and velocity profile in Figure 6 are constructed along the profile P1, P2, and P2 lines. The different river valleys are marked by blue lines.

To understand, the variations of PSI-derived mean surface velocity along the major thrusts, we constructed three orogen perpendicular topographic swath profiles with a swath width of 10 km (Figures 5b and 6). In all three profiles, it is observed that the trend at the MT increases abruptly and then it shows a sudden decrease from the MCT (Figure 5b). Further, the MBT, RT, and SAT zones exhibit an uplift. The zone between SAT and NAT displays subsidence, except in profile P1, as it shows an uplift trend which might be towing to the activity related to the Saryu River fault. From the BT, the trend starts increasing indicating an uplift until MCT in profiles P1 and P2 (Figure 6). However, this is not the same in P3 as it displays subsidence which is assumed to be caused by ground deformation along the Goriganga River. Profiles P2 and P3 show an uplift in between the MCT and PT zones. From the PT, the trend abruptly decreases and rises suddenly again from the THF/STD zone (Figure 6). From these profiles, it can be observed that the uplift is dominant in the Higher Himalayan region and also in the MCT zone. The outer lesser Himalayan region is dominated by subsidence whereas, the inner lesser Himalaya is dominated by uplift in the central part and subsidence by eastern and western parts. The Higher Himalayan region does not show any distinct pattern of uplift or subsidence (Figure 5b).

Figure 6. The integrated mean velocity (mm/yr) and Elevation (m) swath profiles constructed along P1, P2, and P3 lines show N-S variations in topography. The uplift and subsidence along the topographic swath is represented by positive and negative velocities. The velocity profile is marked by a red line against the swath topography.

6. Discussion

The surface deformation can also be affected by climatic factors like erosion, precipitation etc. Therefore, we correlated our PSI results with several parameters where the role of climate is ineffective e.g. GPS results and PGA anomaly. In the Kumaun area, the monsoon season lasts from July to September. According to Figure 4, the monsoon season did not bring about any significant changes, and the places that were already deforming during this period experienced deformation.

6.1. Spatial correlation of PSI results with GPS-derived deformation

Dumka et al. (2018) estimated the average deformation rate of the Lesser Himalaya to be 2.8 mm/yr, whereas the mean deformation rate of Inner Lesser Himalaya is 5.2 ± 1.2 mm/yr which is significantly higher than the mean deformation rate of the Outer Lesser Himalayan region i.e. 1.7 ± 1.3 mm/yr. Additionally, Dumka et al. (2018) suggested that the MCT zone is deforming at the rate of 8.7 ± 1.7 mm/yr, which is comparatively higher than the Inner Lesser Himalaya. More specifically, Dumka et al. (2018) concluded that the crustal shortening is mainly concentrated in the vicinity of NAT, MT, and MCT zones. The mean velocity of the Lesser Himalayan region estimated from PSI in the present study is ±3 mm/yr whereas, for the Inner Lesser Himalayan region it is ±5.5 mm/yr and for Outer Lesser Himalaya it is ±1.5 mm/yr which closely correlates with the GPS deformation rates of Dumka et al. 2018. For the MCT zone, the mean deformation velocity rate ranges from ± 6.4 mm/yr (Table 1). Interestingly, the surface deformation activity is majorly concentrated in the NAT, MT, and MCT zones which is in correspondence with the results of Dumka et al. 2018. Looking at the pattern, the deformation is gradually becoming prominent in the NE direction which is consistent with the Indian plate motion (Figure 4). The relative correlation of crustal deformation rates observed from PSI and GPS techniques is provided in Table 1. Table 1 clearly shows the consistency of observation points of GPS data of Dumka et al. (2018) and PSI observations along various thrusts and faults systems in Kumaun Himalaya.

Table 1. Correlation of GPS and PSI base velocities of Kumaun Himalaya (GPS results after Dumka et al. 2018).

GPS SITE	LONG	LAT	GPS base velocities (mm/yr) After Dumka et al. (2018)				PSI mean velocity (mm/yr) Present study		Tectonic structure
			Ve (mm/yr)	σVe (mm/yr)	Vn (mm/yr)	σVn (mm/yr)	Min Velocity (mm/y)	Max Velocity (mm/y)
TOPI	80.2	30.5	−6.8	1.1	−5.5	1.4	−2	5.4	THF-HW
BELJ	80.2	30.4	−5.8	1.3	−6.4	1.5	−5	5	THF
BRFU	80.2	30.4	−5	0.8	−6.6	0.9	−5	5	THF
MILA	80.2	30.4	−7.6	0.9	−5.5	1	−5	5	THF
MART	80.2	30.3	−6.5	1.2	−4.8	1.3	−1	3.4	THF
LILM	80.2	30.2	−7.9	1.1	−7.7	1.3	−1	3.9	HW-MCT
KHAL	80.2	30.1	−7.2	0.8	−4.4	0.9	0.1	3.1	MT
MUNS	80.2	30.1	−8.9	1.1	−4.2	1.3	0.1	3.1	MT
MUNT	80.2	30.1	−9.3	0.8	−4.9	0.9	0.1	3.1	MT
GIRG	80.17	30.03	−11.4	1.6	0.72	0.82	−4	2.9	MT
BALA	80.2	30	−6.6	1.1	−7.2	1.3	−0.5	4.1	HW-BT
BAGE	79.77	29.82	−0.65	1.25	0.64	0.73	−4	−0.5	BT
BERI	80	29.8	−3	0.7	−2.8	0.9	−1.5	1.26	BT
KKOT	79.9	29.78	−3.39	1.64	0.69	0.77	0.8	2.6	BT
ASKT	80.34	29.76	−5.13	0.48	0.74	0.84	−5	−0.41	AT
BAGO	80.3	29.7	0.1	1.4	−2.9	1.5	−4	3.6	BT
KANA	80.3	29.7	1.1	1.4	−7.7	1.5	−4	3.6	BT
PGRH	80.2	29.6	0.1	0.7	−2.3	0.9	0.4	2.2	KF
ALMO	79.65	29.59	0.17	0.85	0.65	0.72	−2	0.5	KT
PITH	80.23	29.59	−2.39	0.3	0.65	0.74	−0.3	0.64	KF
CHAM	80.09	29.33	−2.09	1.74	0.73	0.81	0.54	1.75	SAT
AULI	79.6	30.5	0.9	0.7	0.1	0.9	−2	1.5	MCT
LASP	80.2	30.3	−8.5	0.9	−5.5	1	−1	3.4	THF
BOGD	80.2	30.2	−5.9	1.2	−5.8	1.5	−1	3.9	HW-MCT
RGAD	80.2	30.2	−10.98	1.1	−5.6	1.3	−1	3.9	HW-MCT

6.2. Spatial correlation of PSI results with geodetic strain, Peak Ground Acceleration, and b-value variations

The results obtained from the PSI analysis were correlated with the Geodetic strain (Dumka et al. 2018; Figure 7a), Peak Ground Acceleration (PGA) (Rout et al. 2015; Figure 7b), and b-value variations. The results of the present study suggest a higher amount of deformation in the Inner Lesser Himalayan region and the MCT zone, where the maximum subsidence is 4 mm/yr and the maximum uplift is 7 mm/yr (Figure 5b). From Figure 5b, a stretch of uplifted region is observed near the vicinity of MT, which is maximum in between MT and MCT and extends up to PT where the velocity rate ranges from 5 to 7 mm/yr. This is well supported by the geodetic strain rate where this region in the MCT zone shows 0.4–0.45 μ strain/yr which is the highest in the study area (Figure 7a). The PGA values range from 0.027 to 0.260 g and end exactly at the eastern limit of this uplifted zone (Figure 7b). This trend of the PGA (0.027–0.260g) closely correlates with the NE-SW trend line of the b-value (Figures 2b and 7b). Previous studies in the Kumaun Himalaya region revealed that a mid-crustal ramp of the MHT is situated beneath the MCT (Molnar 1990; Srivastava and Mitra 1994; Avouac 2003; Wobus et al. 2006; Yin 2006; Caldwell et al. 2013; Negi et al. 2017; Dumka et al. 2018). Hence, the deformation in this region can be associated with the crustal ramp structure along the MHT where a disastrous high-magnitude earthquake will eventually occur in this region due to the stress accumulation. The PGA values increase from west to east in the study area. The zone between SAT and NAT is relatively stable and shows a deformation of ∼ ±2 mm/yr which is very well correlated with the seismicity pattern in that area where very few earthquakes are situated (Figures 2 and 5b). In this zone, the geodetic strain lies 0.2–0.35 μ strain/yr (Dumka et al. 2014, 2018) which is increasing towards the north and the PGA values range from 0.227 to 0.292g (Rout et al. 2015) (Figure 7a and b). The zone between MBT and SAT is mostly stable except near the Champawat area and in the Nepal region where it shows an uplift. The seismicity pattern correlates with our observation where only five earthquakes have occurred in this zone (Figures 2 and 5b). This was further supported by low geodetic strain values (0.15–0.25 μ strain/yr) (Dumka et al. 2018) (Figure 7a). The PGA values in this zone remain similar to that of the zone between SAT and NAT i.e. 0.227–0.292g (Rout et al. 2015) (Figure 7b). The AT near Gwaldam shows high subsidence rates i.e. 4–5 mm/yr (Figure 5b). Here, the geodetic strain shows maximum values of 0.45 μ strain/yr (Dumka et al. 2018) which is well corresponding with our results (Figure 7a). Further, the b-value in this zone lies exactly on the NE-SW trend where the values transition from 0.52 to 0.48 (Figure 2). The Goriganga River exhibits moderate uplift in the THF/STD zone, whereas this trend changes when the river crosses the MCT which shows moderate subsidence up to NCT and this trend increases abruptly from here till BT (Figure 5b). The seismicity pattern shows a cluster of earthquakes in the vicinity of THF/STD whereas very few earthquakes occurred in the region from PT to NCT and from NCT, a cluster of earthquakes has been observed up to BT (Figure 2). The geodetic strain 0.3–0.35 μ strain/yr (Dumka et al. 2014, 2018) in the THF/STD zone and it increases towards the downstream of the river where it becomes prominent from NCT to BT (Figure 7a). The PGA values show 0.227–0.260 g (Rout et al. 2015) from the THF to MT along the river, and then the PGA value increased from 0.261 to 0.292g from MT to NAT, which correlating well with our results (Figure 7b).

Figure 7. (a) Spatial distribution of GPS-derived strain map of Kumaun Himalaya (modified after Dumka et al. 2018). (b) PGA distribution map of the study area (Modified after Rout et al. 2015).

7. Conclusion

We applied PSInSAR time series technique to estimate the active surface deformation of the Kumaun Himalayan region from 2017 to 2021 and based on that the following conclusions can be drawn:

• The total cumulative displacement observed during the span of 4 years (Feb 2017–Feb 2021) in the Kumaun Himalayan region is ± 55 mm and the mean deformation rate is ±7 mm/yr.

• The high-velocity uplift zone (5–7 mm/yr) observed at the MCT zone is positioned exactly above the mid-crustal ramp of the MHT. The buildup of stress in this zone makes it the most active zone in the study area which is highly susceptible for future large magnitude earthquakes in the near future.

• Our investigations reveal, that the Inner Lesser Himalaya and the MCT region are the most actively deforming regions in the study area.

• The b-value estimated from 1803 to 2021 in the Kumaun Himalayan region is 0.54 ± 0.03 which is low indicating a buildup of stress in the area. Further, a distinct NE-SW trend of the b-value is observed where a bulge in the trend is observed near the MCT zone.

• The results obtained from the present study may be helpful in estimating the geohazard assessment and mitigation planning.

Credit authorship contribution statement

GCK designed the structure of the manuscript, RSK and GCK – Conceptualization, Investigation, Methodology, Analysis, Writing. RKD, KL AJ, AT, KM, RT – Review and Editing.

Acknowledgment

GCK AP and UB are thankful to the cluster head, Department of Petroleum Engineering and Earth Science for granting necessary facilities to carry out this research work. Kapil Malk is thankful to SARPROZ. We are thankful to all receivers and to the regional Editor, Geocarto International for the constructive comments and encouragement.

Disclosure statement

On behalf of all co-authors, I declare that there is no conflict of interest between co-authors. The funding agencies are highlighted and acknowledged in the manuscript.