Spatial distribution of bedrock level peak ground acceleration in the National Capital Region of India using geographic information system

Abstract The probabilistic seismic hazard analysis is conducted for the seismically active National Capital Region (NCR) of India, one of the most densely populated regions in the world, to quantify the peak ground acceleration (P GA) at bedrock level by considering seismogenic source characteristics, smoothened grid seismic zoning, magnitude recurrence model, uncertainties, and suites of suitable ground motion prediction equations. The spatial distribution of the resulting P GA value is presented using GIS for the entire NCR of India for 10% and 2% probability of exceedance in 50 years and 100 years. The estimated P GA values are found to vary in the range of 0.12 to 0.37 g and 0.15 to 0.40 g for a 2% probability of exceedance in 50 years and 100 years, respectively. And the P GA values range from 0.07 to 0.33 g and 0.08 to 0.35 g for the 10% probability of the exceedance in 50 years and 100 years, respectively. The results reveal that the eastern region of the study area falls under high seismicity due to its proximity to the Himalayan zone, a highly seismically active region. High P GA values are also observed for south-eastern regions of the NCR of India due to the vicinity of the Great-Boundary fault, Mathura fault, Moradabad fault.


Introduction
India has its diversity in terms of geological, geomorphological, and socio-economic conditions that make it vulnerable to a multi-hazard environment. The northern part of India is one of the most seismically active zones of the world and tectonically complex due to the presence of the Himalayan young fold mountains (Gupta and Gahalaut 2014). The dynamics of the Himalayan Mountains result in the formation of major faults and structural discontinuities in the region (Bungum et al. 2017). The entire Himalayan region has faced extensive damages in the past due to the occurrence of devastating earthquakes, namely 1897 Shillong (M 8.7), 1905Kangra (M 8.6),1934Bihar (M 8.4), 1950, 2015 Nepal earthquake (M 8.1), 1991Uttarkashi (M 6.8), 1999Chamoli (M 6.4),2005 Muzaffarabad (M 7.6), 2011 Sikkim earthquake (M 6.9). The consequences of earthquakes are not limited only to structural collapses, but also cause lateral sliding, landslides, soil liquefaction, and fire, etc. . The extent of damages in a defined space and time may vary from local to regional level in terms of structural, socio-economic, or environmental losses (Erdik 2017;Agrawal et al. 2021).
The assessment of seismic hazards is vital in predicting the ground shaking level for the reduction and management of risk as well as for sustainable urban planning and emergency response planning. Seismic hazard analysis (SHA) is an effective tool to present the potential and possible damaging earthquake scenarios. The probabilistic seismic hazard assessment (PSHA) approach is based on the total probability theorem, considering the uncertainties in terms of magnitude, location, and intensity at the site (Cornell 1968;Bommer 2002;Raghukanth et al. 2011;Wang 2011).
Determination of seismic hazard is associated with epistemic or aleatory uncertainties (Stepp et al. 2001). Epistemic uncertainty is primarily modelrelated uncertainty, and its range can be varied based on the alternate models considered for the source description and application of statistical methods (Crowley et al. 2005). The errors in the estimated parameters can be reduced by using a large data set. Aleatory uncertainty is inherent to the observed natural process, namely the source characteristics of the earthquake, propagation of seismic waves through the path, and site amplification factors (Baltay et al. 2017). In such cases, PSHA is very useful to deal with both types of uncertainties. Implementation of different source models and attenuation relationships can be used to reduce the uncertainties. Many seismic hazard studies have been conducted at the national scale (Parvez et al. 2003;NDMA 2010;Nath and Thingbaijam 2012;Sitharam and Kolathayar 2013) and some at regional scales (Iyengar and Ghosh 2004;Raghukanth and Iyengar 2006;Sitharam and Anbazhagan 2007;Raghukanth 2010;Raghukanth and Kavitha 2014 ;Anbazhagan et al. 2015;Lindholm et al. 2016;Sarkar and Shankar 2017;Anbazhagan et al. 2019;Keshri et al. 2020). Several studies have been conducted for seismic hazard assessment using the combination of the traditional approach and geospatial methods (Mohanty et al. 2007;Ahmad et al. 2017;Deligiannakis et al. 2018).
Over the past few decades, the NCR of India has shown explosive growth of population causing increased vulnerability of the region (Agrawal et al. 2021). According to the Seismotectonic Atlas of India, the regions of western Uttar Pradesh, Uttarakhand, the eastern part of Haryana, Delhi, and the NCR of India are in the moderate to high-risk seismic zones (GSI 2008). There is an urgent need for an updated seismic hazard map at the micro-level for the entire NCR of India by incorporating new seismic data and improved methodologies. This will help to examine the preparedness strategies to take timely measures to reduce the impact and consequences of any future earthquakes in and around the region Singh et al. 2013).
In the present study, a detailed probabilistic seismic hazard assessment is performed for the entire National Capital Region (NCR) of India. The seismotectonic model for the study area has been developed by considering Shiv Nadar University, Uttar Pradesh (28 31ʼ 29ʼ ʼ N, 77 34ʼ 28ʼ ʼE) as a center of a circle of 350 km radius (Sitharam and Anbazhagan 2007;. The entire region is divided into two seismic zones: i) the Himalayan region and ii) the NCR of India (Iyengar and Ghosh 2004;Anbazhagan et al. 2015;Anbazhagan et al. 2019). The major steps involved in the estimation of seismic hazard parameters are identifying the seismic sources, preparing the earthquake catalogue, homogenizing earthquake magnitudes, data declustering, scrutinizing the completeness of the data, determining seismic parameters using G-R (Gutenberg and Richter 1944) relationship, developing individual fault recurrence relationships, estimating strong ground motion, and selecting suitable attenuation relationships of strong motion. The peak ground acceleration (PGA) values at the bedrock level have been quantified by dividing the whole region into grids of size 10 km Â 10 km. The results are presented in the form of seismic hazard maps for a 2% probability of exceedance in 50 years and 100 years with a return period of 2475 and 4950 years, respectively using the GIS tool. The maps are also developed for a 10% probability of exceedance in 50 years and 100 years with 475 and 995 years return period, respectively. The estimated PGA values are compared with the reported values with other literature and standard reports.

Description of the study area
The National Capital Region (NCR) of India consists of the entire National Capital Territory (NCT) of Delhi and 24 districts of its surrounding states like Uttar Pradesh, Haryana, and Rajasthan, as shown in Figure 1, covering an approximate total area of 54,984 km 2 and more than 46 million people are residing according to Census of India (2011). The NCR has a mixed infrastructure setting of urban, semi-urban, and rural regions. From north to south the region is drained by two major rivers in India, Ganga, and the Yamuna. The NCR is characterized by Indo-Gangetic alluvial plains of quaternary age (Chhabra et al. 2010). The prominent cities of the NCR are Delhi, Faridabad, Gurugram, Noida, Greater Noida, Meerut, Ghaziabad, Bulandsahr, Karnal, Panipat, and Alwar. According to the seismic zonation map of India, the NCR falls under seismic zone-IV (IS-1893 2016). The ongoing developmental infrastructure projects in the region will further result in creating new challenges during seismic events. The seismic resilience of the region can be achieved by a detailed seismic hazard assessment of the region, which is very important for future and existing developmental projects.

Geological setting
The NCR is characterized by the presence of the Indo-Gangetic plains of quaternary age in the north and east, the extension of the Great Indian Thar desert in the west, and the Aravalli ranges in the south (Mittal et al. 2013). The river Yamuna passes through the eastern part of the Delhi region. The terrain of the region is generally flat except in the central part of the NCR. The formation of a low NNE-SSW trending Delhi ridge is observed. The characteristics of rock formation in the region can be classified from quartzite to quaternary. The geomorphological characteristics of the entire region can be distinguished as (i) the erosional surface which formed the top of denudation hills, (ii) the older alluvial plain (Pleistocene), and (iii) the depositional younger alluvial plain of river Yamuna (Holocene) (Thakur 2013). The rock system of the region has experienced several folding and has undergone multiple processes of metamorphism (Naha et al. 1984). Many faults are found to be trending along NNE-SSW, NE-SW, and WNW-ESE. The stratigraphic succession of the NCR is dated from Precambrian to recent age. Rapid urbanization led to such a dynamic transformation of the geology of the region.

Seismotectonic setting
The high seismicity of the region is due to the presence of many faults and different geological structures. The major faults and lineaments are georeferenced and digitized using ArcMap 10.5 in the ArcGIS environment (ESRI 2016) based on the seismotectonic atlas, prepared by the Geological Survey of India, and are presented in Figures  2 and 3. In the north-eastern part of the study area, the tectonic belt of young fold mountains of the Himalayas is located and in the southern part, the Proterozoic Delhi fold belt and gneissic batholithic complex are predominant. Highly jointed and Alwar quartzites folding are also observed in surrounding areas and regions to the south of Delhi (Mohanty et al. 2007). The region also shows the presence of the Moradabad fault zone trending along NE-SW direction forming a boundary between two tectonic sub-provinces namely Delhi-Moradabad and Kasganj -Ujhani province (Kumar and Narayan 2020). This tectonic feature is described as a tectonic boundary between the Delhi fold belt and the Vindhya. The northern and south-western boundaries of the Delhi-Moradabad tectonic regime are formed by the Proterozoic Delhi fold belt and the Frontal Fold zone of the Himalaya, respectively (Srivastava et al. 2018). The Delhi-Hardwar ridge can be characterized as the extension of the north-  eastward ridge of the Delhi fold belt towards the Himalayas (Chouhan 1975). The Delhi-Haridwar ridge and Monghyr-Saharsa ridge forms the western and eastern margin of the Ganga basin, respectively (Srivastava et al. 2018). In the south-eastern part of Delhi lies a major Rajasthan Great Boundary Fault (GBF) characterized by well-defined Chittaurgarh-Machilpur lineament due to the presence of distinct geomorphic units on either side (Sharma et al. 2003). The presence of the GBF and the Mahendragarh-Dehradun subsurface fault (MDSSF) contributes to the seismicity of the region. The major tectonic features are identified based on previously reported studies (Iyengar and Ghosh 2004;Prakash and Shrivastava 2012) and are presented in Figure 2. Therefore, it can be concluded that the study area is dominated by many seismotectonic features such as faults, structural discontinuities, and lineaments.
The characterization of seismic source and its behavior is dynamic and varies with the surrounding seismotectonic environment. The study area is divided into two main regions Himalayan region (Region 1) and the NCR of India (Region 2). In the Himalayan region, the events are associated with Main Boundary Thrust (MBT), Main Central Thrust (MCT), and Main Frontal Thrust (MFT) (Mohanty et al. 2007). In the NCR of India, the events are concentrated along Mahendragarh-Dehradun Fault (MDF), Moradabad Fault (MF), and Great Boundary Fault (GBF). For the identification of seismic sources, 19 major faults and lineaments are identified within a 350 km radius from the center of the study area, as shown in Figure 3a, based on past seismic events and seismicity.

Seismicity
The NCR lies in the proximity of the Himalayas as well as the seismotectonic setting of the region discussed in the previous subsection which indicates the presence of major faults and ridges. The impact of historical and significant past earthquakes which caused great damage to Delhi and its surroundings also corroborates the high seismicity of the region. The 1720 Sohna earthquake (M 6.5) destroyed the walls of the fortress (Red Fort) and many houses in Delhi, and it was followed by 4 to 5 aftershocks per day for 40 days and occasional tremors for 4 to 5 months (Oldham 1883). The occurrence of the 1803 Mathura earthquake (M 6.8) caused damage to Qutab Minar (Oldham 1883) and the 1960 Gurgaon earthquake (M 6.0) caused minor property damages and injuries to about 50 people in Delhi ). Damages to the minarets of Jama Masjid were caused during the 1994 Delhi earthquake (M 4.0) (Iyengar and Ghosh 2004). The geological setting of the NCR can sustain large and amplified ground shaking. The region has been affected many times from the far seismic sources as well as from local seismic activity. It was observed that the damages caused to the buildings and other man-made structures of NCR were severe although the seismic sources were very far away (Mohanty et al. 2007). The tremors of the 1999 Chamoli earthquake (M 6.8), 2011 Pakistan earthquake (M 7.4), 2011 Sikkim earthquake (M 6.9), 2015 Nepal earthquake (M 7.8), and many other earthquakes were strongly felt in the region. Recently, in the year 2020, the region has experienced more than ten earthquakes (M > 4.0) and its epicenters are in and around the NCR. As the region is densely populated and has many critical infrastructures and world heritage monuments, an increase in seismic activity poses a threat to both life and property.

Preparation of earthquake catalogues and their homogenization
The preparation of the earthquake catalogue is considered the most important step for PSHA. A comprehensive earthquake catalogue is prepared for the study area considering the radius of 350 km centered on Shiv Nadar University (28 31ʼ 29ʼ ʼ N, 77 34ʼ 28ʼ ʼ E). Earthquake data have been collected from national as well as international databases including the National Centre for Seismology (NCS), Bhukosh-Geological Survey of India (GSI), India Meteorological Department (IMD), Institute of Seismological Research (ISR), International Seismological Centre (ISC-GEM Earthquake Catalogue 2021), United States Geological Survey National Earthquake Information Centre (USGS NEIC), and many peer-reviewed research articles (Oldham 1883;Chouhan 1975;Naha et al. 1984;Sharma et al. 2003;Iyengar and Ghosh 2004;Mohanty et al. 2007;Prakash and Shrivastava 2012;Mittal et al. 2013;Thakur 2013;IS-1893Srivastava et al. 2018;Kumar and Narayan 2020). A total number of 2034 events are collected from 1802 to 2020 with detailed information about the date and time of occurrence, magnitude, hypocentral depth, and location with specified coordinates (Figure 3b).
The earthquake catalogue comprises earthquakes recorded on different magnitude scales such as surface-wave magnitude (M s ), moment magnitude (M w ), body-wave magnitude (M b ), and local magnitude (M L ) with a majority of events recorded on the moment and local magnitude scales. Most of the data points are originally recorded in terms of moment magnitude and local magnitude. The magnitudes are homogenized to a unified and common magnitude scale i.e., moment magnitude (M w ) by using the orthogonal regression approach (Orthogonal Distance Regression) (Das et al. 2012;Lindholm et al. 2016). Limited data points are available for the body-wave magnitude scale which are found to be common with either M L or M s scales. There are 297 and 53 data points for the orthogonal regression between M L and M s ( Figure  4a) and M w and M L (Figure 4b), respectively. So, the following relations are obtained based on orthogonal regression and used for conversion of M s and M L to M w (Equations (1) and (2)).

Declustering of earthquake catalogues
The earthquake catalogue obtained for the study region comprises of mainshocks along with their foreshocks and aftershocks. Both foreshock and aftershock depend spatially and temporally on the mainshock. As mutually exclusive events are considered in the probabilistic seismic hazard analysis of a region, it is important to remove the dependent events from the prepared earthquake catalogue. In the present study, the dependent events are separated from the main events by declustering of earthquake data based on Reasenberg declustering algorithm (Reasenberg 1985;Wiemer 2001;Zhuang et al. 2002) using MATLAB based opensource software ZMAP (Wiemer 2001;Vipin et al. 2013;Telesca et al. 2016;Sinha and Sarkar 2020;Borah and Kumar 2021). The parameters associated with the Reasenberg declustering algorithm are minimum time interval (s min ) and maximum time interval (s max ) for the occurrence of the next earthquake event with a certain probability p. Relationship between these three parameters used in the declustering algorithm is expressed in Equation (3).
Where DM is the difference between maximum and minimum magnitude, t is the time interval. A total of 14 clusters of earthquakes are identified and 90 dependent events out of 2034 are found. A map of the control region showing 1944 mutually exclusive seismic events after the declustering of earthquake data is shown in Figure 5.

Analysis of completeness of seismic data
The homogenized earthquake data is further analyzed to check the completeness of the data set in the given time interval based on the methods proposed by Stepp (Stepp 1972;Khan and Kumar 2018). The data is grouped into several magnitude classes and each class is modeled as a point process in time to obtain the variance of the sample mean. In this procedure, the number of earthquakes occurring per decade is classified into magnitude ranges of 3.0 M w < 4.0, 4.0 M w < 5.0, 5.0 M w < 6.0, M w ! 6.0. The earthquake sequence is modeled by assuming Poisson's distribution to evaluate the variance of the sample mean. If k 1 , k 2 , k 3 , … , and k n are the number of earthquakes per unit time interval, then an unbiased estimate of the mean rate per unit time interval of the sample will be, And its variance will be given by, Where n is the number of unit time intervals (1 year). Now the standard deviation of the mean number of annual events will be calculated as Where T is the sample length in years (Stepp 1972). The standard deviation versus time intervals is plotted for every magnitude class as shown in Figure 6a and 6b. If the data for a given magnitude interval is completed, then the plotted points will give a straight line parallel to the slope of 1/ͱT. Fig. 6a and b reveal that completeness of 20, 40, 50, and 80 years for the region I and 20, 60, 90, and 100 years for region II from the year 2020 in the moment magnitude classes of 3-4, 4-5, 5-6 and > 6, respectively (Table 1).

Estimation of seismic parameters a, b and M C
After checking the completeness of the earthquake catalogue of the study area temporally, it is also necessary to evaluate the magnitude of completeness (M C ). M C is defined as the lowest magnitude at which 100% of the seismic events in a space-time volume are identified (Wiemer and Wyss 2000). Different methods are used in the past for the estimation of M C values, such as the maximum curvature method, fixed minimum magnitude observed (M min ), the goodness of fit (M min90 and M min95 ), the best combination of M min90 and maximum curvature, entire magnitude range, and bootstrap method ). In the present study, the values of M C for regions I and II are determined using the maximum curvature method (Figures  7a and 7b) and the values are presented in Table 2. The obtained Mc values for regions I and II are 3.0 and 2.1 respectively. So, taking that into account, 3.0 Mw has been adopted as minimum magnitude (m 0 ) for both the regions for further analysis (Anbazhagan et al., 2015). The seismic activity of the region is estimated in terms of the recurrence rate of earthquakes for the different range of magnitudes using an empirical relationship in Equation (7) proposed by Gutenberg and Richter (1944). It is one of the simplest and widely used approaches for the quantification of seismic parameters a and b of a region (Sitharam and Sil, 2014;Anbazhagan et al. 2015). N denotes the number of earthquakes greater than or equal to magnitude M w per year, a represents seismic activity and b is known as a tectonic parameter. Based on the magnitude of completeness seismic parameters are estimated using GR recurrence law. The total number of earthquakes greater than magnitude 3.0 is 638, i.e., 461 in the region I and 177 in region II. The value of a and b are determined for both regions and it is compared with the other reported past studies as shown in Table 2 and Figure 8.
It is observed that the value of b is lower for the Himalayan region as compared to the NCR of India, and this is because the seismic parameter b is related to the occurrence of large to smaller seismic events and the Himalayan region experiences earthquakes of higher magnitude as compared to the NCR, where lower magnitude earthquakes are predominant.

Fault deaggregation
19 major faults are identified in the study area and any site in the study area will be subjected to ground shaking from any of these seismic sources. The recurrence relationship developed earlier is related to the entire study region and it is not specific for a particular fault. To differentiate seismic activity rates among different faults and to distinguish near seismic sources from distant sources, the development of faultlevel recurrence relation is essential. The approach to developing a frequency magnitude relationship for individual seismic sources can lead to discrepancies in results due to the scarcity of seismic activity rates. Another approach is to calculate the value of b by determining the slip rate of individual faults, but the estimation of slip rate is also difficult as slip values of all the faults are not recorded for historical events. The regional recurrence of the NCR is disaggregated into individual fault recurrence by using the method suggested by Raghukanth and Iyengar (2006), based on the heuristic principle of conservation of seismic activity stating the regional seismicity defined  as the number of earthquakes per year with m > m 0 , should be equal to the sum of the number of earthquakes occurring on individual faults (Raghukanth and Iyengar 2006).
In the determination of N i (m 0 ), two assumptions are important in terms of fault length weighing factor (a i )and earthquake events weighing factor (v i ). The first assumption is that the longer the fault, the higher is the capacity of producing a greater number of earthquakes of smaller magnitude m 0 than the faults of shorter length. Therefore, N i (m 0 ) is proportional to the length of the fault and this observation leads to the weighing factor a i ¼ L i = P L i , L i is equal to the length of the fault i. It is also observed that irrespective of the length of the fault, a shorter fault may be more active to produce smaller magnitudes of earthquakes than a longer fault. The second assumption is that the future seismic activity will remain similar to past activity, at least for the short term. So, seismic activity needs to be related to the number of past events associated with each fault in the earthquake catalogue as depicted in Figure 9. Another weighing factor, v i , is defined as the ratio of the number of earthquakes associated with fault i to the total number of earthquakes in the region. The final recurrence relation for a fault i is given as In the PSHA of NCR, 19 seismic sources have been identified, the details of each fault in terms of length of the fault, shortest hypocentral distances from Shiv Nadar University, Delhi-NCR, India, the maximum magnitude of earthquake associated with each fault, number of earthquake events occurred in the proximity of each fault and weighting factors for each fault are given in Table 3.
To obtain the maximum possible earthquake (M max ) that each fault can produce, the conventional incremental method is adopted (Gupta 2002;Anbazhagan et al. 2015). The observed maximum magnitude (M max obs ) associated with each seismic source is increased by 0.5 units based on the obtained b values (Anbazhagan et al. 2015;Anbazhagan et al. 2019). The details of each fault in terms of the total length and M max with an increment of 0.5 units higher than the magnitude reported in the past are shown in Table 3 (Raghukanth and Iyengar 2006).
The b value of each fault is equal to the b value of the region in which the fault is located. The equation obtained is: m 0 is the value of threshold magnitude, m max is the maximum potential magnitude of the fault i and b ¼ 2.303 b. The final individual fault level recurrence relation is shown in Figure 10a and 10b.

Probabilistic seismic hazard assessment
The entire study area is divided using a grid of 10 km Â 10 km with a total of 524 grid points, as shown in Figure 11, for the estimation of seismic hazard at the bedrock level for the NCR of India. In this study, every grid point is selected as a site of interest. With the help of ArcGIS 10.5, the shortest and longest distances from each grid point to all active faults present in the study area are estimated. MATLAB programming is used for all computational purposes and probabilistic seismic hazard assessment to calculate the ground motion parameters. The number of earthquake events that occur on a fault is assumed to follow a stationary Poisson process, and the probability of a random variable Y exceeding level a, in a period T year at a site can be represented by the annual rate of exceedance m a .
The expression for m a is given as K is the number of faults present in the study region, p M m ð Þ is the probability density function of magnitude, p R jM rm ð Þ is the probability density function of hypocentral distance and the conditional probability of exceedance of the ground motion parameter is denoted by PðY > a j m, rÞ: The return period for the corresponding ground motion value is given by the reciprocal of the annual probability of exceedance.

Uncertainty in magnitude
Seismic sources present in a region are capable of producing earthquakes of different sizes and smaller magnitude earthquakes have a higher frequency than the earthquake of greater magnitude. A seismic source is assumed to cause an earthquake of magnitude as a random variable distributed in the range of M min (m 0 ) to M max (m max ), following Equation (13) for the recurrence relation.

Uncertainty in hypocentral distance
Another uncertainty exists in terms of the source-to-site distance in the PSHA. It is assumed that all points on an active fault have equal potential to get ruptured and to produce earthquakes of significant magnitude. Therefore, the hypocentral distance R is treated as a random variable and it depends on the relative orientation of the fault with respect to the site (Kramer 1996). In this study, the faults in the study region are the linear sources and the occurrence of a seismic event can take place at any point along the fault. The probability distribution function of the hypocentral distance R can be expressed as Equation (15).
Where L is the length of the fault and r min is the minimum hypocentral distance from the site. Equation (15) is applied for the calculation of the probability density function of uncertainty in hypocentral distances for all the seismic sources present within the study region and is applicable when the fault is located completely to one side of the site of interest (Kramer 1996). In most cases, the faults extend on both sides of the seismic source, and the probability for the total fault is estimated by the summation of the product of conditional probabilities of both sides and a fraction of corresponding sides.

Attenuation relationships
An attenuation relationship, also called ground motion prediction equation (GMPE), can be derived by regression analysis of past recorded earthquake data of the region and it depends on several factors like the geometry of the seismic source, the magnitude of earthquake event, corresponding source to site distance, path characteristics, and local site conditions. The knowledge of site-specific attenuation relation is very important for evaluation of ground motion parameter (Singh et al. 1996), but in case of limited data, one can use the previously developed attenuation models for the region or models based on similar tectonic features as shown by various researchers such as Nath and Thingbaijam (2012), Anbazhagan et al. (2015) and Anbazhagan et al. (2019). The development of attenuation relationship for the entire NCR of India has limitations due to lack of strong-motion data, correct epicentral distances, and range of earthquake magnitude. Very few ground motion prediction equations are available for Northern India (Sharma 1998;Iyengar and Ghosh 2004;Rao and Rathod 2014) and Himalayan region (Aman et al. 1995;Singh et al. 1996;NDMA 2010;Joshi et al. 2012;Ahmad and Singh 2016, etc.). However, there is no specific GMPE derived for the entire NCR. In the present study, six different ground motion equations suggested by Bajaj and Anbazhagan (2019), Rao andRathod (2014), Joshi et al. (2012), Iyengar and Ghosh (2004), Ambraseys andBommer (1991), andCornell et al. (1979) are used. The average of the ground motion parameter values, calculated from the following six GMPEs, are used to determine peak ground acceleration (PGA) values of the NCR of India. Figure 12 shows attenuation relations for PGA versus hypocentral distance at magnitude 7.0 represented for all six GMPEs and their average.
For individual ground motion prediction equation (GMPE), seismic hazard curves at Shiv Nadar University of the study region are shown in Figure 13a through Figure  13f. The GMPEs used in this study are expressed below.
The attenuation relation (GMPE 1) developed by Bajaj and Anbazhagan (2019) for the Himalayan region is expressed as ln Y ¼ a 1 þ a 2 ðMÀ6Þ þ a 3 ð9ÀMÞ 2 þ a 4 ln R þ a m ln RðMÀ6Þ þ a 7 R þ r ln Y is the logarithm of ground motion, M is magnitude and R is hypocentral distance. The value of standard deviation, r is equal to 0.817. The regression coefficients are a 1 , a 2 , a 3 , a 4 , a m (for M w < 6.0 and R < 300, a m in Equation (16) is a 5 , otherwise a m is a 6 ) and a 7 and corresponding values are 1.071, À0.257, À0.184, À0.479, 0.078 (a 5 ), 0.076 (a 6 ), and À0.0085, respectively. Rao and Rathod (2014) performed regression analysis using 1,650 data to obtain the attenuation relation (GMPE 2) for rock site condition as a function of magnitude and distance given in Equation (17). log 10 ðPGAÞ ¼ c 1 þ c 2 ðM W À6Þ þ c 3 ðM W À6Þ 2 þ c 4 R þ c 5 log 10 R þ log 10 e The unit of PGA is g, M w is moment magnitude, R is the hypocentral distance expressed in km, and the values of regression coefficients c 1 , c 2 , c 3 , c 4, and c 5 are 0.66845, 0.49908, À0.04665, À0.00257, and À0.85968, respectively. The value of the standard deviation of random error r (log 10e) is 0.1118. Joshi et al. (2012) derived ground motion prediction equation (GMPE 3) by using data from 8 stations installed in the Pithoragarh region of Kumaon Himalaya. The ground motion prediction model is expressed in Equation (18) and the range of hypocentral distance in this study varies between 10 km to 100 km. ln ðPGAÞ ¼ a 1 þ a 2 M W þ a 3 ln r hypo þ a 4 ln ðr epi þ 15Þ The value of PGA is in Gal, and the coefficients a 1 ¼ À5.8, a 2 ¼ 2.62, a 3 ¼ À0.16, a 4 ¼ À1.33 and r ¼ 0.42. Iyengar and Ghosh (2004) proposed attenuation relation (GMPE 4) for Delhi City, which is used in the present study, is as follows in Equation (19). log 10 ðyÞ ¼ C 1 þ C 2 MÀB log 10 ðr þ e C 3 M Þ þ Pr Where C 1 ¼ À1.5232, C 2 ¼ 0.3677, C 3 ¼ 0.41, B ¼ 1.0047 and r ¼ 0.2632. The value of P equals to 0 for the median and 1 for 84 percentile value of PGA.   The ground motion model (GMPE 5) proposed by Ambraseys and Bommer (1991) can be represented as Equation (20).

Determination of temporal uncertainty
Earthquakes occur randomly from time to time, and therefore, temporal uncertainties must be considered in the assessment of seismic hazards. The probability of occurrence of an earthquake event in a given period can be calculated considering the distribution of earthquakes with time. The number of seismic events that occur on a fault and the temporal occurrence of earthquakes can be determined with the help of the stationary Poisson model (Iyengar and Ghosh 2004). The probabilities of exceedance of specific ground motion parameter a for period T is expressed as Equation (22) that represents the relation between the mean annual rate of exceedance (m a ) and period (T) The values of annual mean exceedance m a are calculated, in this study, for 2% and 10% probabilities of exceedance in 50 and 100 years.

Seismic hazard mapping and comparison
The seismic hazard curves are prepared for all the 19 faults and then combined to quantify the value of aggregate hazard at the sites of interest. The entire NCR is divided into 524 grids of a grid size of 10 km x 10 km. After estimating the value of the mean annual rate of exceedance (m a ) at each of these sites, PGA at rock level for 2% and 10% probability of exceedance in 50 and 100 years are computed using the six most suitable attenuation relationships and their average. The seismic hazard map for NCR of India is generated and presented in Figure 14a-d using the average PGA value estimated at each site corresponding to 2% and 10% probability of exceedances for both 50 years and 100 years. The PGA values estimated using the average values of all six GMPEs are 0.12 to 0.37 g, 0.15 to 0.40 g, 0.07 to 0.33 g, and 0.08 to 0.35 g for 2% probability of exceedance in 50 years, 2% probability of exceedance in 100 years, 10% probability of exceedance in 50 and 10% probability of exceedance in 100 years, respectively (Mehta and Thaker 2020). The eastern and south-eastern portions of the study region are exhibiting high PGA values (Figure 14a-d). A similar spatial pattern is also reported by some previous studies like Sharma et al. (2003),   Sharma and Wason (2004), Sarkar and Shankar (2017) for the same study area. The PGA maps reveal that the south-eastern part of the study area is more affected due to the presence of Great-Boundary fault (GBF), Mathura fault, and Moradabad faults. The eastern part of the study area also indicated high PGA values due to its vicinity to the Himalayan zone, which is associated with high seismicity due to the presence of numerous faults and lineaments (Figure 14a-d).
The results obtained from the present study and its comparison with the results reported on the seismic hazard studies carried out for the Delhi region (Iyengar and Ghosh 2004;Rao and Rathod 2014;NCS-MoES 2016;Sarkar and Shankar 2017) are given in Table 4. The comparison with Sharma et al. (2003) cannot be tabulated as they reported the PGA values for 20% probability of exceedance in 50 years. The values obtained in the present study are on the higher end as compared to the previous study. The PGA value estimation depends upon the factors like selection of appropriate GMPE, seismic sources, and associated seismic parameters. The use of under-predictive GMPEs may be the reason for lower PGA values by the previous studies.
The PGA values estimated using the average values of all six GMPEs are 0.12 to 0.37 g, 0.15 to 0.40 g, 0.07 to 0.33 g, and 0.08 to 0.35 g for 2% probability of exceedance in 50 years, 2% probability of exceedance in 100 years, 10% probability of exceedance in 50 and 10% probability of exceedance in 100 years, respectively. The seismic hazard maps for NCR developed in the present study will serve as a significant input in sustainable urban planning, city planning for critical and lifeline infrastructures, preparedness, and mitigation planning, designing, and constructing earthquake-resistant structures as well as in hazard and risk mitigation and management.

Conclusions
A comprehensive probabilistic seismic hazard assessment for the NCR of India at the bedrock level is carried out in this study. The faults and geological discontinuities are identified with the help of Seismotectonic Atlas (SEISAT) maps and the seismotectonic map of the study area is prepared by using ArcGIS 10.5. The regional seismicity of the area is due to the presence of major 19 faults of different lengths and maximum potential magnitude. By using the earthquake data from 1802 to 2020, regional recurrence for both regions is obtained. The entire area is divided into 524 grids and for each grid PGA value at bedrock level is determined with the help of suites of suitable ground motion prediction equations (GMPEs).
For every grid point, hazard curves are also generated. The final seismic hazard mapping is based on the PGA values calculated based on the suite of suitable ground motion prediction equations. In this study, it is found that the resulting PGA value of the region varies in the range of 0.12 to 0.37 g, 0.15 to 0.40 g, 0.07 to 0.33 g, and 0.08 to 0.35 g for 2% probability of exceedance in 50 years, 2% probability of exceedance in 100 years, 10% probability of exceedance in 50 and 10% probability of exceedance in 100 years, respectively. This probabilistic seismic hazard analysis reveals that the eastern and south-eastern regions of the study area are more vulnerable to seismic hazards. The present study can be considered as a comprehensive assessment of seismic hazards for NCR India, but it has its limitations. Seismic hazard assessment itself consists of many inherent uncertainties, all the uncertainties could not be considered in the present study although major uncertainties in terms of magnitude and epicentral distances have been addressed. The study area consists of many major and minor seismic sources like faults, lineaments, thrusts, for the present work only major faults are considered based on literature review and historical earthquake occurrence. Limitations of the study also lie in the lack of geotechnical, geological, seismological, and other seismic hazard-related data. These findings will provide greater insights and input for the design and construction of earthquake-resistant structures and in decision-making for disaster risk reduction.