Crustal velocity and interseismic strain-rate on possible zones for large earthquakes in the Garhwal–Kumaun Himalaya

The possibility of a major earthquake like 2015 Gorkha–Nepal or even greater is anticipated in the Garhwal–Kumaun region in the Central Seismic Gap of the NW Himalaya. The interseismic strain-rate from GPS derived crustal velocities show multifaceted strain-rate pattern in the region and are classified into four different strain-rate zones. Besides compressional, we identified two NE–SW orienting low strain rate (~ 20 nstrain/a) zones; namely, the Ramganga-Baijro and the Nainital-Almora, where large earthquakes can occur. These zones have surface locking widths of ~ 72 and ~ 75 km respectively from the Frontal to the Outer Lesser Himalaya, where no significant surface rupture and associated large earthquakes were observed for the last 100 years. However, strain reducing extensional deformation zone that appears sandwiched between the low strain-rate zones pose uncertainties on the occurences of large earthquakes in the locked zone. Nevertheless, such zone acts as a conduit to transfer strain from the compressional zone (> 100 nstrain/a) to the deforming frontal active fault systems. We also observed a curvilinear surface strain-rate pattern in the Chamoli cluster and explained how asymmetric crustal accommodation processes at the northwest and the southeast edges of the Almora Klippe, cause clockwise rotational couple on the upper crust moving over the MHT.


Results and discussions
Regional crustal velocity and surface convergence. GPS data from Wadia Institute of Himalayan Geology (WIHG) local network were processed along with the surrounding International GNSS Service (IGS) stations data using the Gamit/Globk software and estimated the site velocities. Details on data and the processing are given in the Data and Methodology section. The estimated regional crustal velocities in the International Terrestrial Reference Frame (ITRF) are shown in Fig. 2a and values are given in Table 1. In the Himalayan mobile belt the crustal velocities are oriented towards the northeast as that of the movement of the Indian Plate 13,14 . However, crustal velocities along the Himalayan arc increases from 48.62 ± 0.10 mm/a in the Garhwal Sub-Himalaya to 49.63 ± 1.05 mm/a in the Kumaun Sub-Himalaya. Whereas, across the arc velocities decreases from 48.81 ± 0.14 mm/a in the south of HFT to 41.24 ± 0.02 mm/a in the Higher Himalaya. While in the Indian www.nature.com/scientificreports/ reference frame 32 (Fig. 2b) almost all the stations in the Higher and the Lesser Himalaya situated at the west of river Kali shows relatively larger velocities and oriented towards SW. But stations at the Sub Himalaya (TSW1) and the Frontal Himalaya (HARI) show lower velocities of 1.62 ± 0.10 mm/a and 0.36 ± 0.14 mm/a respectively, which emphasize that the detachment fault MHT is locked towards the Frontal Himalaya. But beyond the locking zone and towards the north, the Lesser and Higher Himalayan crustal velocity scenarios are different. The inner Lesser Himalayan station GHUT (Ghuttu) and the MUNS (Munsyari) station in the Higher Himalayan mobile belt show much higher velocities of 4.02 ± 0.06 mm/a and 9.20 ± 0.02 mm/a respectively. The shortening rate of the Indian crust has been estimated by taking the residual velocity between the Haridwar station situated at the south of HFT and the Hanle station 33    Apart from main thrusts, there are local thrusts; such as Tons Thrust, North Almora Thrust, and the South Almora Thrust. These thrusts are actively regulating the surface strain pattern in the Garhwal-Kumaun Himalaya. The south-dipping Tons thrust separates the inner and the outer Lesser Himalaya along with the North Almora and the South Almora thrusts. And also, the north and the south Almora thrusts demarcate the northern and the southern boundaries of the Almora Klippe 26,27 (supplementary Fig. 1). The Almora Klippe on which the Higher Himalayan crystalline thrusts overlayed on the Lesser Himalayan metamorphic rocks created complex strain patterns 24,25 . The south-dipping NAT and the north-dipping SAT between the MBT and the MCT oppose www.nature.com/scientificreports/ the northward movement and make complex deformation of the Kilppe. Low seismicity is observed across the Almora klippe which indicates that the region is locked and the local thrust system also adjusts accordingly.
Surface strain-rate analysis and the classification of different strain-rate zones. The existence of varied surface accommodation rates in the study area poised us to look into the regional strain-rate pattern, seismicity, and the available fault plane solutions. Crustal strain rates in the study region were calculated from the GPS estimated crustal velocities based on the modified least square approach as explained by Teza et al. 34,35 . Details of the strain-rate calculations are given in the Data and Methodology section. Analysis shows the existence of multifaceted strain-rate zones, where local and regional thrust systems have pivotal roles in regulating the overall regional seismicity distribution. Figure    www.nature.com/scientificreports/ In the HCZ, higher strain rates are observed from the north of TT and the NAT. The HCZ in the Garhwal-Kumaun Himalaya witnessed Light to Strong Magnitude earthquakes (supplementary Fig. 2a) and are mainly dominated by thrust fault mechanisms with their fault plane dipping towards NNE-SSW 36,37 . The presence of steeply dipping (> 16°) 3038 MCR plausibly enhances the slip of overburden over the MHT and thereby release the lithostatic stress as frequent earthquakes. The observed surface strain-rate vectors increases from the north of Ton thrust in the inner Lesser Himalaya, where the underlying MHT dips steeply right below the physiographic transition (PT2). Moreover, we had seen relatively larger crustal accommodation rate across the Kumaun Himalaya compared to the Garhwal Himalaya. Varied crustal accommodation rates are proxies that represent variable dip slips along the MHT which is reflected in the observed asymmetric surface strain vectors. Relatively smaller crustal accommodation is taking place in the HCZ of the Garhwal Himalaya through intermittent small to strong magnitude earthquakes owing to the relatively small seismic slips over a gentle MCR. Hence the accumulated elastic strain energy would be insufficient for a relatively larger earthquake between its Lesser and the Higher Himalaya. However, in the Kumaun region the seismicity rate is relatively less compared to the Garhwal region. Particularly at its inner Lesser Himalaya, where the reduced brittility of the rocks from brittle to semi-brittle enhances greater crustal accommodation which is adjusted as aseismic slips over a relatively steeply dipping MCR.
In general, at the south of HCZ the PT2 demarcates the southern boundary of earthquakes and the temperature at the subsurface of MHT (~ 300-350 °C) corresponds to the PT2 also favors the transition of rocks from brittle to semi-brittle 28 . Besides, there is also a reduction in effective rock stress aided with the increase of subsurface pore fluid pressure 28,29 . These processes enhance earthquake activities in this HCZ and hence validate higher crustal strain-rate vectors. We also observed the Chamoli earthquake cluster in the HCZ (encircled by green colour) with a maximum strain-rate of ~ 150 nstrain/a, and the same will be discussed separately in the subsequent section. Beyond the Chamoli cluster, high compressional strain is continuing towards Badrinath in the NNE direction with a maximum strain-rate of ~ 200 nstrain/a. High compressional strain (~ 150 nstrain/a) is visible throughout the MCT zone and shows good agreement with the seismic activities.
From Fig. 3 it is evident that the direction of compressional strain is in general towards NNE, which is same as that of the present-day movement of the Indian Plate. But it is interesting to note the existence of an extensional deformation zone (EDZ, ~ 100 nstrain/a) in the frontal part of the Kumaun Himalaya where many active faults like Dhikala, Pawalgarh 39 are present along with the NE-SW oriented transverse Moradabad fault. Trench excavation studies suggest that these frontal active fault systems are responding to the strain transfer process through long-term aseismic deformation 40,41 . The southward orientation of extensional strain vectors from the Higher Himalaya towards the frontal part also elucidates that the strain energy is getting transferred towards the frontal active fault systems and subject the faults to aseismic deformation and their southward propagation. In this process the detachment fault MHT acts as a conduit for the strain energy transfer to the frontal Himalaya. However, there is a possibility that the transferred elastic energy might be acting beyond the HFT towards the south and contributing to the local seismicity around the National Capital Region (NCR) of Delhi.
In the Equal Strain-rate zone (ESZ, ~ 100 nstrain/a), the compressional and extensional strain-rates are acting on the rock mass more or less equally. This zone is mainly seen throughout the eastern part of the Kali river near Nepal Himalaya. In general, the extension or along the arc movements are mainly seen in the frontal part of the Garhwal-Kumaun Himalaya. However, in the ESZ case, the extensional phase is present not only all along the HFT but also towards STDS with an equal amount of compressional strain component. This shows strain partitioning is predominant in this region and can be linked with larger across the arc convergence rate as well as along the arc oblique movement at the east of river Kali and adjacent to the Nepal Himalaya.
Interestingly, there are two corridors of low strain-rate zones (LSZ) visible in the vector strain map (Fig. 3). The LSZ (< 20 nstrain/a) is also clearly observable from the second invariant of strain as shown in Fig. 4a. These corridors include one near the Ramganga reservoir in the Sub-Himalaya and extends up to Baijro in the outer Lesser Himalaya and the other is at the Nainital-Almora Region. From Figs. 3 and 4a it is evident that the Tons thrust and the north Almora thrusts which separate the outer and the inner Lesser Himalaya act as a structural barrier between the high compressional zone and the extreme low strain-rate zones. Here the strain-rate increases from ~ 100 nstrain/a to a maximum value of ~ 200 nstrain/a towards Badrinath region at the north of MCT. Feeble strain rates in the Sub and the outer Lesser Himalaya are plausible because of two reasons. (1)There may be no net force acting on the frontal fault systems of the region or (2) High frictional force is acting on the locked detachment fault, which might be opposing the surface deformation. The gently dipping MHT and the Moho with greater frictional force in the Sub and the outer Lesser Himalaya obstruct seismic activity and hence the zone experiences lesser seismic deformation (supplementary Figs. 1b and 2b). In the Ramganga-Baijro window, the nature of strain starts changing from Baijro in the proximity of NAT towards further north, and ultimately it converges towards the Chamoli region. Whereas in the Nainital-Almora corridor, the low strain zone extends up to the south of NAT.
The intercomparison of geodetic strain-rate with seismicity as reported 42 earlier suggests that low strain rate zones present within the locking width are more likely to produce a great earthquake. Figure 4b shows the interseismic geodetic strain-rate plotted against the earthquake size in our study area. Most of the earthquakes in the study region are occurred in the high strain rate zones (> 100 nstrain/a) including the two significant earthquakes, the 1999 Mw 6.6 Chamoli and the 1991 Mw 6.8 Uttarkashi earthquakes 8 . Figure 4b thus supports that the high compressional zone which lies in the Himalayan seismic belt close to the MCT, shows Moderate to Strong Magnitude earthquakes having recurrence intervals of ~ 5 to 7 years (supplementary Fig. 2a). Thus the region away from the locked zone and towards the Himlayan Seismic Belt (HSB) the stored elastic strain energy is sequentially draining and thereby weakening the occurence of large earthquakes as long as the region maintains its current level of seismicity. www.nature.com/scientificreports/ The nearly NE-SW oriented low strain-rate corridors observed in the Lesser and the Sub-Himalaya (Fig. 4a) are lying over the gently dipping MHT and the subsurface duplex structures within the locked portion. This region lies south of the PT2 is considered highly locked and it shows strong interseismic coupling 14,28 . Figure 4b shows that, the region with low strain rates are not witnessed any significant events (Mw > 5) in the last 50-60 years (Also see supplimentary Fig. 2b). Paleo seismological studies also suggest that there are no transverse rupture zones of any great historical earthquakes 43 in the outer Lesser and the Sub-Himalayan sections within the LSZ. Had it been, then there would have dissipation of strain energy as seismic or aseismic slips. Besides, the lowgrade metasedimentary rocks present in the Outer Lesser Himalaya are relatively brittle and could accommodate little elastic strain energy. But the heightened concern is about the possibility of a Major or a Great earthquake as there are evidences of many such earthquakes reported from the frontal part within the locked portion of the great Himalayan arc 4,40,41 . The estimated recurrence interval of such a Mw 7.0 event in the Himalayan region is about 100 years 4 . Thus when the region exceeds its elastic limit then a large earthquake can produce a large rupture extending towards the HFT and the Gangetic plain, where secondary seismic effects like liquefaction 10 or cavity formation can cause great havoc in the densely populated foothills. Thus the regions marked as interseismic low strain rate corridors in Figs. 3 and 4a have the potential for a large earthquake. The estimated surface locking width of Ramganga-Baijro and Nainital-Almora LSZ is ~72 and ~75 km respectively from the HFT. However considering the possible recurrence of 1905 Kangra earthquake, then the Sub-and the Frontal Himalaya are on long over due for such an earthquake. This implies there are strain diffusing factors that could put off or create uncertainties on the occurrences of such events in these regions due to the existence of EDZ.
The second invariant of strain (Fig. 4a) as calculated from the principal strain components also identify the other three strain-rate zones apart from the LSZ. These estimated strain-rates are well-matched with the strain  44 . The second invariant clearly shows the nature of the extensional deformational zone (EDZ) which appears sandwiched between the earlier mentioned NE-SW orienting low strain-rate corridors (< 20 nstrain/a) and predominantly seen in the Sub and the Frontal Himalaya. The change in the surface area also supports extensional zone having positive eigenvalue contours (6 × 10 -08 /a) near the MF in the EDZ, while compressional zones have negative eigenvalue contours near the MCT (Supply. Figure 3). The EDZ is the strain energy diffusion zone, where extensional deformation happens through the transfer of strain energy from the high strain-rate region in the Lesser Himalaya and extends towards the frontal Himalaya. Greater dissipation of strain energy takes place at the frontal part where the active fault systems undergo along the arc aseismic deformation. Interestingly, transverse fault like the Moradabad fault (MF) that is connected with the frontal EDZ plausibly getting a share of the elastic strain energy at the rate of ~ 100 nstrain/a; which in turn get released as Minor to Moderate magnitude earthquakes in the NCR of Delhi. Gaur (1993) 45 explained how such frontal transverse faults segement the overthrusting Himalayan arc and facilitate strain relaxation through staggered slips over a period of time. Thus the significance of sub and frontal EDZ in deferring the much anticipated large earthquakes in the locked portion of the region cannot be ruled out.
Correlation between strain-rate and the topography. The correlation of strain-rate profile with topography would explain the mechanism of strain-rate transfer from the Higher to the Frontal Himalaya. As we traverse north across the thrust systems in the Garhwal-Kumaun Himalaya the topography increases, hence the gravitational potential energy 29 . Thus a positive correlation between the topography and the strain-rate is expected. Similarly, the subsurface duplex systems in the Lesser Himalaya are also influencing the seismicity and the deformation characteristics in the region 28 . Thus we attempt to understand the mechanism of strain transfer towards the frontal, by unraveling (1) how the strain-rate changes as one moves across-the-thrust systems? (2) the role of subsurface Lesser Himalayan duplex system on across-the-thrust strain-rate and (3) how the overburden topography affects strain-rate distribution in terms of Gravitational potential energy change (GPE)? Three strain-rate profiles are taken across the thrust systems in the LSZ and the EDZ as shown in Fig. 4a. Figure 5a,b represent the strain-rate profiles across the Nainital-Almora and the Ramganga-Baijro LSZ corridors respectively; and Fig. 5c shows the strain-rate across the Ramnagar-Badrinath profile. We observe a mismatch between the strain-rate and the topography in both profiles (a) and (b); particularly, in the Lesser Himalaya where the strain-rate changes from low to high values and seen out-of-phase with the topography. However, in the Higher Himalaya reasonable correlation exist between the strain-rate and the first order topography and seismicity. While moving towards the mountain front along the profile (a) the topography falls from PT2 at the south of MCT; but the strain-rate still stands high at PT2 and it falls only at NAT. Whereas, for the profile (b) the high strain-rate further continuing up to the SAT. As we observed that the topography is not well correlated with the strain-rate changes in the Lesser Himalayan part of the profiles (a) and (b). Hence the GPE aided strain transfer from the Higher to the frontal Himalaya is not prominent in these sections. Thus instead of a single mid-crustal sub-surface ramp, an alternative model involving Lesser Himalayan sub-surface duplex system has been considered to explain the strain transfer processes. So that, the duplex systems are significant in keeping the strain rates higher in the inner Lesser Himalaya even in the absence of seismicity, but efficient to hold and transfer the strain. Thus from Fig. 5a,b it is evident that the interseismic strain adjustments are mostly occurring in the Lesser Himalayan region and marked as aseismic deformation zones in the profiles and the same is observed at the south of PT2, where the strain energy transfer is happening at different time scales. It is also observable from Fig. 4a that 100 nstrain/a is the threshold, from where the strain-rate trend falls towards the frontal part. The TT and the NAT situated at this strain-rate threshold act as structural barriers. The NAT regulates the transfer of strain energy from the HCZ towards the Lesser and the Frontal Himalaya along the low strain-rate corridors. Thus the local thrust systems like the TT, NAT, and the SAT and their aseismic deformation modulate the overall strain partitioning in the region.
It is also inferred that, if there were no local thrust systems like the NAT or the SAT and the connected duplex systems, then the strain-rate could have been correlatable with the mean topography. This is observed for the case of profile (c) that extends from the EDZ in the frontal part to the HCZ towards the north. We can observe a relatively better correlation between the mean topography and the strain-rate in the Ramnagar-Badrinath profile. Here gradual decrease in the strain rates from the Higher Himalaya towards the extensional frontal Himalaya can be observed, therefore indicating relatively lesser influences of NAT and the Lesser Himalayan duplex system compared to the profiles (a) and (b) in the low strain-rate corridors. Mean topography and the long-wavelength GPE are highly correlatable here and the GPE aided strain transfer towards the frontal EDZ is the dominant mechanism of strain energy transfer in this section.
Curvilinear strain-rate pattern and the Chamoli earthquake cluster. Chamoli region in the Garhwal Higher Himalaya is situated in the HCZ (encircled by green colour in Fig. 3) and has been known for its high seismicity rate. More than 500 shallow earthquakes are recorded above the MHT depth in this region after the significant Mw 6.6, 1999 Chamoli earthquake 46 30,47,49 contributes to the enhanced seismicity. Fault plane solutions show strike-slip earthquakes close to the north of the MCT zone; however, thrust solutions are also seen between PT2 and the MCT zone. Although, the presence of sub-surface fluid layer is a sufficient reason to explain the frequent occurrences of small or even minor earthquakes, but not enough for the generation of a great or even significant magnitude earthquake like the Chamoli www.nature.com/scientificreports/ earthquake of 1999. Thus it is imperative to know how the strain adjustment of local and other thrust systems around the Chamoli cluster could contribute to larger earthquakes and clustered seismicity. Strain-rate analysis in and around the Chamoli region shows an interesting curvilinear pattern (Encircled by violet dots in Fig. 3). Here the deformation field starts from Baijro in the Inner Lesser Himalaya at the western edge of Almora Klippe, where the TT terminates. The NE trending compressional strain vectors surpass the NAT and converge towards the Chamoli earthquake cluster. As mentioned earlier, in terms of crustal scale, the overall accommodation rate between the Lesser and the Higher Himalayas is asymmetric in the Garhwal and the Kumaun region, which is 7.38 ± 0.24 mm/a and 5.41 ± 0.55 mm/a respectively.
Apart from the crustal-scale accommodation process of MHT, the contribution of local thrust systems toward modulating the Chamoli strain pattern also has to be considered. The linear crustal accommodation rate across the TT in the Garhwal Himalaya is 1.30 ± 1.22 mm/a (Table 2), while the same is 2.12 ± 0.55 mm/a across the Almora Klippe at the eastern Kumaun region. Similarly, across the NAT and the SAT the accommodation rates are 0.68 ± 0.65 mm/a and 1.44 ± 0.48 mm/a respectively. However, larger linear crustal accommodation rate (2.64 ± 1.66 mm/a) is taking place at the western edge of the Almora Klippe between Baijro and Chamoli. Extended level of seismicity from MCT towards Baijro and greater deformation in this region suggest that the local thrusts (NAT and SAT) are asymptotically connected to the gently dipping MHT as compared to the eastern part of the Almora Klippe 30,50 (Fig. 6). This observed asymmetry in the crustal accommodation processes over www.nature.com/scientificreports/ a gently and steeply dipping ramps at the northwest and the southeast edges of the Almora Klippe, respectively; creates a clockwise (NE orienting) rotational couple and variable slips of the upper crust over the MHT (Fig. 6). This rotational couple creates a curvilinear surface strain pattern with their point of convergence is directed towards the Chamoli region. The subsequent release of strain energy owing to the rotational block movement could be the plausible reason for a Mw 6.6 Chamoli type earthquake and the sub-surface fluidity might have aided the rotational slip and the consequent clustered seismicity.

Conclusion
Studies on the inter-seismic crustal deformation and strain-rate pattern in the Garhwal-Kumaun Himalaya, in the context of an anticipated 2015 Gorkha-Nepal like Major magnitude earthquake, reveal four types of strainrate zones; namely, (1) High compressional zone (HCZ), (2) Extensional deformation zone (EDZ), (3) Equal strain-rate zone (ESZ) and (4) Low strain-rate zones (LSZ). The HCZ (~ 150 nstrain/a) corresponds to the zone of active seismicity close to the MCT. A curvilinear strain-rate pattern is observed in the Chamoli cluster which is attributed to the asymmetry in the crustal accommodation processes over a gently and steeply dipping ramps at the northwest and the southeast edges of the Almora Klippe, respectively. This has created a clockwise rotational couple on the upper crust moving over the MHT. The variable slips of the NAT and the SAT regulate the shallow crustal accommodation processes and sustained the rotational slip aided by sub-surface fluids caused the Chamoli earthquake cluster. We identified two NE-SW corridors of LSZ (~ 20 nstrain/a); namely, Ramganga-Baijro and the Nainital-Almora in the Lesser and the Sub-Himalaya, where the surface locking width is ~72 and ~75 km respectively. This highly locked LSZ does not show any surface rupture and had not witnessed any great earthquakes, but the upper crust is strained and has the potential for large earthquakes. However, strain reducing EDZ is seen sandwiched between the LSZ; particularly, in the frontal Himalaya where the GPE aided transferred strain from the HCZ zones gets diffused by the deformation of frontal active fault systems and transverse structures like the Moradabad fault. This causes more ambiguity towards the possible occurrence of large Sub-and Frontal earthquakes in the Garhwal-Kumaun region. The TT and the NAT thrust systems, stand at the threshold of 100nstrain/a, act as structural barriers in the strain transfer processes from the MCT zone to the southern thrust systems. www.nature.com/scientificreports/ Data and methodology. GPS data from 42 stations (details are in Table 1) including data from WIHG local network and IGS (International GNSS Service) for a span of eight years from 2010 to 2017 were used in this study. The data were processed using the GAMIT/GLOBK software version 10.6 51,52 in the ITRF 2008 reference frame and the position coordinates were estimated by stabilizing surrounding IGS stations. Ocean loading effects of the sites in the GPS data were corrected using the ocean tide model FES2004 and the corrections in the displacements associated with the earth's solid tides were also removed by applying the tidal model IERS03.
Atmospheric loading effects and the tropospheric corrections at sites were removed using the atmospheric loading data and the global mapping function model respectively. Generated long time series of each station coordinates are corrected for outliers and velocities were estimated in the International Terrestrial Reference Frame (ITRF) ( Table 1). We also used the data from published works 14,33,53-56 along with our data. All the published data are converted to the ITRF 2008 reference frame and estimated the velocities in the Indian reference frame using the pole of rotation of Ader 32 . The linear surface shortening rates across each fault are calculated from the horizontal velocity residuals. Residual velocities of stations situating north and south of each fault zone are calculated from the Indian reference frame velocities 32 . Crustal strain rates in the Garhwal-Kumaun region were calculated from the GPS estimated crustal velocities based on the modified least square approach explained by Teza et al. 34,35 . We have excluded the stations having velocity uncertainties of more than 2σ from the surface strain estimation. We choose a local grid of spacing 7.5 km and defined the scale factor (15 km). Local strain rates at each grid node are estimated using the modified least square approach and a weight function is used for the error reduction. The principal strain components are estimated and the resulting eigenvalues represent the strain values at each grid node. The second invariant of horizontal strain rate is calculated using the equation where e1 and e2 are the principle strain components. The second invariant of the strain rate represents the total strain as a scalar quantity 44,57 .
The seismicity data of the study region were collected from the ISC catalog and converted to the moment magnitude (Mw) using the regression relations 20 . The earthquake data of the study region were classified into two classes; earthquakes above and below 100 nstrain/a zones. The earthquakes in the both classes are subdivided according to their size (Mw) and the recurrence interval for each class of earthquakes are calculated (Supply. Figure 2a,b). The earthquake and the strain invariant data sets of the study region are brought in to a common grid frame. Then for each grid, we identified the different class of earthquakes present and its corresponding second invariant of strain rates (Fig. 4b).