Earthquake Shocks Around Delhi-NCR and the Adjoining Himalayan Front: A Seismotectonic Perspective

An increase in the number of earthquakes and subsequent clustering in northwest India, particularly around the Delhi-National Capital Region (NCR) and adjacent NW Himalayan front, provides a good opportunity to understand the underpinning tectonic controls and the likelihood of any large earthquake in the future. The 2001 Mw 7.7 Bhuj, 2011 Mw 6.9 Sikkim and 2015 Mw 7.8 and 7.3 Nepal earthquakes (and 2004 Mw 9.2 Sumatra event) are important in this context. We analyzed the seismicity around the Delhi-NCR and the adjoining Himalayan front, including event clustering and the spatio-temporal distribution of b-values, in the context of kinematics and the regional geodynamics. The overall moderate-to-low b-values, both in time and space, since 2016, provide information regarding an increase and subsequent stabilization of the stress field in the study area. The analysis led to the identification of (1) a structurally guided stress field in the region between the Kachchh and the NW Himalaya that coincides with the direction of Indian plate convergence and (2) frequent occurrences of earthquakes particularly in the Delhi, Kangra and Uttarkashi areas. We propose that faults in western Peninsular India, which pass through the margins of the Aravalli Range, the Marwar basin, and the isostatically over-compensated Indo-Gangetic Plains beneath the under-plated Indian lithosphere, act as stress guides; concentrating and increasing stress in regions of lithospheric flexure. This enhanced stress may trigger a large earthquake.


INTRODUCTION
An increase in clustering of seismic shocks of magnitudes (M L ) ranging from 2.0 to 5.2 around the Delhi-NCR and in the adjacent region surrounding large historical damaging events ( Figures 1A,B), has raised alarm in the society and scientific community. The seismicity record between January 2016 and June 2020 includes 414 events of M L ≥ 2.0, of which 135 shocks (M L ≥ 2.0) took place in the first six months of 2020 in the Delhi-NCR and adjacent Himalayan front. The fallout news coverage, claims and counter claims in electronic and print media raised some pertinent scientific questions, e.g., (1) possibility of any great event in the Delhi and adjoining areas in the future, (2) triggering of the recent shocks, and (3) possible role of pre-Himalayan heterogeneities behind the along-strike seismic behaviour of the Himalayan front including the Delhi-NCR region. It is interesting to note that majority of the shocks are concentrated around the regions of historical seismicity in the northern Indian plains, including the Delhi-NCR. Indeed, the catastrophic events of the 2001 M w 7.7 Bhuj ( Figure 1C), 2011 M w 6.9 Sikkim and the 2015 M w 7.8 Nepal ( Figure 1A) earthquakes are still alive in the minds of population. The faultbound 2001 Bhuj earthquake occurred at the shallow-level of the lower crust caused widespread damage with a toll of ∼20,000 human lives. Small-to-moderate magnitude aftershocks were recorded in different parts of northern Gujarat till the end of 2014 Aggarwal et al., 2016). In addition, the incidences of historical major earthquakes (1819 M w 7.8 Allahbund and 1956 M w 6.0 Anjar, Figure 1C) in the Bhuj region had rattled the basic tenet of intraplate origin of damaging earthquakes (Chung and Gao, 1995). Copley et al. (2011) linked the Bhuj event with the subduction and collisional dynamics in the western segment of the Himalayas. It is suggested that the seismogenic Indian crust, extending from Gujarat to the Himalayas, supports the compressive stress linked with the subduction of the Indian plate in the region south of Tibet, and is actively involved in triggering moderate-togreat magnitude earthquakes in these areas. In this backdrop, we have investigated the possible reasons behind the recent incidents of small-to-moderate magnitude earthquakes around the Delhi-NCR and the adjoining Himalayas.
The ∼2500 km long Himalayan orogeny is structurally and tectonically segmented by the trench-orthogonal ridges (Delhi-Hardwar, Faizabad and Monghyr-Saharsa) and depressions (Sarada and Gandak) in the underthrusted Indian plate (Valdiya, 1976;Gahalaut and Kundu, 2012;Hetényi et al., 2016). The segmented orogeny accommodates stress, which triggered several damaging earthquakes (e.g. 1803M 7.7 Garhwal, 1833M 7.6 Nepal-Bihar, 1897M 8.1 Shillong, 1905M 7.8 Kangra, 1945M 6.3 Chamba, 1950M 8.6 Assam, 1975M 6.8 Kinnaur, 1988M 6.8 North Bihar, 1991M 6.8 Uttarkashi, 1999M 6.5 Chamoli, 2011M 6.9 Sikkim, 2015M 7.8 Nepal and 2015. The Delhi-NCR region, close to the foothills of the Himalayas, also experienced a few damaging earthquakes (e.g. 1720 M 6.5, 1956 M 6.7, 1960 M 6.0) in the past (Figures 1-3). Our aim is to document the trends of concentrated seismic activities that have been occurring since 2016, and identify the causative factors, selecting the area encompassing the Delhi-NCR and the adjoining mountain front in western India ( Figure 1A). We attempt our study with the analysis of seismic b-values in space-time domains. Seismic b-value has been used in a variety of seismological studies, including earthquake forecasting, prediction, and seismic hazard analysis (Kagan and Knopoff, 1987;Rydelek and Sacks, 1989;Geller, 1997;Wiemer and Wyss, 2000;Beauval and Scotti, 2004;Woessner and Wiemer, 2005;Ogata and Zhuang 2006;Khan and Chakraborty, 2007;Felzer, 2008;Khan et al., 2011;Sobolev, 2011;Cheng and Sun, 2018;Peresan and Gentili, 2018). Special emphasis was given to simulate the distribution and the confinement of seismicity around the Delhi-NCR, the adjoining Himalayas, in general, and for three main tectonic domains of interest, i.e., Delhi (Zone X), Kangra (Zone Y) and Uttarkashi (Zone Z) regions, in particular ( Figure 1B). We finally integrate all the observations into a holistic seismic model that may account for interpretation of signatures obtained from structurally complicated western Indian plain and the adjoining Himalayan front.

REGIONAL TECTONIC SETUP
The Late Mesozoic subduction of the Neotethyan oceanic crust at the leading edge of the Indian Plate below the Eurasian plate and the southward migration of the convergence zone in the Cenozoic gave rise to the evolution of the Himalayas (Patriat and Achache, 1984;Copley et al., 2010). The suturing migrated eastward along the Indus-Tsangpo zone, with concomitant oroclinal bending of the Himalayas (Klootwijck et al., 1985). The Paleozoic-Mesozoic Tethyan sedimentary cover on the subducting Indian plate got mobilized in the Himalayan orogenic front. The deformation, thus initiated, migrated southwards in the form of a series of crustal-scale thrust systems (Chakraborty et al., 2019). From north to south, these include the Main Central Thrust (MCT, between Tehyan Himalayas and Lesser Himalayas), Main Boundary Thrust (MBT, between Lesser Himalayas and Sub-Himalayas) and Main Frontal Thrust (MFT, between Sub-Himalayas and Indo-Gangetic Plains) (Le Fort, 1986). It is presumed that the Main Frontal Thrust (MFT) accommodates the maximum deformation of the southward verging orogeny (Lavé and Avouac, 2000). The MFT, MBT, and MCT are merged in a northward dipping dećollement, i.e., the Main Himalayan Thrust (MHT) (Makovsky et al., 1996;Hauck et al., 1998;Bollinger et al., 2006), which is modeled as a duplex with two north-dipping sub-horizontal planes connected by a system of thrust faults (Mendoza et al., 2019), act as an impediment to the plate convergence.
The convergence rate of the Indian plate near the Himalayan front varies significantly between ∼5.4 cm/yr and 4.2 cm/yr along the east-west direction (DeMets et al., 1994). The hard collision of the Indian plate during the Late Miocene caused the maximum upliftment of the Himalayas and the Tibet (Royer and Gordon, 1997). The enormity of deformation of the orogenic process has resulted in the development of unconformity in the floor of the Bay of Bengal, folding and faulting in the equatorial Indian Ocean, and the evolution of basins and ridges near the foothills of the Himalayas-all in response to tectonic processes in the orogen that evolved during Miocene-Pliocene time (Cochran, 1990;Curray, 2005;Gordon, 2009). Thus, the active resistance to the under thrusting distributes and perturbs the stress fields up to the Indian Ocean, and accounts for occasional incidents of great intraplate earthquakes Aggarwal et al., 2016;Khan et al., 2016;Khan et al., 2020). Further, the excess weight of the orogeny is supported in parts by the strength of the Indian plate, which distributes the load by flexing down along the ∼1000 km laterally extended and isostatically over-compensated Indo-Gangetic Plains (e.g., Lyon-Caen and Molnar, 1983;Tandon et al., 2014).
Experiments on deformation and numerical modeling involving tectonic units along collisional boundaries have revealed a complex geometry for sheared suture zone (Beaumont et al., 1996;Selzer et al., 2008). Continental plates, involving collision, get stacked and show the presence of double Moho below the overriding plate (Beaumont et al., 1996;Selzer et al., 2008). In these models, the curved geometry of the crustal elements in plan and cross-section act as the sites for stress accumulation and deformation leading to seismicity. The deformation scenario becomes more complex when the subducting plate inherits tectonic elements at high angles to the regional trend of the mountain belt. Converging plate, containing pre-orogenic anisotropies (e.g., ridge, fault systems, lineaments, etc.), leads to lateral variation of crustal thickness during subduction. The processes may further lead to segmentation of the mountain belt having variable stress-strain patterns Hetényi et al., 2016 (Valdiya, 1976;Karunakaran and Ranga Rao, 1979;Valdiya, 2003;Goswami, 2012;Godin and Harris, 2014). Many regional faults and lineaments, viz. the SONA Fault, the Great Boundary Fault (GBF), Bachau-Jodhpur-Saddarshahar-Chandigarh lineament (BJSC), Kori Creek-Jaisalmer-Suratgarh-Ropar lineament (KJSR) and Sukkur-Bahawalpur-Lahore-Kathua lineament (SBLK) parallel to these ridge systems ( Figure 1A) define large-scale horsts, grabens and sedimentary basins with crystalline basement and sedimentary fills of variable thickness. Understandably, such heterogeneities cause segmentation in the subducting Indian Plate (Figure 1).

EARTHQUAKE DATA AND ANALYSIS
The present analysis has been carried out using the earthquake data recorded by local network of the National Centre of Seismology (NCS), Government of India, New Delhi (Srivastav et al., 2005;Dattatrayam et al., 2014;Prakash et al., 2020). The seismic network in Delhi NCR and the surrounding region has more than 30 broadband stations ( Figure 3 of Prakash et al., 2020), and the earthquake catalogue contains information on origin time/date, latitude and longitude of the epicentre, focal depth and magnitude. The network has a detection and location capabilities down to a magnitude 2.0 (M L ) earthquake. Also, the earthquake data in the catalogue have an accuracy of ±2 km in the epicentre. Despite upgradation of basic seismic network in recent time with addition of state-of-the-art broadband sensor having VSAT connectivity (Prakash et al., 2020), the accuracy in focal depth determination is not improved further because of limitation in recording station distribution. This has restricted any analysis based on hypocentre-information available in the Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 598784 catalogue. Taking into consideration the accuracy of locations and magnitudes of the events available in the NCS catalogue, a total of 414 events of M L ≥ 2.0 have been considered ( Figure 4B) for detailed analysis of the earthquakes in the present study. Further, regional clustering of events (M L ≥ 4.0) has also been investigated in view of the major damaging earthquakes in different segments of the Himalayas and the Delhi region ( Figure 2). We have compiled 441 such events for the duration 2016-19 from the catalogue of International Seismological Centre (ISC), and 84 events for the first 6 months of 2020 from NCS catalogue based on the coverage of recording stations.

Regional Seismicity Analysis
The elliptical segment of the Nepal-Bihar-Sikkim, a zone of interaction of MSR with the foothills of the Himalayas (shown in Figure 1A), covers the region of past six major earthquakes. Along with the MSR, many regional transcurrent faults (e.g., east and west Patna faults, Pingla fault and many others, Godin and Harris, 2014;Hetényi et al., 2016) passing through this region, have been accommodating the plate movements since the Holocene times (Valdiya, 1976, Valdiya, 2003. These are particularly involved in concentrating and raising the stress field (Sibson, 1980;Marshak and Paulsen, 1997), and occasionally trigger major earthquakes in this part of the Himalayas (Dasgupta et al., 2000;Paul et al., 2015;Singh et al., 2020). The spatially clustered seismic activities in this area ( Figure 1A) are caused by occasional release of strain energy through occurrences of aftershocks and other main shocks subsequent to 2015 M w 7.8 and 7.3 Nepal events. The great 2001 M w 7.7 Bhuj earthquake occurred along the Kachchh Mainland Fault (KMF), where another two major earthquakes occurred in 1819 and 1956 (block B, Figures  1A,C), and the aftershocks continued for a period of longtime (Aggarwal et al., 2016;Khan et al., 2017). The strain energy, released in form of earthquakes, apparently has migrated toward the western part, particularly in the northwest Himalaya and adjoining regions. Concentrated activities are noticed in the region with epicentres of major historical earthquake events near Kangra (no. 4, 5) and Delhi (no. 16-18) ( Figure 1A).
The regional seismicity has been investigated in some details over four time-domains, viz., 2016 ( Figure 1B displays the distribution of seismicity (M L ≥ 2.0) for the duration January 01, 2016 to June 30, 2020 (414 events, altogether) in the Delhi and adjoining mountain front ( Figures 4A, 5A). Seismicity is found to concentrate within specific areas, namely surrounding the Delhi (120 events, Zone X), Kangra (151 events, Zone Y) and Uttarkashi (87 events, Zone Z) regions (Figures 1, 4, 5). Spatial zone-specific distribution of seismicity, within the Delhi and adjoining Himalayas region, has been further investigated for different time intervals 2016 (59 events, Figure 3A), 2017 (89 events, Figure 3B), 2018 (57 events, Figure 3C), and 2019-2020 (209 events, Figure 3D). The concentration of seismicity in these domains has changed considerably with maximum of 209 in 2019-20; and this number has increased more sharply in 2020 (134 events over first six months). Significant number of events with magnitudes M L ≥ 4.0 has also been recorded in each time-domain. In particular, the Zone Y surrounding the Kangra area records significant number of events with magnitude M L ≥ 4.0 ( Figure 4). Although, low-magnitude earthquakes are recorded in high numbers near the Delhi (Zone X) and Uttarkashi area (Zone Z), events with M L ≥ 4.0 are not rare in these areas. Figure 5A shows that the seismic activities are distinctly divided into two time-domains; one over 2016-18 and another over 2019-20, and apparently the energy release is explained to be discontinuous in nature. More or less similar oscillations are also found in Zones X and Y ( Figures 5B,C). A distinct cluster of seismic activities with two events of magnitudes M L ≥ 5.0 was associated with the Kangra area since 2019 ( Figure 5C). However, such distinct signature of clustering is apparently less prominent in Zone Z ( Figure 5D).

Seismic b-Value Analysis
Earthquake magnitude-frequency domain is approximated by a power-law distribution as observed for other similar selforganized non-equilibrium systems. Gutenberg and Richter (1954) proposed an empirical relationship (the GR law) between frequencies of occurrence of earthquakes as a function of magnitude, which is expressed mathematically as: where N(M) is the cumulative number of shocks with magnitude M ≥ 0, and 'a' is a measure of the background level of seismic activity/efficiency of the concerned area. It depends on the areal extent of the study-region and duration. The seismic b-value is computed from the slope of the Log 10 (N) vs. M regression line. It accounts for largemagnitude earthquakes over the number of smallmagnitude shocks (Khan and Chakraborty, 2007;Cheng and Sun, 2018), and is inversely related with the differential stress-level (σ 1 −σ 3 ) and/or shear strength of the materials (Wyss, 1973;Urbancic et al., 1992;Schorlemmer et al., 2005;Spada et al., 2013;Scholz, 2015). It also determines the decay in the rate of occurrences of earthquake events with increasing magnitude. In natural situations, b-values are found to lie between 0.5 and 1.5, depending on the tectonics of a region (Khan et al., 2011). The temporal and spatial heterogeneity in crust may break the power-law scaling for distribution of earthquakes (Rydelek and Sacks, 2003;Woessner and Wiemer, 2005), and result in anomalous b-value distribution owing to anomalous stress distribution at high and low levels. However, the selforganized and self-similar non-equilibrium condition in the subsurface medium, may still remain valid for earthquake generations (Rydelek and Sacks, 1989;Woessner and Wiemer, 2005;Sobolev, 2011), and stable for eqn. 1 over a definite timewindow (Wiemer and Wyss, 2000). The plot systematically deviates from linearity because of rare incidents of large magnitude events at upper bound, and is constrained by the incompleteness of catalogue at lower bound (Rydelek and Sacks, 1989;Taylor et al., 1990;Wesnousky, 1994;Naylor et al., 2010). The latter may be caused by missing a large number of small magnitude earthquakes due to: detection threshold; aftershock (or other shocks) occurrence within the tail of larger events; aseismic slip during nucleation along a creeping fault (Heimpel and Malin, 1998;Vorobieva et al., 2016), or lack of enough stations to record very small magnitude events (Pacheco et al., 1992). Wiemer and Wyss (2000) and Woessner and Wiemer (2005) have examined such limitations, and advocated the essential requirement of minimum magnitudes of completeness (Mc) for seismicity analysis. Although, the Mc varies with time in most of the catalogues for several regions, namely it decreases with increase in the number of recording stations and advancement in the methodology of analysis, it is widely used for b-value estimation over high seismic areas (Rydelek and Sacks, 1989;Wiemer and Wyss, 2000). Further, the GR law needs to be estimated for different segments over the study area in view of self-similar seismotectonic character (Molchan et al., 1997), or distinct locations of geological units (Peresan and Gentili, 2018). Thus, the zonation in the region of present concern is constrained mainly by the self-similar seismotectonic behaviour, and a specific Mc value is determined for each zone.
With all discussed constraints, we carried out the analysis of seismic b-values following frequency-magnitude distribution (FMD, Eq. 1) of Gutenberg and Richter (1954) for understanding the possibilities of major earthquakes in any of the three zones (X, Y and Z, Figure 1A). The technique is being used ubiquitously in different tectonic settings at different levels of depths for several decades considering earthquake catalogues covering from few months to several years. This makes the GR law very effective for analysis of earthquake hazards and also attractive for its theoretical significance. Assuming the selfsimilarity in the GR law (Mignan and Woessner, 2012), we have performed the analysis for the entire study area and for the three demarcated zones in different time-domains between 2016 and June 30, 2020 (cf. Figures 3, 5A,6). On the basis of spatial clustering of earthquakes epicentres in different time windows, the b-values are clubbed under two blocks A and B ( Figure 5A). The uncertainties in the calculations of b-values are based on the model proposed by Shi and Bolt (1982), which can be written as where σ(b) represents the standard error of b-value calculation, N and M are the number and the average value of the magnitude distribution (M i ) of the events used for computation of b-values for magnitude M c and above. We found few events of larger magnitudes in each domains fall little away from the best-fit line at its upper bound. These few events have been excluded from the analysis of b-values in view of scarcity of such events at the higher magnitude range. This also has improved the value of coefficient of determination (i.e. R 2 ) in this statistical regression analysis. Selection of completeness magnitude (M c ) of a sequence of natural shocks is crucial for any statistical analysis (Mignan and Woessner, 2012 and the references therein). The selection of lower value of Mc will normally contaminate the analysis, largely by threshold characteristics of the instruments as well the detection of events, whereas the high value results in undersampling of the sequence of dataset beyond the Mc value. In the present study, attempt was made for Mc's selection for different zones with critical examination by visual inspection of departure of linearity (apparent break in the gradient) in the log-linear plot of cumulative FMD of earthquake events. However, the gradual change in the curvature of the plots made the selection of Mc difficult (Mignan and Woessner, 2012), particularly when resorting to maximum curvature method for Mc detection. This typical trend in the cumulative FMD is usually caused by the spatio-temporal heterogeneities in the Mc (Wiemer and Wyss, 2000;Mignan, 2011). Thus, to verify deviations from log-linearity of the FMD distribution, the Mc for each dataset has been computed on the basis on non-cumulative (i.e. computed log 10 value of number of events for each magnitude increment) log-linear plot of events (Woessner and Wiemer, 2005;Mignan and Woessner, 2012). The exercise was carried out to avoid the limitation in the selection of Mc threshold from the cumulative log-linear plot of FMD, and to provide a rather conservative specific Mc for each zone ( Table 1).

RESULTS
The maximum computed b-value of 1.04 is recorded in Zone Z (i.e., the Uttarkashi area; Figure 1), and the next low b-value of 1.03 is found for the entire region for the data ranging from 2016 to June 30, 2020 ( Figure 6B; Table 1). The minimum b value of 0.76 is found for Zone X (i.e., Delhi-NCR), and the next higher one of 0.78 is noted during 2018. Other b-values lie at moderateto-low-level with a maximum of 0.96 during 2019 ( Figure 6B; Table 1).
We found clustered seismicity in Delhi (Zone X, Figure 5B) and Kangra (Zone Y, Figure 5C) areas but not so strongly at the Uttarkashi area (Zone Z, Figure 5D). Maximum number of events were found in the Kangra region in 2019-20, with the close next at Delhi. Figure 6B and Table 1 Table 1). The overall two clustered seismicity in block A and B during 2016-2020 ( Figure 5A) also shows moderate b-values of 0.92 and 0.94, respectively ( Table 1).
Although we found a little high standard errors (≥0.09) for 2016, 2017, 2019, and Zone Y and Z, all the estimates are well compatible within the range of errors. The large uncertainties for some domains might have been found due to limited available observations and heterogeneities in the dataset; the b-value estimate for zone X, in particular, though being characterized by a pretty low standard error (∼0.06), it has a lower value of coefficient of determination (96%) compared to other zones. The robustness of b-values estimates was also tested with other methods (i.e., maximum likelihood, Aki, 1965) and was found to be compatible within the errors limits; also it was found that such analysis is less sensitive to magnitude upper bound, although producing lower values of coefficient of determination R 2 .
The persistent moderate (∼1.0) to low b-values (∼0.8-0.76) in all the Zones X,Y and Z, and blocks A and B (Table 1; Figure 6B) invariably account for higher stress accumulation in the entire study-area (cf. Tormann et al., 2012;Spada et al., 2013;Khan et al., 2018). The low b-values, over time and space domains, usually indicate high-strength homogeneous rock-mass, whereas high b-values are found for a large number of small earthquakes expected in regions of low strength and large heterogeneity (Tsapanos, 1990). Low b-value in a region is indicative of high stress accumulation and release of high seismic moment (Scholz, 1968;Wyss, 1973;Cao and Gao, 2002). From stable and unstable friction experiment, Rivière et al. (2018) have shown that the FIGURE 6 | Plots illustrating (A) the computation of seismic b-value (after Gutenberg and Richter, 1954) for earthquakes of 2020 and (B) the variation of seismic b-values for seismic zones (X, Y and Z) during different time domains. Index for time domains along the x-axis: 1-2016 (entire area); 2-2017 (entire area); 3-2018 (entire area); 4-2019 (entire area); 5-2020 (entire area); 6-2016-2020 (entire area); 7-Zone X; 8-Zone Y; 9-Zone Z; 10-Block A of Figure 5A; and 11-Block B of Figure 5A. Mc: Completeness magnitude, σ (b) : Standard error (after Shi and Bolt, 1982), R 2 : Coefficient of determination. reduction of b-value occurs when shear stress systematically rises before stick-slip failure, and results in large events when faults become highly stressed. Incidences of large number of foreshocks are also reported in such situation (Scholz, 1968;Goebel et al., 2013, Goebel et al., 2015. Scholz (1968) and Wyss (1973) suggested the low b-values to be compatible with release of higher seismic moment energy in an active seismic zone. Scholz (1968) also suggested that the b-value decreases with increasing differential crustal stress (Schorlemmer et al., 2005;Spada et al., 2013;Scholz, 2015;Khan et al., 2018), and the declining b-value accounts for an impending main shock. The increased number of events in 2019-20 (block B, Figure 5A) and an overall moderate-to-low b-value indicate the possibility of increasing stress accumulation in the study area and a future great earthquake.
We also found a cluster of seismic events near the regions where ridges (horsts) and basins (grabens) interact with the Himalayan orogenic front. This observation indicates stress concentrations at the intersection points (Sibson, 1980;Marshak and Paulsen, 1997;. Indeed, concentration of most of the events in the re-entrants (e.g. 1803M 7.7 Garhwal, 1905M 7.8 Kangra, 1945M 6.3 Chamba, 1975M 6.8 Kinnaur, 1991M 6.8 Uttarkashi, and 1999, and historical earthquakes along the margins of the regional faults in the Western India (e.g. 1720 M 6.5 Delhi, 1956M 6.7 Delhi, 1960M 6.0 Delhi, 1966M 5.8 Delhi-NCR, 1819M w 7.8 Allahbund, 1956M w 6.0 Anjar, 1969M w 5.3 Mt Abu, 2001M w 7.6 Bhuj, 2009 Barmer) indicate a link between the seismic events with regional faults striking NNE-SSW and cross-faults forming horsts and grabens in the subducting Indian plate and their interaction with the MHT.
We have also noticed increased frequency of earthquakes of magnitudes M L ≥ 4.0 along the western part of India, particularly in 2019-20 ( Figures 1A, 2), with many events associated with the NNE directed lineaments transcend the Himalayan mountain front in the area. The 2012 M w 5.1 Dholavira and 2020 Mw 5.3 Bhachau earthquakes with five other earthquakes of magnitudes 5.0 and above along the western part of India ( Figure 1A) are indicators for increasing stress fields in Delhi and the adjoining regions. The fluctuation and multi-fractality in the sequence of aftershocks of 2001 Bhuj earthquakes (Aggarwal et al., 2015, Aggarwal et al., 2017 bear testimony for the oscillating nature of seismic energy release, and support the chaotic nature of shocks in Delhi and the adjoining regions ( Figure 6). Further, the occurrences of earthquakes in a region follow the non-linear, chaotic and scale-invariant pattern, and the non-constraint of self-similarity in the earthquake sequence limits the predictability of large future earthquake (Sykes et al., 1999).

DISCUSSION
The analysis of regional seismicity indicates that the stress-strain patterns in the Himalayas, generating major earthquakes, are apparently linked with the regional structural patterns of the converging Indian plate. It is noted that the crustal heterogeneities are responsible for the lateral segmentation of the 2500 km long orogen into distinct tectonic domains (Valdiya, 1096;Gahalaut and Kundu, 2012;Godin and Harris, 2014;Hetényi et al., 2016). To resolve specific-controls, we have attempted tracking of regional structures in the study area to find their continuity, if any, from Indian plate interior through the foredeep to the Himalayas and the link between these structures and seismicity in the area. Transfer of stress into the Indian plate interior by the regional fault systems is documented and identified as an earthquake trigger in selected domains of the Peninsular India (Copley et al., 2011). One such zone of high seismicity is found in the western part of India, i.e. the Kachchh region. The 2001 M w 7.6 Bhuj earthquake was considered to be the result of stress propagation from the Himalayas (Copley et al., 2011) or the Suleiman Range (Stein et al., 2002;Khan et al., 2016). The analysis of Khan et al. (2016) indicates that the far-field stress generated by the collision resistance between the Indian and Eurasian plates and anticlockwise rotation of the Indian plate are the main factors behind the stress accumulation in the western India and its periodic release through earthquakes. In addition to the earthquakes in the Kachchh region, two other events (1969 M w 5.3 Mt. Abu and 2009 M w 5.1 Barmer at 15 and 39 km depth, respectively) were also reported along regional fault systems demarcating the area .
With this understanding, we have further attempted an analysis of regional geotectonic backdrop that may have been causing the stress accumulation. Definitely, the northwestern margin of the Indian plate has a crucial role in this. It includes the Aravalli craton, the Indo-Gangetic Plains and the Kachchh-Saurashtra blocks. The eastern boundary of the block is demarcated by the Great Boundary Fault (GBF), and the western margin is defined by the Suleiman Range-Jacobabad Ridge (Figures 1-3). The intervening region is divided into five blocks (from west to east): the Indus basin, Pokhran-Bahawalpur ridge, Jodhpur-Marwar basin, the Aravalli Mountain Range and the Eastern Plains (Mahi-Banas-Chambal basin). The regional faults demarcating the boundaries between the morphotectonic units extend from the Arabian Sea to the Himalayas, and have surface expressions of: Sukkur-Bahawalpur-Lahore-Narowal-Kathua (-Dalhousie) lineament (SBLK; Ravi lineament), Kori Creek-Umarkot-Jaisalmer-Suratgarh-Ropar (-Bilaspur-Mandi) lineament (KJSR), Bhachau-Nagar Parker-Luni-Umednagar (Jodhpur)-Nokha (Nagaur)-Sadarshahar-Sirsa-Pachkula (Chandigarh)-(Kalka-Kufri) lineament (BJSC), Jamnagar-Sirohi-Phulad-Degana-Sikar-Rohtak-Karnal-Jamuna-Khairi lineament (JSPDJ) and East Aravalli dislocation zone (Agucha-Dariba Fault). Though major parts of the region to the west of the Aravalli Mountains is occupied by the eolian sands of Rajasthan desert, the faults and lineaments have been tracked and mapped from surface geomorphic features and shallow surface geophysical signals, viz., straight courses of river systems, regional gravity anomalies and distribution of sedimentary basins on the surface (cf. Balakrishnan, 1997;Sinha-Roy et al., 1998). In the tectonic map of India by Merh (1995), the SBLK lineament is displayed as the "foredeep zone" and the combined block of the KJSR and BJSC as the "hinge zone", bordering the "shelf zone" to the west of the areas of the Aravalli-Delhi Fold Mountains. The interaction of alternate ridge (horst) and basins (grabens) with the MHT has led to the development of re-entrants over the ridges (horsts) and salient over the basins (grabens) in the frontal part of the Himalayas (Figure 7), which are manifested as the Dehra Dun re-entrant over the Delhi-Hardwar ridge, Kangra re-entrant over the Pokhran-Bhawalpur ridge and the Nahan salient against the Jodhpur-Marwar basin (graben). Over an area, delimited by the Rann of Kachchh in the west and fault system at the Shillong plateau margin in the east, the undeformed Indian lithosphere records signatures of breaking, in association with its flexure and underplating below the Eurasian plate along the MHT. The interactions between MHT and seismically active fault systems bounding the horsts and grabens present in the flexed Indian lithosphere are acting as zones of stress concentration in this highly dynamic, seismically active terrain. Since, it is argued that the fault systems, with effective coefficient of friction of ∼0.8 in the seismogenic Indian crust, supports the compressive force that transpires through the Indian lithosphere due to subduction process (Copley et al., 2011), it may be logical to think that concentration of additional stress in these stressequilibrated fault systems is the reason behind recent shallow-focus, low-magnitude events recorded around Delhi-NCR and the adjoining Himalayan front.
Additionally, when hypocentres of major damaging earthquakes are plotted on the reconstructed subduction profile of the Indian plate (Figure 7), they are found to fall in both Upper and Lower Crustal domains, invariably around the flexing zone. This observation raised question on triggering behind high-magnitude deep-focus earthquakes in the area; those, more often than not, are extremely damaging. Crustal heterogeneities including the fault systems may continue up to deep crustal level in the lithosphere, however the role behind generation of great events is still not fully understood (Schulte-Pelkum et al., 2019). Alternatively maximum strength is found to be associated with the bending zone (Kohlstedt et al., 1996;Conrad and Hager, 1999), and can trigger great to megashocks along subduction margin (Khan and Chakraborty, 2009;Khan et al., 2012). Using joint inversion of seismic and geodetic data sets, Copley et al. (2011) described the source model of 2001 M w 7.7 Bhuj earthquake and suggested that the stress associated with event was possibly related to the flexural effects of down-going Indian lithosphere FIGURE 7 | Block diagram illustrating the distribution of hypocentres (+) of 12 great damaging earthquakes within the reconstructed subduction profile of the Indian lithosphere (after Ansari et al., 2014). 4: 1945M 6.27 Chamba, 5: 1905M 7.8 Kangra 6: 1975M 6.8 Kinnaur, 7: 1991M 6.8 Uttarkashi, 8: 1803M 7.7 Garhwal, 9: 1999M 6.5 Chamoli, 10: 2015M 7.9 Nepal, 11: 1833M 7.6 Nepal, 12: 2015M 7.3 Nepal, 13: 1988M 6.8 Nepal-Bihar border, 14: 1934 beneath the Himalaya. Buckling of the Indian lithosphere due to the curved geometry, resulting from the underthrust below the Eurasian plate in the Himalayan Mountain Range-Tibetan Plateau, has been identified as the causative factor for the extension of the upper crust and compression in the lower crust; the finite neutral surface of the buckled layer was estimated to be at a depth of 25 ± 5 km (cf. Figure 1 in Copley et al., 2011). Two alternative possibilities were suggested for the deeper depth of the high-magnitude compressional earthquake events in the Indian lithosphere: (a) increase in strain rate, and/or (b) geometrical hardening of the Indian lithosphere due to along-strike and across-strike curvature related to subduction. Analysis of the distribution of earthquakes in the Indian plate interior is also attributed to the enhanced strain rate, causing a deeper brittle-ductile transition (Mohanty, 2011). Mohanty (2011) has further shown that some of the critically oriented faults (shear zones) in the Indian plate interior are associated with periodic accumulation of high stress/shortening and release of stress/extension. These observations regarding strainhardening due to high strain rate (∼266 × 10 -9 per year along N013°) and alternating sinistral and dextral slip (compression and extension) along the same fault has been corroborated by the field mapping and data analysis of the Bhuj area (Sinha and Mohanty, 2012). The critical fault passes through the epicentres of the 1956 Anjar event and 2001 Bhuj event and continues up to the Himalayas (BJHC in Figure 1A; see also F7.12 in Balakrishnan, 1997).
From above discussion it may be understood that the recent events around Delhi-NCR area cannot be analyzed in isolation without taking into consideration of events in entire northwestern India and associated Himalayan front. Khan et al. (2016) analyzed the frequency of occurrences of aftershocks following the 2001 Bhuj main shock and found it to be decreasing after 2014. The predominant stress direction of the aftershock activities were operative along the NNE direction towards the Himalayan front, passing through the Delhi region ( Figure 1C). The self-organized criticality (SOC) behaviour is occasionally persistent over a chaotic nature of sequence of shocks because of gradual restoration of stresses by tectonic loading. Our integrated study involving b-values and operative crustal-scale kinematics and geodynamics, might be accounting for a process-based phenomenon that indicates an impending earthquake in the seismic cycle ( Figure 6). Indeed, moderate-tolow 'b' values over the study-area and frequent incidents of smallmagnitude earthquakes bear indication for a large failure by triggering a great event. The regional crustal-scale heterogeneities including fractures forming the boundaries between NNE-SSW aligned horst and graben structures are found to be the sites for stress-concentration; their interaction with the MHT is enhancing the stress in the bending zone of the subducting Indian lithosphere. Failure is inevitable when the accumulated stress exceeds the yield strength of the materials. Past failures in the bent portion were recorded in 12 disastrous seismic events in the Himalayan region ( Figure 7). However, the conjecture drawn from this study can be more robust if results obtained from analysis of b-value clusters are backed with the 3D modeling of a rheological state, taking into account all major structural features involved in the complex structural heterogeneity in the region.

CONCLUSION
Frequent occurrences of earthquake shocks in Delhi-NCR, adjoining areas and in the immediate Himalayan front pose a serious alarm for population. Although no models are available to predict the impending earthquake in any region, statistical analysis backed by regional kinematics and overall geodynamic understanding can give insights into the details of recent seismic activities. We found that the west -northwestern part of India has been activated by triggering 134 earthquakes of M L ≥ 2.0 over only first half of 2020. The compilation and analysis of well-constrained dataset lead towards understanding of clustering and stabilization of the general moderate-to-low seismic b-values in the study area. The regional structural trends, tectonics in the western India, the operative stress field and seismicity, extended from the Kachchh region, are apparently found to be compatible up to the northwest Himalaya ( Figure 1B). The present study, thus, is summarized under the following concluding points: ➢Large number of events of magnitude M L ≥ 2.0 with occasional incidents of earthquakes of magnitudes reaching up to more than 5.0 ( Figure 4) indicate high-level stress transmission, accumulation and release in the study area. ➢Spatially clustered seismicity (Figure 1) is cyclic in time as well as in space domains ( Figure 5). The incidents of large number of events of M L ≥ 4.0 in the western India are related to the general regional kinematics and tectonics. ➢Moderate-to-low seismic b-values over time and space ( Figure 6) indicate stable high stress field, which may lead to fracturing and trigger a large earthquake in the future. ➢The fracture system in the down-flexing, over-compensated Indian lithosphere below the Indo-Gangetic Plains is equally accommodating the stress field, and may initiate nucleation for large earthquake. ➢Insufficient coverage of stations/network for the recorded events under the present study limits the determination of accurate focal depths of the earthquake events within reliable limit of uncertainties. Reliable focal depths of the events would have been more impressive for analysis of stress regime through a 3D structural analysis involving rheological understanding of the area. ➢The heterogeneities appear in the log-linear plots of FMD might be related to the complex structural setup of the study area. The used earthquake catalog seems to be incomplete because of poor coverage of recording stations. ➢The homogeneity of the catalogue could have been achieved by recording more number of low magnitude events and the present study would have been more convincing for forecasting the moderate-to-great earthquakes in this part of the Himalaya based on b-value analysis.