Seismic hazard maps of Peshawar District for various return periods

Probabilistic seismic hazard analysis of Peshawar District has been performed for a grid size of 0.01. The seismic sources for the target location are defined as the area polygon with uniform seismicity. The earthquake catalogue was developed based on the earthquake data obtained from different worldwide seismological networks and historical records. The earthquake events obtained at different magnitude scales were converted into moment magnitude using indigenous catalogue-specific regression relationships. The homogenized catalogue was subdivided into shallow crustal and deep-subduction-zone earthquake events. The seismic source parameters were obtained using the bounded Gutenberg–Richter recurrence law. Seismic hazard maps were prepared for peak horizontal acceleration at bedrock level using different ground motion attenuation relationships. The study revealed the selection of an appropriate ground motion prediction equation is crucial for defining the seismic hazard of Peshawar District. The inclusion of deep subduction earthquakes does not add significantly to the seismic hazard for design base ground motions. The seismic hazard map developed for shallow crustal earthquakes, including also the epistemic uncertainty, was in close agreement with the map given in the Building Code of Pakistan Seismic Provisions (2007) for a return period of 475 years on bedrock. The seismic hazard maps for other return periods i.e., 50, 100, 250, 475 and 2500 years, are also presented.


Introduction
Peshawar is the capital city of the Khyber Pakhtunkhwa province of Pakistan that has an important background in the history of the Indian subcontinent. The city provides key access to the Central Asian states through Afghanistan along the western borders of Pakistan. It is located at 710 • 43.4 N, 330 • 93.7 E, in the western Himalayan region.
Peshawar is characterized by high seismicity rates due to its proximity to the active plate boundary between the Indian and Eurasian plates, which are converging at the rate of 37-42 mm yr −1 (Chen et al., 2000). The Main Boundary Thrust (MBT) system, along which the devastating Kashmir earthquake occurred in 2005, is located in the northern parts of the country together with some other active regional fault systems, including the Main Mantle Thrust (MMT) and Main Karakorum Thrust (MKT). These faults, if reactivated, can act as a potential source of seismic hazard for the region including Peshawar (Waseem et al., 2013). This was confirmed also by the recent 2015 Hindu Kush earthquake that caused widespread damage in the province of Khyber Pakhtunkhwa (Ahmad, 2015), including in Peshawar, damaging a number of important structures in the historic city (Fig. 1).
The Pakistan Building Code in 1986 placed Peshawar in Zone 2 which corresponds to intensity V-VI on the Modified Mercalli Intensity Scale. Lisa et al. (2007), based on the probabilistic seismic hazard analysis for the NW Himalayan thrust, recommended a value of 0.15 g for Peshawar. Hashash et al. (2012), using a discrete faults model for northern Pakistan, suggested a peak ground acceleration (PGA) value in the range of 0.20-0.4 g. Rafi et al. (2012), based on the probabilistic seismic hazard analysis and zonation for Pakistan and Azad Jammu and Kashmir, evaluated a value of 0.175 g for Peshawar. Several researchers either regionally or partially have studied the seismic hazard of Peshawar District (Table 1). The Geological Survey of Pakistan (Khan et al., 2006) seismic zoning map suggests a PGA value in the range of 0.03-0.1 g; Zaman and Warnitchai (2010) suggest a range of 0.33-0.40 g, while Zhang et al. (1999) suggest a range of 0.166-0.244 g. The Building Code of Pakistan Seismic Pro- visions (BCP, 2007), which is legally binding for the seismic design of structures in Pakistan, has placed Peshawar in Zone 2B. This zone has peak ground acceleration in the range of 0.16 to 0.24 g for a return period of 475 years. This has revealed that previous seismic hazard studies of Peshawar and northern Pakistan report widely conflicting results (Ahmad et al., 2019;Ambraseys et al., 2005;Khaliq et al., 2019;Seşetyan et al., 2018;Waseem et al., 2018Waseem et al., , 2020. The present study aims to recalculate the seismic hazard of Peshawar, based on the up-to-date earthquake catalogue and ground motion prediction equations, and compare it with that recommended by BCP (2007). The PGA value at bedrock level was calculated using the classical probabilistic seismic hazard analysis procedure. The area sources as suggested by BCP (2007) were those for which the earthquake catalogue was obtained from worldwide seismogram networks and historical records. The modified Gutenberg-Richter empirical model was used to calculate the seismic zone parameters for both shallow crustal and deep-subduction-zone earthquakes. The seismic hazard in terms of PGA at bedrock was calculated and plotted with a GIS tool. Different ground mo- tion attenuation relationships compatible with the geology and seismicity of the local environment were used to quantify the model in terms of variability in seismic hazard of Peshawar District. Furthermore, the logic tree approach was used to take into consideration the epistemic uncertainty. The GIS-based seismic hazard map developed for a return period of 475 years was compared with that given in BCP (2007). Seismic hazard maps were prepared for various other return periods, i.e., 50, 100, 250, 475 and 2500 years.

Probabilistic seismic hazard analysis
The uncertainties in the location, size and rate of recurrence of earthquakes along with the variation in the ground motion intensity and spatial variability can be well considered in the probabilistic seismic hazard analysis procedures (Ornthammarath et al., 2011;Çagnan and Akkar, 2018;Rowshandel, 2018). The probabilistic seismic hazard analysis (PSHA) provides a framework in which these uncertainties can be identified, quantified, and combined in a rational manner to provide a holistic view of the seismic hazard.
According to the modified Gutenberg-Richter law the earthquake exceedance rate λ(M) for an earthquake magnitude M can be defined using Eq. (1): λ 0 is the exceedance rate in the range of the lower limit M 0 and upper limit M u of magnitude; β is the earthquake source parameter. Considering an earthquake as a Poisson process, the probability density of the earthquake magnitude can be obtained using Eq. (2): The strong ground motion parameters, i.e., acceleration, velocity and displacement, are characterized using attenuation relationships that shows the variation in strong-motion amplitude with source-to-site distance and depend on a number of source, path and site parameters (Douglas, 2019;Kramer, 1996;McGuire, 2004;Rupakhety and Sigbjörnsson, 2009). For example, the attenuation relationship for the peak horizontal acceleration has been developed by Campbell (1981) within 50 km of the fault rupture in magnitude 5.0 to 7.7 earthquakes. Campbell and Bozorgnia (1994) developed attenuation relationships using worldwide moment magnitude in the range of 4.7 to 8.1. This relationship is more specific and provides additional terms for source characterization. Toro et al. (1997) have developed an attenuation relation in terms of peak horizontal acceleration on the rock side for the continental portion of northern America. Among others, Boore and Atkinson (2008) and Akkar and Bommer (2010) have developed a site-specific attenuation relationship that can calculate peak acceleration in terms of earthquake magnitude, source-to-site distance, fault mechanism and site condition. The Boore and Atkinson (2008) model was developed based on the empirical regression of the PEER NGA strongmotion database, while that of Akkar and Bommer (2010) was developed for Europe, the Mediterranean and the Middle East region. The mentioned attenuation relationships can be used for ground motion prediction of shallow crustal earthquakes. However, several researchers including Crouse et al. (1988), Crouse (1991), Molas and Yamazaki (1995), and Youngs et al. (1995) have pointed out different conditions of attenuation relationships for shallow and subduction zones. Lin and Lee (2008) and Kanno et al. (2006) have developed attenuation relationships for earthquake records of Taiwan and Japan, respectively. The study of Lin and Lee (2008) showed lower attenuation for subduction zones than for crustal shallow earthquakes. Therefore, the use of shallowcrustal-earthquake attenuation relationships may lead to underestimation of the seismic hazard for subduction earthquakes in probabilistic analysis.
In probabilistic seismic hazard analysis, the peak acceleration at a location is a function of magnitude and distance that is lognormally distributed with standard deviation. In the hazard analysis, the study area is first divided into seismic sources based on tectonics and geotechnical characteristics. The different seismic sources are assumed to occur independently, and the seismic events are considered to occur uniformly over the source. The acceleration exceedance rate v i (a) for the ith single seismic source is calculated using Eq. (3): where M 0 is the smallest and M u is the largest magnitude of the seismic source, Pr(A > a|M, R ij ) is the probability that acceleration A exceeds the value a at distance R ij for an earthquake of magnitude M. The acceleration exceedance v(a) due to all sources N is calculated through combining all sources, as given in Eq. (4):

Seismicity of Peshawar
The collision of the Eurasian and Indian plate has resulted in the formation of an active Himalayan orogenic system that is further classified into the Tethyan Himalayas, Higher Himalayas, Sub-Himalayas and Lesser Himalayas (Gansser, 1964). The divisions are based on the tectonic blocks formed and separated by major fault boundaries. The Microsoft Encarta reference library (Ali and Khan, 2004b) shows that the valley of Peshawar consists of the southern part of the Eurasian plate and northern part of the Indo-Australian plate. This part of the Himalayas is variably interpreted to be the Lesser Himalayas (Tahirkheli et al., 1982) and Tethyan Himalayas (DiPietro and Pogue, 2004). The seismic hazard study of Waseem et al. (2013) identified about 21 seismogenic faults around Peshawar. Most of these faults have a reverse fault mechanism and have a Joyner-Boore distance R JB in the range of 19-100 km. According to Ali and Khan (2004a), most of the significant earthquakes felt at Peshawar have their origin in the Hindu Kush region of Afghanistan and a few in northern areas of Pakistan.

Case study -PSHA of Peshawar
The seismic hazard software CRISIS (2007) was used to calculate the peak acceleration at bedrock level for Peshawar District. Figure 2 shows the geographical location of Peshawar District within the geopolitical boundaries of the KP province of Pakistan. The hazard analysis requires seismic source geometry, the earthquake reoccurrence relationship and the selected ground motion attenuation relationship. In the present study the ground motion attenuation relationships of Boore and Atkinson (2008) and Akkar and Bommer (2010) were used for shallow crustal seismic earthquakes and those of Lin and Lee (2008) and Kanno et al. (2006) for deep-subduction-zone earthquakes. The earth-quake events within a 50 km depth were considered shallow, while earthquake events occurring at depths larger than 50 km were considered deep earthquakes. The seismic hazard maps were prepared in a GIS environment based on a grid size of 0.01 • for various return periods, i.e., 50, 100, 250, 475 and 2500 years.

Seismic sources identification and characterization
The Building Code of Pakistan (BCP, 2007) has defined the potential shallow seismic sources for Pakistan including northern areas. Those sources within 200 km of Peshawar were considered potential sources for earthquake activity impacting Peshawar (Fig. 3). The potential seismic sources (seven seismic sources in present study) for the Peshawar region in a rectangular shape with lat 31.888-36.006 • , long 69.562-73.620 • , as shown in Fig. 4, were considered for the compilation of the earthquake catalogue. The earthquake catalogue was obtained using worldwide seismogram network sources, i.e., the United States Geological Survey (USGS), National Geophysical Data Center (NGDC), Global Centroid Moment Tensor (GCMT) catalogue and International Seismological Centre (ISC), using the time span of 1500 till 2015 CE with a focal depth of up to 1000 m. The catalogue also included historical data from Ambrasey (2000) and Ambrasey and Douglas (2004). Khan et al. (2018) also reported an updated earthquake catalogue for Pakistan; however, the majority of their events relevant to Peshawar were already included in the catalogue of the present study for seismic sources characterization. These different networks already discussed give earthquake magnitude on different scales, e.g., moment magnitude, surface magnitude and low magnitude. According to Kanamori (1977) and Hanks and Kanamori (1979) the moment magnitude is the most accurate scale that does not saturate in highermagnitude events. Therefore, all the magnitudes were converted into moment magnitude (M w ) using regression analysis. Figure 5 shows the empirical relationships established in the present study based on the catalogue obtained for earthquake magnitude conversion. These were used for the catalogue homogenization.
The homogenized catalogue was further subdivided into shallow (depth less than 50 km) and deep (depth more than 50 km) earthquake events. Figure 6 shows the shallow and deep earthquake records along with seismic zones as defined in BCP (2007). Furthermore, Table 2 reports the number of earthquakes in each seismic source along with the maximum and minimum magnitude of each source. Deep earthquakes were found primarily in seismic source 1 and seismic source 2 which included the Hindu Kush seismic region. The deep sources were selected in consultation with the National Centre of Excellence in Geology, Peshawar since deep sources had not been studied before for Peshawar. In seismic hazard analysis the probability of earthquake occurrence is considered to follow a Poisson process, which considers the independent events to occur randomly in time and space. Only the main shocks are considered for hazard analysis. This is to avoid overestimation of the seismic hazard. The dependent events (foreshocks and aftershocks) are temporally and spatially dependent on the main  shocks. For this purpose declustering was performed to remove the dependent events for the catalogue. The Gardner and Kenopoff (1974) declustering algorithm method was used for removing foreshocks and aftershocks. This performs a windowing procedure in time and space on the event magnitude to identify the dependent events. To perform this analysis ZMAP coding developed by ETH Zurich (freely available) was used. The homogenized catalogue was converted into ZMAP-specified format to perform the routine analysis.
A total of 926 independent events remained after declustering.

Completeness analysis
The catalogue also report events from very far in the past, which cannot be considered complete for all the magnitudes and the whole time span. The time window starts from the year 1500; however, since then the catalogue has not been updated on a regular basis. The instrumental observation of seismic data started after 1960, and it now observes and documents complete details of the earthquake events on a regular basis. Due to these reasons the specified time window (1500-2015) cannot be considered in obtaining the activity rate, as this would result in the underestimation of the activity rate. For this purpose completeness analysis was performed using the visual cumulative method (CUVI) proposed by Mulargia and Tinti (1985). It is a simple procedure based on the observation that earthquakes follow a stationary occurrence process. It is used to find the completion point (CP) after which the catalogue is considered to be complete (Tinti and Mulargia, 1985). The procedure is to divide the magnitudes from 4 to 8 into various bands with a 0.5 step size. The selected bands are 4.00 to 4.50, 4.51 to 5.00, 5.01 to 5.50, 5.51 to 6.00, 6.01 to 6.50, 6.51 to 7.00 and 7.01 to 7.7. In each band the cumulative number of total earthquakes is plotted against the year of earthquakes; the period of completeness (T c ) is considered to begin at the earliest time when the slope of the fitting curve can be well approximated by a straight line (Fig. 7). Table 3 reports the completeness points and time periods for each magnitude band.

Seismic source parameters
The modified Gutenberg-Richter (G-R) reoccurrence law, as mentioned earlier, was used in the present seismic hazard analysis to characterize the G-R parameters. The seismic source parameters (i.e., η 0 , β) were calculated from setting a linear trend line to the graph between log λ m and M w as shown in Fig. 8 for both shallow and deep earthquakes for all seismic zones. Table 4 reports the seismic sources' G-R parameters for all seismic sources for both shallow and deep earthquakes.

Attenuation relationships and peak ground acceleration
The attenuation relationships for a site are developed using substantial dataset information (Cotton et al., 2006); however, this is not available for Pakistan because of the scarcity of available strong-motion data. The alternative to this is to use the already-available attenuation relationships of other regions which have similar tectonic and geological conditions to Pakistan. In the case of shallow earthquakes, the candidate attenuation relationships for north Pakistan should be the ones developed for the active tectonic crustal earthquake region. Thus, the ground motion attenuation relationships of Akkar and Bommer (2010) and Boore and Atkinson (2008) were used to calculate the PGA for shallow seismic sources. However, the ground motion attenuation relationships of Lin and Lee (2008) and Kanno et al. (2006) developed for subduction zones were used for deep seismic sources. The seismic hazard in terms of PGA was then calculated at the bedrock site for different return periods, such as 50, 100, 250, 475 and 2500 years, as the cumulative seismic hazard due to both shallow and deep seismic sources. The various ground motion prediction equations (GMPEs) were combined through a logic tree approach and assigning equal weightings to each GMPE. The ground motions calculated were plotted in a GIS environment to obtain the seismic hazard maps for these different ground motion attenuation relationships.

Seismic hazard maps
The seismic hazard levels (Table 5), based on peak acceleration, defined in BCP (2007) were considered as a basis for the zoning of the seismic hazard at bedrock level. The seismic hazard maps for a return period of 475 years in the case of shallow crustal earthquakes on the one hand and deep earthquakes on the other for Peshawar District are reported in Figs. 9 and 10, respectively. Figure 9 shows that for a return period of 475 years, the predictive relationship of Akkar and Bommer (2010) overestimates the PGA value in comparison to that of Boore and Atkinson (2008), especially  in the northern parts of the district. According to Arango et al. (2012), the distance scaling factor of the latter appears to be more adequate then the previous. Furthermore, Table 6 shows a slight comparison of both ground motion prediction equations that suggests that in terms of N R (number of records), T max (longest response period), M w (moment magnitude) and [R] (distance range), the prediction equation of Boore and Atkinson (2008) is more appropriate and reliable than that of Akkar and Bommer (2010) for hazard assessment. Figure 10 shows the seismic hazard maps for deep subduction earthquakes using the Lin and Lee (2008) and Kanno et al. (2006) attenuation equations for a return period of 475 years. According to Fig. 10 both the attenuation equations resulted in roughly similar seismic hazard for Peshawar District. Furthermore, it is also evidenced from Fig. 10 that  the inclusion of deep subduction zones in the seismic hazard does not contribute significantly; i.e., it remains low (0.08-0.16 g) to very low (< 0.08 g). The cumulative seismic hazard may slightly increase the ground motion level, especially in the northern parts.
In probabilistic seismic hazard analysis (PSHA), one of the major sources of uncertainty is the epistemic uncertainty arising from the selection of predictive relationships. Thus, the different ground motion attenuation relationships already discussed were further used to find out the epistemic uncertainty in the seismic hazard analysis. This was accomplished through the logic tree approach, assigning an equal weighting factor to each GMPE (Fig. 11); the seismic hazard was combined from all the GMPEs. Figure 12 shows the seismic hazard maps for shallow and deep events after incorporating the epistemic uncertainty. As can be seen in Fig. 12a, the seismic hazard of Peshawar District becomes balanced when the average of the seismic Figure 10. Seismic hazard maps for deep subduction earthquake using different attenuation equations and for a return period of 475 years.  hazard calculated using the Akkar and Bommer (2010) and Boore and Atkinson (2008) predictive equations was taken. The reason for this is the provision of equal weighting to both the predictive relationships in hazard analysis. The seismic hazard in the case of deep-subduction-zone earthquakes remains roughly the same after incorporating epistemic uncertainty (Fig. 12b). It can also be further concluded that the earthquakes produced by deep subduction zones are not significant in terms of seismic hazard and may be reasonably ignored. Thus, the shallow seismic sources are sufficient for the seismic hazard assessment of Peshawar. The calculated seismic hazard map after incorporating epistemic uncertainty is compared with the hazard map from BCP (2007). For the return period of 475 years, a close agreement between the two seismic hazard maps can be noticed (Fig. 13). After this check the seismic hazard maps for other return periods, i.e., 50, 100, 250, 475 and 2500 years, were prepared (Fig. 14), which may be used for seismic risk assessment. Hazard maps for various cases are reported in Fig. A1 through Fig. A8.

Conclusions and recommendations
The following was concluded on the basis of a literature review of past seismic hazard studies of Peshawar and classical PSHA conducted for Peshawar in the present study: -The selection of an appropriate ground motion prediction equation is crucial in defining the seismic hazard of Peshawar District. In the case of shallow crustal earthquakes, the predictive relationship of Akkar and Bommer (2010) provides a higher estimate of the PGA value in comparison to that of Boore and Atkinson (2008). The distance-scaling factor of the latter appears to be the reason for this disparity between the two models.
-The inclusion of deep subduction earthquakes does not add significantly to hazard and may be neglected in terms of seismic hazard. Therefore, only the shallow crustal earthquakes contribute to the seismic hazard of Peshawar District. However, recent earthquakes in Pe- Table 6. Comparison of predictive equations used for shallow crustal earthquake (after Arango et al., 2012).

Predictive equation Tectonic regime Region
Boore and Atkinson (2008)   shawar from deep sources have caused widespread destruction in various parts of the district. This raises concern for the existing GMPEs and the classical PSHA procedure to simulate such effects.
-The epistemic uncertainty was used by providing equal weighting to the attenuation equation of Akkar and Bommer (2010) and Boore and Atkinson (2008). The mean seismic hazard map thus produced was balanced and was found to be in close agreement with the design base seismic hazard given in BCP (2007) for bedrock hazard. However, the BCP places Peshawar in Zone 2B, which is reasonable for most of the locations, but it  underestimates ground motions especially in northern parts of the district.
-The mean seismic hazard calculated for Peshawar was also compared with previous studies (Table 7). It can be observed that the seismic hazard obtained by independent researchers suggests an average PGA equal to about 0.24 g, which is in agreement with the PGA = 0.24 g given in BCP (2007) for Seismic Zone 2B (0.16 to 0.24 g) for bedrock. The present PSHA study performed using the most up-to-date earthquake catalogue, recent GMPEs and considering both the shallow and deep seismic sources confirmed the validity of seismic hazard given in BCP (2007). It is worth mentioning that the calculated mean hazard may be approximated as the 50th percentile seismic hazard. Table 7 reports that recent hazard studies considering the fault sources have resulted in larger estimate of seismic hazard that place Peshawar in Seismic Zone 3 and Zone 4; however, the idealization of seismic sources as discrete faults for Peshawar is not reliable due to the scarcity of detailed information regarding the fault sources in northern Pakistan. This higher seismic hazard is not justified by the earthquake history of Peshawar.
It is worth mentioning that the focus of the present study was to provide the base maps for seismic hazard in Peshawar. Site-specific soil properties were not known; therefore, they were not addressed in the present study. Alternatively, the code suggests amplification factors for various soils from Type C to Type E as per NEHRP soil classification. This may be considered to amplify or deamplify the seismic hazard provided in the present study.       Data availability. All data, models and code generated or used during the study are available from the corresponding author by request (Naveed Ahmad, naveed.ahmad@uetpeshawar.edu.pk). Items which may be requested are the earthquake catalogue (raw and processed data), Excel sheets used for G-R parameters derivation, CRI-SIS input files, etc.
Author contributions. NA contributed as the advisor and supervisor of the research and in the analysis of the earthquake catalogue for the derivation of seismic parameters, selection of GMPEs, preparation of input files for CRISIS program and paper drafting. UK contributed in the compilation of the earthquake catalogue and processing of data, analysis of hazard through the CRISIS program and developing hazard maps. KM has contributed to the data organization, literature review, integration of work tasks and preparation of the initial paper draft. QI has contributed in the selecting and designing of seismic sources for both shallow and deep earthquakes, assignment of source parameters in the CRISIS program, output data compilation, and result plotting.
Competing interests. The authors declare that they have no conflict of interest.