A holistic seismotectonic model of Delhi region

Delhi region in northern India experiences frequent shaking due to both far-field and near-field earthquakes from the Himalayan and local sources, respectively. The recent M3.5 and M3.4 earthquakes of 12th April 2020 and 10th May 2020 respectively in northeast Delhi and M4.4 earthquake of 29th May 2020 near Rohtak (~ 50 km west of Delhi), followed by more than a dozen aftershocks, created panic in this densely populated habitat. The past seismic history and the current activity emphasize the need to revisit the subsurface structural setting and its association with the seismicity of the region. Fault plane solutions are determined using data collected from a dense network in Delhi region. The strain energy released in the last two decades is also estimated to understand the subsurface structural environment. Based on fault plane solutions, together with information obtained from strain energy estimates and the available geophysical and geological studies, it is inferred that the Delhi region is sitting on two contrasting structural environments: reverse faulting in the west and normal faulting in the east, separated by the NE-SW trending Delhi Hardwar Ridge/Mahendragarh-Dehradun Fault (DHR-MDF). The WNW-ESE trending Delhi Sargoda Ridge (DSR), which intersects DHR-MDF in the west, is inferred as a thrust fault. The transfer of stress from the interaction zone of DHR-MDF and DSR to nearby smaller faults could further contribute to the scattered shallow seismicity in Delhi region.


Geology and tectonic setting of Delhi region
The Ganga basin, with an area of about 250,000 sq km falls within Long. 77° E-88° E and Lat. 24° N-30° N (Fig. 1b). It is located between the northern fringe of the Indian peninsula and the Himalaya and extends from Delhi-Hardwar Ridge (DHR) in the west to Munger-Saharsa ridge in the east (Fig. 1a,b). Delhi is located near the northern fringe of the Proterozoic Aravalli-Delhi fold belt and western edge of the Ganga basin (Fig. 1b). The terrain is generally flat except for a low NNE-SSW trending Delhi Hardwar ridge in the southern and central part of the area which consists of Quartzite while the Quaternary sediments, comprising the older and newer alluvium, cover the rest of the area. The thickness of the alluvium, both on the eastern and western side of the ridge, is variable but west of the ridge it is generally thicker (290 m).
The thick deposits of soft sediments of Yamuna plains plays a dominant role in ground motion amplification as experienced during past earthquakes 7 .
Historically, studies of the Himalayan foot-hills belt were initially conducted by Medlicott 8 , Theobald 9 , Oldham 10 , and Middlemiss 11 . Later, Wadia 12 and Auden 13 and several other officers of the Geological Survey of India mapped different parts of this belt. Agocs 14   www.nature.com/scientificreports/ proposed that the great Vindhyan basin must be extended into the Lesser Himalayan region. Sengupta 16 using the aeromagnetic data subdivided the Ganga basin into four parts separated by basement ridges or faults or both (Fig. 1b). Based on a re-interpretation of earlier gravity data, Sengupta 17 (1964) correlated the evolution of the Himalaya with the subcrustal movements below the Gangetic plains. The Oil and Natural Gas Commission, based on geophysical surveys (aeromagnetic, gravity, and seismic) and drilling data, in 1968 identified the three ridges in the Ganga basin named (from east to west) as Munger-Saharsa ridge, Faizabad ridge and Delhi-Hardwar Ridge (DHR) (Fig. 1b). Valdiya 18 correlated the transverse structures in the Himalaya to these three hidden basement ridges. The western boundary of the Ganga basin is delineated by DHR and the eastern margin by the northeastward continuation of the buried basement ridge (Munger-Saharsa ridge) 2 . The DHR was proposed with the least areal extent (6000 sq km) among all three ridges. Sastri et al. 2 and Karunakaran and Ranga Rao 19 described the shallow character of the DHR. The ridge was not traced with seismic survey beyond Meerut; 'the trend is probably obscured by a thick Neogene cover' 6 . Based on magnetic survey, Arora et al. 20 proposed a major conductive structure, namely, the Trans Himalayan conductor (THC), that strikes perpendicular to the Ganga basin into the foothills of the Himalaya and located east of Delhi (Fig. 2). Later, through a magnetic survey, Arora and Mahashabde 21 characterized the THC as a major electrical conductive structure (having a resistivity of 2 Ohm.m) with a width of 45 km and depth of 15 km following the strike of the Aravalli range and running into the Himalaya (Fig. 2).
Mallick et al. 22 , following the study of Raiverman et al. 23 , have suggested a deep-seated fault along the course of the Yamuna River formed by the flexure of the Indian Plate due to subduction beneath the Himalaya. Valdiya 24 and Chandra 25 have also indicated a fault zone along the strike of Aravallis in this area. By correlating seismicity with the changes in the Coulomb stress, Arora et al. 26 proposed along-strike segmentation of NW Himalaya, controlled by the subsurface ridges (underthrusting the Indian Plate) and by rift and nappe structures. They www.nature.com/scientificreports/ suggested the episodic reactivation of Delhi-Hardwar Ridge due to the strains resulting from the locking of Indian-Eurasian Plates as proposed by Arora 27 . Dubey et al. 28 inferred three NW-SE trending reverse faults in the Delhi region using Remote Sensing, Ground-Penetrating Radar (GPR), and Bouguer gravity anomaly data. However, due to very limited depth of penetration of GPR survey (a few meters), modern geophysical surveys with a higher depth of penetration (e.g., Magnetotellurics / Seismic) are imminent to verify and precise characterization of these faults. Dubey et al. 28 have also suggested that earthquakes that occurred near Rohtak and have orientation other than MDF (i.e. NE-SW) might be related to lithospheric crustal loading of the Himalaya orogeny on the Delhi-Sargoda Ridge. Based on the gravity and aeromagnetic investigations, GSI 29 proposed a NE-SW trending, 295 km long fault linking Indian peninsular craton in the south to Himalayan Frontal Thrust (HFT) in the north along the DHR and named it as MDF. At the junction of MDF and HFT, Jade 30 estimated the convergence rate of 10-18 mm/year between India and Tibet. The information on slip rates along major faults of Delhi region is not available. Patel et al. 31 delineated the shallow steep vertical faults near MDF (though MDF was not traced) using GPR survey and suggested MDF as a normal fault system at shallow subsurface and showing normal with oblique-slip motion.
Ravi Kumar et al. 32 published the Bouguer anomaly map of North India including the Ganga basin using well-controlled ground data and inferred that the Aravalli Delhi Mobile Belt (ADMB) and its marginal faults extend to the Western Himalayan front via Delhi where it interacts with the Delhi-Lahore ridge and further north with the Himalayan front causing seismic activity. Godin and Harris 33 , using Bouger gravity data of Delhi region derived from Earth Gravitational Model (EGM) 2008 have suggested that NE Delhi-Hardwar trend continues northeastward across the surface trace of the Main Frontal Thrust to the Karakoram fault. They further suggested that the DHR is delimited by the Shimla and Dehradun lineaments and proposed it as a horst with steeply-dipping normal faults on either side (Fig. 3). The Dehradun lineament connects the eastern edge of the Delhi-Hardwar Ridge to the Burang graben north of Shimla, the westernmost N-S graben of southern Tibet.
Dwivedi et al. 34 , through 3D structural inversion of gravity data (from Gravity Map (WGM)-2012 and gravity map series of India-2006 (GSI-NGRI, 2006)), speculated that NE trending Delhi Fold Belt deflected westward towards the shallower DSR and produce clustered seismicity in the hinge zone of this crustal bending near the Delhi region. They suggested that it is happening due to development of high strain resulting from crustal buckling of Delhi Fold Belt and DSR. They opined that the structural setup possibly developed after NW corner indentation and anti-clockwise rotation of Indian plate (post-Eocene collision) (as proposed by Voo et al. 35 ) led to the westward deflection of NE trending Delhi Fold Belt. In addition to geophysical and geological studies, the seismological studies have a special contribution in understanding the subsurface structures and seismotectonics of the region. www.nature.com/scientificreports/

Seismic monitoring in Delhi region
Among the far-field moderate to large earthquakes (1803 GK, 1991 UKS, 1999 CHM, 2015 GRK) experienced in the Delhi region from Himalayan sources, the earthquake of 1st September 1803 (1803 GH) is considered to be important as damage was observed in Delhi and its surrounding region. Different locations were proposed for this earthquake. Initially, this earthquake was considered as the 1803 Mathura earthquake (M 6.8) 10,36,37 . Later, it was studied in detail 38,39 and suggested renaming the event as the 1803 Garhwal earthquake. As mentioned in the preceding section, the Delhi region has also experienced near-field earthquakes from the local sources [Historical earthquake of Delhi (M6.5, 1720); Bulandshahar earthquake (M 6.7, 1956); and Gurgaon earthquake (M 4.8, 1960)] (locations given in Fig. 1b and Supplementary Table 2). The intensity of the 1720 Delhi earthquake was assessed as IX in the Old Delhi area. Though the exact epicenter of this event is uncertain; it was in the vicinity of Delhi 40 . The 1956 Bulandshahar earthquake was felt over a larger area and deaths as well as destruction to property were reported. The 1960 Gurgaon earthquake (M4.8) is the closest instrumentally located event to the Delhi region, though the location and magnitude were debated (Supplementary Table 2 (Fig. 2) with one of the nodal planes in NE-SW direction. Singh et al. 6 also analyzed the 25th November 2007 (Mw 4.1) earthquake in detail and given the strike-slip faulting with some normal component mechanism (having dip of 55° to 86°) (Fig. 2). Shukla et al. 46 used 6 to 10 first motions for estimating the focal mechanism of small-magnitude earthquakes which are insufficient/ and are often difficult to read and focal mechanisms may not be well-constrained 5 . Though Singh et al. 6 emphasized that the focal mechanism estimated by Bansal

Fault plane solutions and structural trends
Recently, three earthquakes occurred on 12th April, 10th May and 29th May 2020 were recorded by the more  Table 1 along with the FPS of 12th April 2020 determined by Pandey et al. 48 . Only those stations with cut-off signal to noise ratio >2 in the frequency range of interest are used for estimation of fault plane solutions of these events (Fig. 2). The computational details are given in Supplementary data. The FPS of 12th April 2020 shows two nodal planes striking at 13° and 253° with a dip of 55° each. The two nodal planes show rake of -135 (normal right lateral oblique) and -45 (normal left-lateral oblique) with dominant normal fault mechanism ( Table 1). The FPS of 10th May 2020 shows two nodal planes striking 32° and 275° with a dip of 75° and 31 o each. The two nodal planes show rake of -117 (normal right lateral oblique) and -31 (normal left-lateral oblique) with dominant normal fault mechanism ( Table 1). The second nodal planes of both these earthquakes have suggested the strike of 253° and 275° (Table 1), respectively, which are not consistent with either the trend of Aravalli belt (NNE-SSW) or with the trend of major fault lines (NNE-SSW to N-S) in the region (Fig. 2). Therefore, the nodal planes with strikes of 13° and 32° that are consistent with the Aravalli and major tectonic trends and are considered.
The FPS of 29th May 2020 shows two nodal planes striking at 10° and 124° with a dip of 37° and 73 o each. The two nodal planes show rake of -123 (normal right lateral oblique) and -29 (normal left-lateral oblique) with dominant normal fault mechanism ( Table 1) (Table 1). Both these earthquakes (29th May 2011 and 01st June 2017) were located close to DSR (~18 km and ~07 km west of MDF, respectively).
An earthquake of magnitude M4.9 was occurred on 5th March 2012 about 20 km south-west of this earthquake. Bansal and Verma 50 proposed a strike of N348 o and a dip of 48° with a rake of 131° (strike slip motion with reverse component) for this earthquake. The earthquake was located ~ 4 km west of the MDF. The estimated fault plane solution has suggested a NNW-SSE trending reverse fault; therefore, it may not be strictly associated with the MDF. Therefore, it is inferred that the earthquakes occurring to the east of DHR/MDF are following the normal with strike slip mechanism and the earthquakes located to the west of DHR/MDF follows the reverse with strike-slip mechanism.
Further, the earthquakes of 12th April 2020 and 10th May 2020 are located to the east of the DHR and southwestern edge of the Trans Himalayan Conductor (THC) proposed by Arora et al. 19 (Fig. 2). FPS of both the recent earthquakes have shown strike of NNE and steep dip of 55°-75° that are commensurate with the geometry of the edge of THC (Fig. 2). The depth of ~ 15 km has been estimated for both these events. Therefore, both the recent earthquakes (12th April 2020 and 10th May 2020) might be nucleated at the southwestern edge of the THC.

Distribution of strain energy in Delhi region
The energy is a direct indicator of the size of an earthquake and estimation of the accumulated energy in a region can provide valuable information regarding the potential seismic hazard of the region. The spatial variation of energy release can provide information on the potential locales of stress accumulation in a region in the absence  www.nature.com/scientificreports/ of geodetic data. Similarly, the temporal variation of seismic energy of a region can provide the different stages of energy release process 42,51 and could be used as a long-term earthquake precursor. In the present study, the spatial distribution of strain energy has been estimated for Delhi and the surrounding areas based on the earthquake catalog of the region for the period 1998-2020 taken from the website of National Centre for Seismology (https:// seismo. gov. in/ conte nt/ seism ologi cal-data). Earthquakes with magnitude range between M0.8 and M5.1 have been considered. The computational details are given in Supplementary material. The spatial distribution of strain energy estimated in the Delhi region has been shown in Fig. 4. The maximum energy, in the range of 08*10 11 Joule has been released in this period in the Delhi region. The energy is released mainly in two areas, (i) the area west of DHR/MDF, and (ii) the area east of the DHR (central and NW Delhi). Almost twice the amount of energy has been released in the western part (at the contact zone of DHR and DSR), compared to the eastern part (east of DHR).
During an earthquake, the normal faulting is caused by the gravitational potential, however in case of reverse and strike slip fault, the energy is accumulated as elastic potential. The rocks get deformed under compression, are characterized by yield stresses about 10 times larger than yield stresses in tensional stress fields 52 . Additionally, reverse faults need more energy to move the rocks as compared to thrust in case of reverse fault, the hanging wall moves against gravity. Therefore, the energy dissipation in reverse fault is always more than the thrust and normal faulting. In Delhi region we infer the reverse faulting in the western part of DHR-MDF and normal faulting in the eastern part.

Proposed seismotectonic model
In general understanding, a seismotectonic model suggests the correlation of seismicity with the fault lines of the area and the slip rates along these faults. Delhi area is, however, under-represented in geodetic studies as the GPS network is to be established and no study on Active fault mapping is available along the major faults to confirm slip/slip rates. Taking advantage of the quality data generated by the local seismological network, we propose a seismotectonic model of the Delhi region through the integration of seismicity characteristics, the associated structural features and the estimates of strain energy released in the region. From subsurface structural appraisal, it is inferred that the MDF/DHR follows the trend of Aravalli Delhi Mobile Belt and is a NE-SW trending horst structure with steep normal faulting that continues northeastward across the surface trace of the Main Frontal Thrust to the Karakoram fault. The past seismological studies [4][5][6]5 have also suggested normal with a strike-slip focal mechanism in the area east of DHR/ MDF. The focal mechanisms of three recent earthquakes (M 3.5, 12th April 2020 and M3.4, 10th May 2020 and M 4.4, 29th May 2020) to the east of DHR/MDF have also shown normal with strike-slip focal mechanism. The energy released in the last two decades also indicates a similar focal mechanism. Shukla 46 through FPS of 10 earthquakes (of 2001-2004, with single component seismographs) and Bansal and Verma 51 through FPS of 01 earthquake recorded with digital strong motion data have suggested DSR, bounded by reverse fault/thrust is an NW-SE trending structure. The fault plane solutions of the two earthquakes (29th May 2011 and 01st June 2017) fallen close to DSR have also shown reverse/thrust faulting with strike slip mechanism. The DSR appears to interact with the DHR in the west of Delhi (Fig. 5a,c). The interaction zone of DHR/MDF and DSR is seismically most active in last two decades (Fig. 2). The scattered and mostly shallow focus seismic activity in the region is inferred to be associated with the presence of faults that got activated time and again due to transfer of stress from the interaction zone of DHR-MDF and DSR. The seismotectonic constraints are presented in the form of a model in Fig. 5. www.nature.com/scientificreports/ We have also examined the possibility of THC as the causative structure of two recent earthquakes, M3.5, 12th April 2020 and M3.4, 10th May 2020. Due to high conductivity, the THC might be filled with mantle derived fluid and is ductile in nature. Hence it may not allow large accumulation of strain. Therefore, the strain can be accumulated on both lateral edges of the THC (Fig. 2). Since the FPS of two recent earthquakes (of M3.5, 12th April 2020 and M3.4, 10th May 2020) has shown strike of 13-32°, i.e., NNE, and steep dip of 55°-75 o to the NE, we reject the possibility of THC as a causative structure in those cases.

Discussion and conclusions
The densely populated, socio-economically important, housing national capital, the Delhi region has been experiencing earthquakes from both regional (Himalayan) and local sources. The seismicity due to local earthquakes, although not associated with significant damage to property, has created panic among the public in the epicentral areas. The estimation of focal mechanism of recently occurred earthquakes of 12th April 2020 (M3.5) and 10th May 2020, (M3.4) and 29th May 2020 (M4.4), the spatial distribution of energy in the Delhi region in last two decades and appraisal of local sources (geological, geophysical, and seismological) have been integrated to understand the causative structural sources and mechanism of seismicity in the Delhi region. The fault plane solutions of these three earthquakes have suggested a nodal plane trending NNE-SSW direction with a steep dip of 37° to 75° in the NE direction and normal with strike-slip mechanism.
Our analysis on appraisal of subsurface structures of the region suggests three probable major earthquake mechanisms in the Delhi region: (i) episodic reactivation of DHR due to the strains resulting from the locking of Indian-Eurasian plate, (ii) lithospheric crustal loading of the Himalayan orogen on the Delhi-Sargoda Ridge, (iii) interaction of Aravalli Delhi Mobile Belt with Delhi-Lahore Ridge and further north with the Himalayan front.
Since the year 2000, low level seismicity has been observed in the northern part of the Delhi within the zone of interaction of DHR and Himalayan Frontal Thrust. Therefore, the seismic potential of minor to moderate earthquakes cannot be ruled out due to transfer of stress from the intersection zone (of DHR and Himalayan Frontal Thrust) through DHR. A simplified subsurface mechanism of seismicity in the Delhi region has also been proposed in the current study. We believe that the proposed seismotectonic model will serve as a starting model for future studies. NCS has already planned to enhance the local seismic network in the region to help in precise location of smaller events and constraining their focal depths accurately. Also, other investigations such as active fault mapping and subsurface imaging, which are in pipeline would help in refining and strengthening the proposed model.
The seismological investigations, energy estimation and appraisal of subsurface structures have led to the following major results: (i) Two major ridges (DHR and DSR) are interacting in the west of Delhi NCT; DSR www.nature.com/scientificreports/ is bounded by thrust/ reverse faulting and DHR is a steep horst structure associated with normal faulting, (ii) The earthquake in the Delhi region have occurred in two major clusters in last two decades; each is located west and east to DHR/ MDF, (iii) Higher seismic energy (almost two times) has been released in the western cluster located west to Delhi as compared to eastern cluster, (iv) two major structural/ fault mechanisms in the region have been recognized-thrust/ reverse with strike slip mechanism in the western part and normal with strike slip mechanism in the eastern part due to the presence of DHR/ MDF associated with steep normal faulting in the east and presence of DSR bounded by thrust in the west, respectively, and (v) The scattered shallow seismicity in the region is inferred to be due to stress transfer from interaction zone of two ridges (DHR and DSR) to the nearby faults.

Methods
The earthquakes occurred on 10th May (M3.4), 29th May 2020 (M 4.4), 29th May 2011 and 01st June 2017 are recorded up to 32 broad band seismometers. The events were located using the Seisan software package 53 maintaining a small root mean square error (rms) of 0.27 to 0.5, respectively using the data of the stations having good signal to noise ratio. The fault plane solutions (FPS) of these events have been determined using 12 local stations waveform data (NPL, AYAN, JHJR, KUDL, SONA, UJWA, GNR, JMIU, KKR, NRLA, NDI and BISR) (Fig. 2). The ISOLA software package has been used to determine the FPS. A correlation coefficient (between the observed and synthetic data) of > 0.5 and the double couple percentage (DC %) of > 50% in the moment tensor decomposition have been considered for finalizing the FPS (For details see supplementary data). Addionally, we also determined the Fault Plane Solutions of the events using first motion data (supplementary Figs. 6 and 7). In this exercise, FOCMAC subroutine in SEISAN software package has been used. In all cases FPS obtained by first motion polarity are consistent (in terms of fault mechanism) with the FPS obtained by waveform inversion technique. However, strike, dip and rake values are found to be different but agreed closely with the data derived by waveform inversion. The strain energy has been estimated for Delhi region using the standard formulation suggested by Kanamori 54 to examine its spatial distribution. The earthquake catalog of the Delhi region for the period of 1998-2020 was used for the purpose (For details see supplementary data). The FPS computed in the present study have been compared with the earlier FPS in Delhi region to help in refining the current understanding of subsurface structural trends. The past works on geophysical, geological and seismological studies have been reviewed and a link between trends suggested by the FPS, energy estimation and the subsurface structures has been established. www.nature.com/scientificreports/