On 03 July 2020, an earthquake M4.6 occurred at local time 14:00 h IST (08:30:00 UTC), located ∼50 km southwest of the National Capital Region of Delhi. Although the people living in the city of Delhi and in other nearby cities were engaged in their day time activity, the event was widely felt across an area occupying more than ∼100 sq km from the epicenter. The maximum intensity for the event, close to the epicenter, was obtained to be V on the Modified Mercalli Intensity scale based on the felt experiences shared by people on the official website of the National Center for Seismology (NCS). As reported in the social media and local newspaper as well as based on the questionnaire responses, it was found that the event was felt prominently by many indoor people as compared to the few outdoor people, and it was not unusual. Nevertheless, the event, though small in magnitude, remained frightening to a few, as reported in the media. Five months later, another event of slightly smaller magnitude M4.2 struck the region on 17 December 2020, the epicenter of which was located very close to the previous event of 03 July 2020 (M4.6), ∼10 km away from it to the NE with a focal depth of 21 km (Figure 1).

The epicentral regions of both the events are characterized by low as well as sparse seismicity (nearly quiescent-like situation) in an area of ∼100 sq km. Furthermore, analysis of the InSAR (Sentinel-1B) data also showed no significant ground surface deformation resulting from the 2020 events. Conspicuously, identification of the causative source for these events in the field remained challenging. The adjacent regions to the north of it, close to Delhi and Haryana, had witnessed considerable seismic activity. Several moderate-magnitude earthquakes have occurred in the recent past (Figure 1), and the area is also characterized to be seismically very susceptible and vulnerable to great earthquakes (Valdiya, 1976; Chandra, 1992). Many studies have been carried out in the region to characterize seismic sources and regional tectonics (Chouhan, 1975; Verma and Greiling, 1995; Singh et al., 2002; Shukla et al., 2007; Bansal et al., 2009; Pandey et al., 2020). Furthermore, it is undeniable that the Himalayan Thrust system and reactivation of the fault systems of Delhi Fold Belt are the prime causative sources of the seismic hazard in Delhi and adjoining regions (Chandra, 1992; Bilham et al., 2001; Gupta et al., 2013; Shukla et al., 2016). Based on the examination of seismic data (1960–2022), it has been observed that there are sporadic instances of seismic activity dispersed throughout significant faults like Mathura fault, Sohna fault, and various mapped lineaments in the region (Figure 1). The dispersed distribution of seismic activity around these geological structures necessitated additional exploration of the fundamental geological and tectonic mechanisms responsible for such activity. The continuous monitoring and meticulous analysis of seismic events occurring in close proximity to these faults and lineaments contribute significantly toward a comprehensive assessment of the seismic hazards in the region. The ongoing investigation and diligent observation of seismic activity are imperative in order to gain a deeper comprehension of the intricate nature of these geological formations and bolster our capacity to evaluate and alleviate seismic hazards in the region.

In the present study, we analyze the best-recorded broadband data from 26 seismic stations of the national network corresponding to the recent earthquakes on 03 July 2020 and 17 December 2020 for source study and to understand the present-day stress fields responsible for the occurrence of these events in the relatively quiescent region. We also analyze the data for understanding the possible connection between the July 2020 and December 2020 earthquakes by computing the static Coulomb stress. This is carried out by using a uniform slip model as well as the source parameters acquired from the waveform inversion analysis. We also examined the stress consequences of the July 2020 event and its relationship with the December 2020 earthquake. The waveform data recorded at seismic stations located up to ∼500 km from the epicentral zone are used in the source parameter estimation and full-waveform inversion analysis for the fault plane solutions. Furthermore, a stress inversion modeling is carried out for understanding the seismogenesis of the 2020 events, and, subsequently, the causative fault of the events is ascertained based on the orientation of the extensional principal stress field.

The study area is located in the northernmost part of the Aravalli Range, bordering the regions of Haryana and northeastern parts of Rajasthan. The area is overlaid with sedimentary deposits on both sides, forming alluvium plains to the east and west of the Aravalli. The folded Aravalli range is the main tectonomorphic feature in eastern Rajasthan. The ensemble consists of Proterozoic age Aravalli and Delhi Supergroup rocks, which are on top of the Archean basement. Most of the Aravalli range consists of Delhi Supergroup rocks (Verma and Greiling, 1995; Bansal et al., 2021). In Delhi, both intrusive and extrusive rocks are found in a large area.

The Aravalli–Delhi orogenic belt abundantly occupies faults, fractures, and lineaments oriented criss-cross along NE–SW, NNE–SSW to NE–SW, and NW–SE to WNW–ESE directions (Singh, 1988). The interrelationship of sedimentary facies and paleo-environments of the Delhi Supergroup can best be explained by an evolving intra-cratonic rift-basin model that is found to be comparable to several well-studied continental rift basins. Ravindra and Bakliwal (1983) discussed some of the faults in their studies and found them related to the Delhi orogeny. However, a few of them were characterized by some unique sedimentary and volcanic features, which indicate syn-sedimentary origins (Singh, 1982), and they are found to be seismically active. A prominent fault/ridge system to the NW of the study region, called the Delhi–Haridwar ridge (DHR) and Mahendragarh–Dehradun subsurface fault (MDSSF), is oriented NE–SW, and prominent small-to-moderate earthquakes have been reported along it (Figure 1).

The earthquakes originating from nearby sources are typically intraplate, with shallow focal depth and intermediate magnitudes. The seismic activity is observed to be clustered in three distinct regions: (a) to the west of Delhi, (b) in close proximity to Sonipat, and (c) surrounding Rohtak (Prakash and Shrivastava, 2012). Over time, it has been noted that the majority of earthquakes in the western portion of Delhi are linked to the Sohna fault, which trends in a north–south direction. These earthquakes also occur along the tri-junction of the Delhi–Haridwar ridge and the Delhi–Sargodha fault as well as along the axis of the DFB (Das et al., 2018). In addition, there are other earthquake occurrences scattered around this region. No mapped faults or features are specifically related with activities around Rohtak and Sonipat. The activity in this area is likely caused by the extension and intersection of multiple unmapped faults that are buried behind a thick layer of alluvium and the Delhi–Hardwar ridge. Delhi is experiencing a significant number of scattered earthquakes, which may be attributed to the presence of unexplored active faults. In recent years, a fault near Khanpur in South Delhi was detected (Bansal et al., 2021).

Earthquake activity in Delhi and its surrounding regions was monitored using 17 permanent broadband seismic stations that are maintained by the National Center for Seismology (NCS) under the aegis of the Ministry of Earth Sciences, India. Later, eight additional seismic stations were deployed in the region soon after the recent earthquake sequence of 2020 (Figure 1). Conspicuously, the present setup of the seismological observatories has improved the detection capability of the NCS, and smaller magnitude earthquakes (i.e., M > 2.0) are being precisely located in and around the Delhi region in near real-time. The earthquakes in the region are being recorded at these seismic stations by using REFTEK 151A-120 tri-axial orthogonal sensors with a flat velocity response in the frequency range 50 Hz–120 s. The VSAT linked network accurately records GPS-synchronized seismological data with a sampling frequency of 100 Hz. In the present study, we analyzed the waveform data of the recent events of 2020 (M_w 4.6 and 4.2) that were recorded by these local seismic stations.

Well-recorded waveform data at about 14 seismic stations were analyzed using SEISAN software (Havskov and Ottemöller, 2003) to relocate the two recent earthquakes, viz., Mw 4.6 and 4.2 (Figure 1). The P-wave velocity model of the Himalayan foreland basin (Mitra et al., 2011) is considered in this analysis. The arrival times of P- and S-waves were picked on the vertical and horizontal seismograms, respectively, in the case of both the events. The average root mean square (rms) errors, which indicate a measure of the time difference between the observed and estimated travel times, were found to be ∼0.40. Upon analysis, it was determined that the uncertainty associated with identifying the epicenters and determining the focal depth of both events fell within the acceptable range, as the error margin was found to be within the permissible limits of ±2.0 km and ±2.4 km, respectively. (Table 1).

Furthermore, to characterize the epicentral region, we analyzed earthquake data (NCS catalog) for the period 2000–2022 for a rectangular area (∼100 sq km) covering the study region (Figure 2A). Prior to the analysis, the catalog was de-clustered for duplicate events, foreshocks, and aftershocks, and further it was homogenized by converting different magnitude scales into M_w. A total of 62 events were homogenized in the magnitude range M_w 3.0–4.6. The homogenized data are used in the analysis for estimation of various statistical parameters, viz., cumulative moment release, changes in the magnitude–frequency distribution of earthquakes (b-value), and temporal (year wise) variation of events versus magnitudes (Figures 2B–D), that characterize the sources in any specific region. Figures 2A, D clearly show that most of the activity is concentrated northwest of the epicentral region and during 2004–2012. Throughout the sequence in years 2000–2022, the distribution of magnitudes is similar, except for the 2020 event, which is largest magnitude event recorded in the area.

The ISOLA software (Sokos and Zahradnik, 2008) is used in the present study to estimate the fault plane solutions of the events that occurred on 3rd July and 17th December 2020, after constraining the hypocentral location in the previous sub-section. The ISOLA code is developed to invert local and regional full-wave seismograms for single- and multiple-point source models following the iterative deconvolution approach (Kikuchi and Kanamori, 1991). The decomposition in ISOLA as part of the inversion process, namely, volumetric (ISO), compensated linear vector dipole (CLVD), and double couple (DC) stipulates as ISO% + CLVD% + DC% = 100% (Vavryčuk, 2001; Benetatos et al., 2012; Kühn and Vavryčuk, 2013). We considered the waveform data of 14 broadband seismic stations (Figure 1, inset), which are spread across the region with good azimuthal coverage, in the analysis to retrieve the source mechanism, focal depth, and origin time of the events. However, the optimum solution is inferred from the variance reduction function between the observed and synthetic seismograms considering nine seismic stations only, which were matching well. Prior to the analysis, the signal-to-noise ratios at each station were analyzed to fix an appropriate frequency band in which the signal could be inverted. In addition, the instrument effects were removed from the waveforms to obtain the actual ground velocity.

The waveform data were band-pass-filtered in the frequency range (0.2–0.35 Hz) to satisfy the assumption of a point-source model. Green’s functions were computed subsequently along the source-to-receiver paths, over a range of focal depths, using a generic velocity model that is adequate to model a seismogram. In the present study, we considered two different velocity models for the region, viz., Mitra et al. (2011) and Rajat et al. (2023 communicated). After a thorough analysis, it was found that the computed output derived from the Rajat et al. (2023) model exhibited the highest degree of conformity with the observation as compared to other models that were evaluated. Figures 3A, B show a correlation plot between the observed and synthetic waveforms as well as the focal mechanism for a series of focal depths ranging from 2 to 50 km with 2-km intervals for the 03 July 2020 event and from 10 to 50 km with 2-km intervals for the 17 December 2020 event, respectively. In order to avoid the potential trade-offs between the origin time and focal depth, a grid search technique was adopted, and a series of moment–tensor inversions were performed for obtaining an optimum solution over a range of focal depths and origin times to maximize the variance reduction, i.e., correlation. The origin times in the range between −3 s and +3 s were considered in the analyses with a time shift of 1 s from the origin time. Figures 4A, B show the variance reductions that quantify the fitness error between the observed and synthetic seismograms of the events that were computed for the optimum deviatoric and full moment–tensor solutions. The best-fit solutions of the events are eventually achieved for the origin times 13:24:59.0 UTC and 18:16:46.00 UTC of the events on 3rd July 2020 and 17th December 2020, respectively (Table 1). Evidently, the seismic moments of the events M_w 4.6 and 4.2 were obtained to be 8.287E+15 Nm and 2.29E+15 Nm, respectively. The moment magnitude of the events was eventually computed using an empirical relationship dependent on the seismic moment, as proposed by Hanks and Kanamori (1979). The moment–tensor solutions have been obtained at different values of decompositions, say, 0.0% VOL, 31.3% CLVD, and 68.7% DC components.

The focal mechanism solutions show normal faulting mechanisms with a minor component of strike–slip with two nodal planes orienting almost NE–SW and NW–SE (Figure 1). The fault parameters of the events that were retrieved from double-couple focal mechanism solutions are listed in Table 1. The uncertainty in the fault plane solution was also performed considering all possible solutions that were matching with the polarities obtained using the stress inversion.

Stress tensor inversion was carried out to calculate the stress tensors that usually best explain the seismological observations. In the present study, a technique developed by Vavryčuk (2011) and Vavryčuk (2014) is used to determine the orientation of local stress axes along the activated fault/lineament as well as to estimate the stress ratio ‘R’ (also called shape ratio) in the dislocation zone. R is defined as the ratio of the intermediate stresses relative to the maximum horizontal stress and minimum horizontal stress as R = σ 2 − σ 3 σ 1 − σ 3 .

In other words, the R value defines the shape of the stress ellipsoid. This technique, which is widely accepted as a powerful tool, is well-known for understanding the stress and tectonic regime of a region. The inputs used in this approach are primarily orientation and slip of the fault nodal planes. In the present study, these input parameters were determined using the fault plane solutions obtained from the moment–tensor inversion of the well-recorded recent events of 2020 (Table 1) as well as P-wave first motion analysis of local events (M > 3.5) that occurred during 2010–2020 (Table 2).

Any variation in the estimation of the focal mechanism solutions may result in errors in the estimation of the stress ratio and identification of the causative fault plane out of the two nodal planes (Vavrycuk, 2011). However, such limitation could be resolved by incorporating an algorithm useful to identify the causative fault plane and is based on evaluation of the fault instability (I) following the formulation proposed by Vavryčuk (2011; 2014) as I = τ ‐ μ σ ‐ 1 / μ + 1 + μ 2 , where μ is the friction coefficient and τ and σ are the normalized shear and normal tractions on a fault, respectively. The method consists of maximizing the slip shear stress component using an iterative algorithm based on robust grid-search inversion over the frictional coefficients ranging from 0.2 to 1.0. The stress axes and fault orientation were determined using the fault instability formulation after a few iterations. In the first iteration, the stresses were obtained by introducing a random perturbation in the focal mechanism. The obtained stress axes later identified the activated causative fault plane, which subsequently was used in the stress inversion analysis. The process was repeated for maximum 100 iterations to obtain the final stress axes. An open-access STRESSINVERSE code written in MATLAB (Vavryčuk, 2014) is used for the stress inversion. In the analysis, four parameters of the stress tensor were recovered that include orientation of three principal stress axes, viz., maximum compressive stress (σ1), intermediate stress (σ2), and minimum compressive stress (σ3), and the shape ratio (R). Furthermore, uncertainties were also ascertained for random perturbation of the focal mechanism by contaminating the signal with an artificial noise. A random noise of 100 realizations was used in the inversion introducing a level of noise of the order 10° in the estimated accuracy of input focal mechanisms. The focal mechanism solutions of total eight events including recent two events of 2020 were inverted to recover the tectonic stress regime in the seismogenic zone (Table 2). Figures 5A–D show the results of the stress inversion analysis indicating an extensional stress field orienting NW–SE with a shape ratio of 0.70 (Table 3).

Stress drop is a measure of the difference between the stress before and after an earthquake, which may potentially be used to distinguish the geodynamic processes causing earthquakes. The strength of an earthquake in the source is represented by seismic moment (Mo). In the present study, we estimated the seismic moment, source radius, and stress drop of the events (Tables 4a, 4b) using amplitude spectra of S-waves following Brune’s model (Brune, 1970), which basically belongs to the circular shear dislocation model. The seismic moment is obtained using the formulation M 0 = 4 π β 3 A 0 R θ ϕ , where A₀ is the low-frequency spectral amplitude of S-waves, R is the hypocentral distance, R_θϕ is the radiation pattern of S-waves, β is the S-wave velocity at the source, and ρ is the density at the source. We mention that the correction for free surface amplification and other effects were appropriately accommodated in the radiation factor of S-waves. The corner frequency fc, at which the low frequency and high-frequency asymptotes of amplitude spectra intersect, is used to calculate the radius of a circular fault r using the formulation r = 2.34 β 2 π f C .

Subsequently, the rupture area could be obtained using πr². Furthermore, the stress drops of the events were estimated using the formulation Δ σ = 7 16 M 0 r 3 .

The moment magnitudes of the events were eventually computed using the Hanks and Kanamori (1979) formulation M w = 2 / 3 * ⁡ log 10 ⁡ M o   ‐ 10.73 .

The S-wave spectra at few seismic stations corresponding to the 03 July 2020 and 17 December 2020 events are demonstrated in Figures 6. The average seismic moment, moment magnitude, corner frequency, and stress drop of both the earthquakes are listed in Tables 4a, 4b. As the events were deep-seated (Table 1), the seismic parameters were estimated considering the S-wave velocity of the granitic layer, i.e., 4.53 km/sec. The moment magnitude is obtained to be Mw 4.6 for the 03 July 2020 event and Mw 4.2 for the 17 December 2020 event. The corresponding average seismic moments were found to be 8.28E+15 Nm and 2.29E+15 Nm, and the respective stress drops were estimated to be ∼164 and 137 bars. Such higher stress drops for the intraplate events corroborate well with the results obtained earlier for small-to-moderate magnitude earthquakes (Allmann and Shearer, 2009; Kumar et al., 2014; Sairam et al., 2018). In a recent study, ∼26 bar stress drop was estimated for a relatively small magnitude earthquake (M_w 3.5) that had occurred close to the northeastern part of the National Capital Territory of Delhi. Shearer et al., (2006) and Pandey et al., (2020) emphasized even higher stress drops for the events of intraplate origin with focal mechanism showing relatively a high fraction of normal faulting.

The seismic fault rupture results in accumulation of the forces that cause a lasting alteration in the stress field in a region nearby the earthquake location, and it is termed as static stress change. Such a change is found to be influenced by different parameters, namely, earthquake location, rupture shape, and type of fault motion (Mitsakaki et al., 2006). One way to quantify this change is through the Coulomb failure function (ΔCFF) analysis, which is influenced by variations in both shear stress (Δτs) and normal stress (Δσn). The change in failure stress, ΔCFF, can be determined using the following equation: Δ σ f = Δ τ s + µ ’   Δ σ n , where µ' is effective coefficient of friction on the fault that ranges between 0.2 and 0.8, depending on the fault mechanism, and it accounts for the unknown impact of changes in pore pressure (King et al., 1994; Mouyen et al., 2010). The parameters Δσ_n, Δσf, and Δτ_s are normal stress, shear stress, and shear strength, respectively.

Coseismic static Coulomb stress changes caused due to an event play a crucial role in earthquake-triggering and clustering processes and have a significant impact on the occurrence of aftershocks and subsequent mainshocks in the region (Deng and Sykes, 1997; Harris, 1998; Stein, 1999; Freed and Lin, 2001; King et al., 2001; Freed, 2005; Steacy et al., 2005; Parsons et al., 2008). A small and rapid increase in static stress of less than 0.01 MPa can shorten the time taken for an earthquake to happen; however, a similar drop in static stress can postpone the occurrence of an earthquake (Marsan, 2006; Parsons et al., 2008; Toda et al., 2012; Kariche et al., 2018). When the shear stress loading on the fault surpasses the shear strength, resulting in an accelerated rupture, a positive Δσf is observed. In this particular scenario, the initial seismic event has the potential to expedite the deterioration of a neighboring fault in close proximity. If the change in the normal stress exhibits a negative value, it implies that the shear strength is enhanced and the loading is reduced. Consequently, it can delay the second event from the impending failure and leads it into a region commonly referred to as “stress shadow.” However, a positive change in the Coulomb failure stress transfer within the specific area surrounding the mainshock leads the existing sources into a state of extreme proximity to failure.

The Coulomb stress change analysis computes stress variations transpiring on a designated recipient fault based on the fault slip model (King et al., 1994; Lin and Stein, 2004; Ma et al., 2005; Toda, 2008; Toda, 2008; Ishibe et al., 2011). In addition, various researchers have assimilated the fault plane solutions in computing the Coulomb stress variation (Hardebeck et al., 1998; Meier et al., 2014; Asayesh et al., 2020). Furthermore, using the Coulomb stress analysis, several researchers have inferred that the fault planes may undergo induced stress along the nodal planes ( Wu et al., 2017; Asayesh et al., 2020). The present study is focused on two consequent significant earthquakes (M4.6 of July 2020 and M4.2 of December 2020) that had occurred in the northern part of the Alwar district and to the southwest of the Delhi NCR, within a very short span of time of about 5 months. Considering the proximity of both events, a potential connection between the fault ruptures, stress transfer, and fault interaction is expected. Consequently, the Coulomb failure analysis was carried out to understand whether the first event (M4.6, July 2020) has potentially led to the advancement of the triggering process for the second event (M4.2, December 2020). In order to analyze the potential correlation between both the seismic events, we used Coulomb 3.3 software, which uses elastic half-space modeling techniques.

In the analysis, we assumed an elastic half-space, an apparent coefficient of friction of 0.4, a shear modulus of 32 GPa, and a Poisson’s ratio of 0.25, which allowed us to compute the static Coulomb stress variations as described by Toda et al. (2011a); Toda et al. (2011b). The present investigation involves computation of the static stress alteration on stationary receiver earthquake ruptures and optimally aligned faults. This is accomplished by employing the effective coefficient of friction (µ') of the fault, in conjunction with the strike, dip, and rake information of the fault as well as its hypocentral parameters. The specific receiver fault approach was used to assess the impact of static stress variations induced by the earthquake that occurred in July 2020 on the fault responsible for the earthquake that occurred in December 2020. The inferred model (Figure 3) indicated that the fault is oriented in a NE–SW direction, exhibiting characteristics of a normal fault with a minor strike–slip component (Table 2). This plane was further used as a source plane for the Coulomb stress changes caused by the July 2020 earthquake. Figure 7 demonstrates the Coulomb stress variations resulting from the July 2020 event that is considered a source fault, which refers to the fault plane that experienced slip during the earthquake. However, the source of the event of December 2020 is considered a receiver fault in the analysis. Evidently, four stress lobes were obtained that exhibited positive values, and three stress lobes demonstrated negative values. The initial classification, however, showed a pair of lobes exhibiting pulses orienting southwest and southeast, in addition to another pair of lobes displaying pulses orienting north–northeast and east–southeast. It is to be mentioned that the selection of the specific nodal plane is based on the mapped fault by GSI (source: BHUKOSH web portal). The analysis indicated a positive Coulomb stress change in proximity to the northeast segment of the fault (Figure 7A), due to the July 2020 event which characterized this segment as a potential zone for future failure.

In the present study, we investigated recent two significant earthquakes that occurred in the Delhi NCR during 2020 to characterize the sources and seismogenesis of the study region. The events were relocated precisely with rms error ∼0.40, considering a P-wave velocity model of the Himalayan foreland basin (Mitra et al., 2011). The events were found to be in proximity, with a separation of ∼10 km, and they are located close to a NE–SW-trending lineament mapped in the study region (Figure 2). Full-waveform inversion using ISOLA software indicates deep-seated focal depths of the events in the mid-crustal range, i.e., >20 km (Table 2). The focal mechanisms of the events demonstrate two distinct nodal planes orienting in the NE–SW and NW–SE directions, which appear to be plausible according to the tectonic settings in the study region. Focal mechanisms of both the events, indicating normal faults with slightly dextral strike–slip movement along NE, however, corroborate well with the lineament mapped in the vicinity of the events. Interestingly, the past earthquakes in the adjacent regions close to the Delhi NCR have mostly witnessed normal faulting with minor strike–slip components, orienting in the NNE direction parallel to the mapped major faults such as DHR and MDDSF.

Stress inversion analysis yields stress parameters, viz., focal sphere with compression (P) and tension (T) axes together with best-fit principal stresses, stereogram of the principal stresses with uncertainties, Mohr circle, and the histogram of the shape ratio (Figures 5A–D). The predominant orientation of the P-axes (red circles) suggests a maximum shortening in the NE direction, in the source region, as depicted in the stereogram (Figure 5B). Evidently, the azimuth of the principal stress, σ1 (red dots), is found to be approximately N25°. Interestingly, the azimuth of σ1 correlates well with the P-axis orientation. The dominant tectonic style is usually determined on the basis of resulting plunge values of σ1 (red dots), σ2 (green dots), and σ3 (blue dots), which are obtained to be around 52°, 38°, and 10°, respectively. The σ1 axis is found closest to the vertical position (highest plunge), indicating the dominance of the extensional stress regime in the region.

The principal stress ratio “R” is determined to be 0.7 (Figure 5D; Table 3). Using the R value, the dominant tectonic regime of the study region can be ascertained based on the stress regime index (R′), which purely depends on the plunge of the three principal stress axes (Kassaras and Kapetanidis, 2018). The factor ‘R′’ may be calculated using the shape ratio ‘R’ as well as the direction of the principal axes (Delvaux et al., 1997; Czirok, 2016). Accordingly, in case of (i) R′ = R; σ1 is found to be vertical and the tectonic stress regime would generally be considered extensional, i.e., of the normal slip type, (ii) R′ = 2 - R; σ2 is found to be vertical in such a case, and the tectonic stress regime would generally be of strike–slip nature, and (iii) R′ = 2 + R; σ3 is found to be vertical, and the tectonic stress would generally be of the thrust regime in this case. Conspicuously, the value of R′ would always be ranging between 0 and 3 as R normally varies from 0 to 1, indicating a clear extensional regime to radial compressional regime following a linear progression (Czirok, 2016). In the present study, σ1 is found to be almost vertical, suggesting a situation where the value of R′ equals to that of R, indicating the presence of normal faults in the region. According to the formula proposed by Zoback (1992), the direction of maximum horizontal stresses (Sh_max) can be determined based on the azimuth and the plunge of the principal stress axes. The Sh_max azimuth is calculated to be N15°E (azimuth of T + 90°) for the study region (Table 3), suggesting the region to be subjected to normal faults with strike–slip faulting and Sh_max direction almost to the north. The stress regime for the present study region is obtained using very limited data, suggesting the region to be heterogeneous; however, the result could be improved based on more available fault plane solutions for the nearby regions. Yadav et al. (2022) carried out a study for the Delhi–Aravalli region and suggested that the maximum principal stress in the region is oriented to N–S to NE–SW and is parallel to the Delhi–Aravalli fold belt, which, however, corroborates well with the findings of the present study.

The moment magnitude of both the events, that occurred on 03 July 2020 and 17 December 2020, was found to be Mw 4.6 and 4.2, respectively. Furthermore, the average seismic moments of the respective events were obtained at 8.28E+15 Nm and 2.29E+15 Nm with stress drops 166 and 138 bars, respectively. Such high stress drops for the intraplate earthquakes corroborate well with the earlier findings for the small- to moderate-sized earthquakes (Allmann and Shearer, 2009; Kumar et al., 2014; Sairam et al., 2018). Usually the stress drop, on occurrence of an earthquake, is influenced by several factors including source geometry, faulting style, tectonic deformation rate, stress regime, rock composition, and focal depth (Kaneko and Shearer, 2014; Goebel et al., 2015). The variation in the stress drop is, however, attributed to the source model assumptions (Huang et al., 2016). In a study, ∼26 bar stress drop was estimated for a relatively small magnitude earthquake (M_w 3.5) that had occurred close to the northeastern part of the National Capital Territory of Delhi. Shearer et al. (2006) and Pandey et al. (2020) suggested even higher stress drops for the events of the intraplate origin with the focal mechanism showing relatively a dominant normal faulting.

Coulomb stress modeling with a fixed receiver fault plane (Figure 7), that is a correlation between the loading lobes, indicated a close interaction between the earthquakes. The Coulomb stress variation map computed using the July 2020 event as a source fault and the December 2020 event as a receiver fault showed that the previous event transferred stress which promoted the failure of the later event (Figure 7A). The December 2020 fault rupture is located in a loading zone which is clearly identified by a positive stress change in this portion. The slip on the July 2020 fault probably increased the stress at the receiver fault and triggered the December 2020 earthquake. Evidently, the later event was potentially triggered by a stress change of 0.002–0.004 bars. The observed alteration in the transferred Coulomb stress resulting from the July 2020 event is found to be causative for triggering the failure of the December 2020 event. The link between the two earthquakes is mainly related to the existence of faults with different types of motion in Aravalli belts, including the N15E-trending normal fault.

The Aravalli range that witnessed the two recent significant tremors is largely composed of an Archean–Proterozoic crystalline basement overlain with fold belts that witnessed multiple episodes of metamorphism and deformation during the Precambrian (Roy et al., 1995; Gupta et al., 1997; Sinha-Roy et al., 1998). The major tectonic features in the region show a predominant NE–SW trend that probably controls the regional tectonics (Gupta et al., 1997). Many of these stratigraphic features that are tectonized (Chetty, 2017) signify the mapped NE–SW lineaments in the region. In the vicinity of the source zone, the prime deformation is observed along the NE-trending ridge fault system, i.e., DHR and MDSSF, characterized by a normal fault with the strike–slip component (Prakash and Shrivastava, 2012). The Sohna and Mathura faults in the region, located to the east of the source zone, were also found to be quite active (Chouhan, 1975; Pandey et al., 2020). The activity might have upsurged the transfer of stresses and their accumulation along the minor active faults/lineaments mapped in the region. It corroborates well with the Coulomb stress pattern (Figure 7A). The seismicity analysis performed using the earthquake catalog complete for magnitudes of 3.0 and above revealed that the epicentral region displayed a reduced rate of seismic activity in comparison to the adjacent regions. The b-value is a measure of the frequency–magnitude distribution of earthquakes in a given region, with higher values indicating a higher likelihood of smaller magnitude earthquakes. The b-value of 1.1 ± 0.12 from the present study suggests that the region is less likely to experience significant seismic activity and is rather more prone to the smaller-magnitude earthquakes (Figures 2C, D). It is worth noting that the seismic activity in the study region during 2000–2022 was recorded within the magnitude range <4.0.

The epicentral zone of the recent activity falls largely in the seismic zone IV (seismically a hazardous zone) of the seismic zonation map of India (BIS, 2002). It is indicative of frequent occurrence of low- to moderate-magnitude earthquakes (two or more events per year) in the region, which is quite intriguing as the region is characterized by cratonized terrain with no major tectonic activity in the past. Kumar and Pandit (2020), however, showed earthquake activity to the southwest of the epicentral zone in recent years and suggested the genesis of such small events may possibly cause stress build-up due to northward drifting and collision of the Indian plate. The optimum solution of the stress modeling indicates a dominant normal stress regime in the source region (Figure 5). The confidence level of principle stresses (Figure 5B) further demonstrates the accuracy of the inverted focal mechanisms with a high stress ratio at 0.7 ± 0.05 (Figure 5C). The distribution of fault planes in Mohr’s circle implies a low shear stress (Figure 5D). About 80% of the fault plane solutions located in the lower Mohr’s semi-circle indicate the NE orientation of the source having a low friction coefficient of ∼0.5. It is emphasized that only a few seismic events in the source region with low activity limited our analysis in achieving accurate inversion to estimate the most likely orientation of the source and stress regime.

Several researchers have suggested that a low friction coefficient in a weak zone along the fault may be causative for seismic activity (Pınar et al., 2010; Fojtíková and Vavryčuk, 2018; Abdelfattah et al., 2020). Furthermore, a low friction coefficient indicates low stress drops in the source, as high shear stress may not be accumulated in extensional tectonics (Abdelfattah et al., 2020). A low friction coefficient, as obtained from the present study, hence characterizes the seismicity in the region. However, the friction coefficient is rather an unstable parameter, and therefore an extensive dataset of accurate focal mechanisms would be required to achieve a better result. Geologically, the epicentral area is composed of gneissic granites and alluvium, which indicate a strong foliation formed on account of shearing stress (Roy, 1976; Chetty, 2017). The shearing of platy minerals like, mica schist, biotite, and chlorite in the vicinity region along the foliation direction might have reactivated the faults/lineaments in the region.

The source parameters of the 03rd July 2020 and 17th December 2020 earthquakes that occurred in the vicinity of Delhi have been estimated using well-recorded data by the seismic network of the Delhi region. The focal mechanism solutions of the events were constructed using the waveform inversion approach, which suggested a normal faulting mechanism with a slightly strike–slip component. The optimum solution reflects a contemporary shear stress regime characterized by major compressive stress (σ1) with 52° plunge and trending 27° NE, intermediate compressive stress (σ2) with 38° plunge and trending 187° SW, and minor compressive stress (σ3) with 10° plunge and trending 285° SE, which corresponds to a dominant normal stress regime that was ascertained using the stress field inversion. The stress regime to the northern and eastern parts of the source zone of the M4.6 earthquake was found to be increased by 0.004–0.008 bars, surpassing the threshold value for triggering an earthquake. Such enhancement in the pressure changes due to the July 2020 event probably triggered the earthquake that occurred on December 2020. However, the stress variations in the off-fault region that is in the peripheral areas adjacent to the main fault, namely, along the Sohna fault and MDSSF fault, due to July event of 2020 were found insignificant. We suggest that the events had occurred along the pre-existing NE–SW lineament mapped in the vicinity and the stress drops for both the events were found to be 166 bar and 138 bars, respectively, that fall well within the range accepted for the intraplate earthquakes. Furthermore, these values of stress drops were also found to be in accordance with the low static friction (∼0.5) for the study region. The moderate stress drops clearly indicate a potential relationship between the local causative source and large-scale regional tectonics in the vicinity region, viz., MDSSF and DHR. The NE–SW-trending lineaments mapped in the region, therefore, appear to be the causative source of the 2020 earthquakes. We emphasize that the seismogenesis of small-to-moderate magnitude earthquakes in the study region may be due to the reactivation of hidden old sutures/weak zones (mapped as lineaments on the surface) owing to stress transfer from the active regional tectonics.