Probabilistic seismic hazard analysis using the logic tree approach – Patna district (India)

Peak ground acceleration (PGA) and study area (SA) distribution for the Patna district are presented considering both the classical and zoneless approaches through a logic tree framework to capture the epistemic uncertainty. Seismicity parameters are calculated by considering completed and mixed earthquake data. Maximum magnitude is calculated using three methods, namely the incremental method, Kijko method, and regional rupture characteristics approach. The best suitable ground motion prediction equations (GMPEs) are selected by carrying out an “efficacy test” using log likelihood. Uniform hazard response spectra have been compared with Indian standard BIS 1893. PGA varies from 0.38 to 0.30 g from the southern to northern periphery considering 2 % probability of exceedance in 50 years.

1 Introduction

Seismic hazard analysis is effective in presenting the potentially damaging phenomenon associated with earthquake. Earthquake disaster is not only associated with collapsing of structures due to ground shaking but also triggers fire, liquefaction, and landslide. Thus, it is indispensable to forecast the ground shaking level to serve the engineering needs in mitigating the risk associated with earthquakes. In India, moderate earthquakes (M_w<7) including Anjar 1956, Koyna 1967, Udaypur 1988, Uttarkashi 1991, and Chamoli 1999 have caused significant damage in last 10 decades (Nath and Thingbaijam, 2012). In addition, many great events (2015, Nepal earthquake) have originated from continent-to-continent collision. The Himalayan seismic gap (Bilham and Wallace, 2005) and thick soft soil sediments make the scenario more dangerous for cities close to the Himalayan region. Apart from this, improper planning, increase in population density, poor land use, and substandard construction practices in these cities magnify the prevailing seismic risk. Most of the existing seismic hazard maps are mainly on a macro-level for different parts of the Indian subcontinent and are not up to state-of-the-art knowledge in engineering seismology. For example, Khattri et al. (1984) developed a hazard map representing peak ground acceleration (PGA) for all of India with 10 % probability of exceedance in 50 years. Under the Global Seismic Hazard Assessment Program (GSHAP), Bhatia et al. (1999) presented a probabilistic seismic hazard analysis (PSHA) of India. Mahajan et al. (2010) delivered PSHA for the northwestern Himalayas. Recently, the National Disaster Management Authority (NDMA, 2010) and Nath and Thingbaijam (2012) have presented the PSHA map for all of India. The deterministic seismic hazard microzonation was proposed by Kolathayar et al. (2012). Further, Sitharam and Kolathayar (2013) developed the seismic hazard map of India using an areal source. In addition, Kumar et al. (2013) developed a deterministic seismic hazard analysis (DSHA) and PSHA map for the Lucknow region considering local and active seismic gaps. Additionally, the current Indian Standard (IS 1893, 2016) code consists of many constraints such as poor delamination of active seismic sources, lack of vulnerable sources to study, improper seismic hazard parameters which are not region-specific, and limited soil amplification consideration (Anbazhagan et al., 2014). Subsequently an updated seismic hazard map at the micro-level is essential for the cities near the Himalayan region, by considering new data, updated knowledge, and improvement in previous methodologies.

There are two types of uncertainties associated with hazard analysis. One is due to randomness of the nature of earthquake, wave propagation, and site amplification, named “aleatory uncertainty”, while the other is due to incomplete knowledge of earthquake process, named “epistemic uncertainty”. The former can be easily reduced by integrating the distribution of ground motion about the median (Bommer and Abrahamson, 2006) and the latter can be assessed using the logic tree approach. Gullu and Iyisan (2016) selected the ground motion prediction equations (GMPEs) for the logic tree based on the weighting factors that were incorporated with a Venn diagram of attenuation models regarding the experimenter's concern and expert's knowledge. Epistemic uncertainty is due to improper knowledge about the processes involved in earthquake events and algorithms used to model them. Hence, in this study, the logic tree framework has been used to reduce the epistemic uncertainty in the final hazard value calculation. In the absence of the appropriate region-specific models of wave propagation, ground motion prediction models are generally used to determine the hazard value. The uncertainty in GMPEs can be reduced by incorporating a logic tree in the hazard analysis study. Logic trees represent the various nodes that define the alternative input choices and each branch is assigned a weight factor that signifies the quantitative degree of likelihood assigned. To quantify the epistemic uncertainty, different branches of logic tree need to be considered, which is based on source models, regionalization of b value, determination of magnitude of completeness and maximum magnitude, and epistemic uncertainty in GMPE using the suitable representative approach.

In the present study, PSHA of Patna district (India) at the micro-level has been prepared along with the response spectrum by reducing the epistemic uncertainty. Patna lies at 250 km from the central seismic gap (Khattri, 1987) in the Himalayan region where huge devastation and destruction due to the 1803, 1934 Bihar–Nepal and 2015 Nepal earthquakes were reported. Similar to Bilham (2015), a large earthquake appears to be imminent in future due to rupturing of the main fault beneath the Himalayas because of the Nepal 2015 earthquake. Hence such studies need to be done for the cities that lie within the vicinity of the Himalayan region and the Indo-Gangetic Basin. Seismic sources and seismic events have been measured for a 500 km radius around the district centre as per Anbazhagan et al. (2015a). The a and b parameters have been obtained by taking into consideration the completed earthquake data using the Gutenberg–Richter (G–R) relationship and mixed data using methods proposed by Woessner and Wiemer (2005). The magnitude of completeness (M_c) is also calculated using nine methods proposed by Woessner and Wiemer (2005). Maximum magnitude has been determined using weighted mean considering three methods: increment factor on maximum observed magnitude, the method of Kijko and Sellevoll (1989), and regional rupture characteristics (Anbazhagan et al., 2015b). Ground motion prediction equations (GMPEs) have been selected from the 28 numbers of applicable GMPEs for the region. The seismic hazard map for Patna district has been developed using PSHA applying probabilistic methods, namely the classical method proposed by Cornell (1968), which was later upgraded by Algermissen et al. (1982), and the smoothed–gridded seismicity models using areal source and four models proposed by Frankel (1995). For the development of a hazard map using the areal approach, delineation of seismic zones has been performed based on the seismicity parameters, i.e. a, b, and M_c. The hazard curves between mean annual rate of exceedance versus PGA and spectral acceleration (S_a) are developed at the rock levels by both models. The final hazard map in terms of the rock level peak ground acceleration values is for 2 % and 10 % probability of exceedance in 50 years, i.e. a return period of 2475 and 475 years based on a logic tree. Additionally, a hazard map for S_a at 0.2 and 1 s for a return period of 2475 and 475 years is also given. Furthermore, uniform hazard spectrum for Patna district at the rock level for a return period of 2475 and 475 years based on a logic tree has been estimated and compared with Indian standard IS 1893.

2 Geology, seismotectonics, and seismicity of study area (SA)

Regional seismicity, geological, seismological, and seismotectonics information of the seismic study area (SSA) have been assembled and evaluated for a desirable radius for seismic hazard analysis. The present study area covers the longitude 84.6–85.65^∘ E and latitude 25.2–25.8^∘ N and is near various rivers such as Gandak in the west, Ganga in the south, and Kosi and Bhagmati in the north (see Fig. 1). Patna lies in seismic zone IV with a zone factor of 0.24 as per IS 1893 (2016). To carry out a seismic hazard analysis, details and documentation about seismic features such as faults, shear zones, and lineaments along with all earthquake events (M_w>4) that have occurred in the SSA are mandatory. Based on a damage distribution map, i.e. isoseismal map (1833 Nepal earthquake and 1934 Bihar–Nepal earthquake), and location of the main boundary thrust, main central thrust, and Himalayan Frontal Thrust (HFT), a radius of 500 km has been selected for the present SSA. A detailed study about selecting the SA of 500 km is given in Anbazhagan et al. (2015a). Geographical information of India demonstrates that approximately 60 % of the land is highly susceptible to earthquakes (NDMA, 2010). The tectonic features of the SA have been compiled from the Seismotectonic Atlas (SEISAT, 2010) published by the Geological Survey of India (GSI, 2000). The seismotectonic map was developed by considering a 500 km radius from the Patna district boundary by considering linear sources (faults and lineaments) from SEISAT and published literature (e.g. NDMA, 2010; Nath and Thingbaijam, 2012; Kumar et al., 2013). Separation of the main boundary thrust (MBT) and the main central thrust (MCT) has been performed and all the faults along with MBT and MCT have also been numbered. The seismotectonic map for Patna district is shown in Fig. 1. A brief description of seismicity and seismotectonics of SSA is given below.

Figure 1 Seismotectonic map of Patna SSA.

Patna district lies near the seismically active Himalayan belt and on the deep deposits of the Indo-Gangetic Basin (IGB). It is also surrounded by various active ridges such as the Monghyr–Saharsa Ridge Fault, many active tectonic features such as the Munger–Saharsa Ridge Fault, and active faults such as the East Patna Fault or West Patna Fault. These faults are acknowledged as transverse faults, and the occurrence of seismic events is due to stimulus of fluvial dynamics in the North Patna plains transverse faults (Valdiya, 1976; Dasgupta et al., 1987). According to Banghar (1991) the East Patna Fault is one of the active faults in the study area and its interaction with the Himalayan Frontal Thrust is characterized by a cluster of earthquakes. Dasgupta et al. (1993) state that all other faults between Motihari and Kishanganj cities have the same possibility of seismic hazard as they form part of a related fault system. Historic earthquakes such as 1833 Bihar, 1934 Bihar–Nepal, and 1988 Bihar–Nepal have affected Patna city as far as economic loss and loss of lives is concerned. Many other earthquakes that have occurred near the Bihar–Nepal border also prove to be devastating for Patna district. In addition, the north side of Patna is near the East and West Patna faults. The frequency of seismic events on these faults is high (Valdiya, 1976; Dasgupta et al., 1987). In addition SSA is also 250 km from the Himalayan plate boundary. These plate boundaries were the source of major historic earthquakes. Considering the above seismic aspects, Patna district can be acknowledged to be under a high seismic risk. Thus, in the present work, PSHA of Patna district has been carried out by considering all seismic sources and earthquake events by reducing epistemic uncertainty using the logic tree approach.

The earthquake data are collected from various agencies such as the National Earthquake Information Centre (NEIC), International Seismological Centre, Indian Meteorological Department (IMD), United States Geological Survey (USGS), Northern California Earthquake Data Centre (NCEDC), and GSI. The events have been selected from all the mentioned agencies. Duplicate events have been deleted, and the magnitude has been homogenized to moment magnitude scale. A total of 2325 events have been compiled which are on different magnitude scales such as local magnitude, surface wave magnitude, and body wave magnitudes. To attain uniformity, all the reported events are converted to moment magnitude (M_w) using relations given by Scordilis (2006) considering worldwide data. Furthermore, the declustering algorithm proposed by Gardner and Knopoff (1974), and modified by Uhrhammer (1986), was used for the separation of the main event from dependent events. Out of 2325 events, 54 % were noticed as dependent events, i.e. 1272 events were documented as main shocks for the Patna region. The complete catalogue contains 454 events having a moment magnitude less than 4 and 1127 events with M_w≥4. To develop the seismotectonic map, the linear source map was superimposed with the declustered earthquake events and given as Fig. 1. Near MBT and MCT, earthquake events are densely located (See Fig. 1) compared to other parts of the seismotectonic map. As per Cornell (1968) and Frankel (1995), the seismic study area needs to be divided based on the seismicity or tectonic provision for calculating the significant hazard value from any potential source. Based on the event distribution SSA is divided into region I (which belongs to MBT and MCT) and region II. These regions were separated using a polygon, as shown in Fig. 1; region I fit into events inside the polygon and region II belongs to events outside the polygon. Region I contained 280 events with M_w 4 to 5, 197 events with M_w 5.1 to 6, 26 events with M_w 6.1 to 7, and 4 events with M_w greater than 7, whereas Region II contained a total of 310 significant events, viz. 168 events with M_w 4 to 5, 121 events with M_w 5.1 to 6, and 21 events with M_w 6.1 to 7. Both the regions were separately analysed for the seismic hazard estimation.

3 Seismicity parameters

3.1 a and b parameters considering period of completeness

The most widely known Gutenberg–Richter (G–R) relationship (Gutenberg and Richter, 1956) is usually used for the determination of a and b parameters for any SSA. The seismic recurrence rate can be precisely calculated only for the complete seismic event data. Stepp (1972) is used for examining the completeness of both the regions. Based on the analysis, it has been observed that for M_w>5, the catalogue is completed for 110 years for both the regions. However, for M_w<5, the catalogue is completed for the last 80 and 70 years respectively for region I and region II. After determining the completeness of catalogue, G–R recurrence law for both the regions has been estimated. The b value for region I and region II respectively was found as 0.87 and 0.97. Whereas the a value for region I and region II respectively for the present study was determined as 5.32 and 4.98. More details about period of completeness and G–R recurrence law are described in Anbazhagan et al. (2015a).

3.2 Magnitude of completeness (M_c)

Magnitude of completeness is defined as the lowest magnitude at which 100 % of the events in a space–time volume are detected (Rydelek and Sacks, 1989; Taylor et al., 1990; Wiemer and Wyss, 2000). M_c is also important for mapping out seismicity parameters such as b value of the Gutenberg–Richter relationship. The magnitude of completeness was calculated using nine different methods defined by Woessner and Wiemer (2005). In addition to magnitude of completeness, these methods also estimate G–R a and b parameters. These methods are maximum curvature method (M1), fixed minimum magnitude observed (M_min) (M2), goodness of fit M_min90 (M3) and M_min95 (M4), best combination of M_min90 and M_min90 and maximum curvature (M5), entire magnitude range (M6), Shi and Bolt (1982) method (M7), bootstrap method (M8), and Cao and Gao (2002) method (M9). Magnitude of completeness for Patna site for region I and region II (shown in Fig. 1) was estimated using software package ZMAP (Wiemer, 2001), a MATLAB-based programme. The a, b, and M_c from each method are represented as Fig. 2 for M1, M2, M3, M4, M5, M6, M7, M8, and M9 for both the regions. It has been observed that M_c varies from 1.7 to 5.0 M_w for region I and 1.9 to 4.9 M_w for region II. It is also observed that an R value of 95 % for the observed magnitude–frequency distribution cannot be modelled by a straight line for region II due to lack of data. The Gutenberg–Richter a and b parameters calculated using these nine methods are different from those calculated using completed data with G–R relationship values for both the regions. Calculated values of G–R a and b parameters for both the regions are given in Table 1. The value of the a parameter calculated from the above methods varies from 3.11 to 6.57 for region I and 3.07 to 6.4 for region II. However, the b parameter calculated from the above methods varies from 0.149 to 0.843 for region I and 0.176 to 0.848 for region II. The lower b value is observed as it is calculated based on the magnitude of completeness, which may be due to the change in the algorithm as it selected the completed magnitude as minimum observed magnitude. This is not used further in the hazard calculation. The difference in a and b parameters is determined using the above methods, as it is calculated based on magnitude of completeness using mixed data (Woessner and Wiemer, 2005) instead of period of completeness for completed data of earthquake events. It can be seen from Table 1 that the average value of the a parameter is 4.95 for region I, which is low compared with the number of earthquakes in the region. Similarly, the average b values of 0.522 and 0.661 for region I and region II are also low when compared to the number of earthquake events having a larger magnitude. According to Boomer et al. (2005) calculation effort increases dramatically with the inclusion of more branches in the logic tree. Therefore, Bommer et al. (2005) suggested avoiding using branches with slight differences between the options, in cases when those options result in very similar nodes. Hence, only M6 has been used because this method synthetically maximizes the available data and stabilizes the M_c value. Woessner and Wiemer (2005) suggested that the M6 method is capable for M_c calculation as it synthetically maximizes the available data and stabilizes the M_c value. Therefore, for further analysis, a and b values of 6.57 and 0.843 and 6.22 and 0.815 respectively had been considered for region I and II. For further study, a weight factor of 0.5 was given to each of the methods (i.e. period of completeness and magnitude of completeness, viz. M6) used to determine the a and b values for both the regions. The final values of 5.0 M_w and 4.8 M_w are adopted as magnitude of completeness for region I and II respectively for further study.

Figure 2 M_c value along with a and b parameters for M1, M2, M3, M4, M5, M6, M7, M8, and M9 for region I (open square) and region II (open triangle).

Table 1 Variation in magnitude of completeness (M_c) and a and b parameters of G–R relationship.

3.3 Maximum magnitude estimation (M_max)

The maximum probable earthquake magnitude has been calculated using both deterministic and probabilistic approaches. Three methods, viz. conventional methods of increment of 0.5 in maximum observed magnitude ($M_{\mathrm{obs}}^{max}$) based on b values, the Kijko method (Kijko and Sellevoll, 1989), and regional rupture characteristics (Anbazhagan et al., 2015b), have been used in M_max calculation. For the estimation of M_max using Kijko and Sellevoll (1989), calculation of M_c is already discussed above. Secondly, M_max magnitude has been calculated by adding a constant value of 0.5 to the $M_{max}^{\mathrm{obs}}$ value at each fault (see Fig. 1) like the NDMA (2010) report. M_max is also estimated using regional rupture characteristics by considering the $M_{max}^{\mathrm{obs}}$ and possible seismic source. The whole procedure to calculate region-specific rupture characteristics was presented in Anbazhagan et al. (2015a). As per Risk Engineering Inc (1988) and others, the increment varies from source zone to source zone, and as per Wheeler (2009) short historical records produce samples of seismicity that are too small to constrain M_max. As per Anbazhagan et al. (2015b), M_max estimated from probabilistic methods, i.e. Kijko and Sellevoll (1989), is sensitive to SSA and seismicity parameters of a region. However, M_max determined using regional rupture characteristics is more reliable as it depends upon the seismic source and rupture length. Taking these points into consideration qualitative weight factors of 0.3, 0.3, and 0.4 has been assigned to the incremental method, Kijko method, and regional rupture method respectively. More weight is given to the regional rupture approach as it accounts for rupture of seismic source, which in turn depends upon the energy released for an event. Maximum magnitude calculated corresponding to each fault is available in Table S1 in the Supplement and Anbazhagan et al. (2015a).

4 Selection of ground motion prediction equation (GMPE)

GMPEs have been selected based on the efficacy test recommended by Scherbaum et al. (2009) and Delavaud et al. (2009). There are various GMPEs available for the active crustal region and basin. Out of various GMPEs, 27 GMPEs are applicable for the present SA. The details of the efficacy test have been given in Anbazhagan et al. (2015c). Details of these GMPEs are given in Anbazhagan et al. (2015a). Similar to Anbazhagan et al. (2015a), the hypocentral distance is divided into three length bins, viz. 0–100, 100–300, and 300–500 km. The applicable GMPEs with the abbreviation used in the present study are given as Table 2. The determined PGA values are used to estimate the log-likelihood (LLH) values. Further, the data support index (DSI) given by Delavuad et al. (2012) is used to rank the best suitable GMPEs. Positive DSI values have been identified for each segment and ranked based on high to low values. Positive DSI values for the Patna earthquake are marked as bold in Table 3. It can be seen from Table 3 that the three GMPEs ANBU-13, NDMA-10, and KANO-06 can be used for up to 100 km of hypocentral distance. For a 100–300 km distance, ANBU-13, NDMA-10, KANO-06, and BOAT-08 are used and for hypocentral distance greater than 300 km NDMA-10 will be used for further hazard analysis. NDMA-10 is used for a distance more than 300 km, as it is the only available equation for the larger distances. Seismic hazard values in terms of PGA and SA can be calculated considering these equations for each seismic source. The variation in PGA with distance for the selected GMPEs is given as Fig. S1. In addition, LLH-based weight as per Delavaud et al. (2012) for selected GMPEs was also calculated. Scherbaum and Kühn (2011) showed the importance of weight treatments through the logic tree approach as probabilities instead of generic quality measures of attenuation equations, which are subsequently normalized. They also indicated the risk of independently assigning grades by different quality criteria, which could result in an apparent insensitivity to the weights. In order to provide the consistency with a probabilistic framework, they proposed assigning the weight factors in a sequential manner, which is used in the present study. The weight factors of 0.72, 0.17, and 0.11 are calculated with ANBU-13, NDMA-10, and KANO-06 at up to 100 km of hypocentral distance according to Delavaud et al. (2012). For a 100–300 km distance, KANO-06, ANBU-13, NDMA-10, and BOAT-08 with weight factors of 0.32, 0.28, 0.26, and 0.14 are calculated, and for a hypocentral distance greater than 300 km a weight factor of 1 has been associated with NDMA-10. It can be noted here that only one GMPE is surfaced with positive DSI for a distance segment of 300 to 500 km and required additional GMPEs in this range, which is important for the far-field damage scenario in the region. These GMPEs with associated weight factors were further used in probabilistic seismic hazard analysis of Patna SSA. These weight factors would be further useful in forming the logic tree to reduce the epistemic uncertainty in final hazard value. Detailed analysis of determination of LLH and weight factor corresponding to each GMPE is given in Anbazhagan et al. (2015a).

Table 2 Available GMPEs with their Abbreviations considered for the seismic study area.

Table 3 Segmented ranking of GMPEs for the Patna region. Positive DSI values for the Patna earthquake are marked as bold.

NA: not available.

5 Delineation and spatial smoothening of seismic source model

Various researchers have delineated the seismic source for various parts of India. Considering tectonic features and past earthquake events, Gupta (2006) delineated the seismic sources for India. Kiran et al. (2008) and NDMA (2010) have done the same on the basis of the seismicity parameters. Furthermore, Nath and Thingbaijam (2011) have delineated based on focal mechanism data from the Global Centroid Moment Tensor database. Kolathayar and Sitharam (2012) identified and delineated India into 104 seismic zones based on similar seismicity characteristics. Vipin and Sitharam (2013) determined the seismic sources on the peninsula considering the seismicity parameters. In the present study, delineation of the seismic sources has been performed based on the seismicity parameters, viz. a, b, and magnitude of completeness (M_c). For delineation of different zones, Patna SSA has been divided into a grid size of $\mathrm{0.02}{}^{\circ }\times \mathrm{0.02}{}^{\circ }$ and from the centre of each grid a radius of 500 km is considered. The number of earthquake events within 500 km of each radius were considered to determine the seismicity parameters. The reason for selection of a 500 km radius was discussed above and given in detail in Anbazhagan et al. (2015a, 2013a). Considering the seismicity parameters (a value, b value, and M_c), the whole study area has been divided into seven areal seismic zones and shown in Fig. 3 (variation in only b value is shown in background). These seven zones are considered areal seismic sources as these are spread over a large area. The seismicity parameter has been calculated for each of these zones considering the frequency magnitude distribution (FMD) at the 90 % confidence level. M_max for each seismic zone has been calculated as per the method discussed earlier. The average values of a, b, M_c, and M_max have been given in Table 4.

Figure 3 New seismic source zones identified for Patna based on seismicity parameters (variation in b value is shown in background).

Table 4 Seismic parameters for adopted source models (uncertainties with bootstrapping).

^∗ Average value of seismicity parameters.

For spatial smoothening of the seismic source model, a grid size of $\mathrm{0.02}{}^{\circ }\times \mathrm{0.02}{}^{\circ }$ along longitude and latitude was selected for representing different kinds of seismic sources and to count the number of earthquakes with magnitude less than or equal to M_c for each grid. To account for the seismicity of the Patna SSA, the maximum likelihood estimates of 10^a for that grid cell have been determined, which correspond to the number of earthquakes per year. Using a maximum likelihood estimate of 10^a, the recurrence rate for different magnitude intervals has been estimated using an algorithm recommended by McGuire and Arabasz (1990). The value 10^a for each grid has been smoothed by applying a Gaussian function, given as Eq. (1), to find the final modified values corresponding to each grid. This smoothing is made to account for the uncertainty related to the location of earthquake events.

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{̃}}{n}}_i={\displaystyle \frac{{\sum }_j{n}_j{e}^{-{\mathrm{\Delta }}_{ij}^{\mathrm{2}}/c^{\mathrm{2}}}}{{\sum }_j{e}^{-{\mathrm{\Delta }}_{ij}^{\mathrm{2}}/c^{\mathrm{2}}}}},\end{array}$$

where n_j is the number of earthquakes in the jth grid, ${\stackrel{\mathrm{̃}}{n}}_i$ is the smoothed number of earthquakes in the ith cell, c is the correlation distance to account for the location uncertainties, and Δij is the distance between the ith and the jth cells. The sum taken over the jth cell should be within the distance of 3c of the ith cell.

6 Computation models for determining hazard value

Probability of exceedance of a ground motion for a spectral period can be determined once the probability of its size, location, and level of ground shaking are identified cumulatively. Seismic hazard map for the Patna district has been developed by applying probabilistic methods, namely the classical method proposed by Cornell (1968), which was later improved by Algermissen et al. (1982), and smoothed–gridded seismicity models (Frankel, 1995).

A tool of 178 seismic sources (shown in Fig. 1 and given as Table ET1) have been used for determining the probability of occurrence of a specific magnitude, probability of hypocentral distance and probability of ground motion exceeding a specific value as per Cornell (1968). Probability of rupture to occur at different hypocentral distances has been determined as per Kiureghian and Ang (1977). The condition probability of exceedance for GMPEs was determined using a lognormal distribution as given by EM-1110 (1999). The ground motion at a site for a known probability of exceedance in a desired period has been calculated by amalgamating all the above probabilities. As a result of PSHA, the hazard curve is determined and shows the variation in PGA or SA with the frequency of exceedance of the different levels of seismic ground motion. A detailed explanation is given in Anbazhagan et al. (2015a). The deaggregation based on the principle of superposition proposed by Iyenger and Ghosh (2004) has been used. The probability of exceedance of ground motion for each seismic source has been computed by merging these uncertainties. Detailed discussion on the methodology of PSHA can be found in Anbazhagan et al. (2009).

It can be noted that in the SSA, the northwest and central parts of Patna are not fully covered by well-identified seismic sources and many sources given in the Fig. 1 are not well studied to prove its seismic activity. Moreover, there are many places where the linear source has not been identified. Thus, to overcome the limitation, a zoneless approach proposed by Frankel (1995) has been used for developing the PSHA map for Patna SSA. This method accounts for the spatial smoothing of historic seismicity to directly calculate the probabilistic hazard. The annual rate of exceedance for a given ground acceleration level is given by Eq. (2):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathit{\lambda }\left(Z>z\right)=\\ & {\sum }_d{\sum }_{\mathrm{i}}{\mathrm{10}}^{\left[{\log }_{\mathrm{10}}\right]\left(N_{\mathrm{d}}/T\right)-b(m_{\mathrm{i}}-m_{\mathrm{cut}})}P(Z>z/D_{\mathrm{d}}{M}_{\mathrm{i}})\end{array},\end{array}$$

where d and i are indices for distance and magnitude bins. N_d is the total of ${\stackrel{\mathrm{̃}}{n}}_{\mathrm{i}}$ values over a given hypocentral distance increment (calculated using Eq. 1), $P(Z>z/D_{\mathrm{d}}{M}_{\mathrm{i}}$) will give the probability that a PGA of Z will exceed z when an earthquake of magnitude M_i occurs at a distance of D_d, and T is the time in years of earthquake catalogue used to determine N_d. The probability that a PGA of Z will exceed z can be determined using EM-1110 (1999). The hazard map has been determined by the four models proposed by Frankel (1995). Model 1, model 2, and model 3 are used for a magnitude less than 7; however model 4 can be used for a magnitude greater than 7. In model 1, the earthquake events having M_w between 3 and 5 are assumed to illuminate areas of faulting which can produce destructive events. Model 2 also ensures that the hazard map reflects the local historic rate of magnitude moment of 5 and larger events. As this model cannot explain the cause of major earthquakes in the active region with certainty, it is prudent to address the possibility of near repeats, i.e. within about 100 km of a historic moderate earthquake. Model 3 is based on a uniform source zone encompassing the active seismicity zone, which is opposite to model 2. Model 4 is associated with hazard from the larger events that is M_w>7. As these events are less in the active seismic region and limited to a few areas only, sources associated with them have been considered for determining hazard. These models are shown in Fig. 4, which is used for the development of the PSHA map using the method proposed by Frankel (1995).

Figure 4 Four models used in the development of the PSHA map of Patna based on the zoneless approach.

7 Modelling of logic tree for hazard analysis

Seismic hazard can be assessed more practically using a logic tree (Kulkarni et al., 1984) as it includes the accounted for epistemic errors, components of seismic models, and ground motion predictions (Fig. 5). For determining the consistent model with different degrees of confidence, each branch of logic tree is to be investigated for implementing the uncertainties in probability models. Important consideration has been given to each branch of logic tree by incorporating the weights for assessing the final hazard of the Patna district. After declustering the catalogue and developing the seismotectonic map, two models have been used with an equal weight of 50 % for both the classical and zoneless approaches. The zoneless approach has been further divided into the areal approach and Frankel approach of equal weight of 50 % each. For the Frankel approach, SSA has been considered for four models (discussed above) with weight factors of 30 %, 30 %, 20 %, and 20 % for model 1, model 2, model 3, and model 4 respectively. These weights have been adopted based on the reliability of the source model. Larger weights are assigned to model 1 and model 2 because they are based on more reliable data and assumedly have better representation of seismicity of SSA. Model 3 deals with the weak assumption that earthquakes with magnitude 3.0–7.0 are equally probable everywhere in Patna SSA whereas there is a great uncertainty in the data used for model 4. In addition, b values were calculated for each of the models using Gutenberg and Richter (1956) and Woessner and Wiemer (2005) (using the entire magnitude range method) by assigning an equal weight factor of 0.5. Furthermore, M_max has been calculated using three methods, namely increment to $M_{max}^{\mathrm{obs}}$, Kijko and Sellevoll (1989), and regional rupture characteristics with weight factors of 30 %, 30 %, and 40 % respectively for each model as shown in Fig. 5. Segmented-based analysis of GMPE was performed and weight was assigned to each GMPE based on the efficacy test. Based on the above discussion, a final hazard map for Patna SSA has been produced for 2 % and 10 % probability of exceedance in 50 years.

Figure 5 Formulated logic tree used in PSHA of Patna SSA.

8 Mapping of probability of exceedance using a different approach considering epistemic uncertainty

8.1 Classical approach (Cornell, 1968)

For determining the hazard value, a different weight has been considered with respect to b value, maximum magnitude and GMPE (see Fig. 5). The seismic hazard using the classical approach (Cornell, 1968) has been estimated using 178 seismic sources. SSA is divided into 1725 grids of size $\mathrm{0.02}{}^{\circ }\times \mathrm{0.02}{}^{\circ }$. The whole procedure can be referred to in Anbazhagan et al. (2015a). Hazard curves from the 10 most venerable sources are given as Fig. 6a and S60 is determined as the most reliable for the Patna district (7.5 M_w and hypocentral distance 55.11 km). Figure 6b shows a cumulative hazard curve obtained at the Patna district centre for 0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 1.6, and 2 s. It can be observed from Fig. 6b that the frequency of exceedance for 0.075 g at 0 s is 0.001, which will give the return period of 834 years. This indicates that PGA of 0.075 g has 5.03 % probability of exceedance in 50 years at Patna. Further explanation can be referred to in Anbazhagan et al. (2015a). The mean deaggregation plot for Patna for return periods of 2475 and 475 years is given as Fig. 7a and b. PGA for 6.0 M_w at 40 km hypocentral distance is notable for 2 % probability of exceedance at 50 years. Likewise, for 10 % probability of exceedance at 50 years the motion for 5.5 M_w at 50 km hypocentral distance contributes most. Hazard curves have been generated at each grid for Patna, and the level of ground motion for ν(z) can be estimated from it. Figure 8a and b show the PSHA maps for Patna district for return periods of 2475 and 475 years respectively. PGA varies from 0.35 g in the northwestern and 0.43 g in the northeastern peripheries to 0.08 g towards the central part (see Fig. 8a). Similarly, the PGA value at the northeastern periphery is 5.3 times more than in the central part of Patna considering 10 % probability in 50 years (see Fig. 8b). These results are similar to those of the previous study carried out by Anbazhagan et al. (2015a).

Figure 6 (a) Hazard curve for 10 most contributing seismic sources at Patna city centre. (b) Hazard curve at Patna district centre for different periods using the classical approach.

Figure 7 (a) Deaggregation of hazard value at Patna at bedrock at PGA for 2 % probability of exceedance in 50 years using the classical approach. (b) Deaggregation of hazard value at Patna at bedrock at PGA for 10 % probability of exceedance in 50 years with the classical approach.

Figure 8 (a) PSHA map for Patna urban centre for 2 % probability of exceedance in 50 years using the classical approach. (b) PSHA map for Patna urban centre for 10 % probability of exceedance in 50 years using the classical approach.

8.2 Zoneless approach

Likewise, in classical approaches, epistemic uncertainty has been considered and weight factors are considered as shown in Fig. 5. The PGA map of Patna has been developed using the zoneless approach by dividing it into seven areal zones based on seismicity parameters (Fig. 3), similar to Kolathayar and Sitharam (2012) and Vipin and Sitharam (2013). For the development of the PSHA map using a simplified areal zonal model, the seven zones along with the seismic parameters (Fig. 3 and Table 4) are used. These seven areal seismic sources are smoothed using the smoothed historic seismicity approach recommended by Frankel (1995). For development of the seismic hazard map, each zone is considered to be of seismic source with a constant seismicity value. The activity rate was calculated for each zone and it was obtained by counting the earthquakes having magnitude greater than or equal to M_c (Table 4) for the whole earthquake catalogue using MATLAB. The calculated activity rate was then spatially smoothed according to Eq. (1), and the chosen correlation distance c=50 km. The annual rate of exceedance at the centre of each grid for the seven zones has been calculated using Eq. (2). The cumulative hazard curves for different periods at the Patna district centre are given as Fig. S2. At the zero period, frequency of exceedance for 0.075 g is 0.012 and estimated return period is 84 years, which means 0.075 g has 44.96 % probability of exceedance in 50 years. Similarly, for 0.5 g, the return period is 24.4 thousand years and probability of exceedance is $\mathrm{2.05}\times {\mathrm{10}}^{-\mathrm{1}}$ % in 50 years at the Patna district centre. As the period of interest rises from 0 s to 0.8 seconds, a huge change in return period has been noticed (see Fig. S2, submitted as Supplement). Primarily the return period decreases from 84 years at zero periods to 13 years at 1.0 s, which has further increased to 28 years at 0.2 s and again to 1.97×10⁵ years for 2 s. The mean deaggregation plot for Patna SSA for return periods of 2745 and 475 years is given as Fig. S3a and b. Figure S3a shows that the motion for 6.0 M_w at 15 km hypocentral distance is dominant for 2 % probability of exceedance at 50 years. It changed to 5.75 M_w at 20 km hypocentral distance for 10 % probability of exceedance at 50 years. Figure 9a and b are the PSHA maps for the Patna urban centre for 2 % and 10 % probabilities of exceedance in 50 years respectively considering the zoneless approach. PGA varies from 0.41 g in the southeastern periphery to 0.34 g towards the central part (See Fig. 9a). However, the southwest part of the district encounters PGA of 1.4 times that of the northwest part of the district. Similar PGA in the southwest part increases 1.57-fold compared to the northwestern part while considering 10 % probability of exceedance in 50 years (Fig. 9b).

Figure 9 (a) PSHA map for Patna urban centre for 2 % probability of exceedance in 50 years using the areal seismic zone. (b) PSHA map for Patna urban centre for 10 % probability of exceedance in 50 years using the areal seismic zone.

8.3 Four models using the zoneless approach (Frankel, 1995)

The hazard value for the Patna district has also been determined with the four models proposed by Frankel (1995). Each of these four models (Fig. 4) has a different spatial distribution of seismic activity. However, the present SSA has five characteristic earthquakes (M_w≥7) so models 1, 2, and 3 have been analysed separately by considering earthquake events, and a PGA map using model 4 has been developed based on seismic sources associated with characteristic earthquakes. The seismic hazard map is generated considering a grid size of $\mathrm{0.02}{}^{\circ }\times \mathrm{0.02}{}^{\circ }$. The activity rate was calculated in every grid cell and it has been obtained by counting the earthquakes having a magnitude greater than or equal to M_c=3.0 and 5.0 for model 1 and models 2 and 3 for different periods of the earthquake catalogue (Fig. 4) using MATLAB. The calculated activity rate was then spatially smoothed according to Eq. (1), and the chosen correlation distance c=50 and 75 km for model 1 and models 2 and 3. The annual rate of exceedance at the centre of each grid for the seven zones has been calculated using Eq. (2). The cumulative hazard curve has been obtained from models 1, 2, 3, and 4 at the Patna district centre for 0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 1.6, and 2 s and shown in Fig. S4. At the zero period, the return period is 85 years and 0.075 g with 43.96 % probability of exceedance in 50 years at the Patna district centre and 0.5 g with a return period increased by 24.4 thousand years in the case of PGA. Primarily the frequency of exceedance declines from 85 years at zero periods to 14 years at 1.0 s, which further increases to 29 years at 0.2 s and again till 2.0×10⁵ years for 2 s. Figure S5a and b show the mean deaggregation plot for Patna for 2 % and 10 % probability of exceedance at 50 years. The deaggregation has been calculated by considering the weighted mean from all four models. PGA for 6.0 M_w at 25.25 km hypocentral distance and 5.75 M_w at 30.3 km hypocentral distance is predominant for 2 and 10 % probability of exceedance at 50 years. With the four models described in Fig. 4, a PGA map has been developed for Patna SSA and given in Fig. S6a–d considering 2 % probability of exceedance in 50 years and Fig. S7a–d considering 10 % probability of exceedance in 50 years. It can be noted from model 1 that the southwestern part of Patna has a high hazard value and a similar trend has been seen from model 2. Model 3 is a map of uniform hazard whereas as far as model 4 is concerned, the northeastern part and central part have high hazard because that portion of SSA is associated with characteristic earthquakes. The weighted-mean PGA map for Patna has been developed by assigning different weight to these four models as 0.3, 0.3, 0.2, and 0.2 for models 1, 2, 3, and 4 respectively. A larger weight is given to models 1 and 2 as they represent real seismic activity because they are based on more reliable data. However, model 3 deals with weak conjecture that earthquake events between 3 to 7 are equally likely everywhere in Patna and model 4 has great uncertainty in occurrence of characteristic earthquakes. Figure 10a and b are the PSHA maps for the Patna district for return periods of 2475 and 475 years respectively. PGA varies from 0.34 g in the eastern periphery to 0.26 g towards the northwestern periphery, while it increases 1.38-fold for the southwest part of the district (see Fig. 10a). Similarly, considering 10 % probability of exceedance in 50 years, the PGA value in the southwestern part of Patna is 1.5 times larger (see Fig. 10b).

Figure 10 (a) Weighted-mean PGA map of Patna SSA for 2 % probability of exceedance in 50 years using the Frankel approach (1995). (b) Weighted-mean PGA map of Patna SSA for 10 % probability of exceedance in 50 years using the Frankel approach (1995).

It can be seen from the mean deaggregation plot that the motion for 6.0 M_w at 40 km hypocentral distance, 6.0 M_w at 15 km hypocentral distance, and 6.0 M_w at 25.25 km hypocentral distance is predominant in the case of the Cornell, areal, and Frankel approaches respectively considering 2 % probability in 50 years. However, the motion for 5.5 M_w at 50 km hypocentral distance, 5.75 M_w at 20 km hypocentral distance, and 5.75 M_w at 30.3 km hypocentral distance respectively predominate in the case of the Cornell, areal, and Frankel approaches. The PGA values vary from 0.08 to 0.43, 0.29 to 0.41, and 0.26 to 0.36 g in the case of the Cornell, areal, and Frankel approaches respectively considering 2 % probability in 50 years. However, it ranges from 0.04 to 0.18, 0.09 to 0.16, and 0.09 to 0.16 g respectively considering 10 % probability of exceedance in 50 years in the case of the Cornell, areal, and Frankel approaches. Comparing hazard maps developed using a classical approach and zoneless approach, it has been seen that the northeastern part of Patna SSA has experienced the maximum PGA value. As per the classical approach (Cornell, 1968), predicted PGA value for the central part of the Patna district is 0.08 g, whereas per Frankel's approach (Frankel, 1995) it is 0.32 g; however as per the areal approach it is 0.31 g. Similarly, PGA values of 0.15, 0.39, and 0.39 g have been observed in the case of the Cornell, Frankel, and areal approaches in the southwestern part of Patna SSA. It is because of absence of well-defined seismic source in that area, whereas earthquake events of moment magnitude of 6 and above have occurred. However, in the northwestern part the PGA value is almost equal when calculated using these approaches. This is the reason why both zoneless and classical approaches have been considered in this study to counter the epistemic uncertainty. Thus, both the seismic sources and earthquake events can be accounted for properly.

Figure 11 (a) Final seismic hazard map of Patna SSA for 2 % probability of exceedance in 50 years using the logic tree approach. (b) Final seismic hazard map of Patna SSA for 10 % probability of exceedance in 50 years using the logic tree approach.

Figure 12 Final seismic hazard map of Patna SSA for (a) 2 % and (b) 10 % probability of exceedance in 50 years at 0.2 s and (c) 2 % and (d) 10 % probability of exceedance in 50 years at 1 s using the logic tree approach.

9 Final hazard map using logic tree approach

The final hazard value has been developed by assigning the weight factor of 0.5 to both PGA values calculated with the classical and zoneless approaches. In the present study the average of the hazard values from both the methods have been considered. Hazard integration has been performed as the SSA seismic sources are not identified fully (e.g. northwest and central parts of Patna); hence, to overcome the limitation, the zoneless approach has also been used for developing the PSHA map. It is necessary here to note that the experimenters performing for the seismic hazard assessment using a weighting factor may lead to complication in the calculations with the inclusion of different branches. To prevent this trouble, Bommer et al. (2005) suggested avoiding using the branches having slight differences between the options that it carries, in cases when those options result in very similar nodes. Therefore, when selecting the weighting factors in the logic tree in this study, the cases contrasting (or different) with each other as much as possible have been taken into consideration. In the zoneless approach, 0.5 weight factors were given to both PGA maps developed using areal and Frankel (1995) approaches as explained earlier. Thus, both the hazard maps were compiled and finally a 0.5 weight factor is given to the zoneless approach. The final PGA variation corresponds to 2 % and 10 % probability of exceedance in 50 years, was shown in Fig. 11a and b. In addition, SA at 0.2 and 1 s considering epistemic uncertainty has been given as Fig. 12a, b, c, and d for 2 % and 10 % probability of exceedance in 50 years. PGA varies from 0.37 g in the southeastern periphery to 0.30 g towards the northwest periphery, whereas the southwest part of the district encounters PGA of 0.31 g (See Fig. 11a). Similarly, PGA corresponding to a 475-year return period is about 0.12 g in the northwestern periphery and 0.15 g in the southeastern periphery (Fig. 11b). The reason for having a high PGA value in the southeastern periphery is due to the location of the East Patna and West Patna faults, and PGA value of 0.35 g in the southwestern part is due to the presence of earthquake events of a magnitude moment more than 6. It can be seen from the mean deaggregation plot that the motion for 6.0 M_w at a 40 km hypocentral distance, 6.0 M_w at a 15 km hypocentral distance, and 6.0 M_w at a 25.25 km hypocentral distance are predominant in the case of the Cornell, areal, and Frankel approaches respectively considering 2 % probability in 50 years. However, the motion for 5.5 M_w at 50 km hypocentral distance, 5.75 M_w at 20 km hypocentral distance, and 5.75 M_w at 30.3 km hypocentral distance respectively predominate in the case of the Cornell, areal, and Frankel approaches. The PGA values vary from 0.08 to 0.43, 0.29 to 0.41, and 0.26 to 0.36 g in the case of the Cornell, areal, and Frankel approaches respectively considering 2 % probability in 50 years. However, it varies from 0.04 to 0.18, 0.09 to 0.16, and 0.09 to 0.16 g considering 10 % probability of exceedance in 50 years in the case of the Cornell, areal, and Frankel approaches. The PGA value varies from 0.12 to 0.15 g for a return period of 2475 years, which is comparable with the PSHA map of India developed by Nath and Thingbaijam (2012). Recently, a major thrust faulting earthquake of magnitude 7.8 on 25 April 2015 occurred in Nepal, which affected various places in India including the Patna district. We had completed our mapping before this earthquake and compared our results with a shake map published by USGS (2019). It is noticed that PGA values for 10 % probability of exceedance in 50 years matches with the USGS (2019) shake map for the recent Nepal earthquake.

In addition, a uniform hazard response spectrum (UHRS) has been developed considering all three approaches and compared with IS 1893 (2002). For developing UHRS, seismic hazard curves of spectral accelerations at a different spectral period for the same probability of exceedance have been developed. The UHRS at 2 % and 10 % probability of exceedance for 50 years at the centre of the district using classical and zoneless approaches, viz. Frankel and areal approaches, has been drawn and given as Fig. 13a (marked as star in Fig. 11a). Similarly, the UHRS has been developed at the northeastern part of Patna considering 2 % and 10 % probability of exceedance, shown as Fig. 13b (marked as plus in Fig. 11a). It can be seen from Fig. 13 that the hazard value at 2 % probability is more for the same return period when compared to 10 % probability of exceedance in 50 years. It has also been observed that spectral acceleration at the zero period, i.e. PGA, is less in the case of the Cornell approach when compared to the Frankel and areal approaches at the centre of the district where as it is more when compared to the northeastern part of SSA. The developed UHRS has been compared with IS 1893 (2002) and it has been observed that the SA predicted is lower at the centre of the district at 2 % and 10 % probability of exceedance in 50 years except for Frankel's approach. However, in the case of the northeastern parts of SSA, the predicted SA values are more compared to IS 1893 (2002) (Fig. 13b). Hence, UHRS should be developed based on the regional characteristics so that it could be effectively used in infrastructural development of a district.

Figure 13 (a) Design spectrum for Patna for 5 % damping from 2 % and 10 % probability of exceedance in 50 years and IS 1893 (2002) at centre of the city (marked in Fig. 11a). (b) Design spectrum for Patna for 5 % damping from 2 % and 10 % probability of exceedance in 50 years and IS 1893 (2002) at the northeastern part of the city (marked in Fig. 11a).

10 Conclusion

A new seismic hazard map for the Patna district was developed considering the earthquake events and seismic sources through the logic tree approach. Based on past earthquake damage distribution, the seismic study area of 500 km was arrived at and the seismotectonic map was generated. The maximum magnitude has been estimated by considering three weighted-mean methods, i.e. incremental method, Kijko method and regional rupture-based characteristics. From 28 applicable GMPEs, GMPEs ANBU-13, NDMA-10, and KANO-06 were selected up to 100 km epicentral distance; however ANBU-13, NDMA-10, BOAT-10, and KANO-06 were selected up to 300 km and NDMA-10 for more than 300 km. These GMPEs were ranked and weights were found based on the log-likelihood method. A new hazard map for the Patna district has been developed using both classical and zoneless approaches considering different weight factor corresponding to b value, maximum magnitude, and GMPE to reduce the uncertainty values. The logic tree has been accounted for to capture this epistemic uncertainty in the seismicity models. The final seismic hazard map corresponding to 2 % and 10 % probability of exceedance in 50 years has been developed by giving a weight factor to the seismicity models, maximum magnitude, and GMPEs. The PGA values vary from 0.08 to 0.43, 0.29 to 0.41, and 0.26 to 0.36 g in the case of Cornell, areal, and Frankel approaches respectively considering 2 % probability in 50 years. However it varies from 0.04 to 0.18, 0.09 to 0.16, and 0.09 to 0.16 g respectively considering 10 % probability of exceedance in 50 years in the case of the Cornell, areal, and Frankel approaches. However, hazard values in terms of PGA at bedrock level after considering logic tree vary from 0.30 to 0.37 and 0.11 to 0.15 g respectively considering 2 % and 10 % probability of exceedance in 50 years. In addition, a spectral acceleration hazard map has been developed at a period of 0.2 and 1 s corresponding to 2 % and 10 % probability of exceedance in 50 years. Hence the logic tree should be used to reduce the epistemic uncertainty in determining the hazard value for any seismic study area. It has been also concluded that uniform hazard response spectra should be developed considering region-specific parameters.