Seismic Hazard Assessment of Shigo Kas Hydro-Power Project (Khyber Pakhtunkhwa, Pakistan)

Abstract

In this paper, a seismic hazard assessment (SHA) of the Shigo Kas hydropower project has been performed by deterministic and probabilistic approaches. The previously developed MATLAB-based code has been used for deterministic SHA, incorporating local site effects through deep soil analysis. On the other hand, for probabilistic SHA, CRISIS 2007 has been used through diffuse areal source zones. The latest updated earthquake instrumental and historical catalogs have been developed. Based on the recommendations of the International Commission on Large Dams, peak ground acceleration (PGA) values for the maximum credible earthquake (MCE), safety evaluation earthquake (SEE), design basis earthquake (DBE) and operating basis earthquake (OBE) have been assessed, which are 0.50 g, 0.68 g, 0.35 g and 0.24 g, respectively, at the intake location, and 0.50 g, 0.61 g, 0.30 g and 0.22 g, respectively, at the powerhouse location. Hazard maps have been developed for scenario-based earthquakes (MCE) and for the peak ground acceleration of 145-, 475- and 2500-year return periods. The de-aggregation process has evaluated the combined effects of magnitude and distance. At a distance of 30 to 70 km from the earthquake source, earthquakes of magnitude 5 Mw to 5.6 Mw and 5.9 Mw to 6.5 Mw are more hazardous for the current project.

1. Introduction

Pakistan faces an acute shortage of energy, and with the passage of time, the gap between demand and deliverable capacity at national power stations is widening. The shortage of energy can be minimized by exploring small and medium hydroelectric projects in the country. There is abundant untapped potential for hydroelectric energy in Pakistan, which needs to be harnessed. Pakistan is endowed with a hydroelectric potential of more than 40,000 MW, most of which lies in the province Khyber Pakhtunkhwa [1]. The Shigo Kas hydropower project (HPP) is a feasible and economical source of energy [2], having an installed capacity of 102 MW. This project was identified by the German Technical Cooperation (GTZ) in the 1990s. However, the study stopped at a preliminary level. This study is undertaken to perform the seismic hazard analysis (SHA) of the Shigo Kas HPP area based on diffuse areal zones (the zones have similar seismic potentials within their region). This technique is beneficial when faults are not well-characterized and it is not possible to correlate them with earthquake hypocenters in areas [3,4].

The damages due to earthquakes can be minimized by the characterization of source, geology and geotechnical conditions, and structural resistance. SHA is used for the estimation of strong ground motion parameters, such as peak ground acceleration (PGA), peak ground velocity (PGV), and peak ground displacement (PGD) [5,6,7,8,9]. These parameters are used for the seismic-resistant design of structures [10]. Earthquakes are a natural phenomenon that cannot be controlled or predicted; however, the vulnerability of structures can be significantly minimized by constructing seismic-resistant structures [11,12,13]. Seismic hazard evaluation enables us to predict future ground motion in a potential area, which can be used to design new structures and retrofit existing structures [14,15].

Pakistan is located in a seismically active region of the world [16]. It has been struck by disastrous earthquakes in the past and in recent times. According to Norsar [17], the ancient city Taxila was destroyed by an earthquake of intensity X. In 1555 AD, an earthquake struck Srinagar (Indian-occupied Kashmir) and its surrounding areas, and caused extensive damage that was greater than the 2005 Kashmir earthquake [18]. Other disastrous earthquakes that affected the study region include the Kunar (Afghanistan) earthquake of magnitude 7.5 in 1842, the Abbottabad (Pakistan) earthquake of magnitude 7.7 in 1878, the Srinagar earthquake of magnitude 6.3 in 1885, the Kangra (India) earthquake of magnitude 7.8 in 1905, the Quetta earthquake of magnitude 7.5 in 1935, the Nepal earthquake of magnitude 7.8 in 2015, and the Kohistan earthquake of magnitude 6.2 in 1974 [18,19].

An extremely destructive earthquake of magnitude 7.6 occurred in Pakistan on 8 October 2005 [16]. The earthquake’s epicenter was about 19 km northeast of Muzaffarabad near Garhee Habib Ullah (34°29′35” N and 73°37′44” E). In total, 73,000 people were killed, 80,000 were injured, and 3.5 million people were left homeless [20]. Moderate earthquakes (i.e., magnitude 5.0 to 6.0) frequently occur in Pakistan due to the seismically active Hindu Kush region [21]. The Himalayan region has the ability to produce a magnitude 8.0 earthquake once every hundred years [22]. Past earthquakes, like the 2005 Kashmir earthquake and the 2011 Quetta earthquake, have caused great life and financial losses to the country. The Building Code of Pakistan [23] has placed the study area in seismic zone 3, a high seismicity zone with PGA values ranging from 0.24 g to 0.32 g.

Earthquakes are dangerous for buildings and can also destroy the entire infrastructure of an area, such as dams, bridges, pipeline systems, roads, electricity poles, and hydropower stations [24,25,26]. Historically, a few dams have failed due to earthquakes, causing disasters and loss of life and property due to flooding [27,28]. Some examples of the earthquakes that have greatly damaged dams and other infrastructure include the 1999 Chi-Chi (Taiwan), the 2001 Bhuj (India), the 2008 Wenchuan (China), and the 2011 Tohoku (Japan) earthquakes, among others [29]. On 11 March 2011, an extremely destructive tsunami was generated due to the Tohoku earthquake. A magnitude 9.0–9.1 (Mw) earthquake struck the Tōhoku region on the northeast coast of Japan, resulting in a devastating tsunami that killed over 15,000 people and caused damage to more than 1 million buildings destroyed an 18.5 m high Fujinuma earth dam, constructed in Sukagawa City, Japan. The Fujinuma dam retained a small reservoir of 1.5 million cubic meters, which was used for irrigation. The main embankment of the earth-filled dam breached within 20 min of the earthquake. The entire length of the dam toppled and the resultant flood destroyed a small village downstream of the dam, killing eight people [30]. The failure of the dam was a geotechnical fault due to extreme ground shaking. Therefore, any new hydroelectric project in a potentially seismically active region of the world should undergo seismic hazard assessment studies, such as the northern Pakistani areas. Earthquakes can also cause crest settlement in dams, which may reduce the freeboard kept to avoid the overflow of water, or even cause the dam to overflow and, consequently, overturn. The 113.5 m high Tikves dam in Northern Macedonia, built without the compaction of the rock fill, has settled by 2.5 m without any earthquake, and it may settle more if there is seismic activity nearby [31].

Numerous earthquake engineering assessments seek to verify that a structure can survive a specified level of ground shaking while still performing as intended. However, at what level of ground shaking should this analysis be conducted? There is considerable uncertainty regarding the location, magnitude, and associated shaking intensity of future earthquakes [32]. In general, two fundamental methodologies are used for seismic hazard analysis: one is referred to as “deterministic seismic hazard analysis (DSHA)” and the other as “probabilistic seismic hazard analysis (PSHA)”. The deterministic approach estimates strong-motion parameters for the maximum credible earthquake that occurs closest to the site of interest, without regard for its likelihood of occurring during a specified exposure period. On the other hand, the probabilistic approach considers the effects of all earthquakes that are expected to occur at various locations over a specified time period [33].

The purpose of this study is to recalculate the project area seismic hazard using the most recent earthquake catalog and ground motion prediction equations, and to compare it to the value recommended by earlier studies, such as [3,4,21,23,34,35,36]. At the bedrock level, the PGA value was determined using the DSHA and PSHA procedures. Different ground motion attenuation relationships that were compatible with the geology and seismicity of the area were used to quantify the model’s variability in terms of seismic hazard in the project area. Seismic hazard maps for various other return periods, e.g., 145, 475, and 2500 years, were created. In view of the above historical examples, it is evident that a comprehensive SHA is absolutely vital for any hydroelectric project in this region. This article presents the latest guidelines and procedures for the SHA of hydroelectric projects, based on which the SHA of Shigo Kas HPP has been performed. This study will be very useful for the detailed design of the project and the retrofitting of the existing structures in the surrounding areas.

The project is located in the vicinity of active faults, as shown in Figure 1. The study focuses on the area denoted by the red rectangle in Figure 2, which includes Pamir, Hindukush, Northwest Himalaya, the Axial Belt, and the northeastern Karakorum block. Pakistan is located near the junction of the Indian and the Eurasian plates. The repeated earthquakes in Pakistan and the surrounding area are consequences of the collision of these two plates, resulting in the formation of the Himalaya, Karakoram, and Hindukush Mountain ranges due to compression and uplift. The Indian plate is being subducted under the Eurasian plate. The movement of the Indian plate towards the north is occurring at 4–5 cm/year [37]. The resulting collision has splintered the Indian plate into several segments beneath the Kashmir basin [38]. The folding and thrusting of the upper plate produced some major faults in the surrounding areas, and include the main Karakoram thrust (MKT), the main boundary thrust (MBT), and the main mantle thrust (MMT), which are shown in Figure 1 and Figure 2 along with major belts. The western branch of the Hazara Kashmir Syntax in the MBT zone was the source of the October 2005 Kashmir earthquake [39].

2. Earthquake Catalogues

For the probabilistic SHA, a catalog has been compiled from both instrumental and non-instrumental (historical) data. The historical data mainly depend on human observations; for the quantitative analysis of the raw data, intensity scales have been expressed. The historical earthquake data were collected from the [23,40,41] catalogs, and summarized in

Table A1. Major historical earthquakes in the study region.

Earthquake	Year	Estimated Mw	Estimated Max. MMI	Details	Source
Paghman, Afghanistan (34.60° N, 68.93° E)	1505	7.0–7.8	IX-X	Region north of Kabul totally destroyed. In Kabul, many houses collapsed and fortress damaged.	[41]
Kashmir, Pakistan (33.50° N, 75.50° E)	1555	7.6	n/a	The intensity report is limited	[40]
Attock, Pakistan	1669	6.3	VIII-IX	The earthquake felt at Attock caused a fissure of 50 m long in the ground.	[41]
Alingar Valley, Afghanistan (34.83° N, 70.37° E)	1842	n/a	VIII-IX	Felt at Peshawar and Kabul, lasted for 3 min, several shocks, rocked the fort in frightful manner.	[41]
Kashmir, Pakistan (34.60° N, 74.38° E)	1885	6.3	VIII-IX	Destructive shocks in Kashmir (Sopor, Gulmarg, Srinagar almost totally ruined and 3000 people killed), Muzaffarabad heavily damaged. Felt at Peshawar, Lahore, Simla, Leh, Kapalu and Gilgit in a radius of 650 km.	[41]
Chitral, Pakistan	1927	6.3	VIII-IX	Destructive earthquake, ruined Skorzor and Drosh. Damage was equally heavy in USSR at Kulyab.	[23]

Note: The modified Mercalli intensity scale (MMI), developed from Giuseppe Mercalli’s Mercalli intensity scale of 1902, is a seismic intensity scale used for measuring the intensity of shaking produced by an earthquake.

in the Appendix A.

In the instrumental data, the time span from 1927 to May 2019 is covered reliably. Before this period, no instrumental earthquake records were available for this project. To assess the seismicity of the area, various catalogs have been evaluated, including the Pakistan Meteorological Department (PMD) historical database, the United States Geological Survey (USGS) database, the Water and Power Development Authority (WAPDA) Tarbela-Mangla seismic network, and the International Seismological Centre (ISC) of England database. The instrumental catalog contains 618 events within a region bounded between 33.50–36.00° N and 70.50–73.00° E. The major instrumented earthquakes records from the composite earthquake catalog used for this study is presented in

Table A2. Major instrumented earthquakes from the composite earthquake catalog used for this study (Mw ≥ 5.5).

Date	Time	Latitude	Longitude	Depth	Magnitude
5 September 1927	20:10:00	34.2	72	NA	5.6
30 September 1927	18:40:00	34.2	72	NA	5.9
14 November 1928	4:33:09	35	72.5	110	6.1
15 May 1969	20:39:49	34.62	70.82	49	5.6
27 August 1969	22:35:55	35.36	71.13	76	5.5
9 May 1971	26:34.20	35.54	71.06	81.7	5.7
27 August 1969	n/a	35.33	71.59	33	5.5
26 February 2005	n/a	35.5	71.5	NA	6.2
6 May 2008	23:38:02	35.59	72.23	33	5.5
11 February 2008	20:56:59	35.57	71.96	217.2	5.8
24 April 2018	14:25:30	34.50	70.18	71	6.2
27 February 2015	02:59:07	34.76	73.36	32	5.8

in Appendix A.

The data collected from international and national institutes were sorted to eliminate multiple entries of events. The catalog has been evaluated with respect to reliability and completeness. After 2004, only USGS and PMD data were available, so these databases were selected as representative for the most recent years. The conversion of magnitude was performed, and all the body wave magnitudes (m_b) and surface wave magnitudes (M_S) were converted to moment magnitude (M_W) using the relationships given by [42] (Equations (1)–(3)) and [43] (Equation (4)), as given below:

For M_S,

$$M_W=0.67M_S+2.07 \mathrm{for} 3.0 \le M_S\le 6.1$$

(1)

$$M_W=0.99M_S+0.08 \mathrm{for} 6.2 \le M_S \le 8.2$$

(2)

For m_b,

$$M_W=0.85m_b+1.03 \mathrm{for} 3.5 \le m_b \le 5.5$$

(3)

$$M_W=1.46m_b-2.42 \mathrm{for} 5.6 \le m_b \le 7.3$$

(4)

Whereas local magnitude (M_L) was converted to M_W using relationships given by [44], as given below:

For M_L

$$M_W=M_L \mathrm{for} 3.5 \le M_L \le 5.5$$

(5)

Figure 3 shows the spatial distribution of the epicenters of all earthquakes occurring within 150 km of the area, and the Shigo Kas HPP site is at the center. Seismic hazard analysis assumes that the probability of an earthquake occurring follows a Poisson process, in which independent events occur randomly in time and space. It has been observed that most of the catalogs have fore and aftershocks, thus resulting in various discrepancies [34]. The fore and aftershocks were removed, and the main events were extracted via a process known as declustering [45]. The hazard analysis considers only the major shocks. This is to avoid an exaggerated assessment of the seismic hazard. Dependent events (foreshocks and aftershocks) are temporally and spatially dependent on the primary shocks. To accomplish this, declustering was used to eliminate the catalog’s dependent events. Foreshocks and aftershocks were removed using the Gardner and Kenopoff [45] declustering algorithm. This applies a time and space windowing procedure to the event magnitude to identify dependent events.

The catalog prepared for this study was studied for completeness in various magnitude ranges. The cumulative visual method (CUVI) of [46] has been used to assess earthquakes of various magnitudes. It is a straightforward approach based on the observation that earthquakes occur in a stationary manner. It is used to determine the catalog’s completion point (CP)—the point at which the catalog is considered complete [46]. The procedure is to divide magnitudes between 4 and 8 into various bands with a step size of 0.5. The bands are as follows: 4.00–4.50, 4.51–5.00, 5.01–5.50, 5.51–6.00, 6.01–6.50, 6.51–7.00, and 7.01–7.7. The cumulative number of total earthquakes is plotted against the year of the earthquakes in each band; the period of completeness (Tc) is considered to begin at the earliest point at which a straight line can effectively approximate the fitting curve’s slope. The completeness points and time periods for each magnitude band are listed in

Table 1. Completeness of various magnitude earthquakes.

S.No	EQ Magnitude Range	Completeness Year
1	Mw ≥ 4.0	1997–2019
3	Mw ≥ 4.5	1991–2019
4	Mw ≥ 5.0	1968–2019
5	Mw ≥ 5.5	1959–2019
6	Mw ≥ 6.0	1928–2019
7	Mw ≥ 6.5	1927–2019

. The completeness period for various magnitude ranges is then used to establish the Gunterberg–Richter recurrence relationship, discussed in Section 4. The resulting completeness periods for various earthquakes are given in Table 1.

3. Modeling of Earthquake Sources

3.1. Crustal Faults Source Model

The project area is situated in a seismically active region. All the faults within a radius of 150 km of the project area have been taken into consideration. Fault lines in the surrounding area identified by different researchers were taken into consideration; for example, the fault map by [23,34] for the Kashmir 2005 earthquake. Wells et al. [47] developed a relationship between 421 historical earthquakes, which has been used for the maximum magnitude potential evaluation of the faults. Due to the lack of sufficient data, 50% of the fault length is considered the rupture length [48]. The study indicates that an MMT located 1.5 km away from the intake site, with a potential magnitude of 7.8, is the most seismogenic source. The names of potential seismogenic faults in the surrounding area of the project, as well as their lengths and types, are presented in

Table 2. Faults within a 150 km radius of Shigo Kas HPP.

S.No	Fault Name	Type	Joyner and Boore Distance from Shigo Kas HPP (RJB)	Mw	Length (km)	Dip Value (Degree)
1	Darband Fault	Normal	121.00	7.1	47	60
2	Himalayan Frontal Thrust	Reverse	149.43	7.8	225	30
3	Hissaratang Thrust	Reverse	95.44	7.7	160	60
4	Khairabad Fault	Reverse	86.17	7.8	205	60
5	Manshera Thrust	Reverse	111.40	7.0	48	30
6	MBT	Reverse	149.70	8.1	353	30
7	MMT	Reverse	1.50	7.8	340	30
8	MKT	Reverse	95.41	8.1	450	35
9	Nowshera Fault	Strike Slip	85.64	7.2	80	60
10	Punjal Thrust	Reverse	133.23	7.4	98	30
11	Puran Fault	Strike Slip	97.42	7.3	99	60
12	Raikot Fault	Reverse	104.34	7.1	63	60

Note: The Joyner Boore distance (RJB) is the shortest distance to an extended fault’s surface projection.

.

3.2. Area Source Models

An area source model is used when it is difficult to elucidate active fault structures from the deforming background. Area source zonation is based on the seismicity of a region. The basic principle is that the seismicity of the area source must be uniform and homogeneous. Different organizations such as NESPAK and PMD have developed source zones for Pakistan [23,35]. The area zonation of PMD has been used for this study. PMD has divided the whole country into 19 zones, as shown in Figure 4. Zones 1, 8, 14, 15 and 16 are in our area of interest, and their seismicity has been discussed in

Table 3. Seismicity of source zones in the surrounding area of the Shigo Kas HPP.

Zone Number	Regions
1	Afghanistan border; having low seismicity.
8	Pattan Kohistan area; the seismicity of this region is high.
14	Malakand and Islamabad Lahore, etc., a very active region.
15	Chitral and Drosh active region.
16	Koh-e-Sulaiman low seismicity.

. The magnitude of the threshold (m_o) indicates the point at which the distribution is truncated. Due to the fact that smaller occurrences are insignificant and poorly limited by the ground motion prediction model, m_o is set to magnitude 4. The upper bound, mmax, is the maximum observed magnitude in the region (within a 150 km radius of Shigo Kas HPP), and is increased by adding a factor of 0.5 for uncertainties and discrepancies in measuring the magnitude of the earthquake.

Table 4. Seismicity model of the study area.

S.No	Zone No.	Total Events (Mw ≥ 4.0)	Zone Name	G–R Relations	a	b	mo	mmmax
1	zone 1	49	Kohistan–Kashmir	log λm = 3.92 − 0.98 m	3.92	0.98	4	6.6
2	zone 8	80	Eastern Afghanistan	log λm = 4.37 − 1.02 m	4.37	1.02	4	6.2
3	zone 14	267	Upper Punjab–KPK	log λm = 5.66 − 1.23 m	5.66	1.23	4	6.4
4	zone 15	192	Chitral	log λm = 5.17 − 1.13 m	5.17	1.13	4	6.7
5	zone 16	25	Koh-e-Sulaiman	log λm = 5.8 − 1.48 m	5.80	1.48	4	5.5

shows the mmax and m_o values for the source regions. The hypo-central distribution in a 150 km quadrangle (Figure 5) shows that the area has a cluster of shallow earthquakes (focal depth 33 km). Still, some medium-depth earthquakes (70–210 km depth) are also visible in the area. In the database, the depth of focus of a few earthquakes is not given, and the plot has taken them as an event with a zero km depth.

4. Development of Earthquake Recurrence Relationship

To compute the frequency–magnitude relationship for each seismic source zone, the Gunterberg–Richter law has been used. The Gunterberg–Richter recurrence (G–R) relationship between magnitude and frequency is given in Equation (6).

$$Lo{g}_{10}\left(N(M_W)\right)=a-b{M}_W$$

(6)

where $N(M_W)$ is the number of events with a magnitude equal to or greater than $M_W$, and a and b are seismic coefficients. The seismicity of a zone is measured by parameter a, which varies from zone to zone, whereas the b value shows the relationship between the number of small and large earthquakes in a zone. A high b value is related to a high number of smaller earthquakes in relation to larger ones.

Once the area sources are identified and the catalog has been developed, the next step is to prepare it for the development of the G–R relationship. Therefore, the catalog has been de-clustered and further subjected to completeness analysis, as explained in Section 2 above. Their magnitude arranges all the events in ascending order. For each event with a magnitude greater than 4, the number of earthquakes is counted and then divided by the time span of the catalog to obtain the number of earthquakes accumulated per year. For the identified regions, the G-R relationship is plotted between the magnitude and number of earthquakes per year, and this is given in Figure 6, Figure 7, Figure 8, Figure 9, Figure 10. The seismicity of all seismic zones is tabulated in Table 4.

5. Seismic Design Criteria for Hydropower Structures

The seismic hazard analysis and design of dams were performed via the guidelines of the International Commission on Large Dams (ICOLD) [48]. The following design parameters are needed for the earthquake-resistant design of large dam structures.

5.1. The Maximum Credible Earthquake (MCE)

The MCE of an earthquake is based on evaluating the maximum earthquake potential, considering the regional and local geology and the seismicity of the area. It is the earthquake that produces the maximum vibratory ground motion against which the water retaining structures and certain other vital components have to be designed, such that they may undergo some displacements, but should continue to retain water, while several structural components may suffer damages and be rendered non-operational. The MCE for each causative fault has to be determined by applying the deterministic seismic hazard analysis (DSHA) approach [48].

5.2. Safety Evaluation Earthquake (SEE)

To prevent the uncontrolled release of water from the dam reservoir, the dam’s safety-critical structures, such as the spillway and bottom outlets, must be designed on SEE. The recommended return period is 10,000 years; however, it may be less than 10,000 years, depending upon the importance of the project [48]. The return period regarding the SEE for the Shigo Kas HPP has been taken as 2500 years, taking into account its smaller reservoir volume.

5.3. The Design Basis Earthquake (DBE)

The DBE for the appurtenant structures has been developed for a 475-year return period [48]. DBE is based on the reference design earthquake for appurtenant structures, such as intake, powerhouse, penstock, etc.

5.4. The Operating Basis Earthquake (OBE)

The OBE is an earthquake event that could reasonably be expected to affect the dam site during the life of the structure. These are the earthquake events that produce vibratory ground motions, in the context of which the project’s safety and continued operational features must remain functional with some repairable damage. This establishes a benchmark against damage and loss of service. The PGA value for a 145-year return period in PSHA is taken as OBE [48].

6. Seismic Hazard Analysis

The SHA of the Shigo Kas HPP has been carried out using the following procedure:

6.1. Ground Motion Prediction Equations

The basic purpose of SHA is to predict ground motion parameters such as PGA and spectral acceleration (SA). Ground motion attenuation equations have been developed for various world regions based on a mathematical process known as the data fit technique from ground motion data. These attenuation models basically depend on earthquake size, the location from the source, the depth of the earthquake, soil type, focal mechanism type, and travel path. Since no attenuation model is available yet for Pakistan owing to the lack of ground motion data, the next generation of ground motion attenuation models presented by Boore and Atkinson (2008) (BA08) [49] and Campbell and Bozorgnia (2008) (CB08) [50] have been selected based on the guidelines of [51]. These attenuation models have been utilized with equal weightage.

6.2. Geophysical and Geotechnical Investigations

Geophysical and geotechnical investigations are very important for the evaluation of local site effects for any hydroelectric structure. A seismic refraction survey (SRS) was used for geophysical investigations. The seismic waves’ velocity, layer thickness, and layer depth were evaluated through the seismic refraction survey. The seismic refraction survey was carried out at important hydraulic structure locations. A total of 2500 m of surveys was conducted, consisting of 10 lines, each of 250 m length.

The geophysical data were processed with a commercially available software, IXRefrax (from Interpex Limited of the USA) to estimate the layer thickness, the depth from the surface, and the primary wave velocities through each layer. The analysis shows that two lithology layers were found at the intake and three layers at the powerhouse location. For the site response analysis performed in the subsequent section, the shear wave velocity is required instead of the primary wave velocity. The primary wave velocities have been converted to shear wave velocities by the formula in [52]. The relationship is valid for rock sites having compressional wave velocities up to 6000 m/s. The Carroll (1969) formula is:

$$V_S=0.756090V_p^{0.81846}$$

(7)

The geophysical data were verified at the intake, the powerhouse, and other structure locations through rotary diamond drilling. A total of 400 m of drilling was done at 10 locations. In-situ testing in the boreholes was carried out at regular intervals. The waxed core samples from the rock strata were preserved as per the requirement of the study, and have been evaluated for percentage extraction, coefficient of rock quality (RQD), and natural density of rock and soil samples.

Based on the above-mentioned experimental tests, the overall project area has been divided into two lithology layers. The uppermost has an average depth of 2 to 3 m and is composed of medium to coarse gravels of igneous and metamorphic origin, which are milky white, greenish, and grayish in color. The second lithology layer was sound bedrock, found beneath the top overburden layer, and quartz with mica and schistose gneiss were the major components found.

6.3. Site Response Analysis

Site response analysis evaluates the effect of ground motion in the overlying soil layers above bedrock. The DEEPSOIL Version 5.1 software [53] has been used for site-specific response analysis. A 1-D equivalent linear approach has been utilized. A seismic refraction survey was performed at different locations of the project site, and shear wave velocity and other soil profile characteristics were obtained. The results were then confirmed by drilling at those sites. The modulus and damping ratio curves from [54] for gravels and rocks were used.

There are limited data available, and all other records consist of low-intensity data, except for the Kashmir earthquake. The PGA value in all of the records was less than 0.02 g; the results of these analyses were discarded. Records of the Kashmir earthquake are available at the Abbottabad, Murree, Nilore, Tarbela, and Mangla recording stations. These recordings are not taken into account because they predict very low-intensity values, most likely due to the effects of the foundations of the building wherein they were taken, or the greater distance from the source. Similarly, recording at the Mangla Dam is not considered due to a lack of information on the recording station. For the devastating 2005 earthquake, there are only a few recordings available. The records for the Taiwan Chi-Chi earthquake in 1999 have been used because they are in close agreement with the 2005 Kashmir earthquake.

Furthermore, the records from the TCU078 (N-S and E-W) station are best suited for our case, as this station is just 5 Km away from the source, and our project is also a few Km away from the potential fault line of MMT. Therefore, the Chi-Chi 1999 (TCU 078 NS and EW) earthquake record form the PEER database [55] has been selected based on the guidelines of [56]. The acceleration time history of Chi-Chi (TCU 078) is shown in Figure 11. Due to the effects of the rocky strata, the overall amplification was not so considerable, and therefore was ignored.

6.4. Deterministic Seismic Hazard Analysis

The NEES Integrated Seismic Risk Assessment Framework (NISRAF) has been used to develop Hazard Maps. The software was developed at the University of Illinois at Urbana-Champaign, USA [57]. NISRAF is a MATLAB-based software with multiple applications, such as SHA, synthetic time history, hazard map development, risk assessment of buildings, etc. In this study, the hazard map development tool of NISRAF has been used for the Shigo Kas HPP area. The shallow crustal earthquake models of BA08 and CB08 have been used with equal weights. For all relations, the site class B boundary condition set by the National Earthquake Hazards Reduction Program (NEHRP, Building Seismic Safety Council, 2001) was assumed as $V_S^{30}=1200 m/s$ [58]. The range of $V_S^{30}$ for the site class B boundary condition is 760–1500 m/s [59], where $V_S^{30 }$ is the averaged shear wave velocity in the upper 30 m section. However, if rock is found under the 30 m depth, then the effective shear wave velocity of the overburden thickness is considered. The $Z_{2.5}$ represents the depth at which the shear wave velocity reaches up to 2.5 km/s. In [51], $Z_{2.5}=0.63$ km, which is calculated using their recommended correlation. $Z_{tor }$ (depth to the top of the rupture) is considered as zero.

The main main mantle thrust (MMT) in our case, with a potential of 7.8 Mw, has been determined for the Shigo Kas HPP site. The coordinates of the nearest point on the projection of the fault were used as the epicenter, which are (34.77° N, 71.96° E). The shear wave velocity at the bedrock was taken as 1200 m/s. The PGA value obtained at the intake and powerhouse locations is 0.50 g. The scenario hazard map developed for the study area is shown in Figure 12.

6.5. Probabilistic Seismic Hazard Analysis

The uncertainties associated with the location, size, and frequency of occurrence of earthquakes, as well as the variation in the intensity of ground motion and spatial variability, can be well accounted for in probabilistic seismic hazard analysis procedures [60]. Probabilistic seismic hazard analysis (PSHA) provides a framework for identifying, quantifying, and combining these uncertainties rationally in order to provide a holistic view of the seismic hazard. In probabilistic seismic hazard analysis, the peak acceleration at a location is a lognormally distributed function of magnitude and distance with a standard deviation. The study area is first divided into seismic sources based on its tectonic and geotechnical characteristics for hazard analysis. The various seismic sources are assumed to occur independently of one another, and the seismic events are assumed to occur uniformly across the source.

In the current study, PSHA has been performed using the CRISIS 2007 software developed by [61]. The software requires the study area maps in a shapefile format, the seismic sources and their seismicity, the attenuation model used, and the return periods for which the analysis is required. The BA08 and the CB08 attenuation models were prepared in a tabular format in CRISIS 2007. The analysis was performed, and PGA values for return periods of 145, 475, 1000, 2500 and 10,000 years were determined. The PGA value for the 475-year return period (10% probability of exceedance in 50 years) is 0.35 g and 0.30 g at the intake and powerhouse locations. This is slightly higher than the previous studies of [23,34,35]. The relatively higher value is due to the influence of MMT, which is 1.5 km away from the intake site.

To understand exactly how distance and magnitude range contribute to the PSHA, it is the most important step in PSHA to perform the de-aggregation. PSHA does not consider specified individual earthquakes, but rather takes into account all the possible earthquakes that could occur and determines the rate at which different ground motion parameters can be expected to be generated at the site as a result. De-aggregation is the decomposition of the total hazard of a site into many single scenarios, in order to determine which scenario contributes significantly to the hazard of a particular location. As in PSHA, the total hazard at a site gives the combined effect of all magnitudes and distances on the probability of exceedance of a given ground motion level due to different seismic sources, and it becomes very difficult to know whether source, magnitude or distance contributes the most to the hazard, and hence de-aggregation is conducted. The de-aggregation is performed in CRISIS 2007 software for PGA at the rock sites, and a return period of 475 years for both BA08 and CB08 attenuation relationships. The de-aggregation results are shown in Figure 13 and Figure 14. According to the BA08 model, earthquakes with a magnitude of 5.9 Mw to 6.5 Mw and at a distance of 30 to 70 Km from the source are considered as a governing scenario of maximum hazard. On the other hand, the CB08 attenuation de-aggregation results suggested earthquake ranges from 5 Mw to 5.6 Mw for a distance from the source of 30–50 Km.

The PGA versus return period relationship for the intake and the powerhouse is shown in Figure 15. The uniform spectral hazard (USH) curve was prepared for the intake and powerhouse sites, and is shown in Figure 16. The seismic hazard maps for 145-, 475- and 2500-year return periods are given in Figure 17, Figure 18, Figure 19, respectively.

7. Conclusions and Discussions

The area of the Shigo Kas HPP is characterized by high seismic activity. Deterministic and probabilistic approaches have been used for seismic hazard analysis of the Shigo Kas HPP. The scenario hazard map has been prepared using the NISRAF software considering local site effects. The nearest point on the projection of the main mantle thrust has been used as the epicenter with a maximum magnitude potential of Mw 7.8. The CRISIS 2007 software has been used for probabilistic seismic hazard analysis. PGA values for the operating basis earthquake, design basis earthquake, and maximum credible earthquake have been obtained at the intake and powerhouse locations. Deterministic as well as probabilistic seismic hazard maps have been developed for the project area. DHSA has been performed, and the PGA value for MCE has been established, which is 0.50 g for the intake structure and powerhouse structure. The PGA value for the 2500-years return period is fixed as SEE, which is 0.68 g at the intake and 0.61 g at the powerhouse. The DBE value corresponds to a return period of 475 years, and a PGA value of 0.35 g has been obtained at the intake location and 0.30 g at the powerhouse location by PSHA. The PGA value for OBE has been estimated as 0.24 g at the intake and 0.22 g at the powerhouse location for a return period of 145 years. De-aggregation analysis has been performed, and it was found that earthquakes of magnitude 5 Mw to 5.6 Mw and 5.9 Mw to 6.5 Mw are more destructive at distances of 30 to 80 km away from the earthquake source.

Recently, Waseem et al. [36] performed SHA for northern Pakistan for shallow and deep seismic events. The resulting study indicated that the northwestern part of Pakistan faces a greater risk of earthquakes than the south and western regions. The maximum observed value of PGA is 0.38 g for deep earthquakes and 0.41 g for shallow earthquakes (475-year return period). The ground motion values associated with shallow earthquakes are also more significant than those related to deep earthquakes. Hashah et al. [3,4] predicted that the most hazardous cities in northern Pakistan are Kaghan and Muzaffarabad, with PGA ~0.6 g for a 475-year return period. The majority of the studies for our project region are characterized by ground motion (PGA) values in the 0.32–0.40 g range for a 475-year return period. On the other hand, the new macro-seismic zoning map of the Building Code of Pakistan [23] calculated ground motion values for the region at 0.24 g to 0.32 g PGA for a 475-year return period. The ground acceleration value for our study area, recommended by PMD [35], is 0.2 g to 0.25 g. Similarly, Monalisa et al. [34] computed the PGA value for a 475-year return period to be 0.20 g for the project area, despite the city being surrounded by several active faults with high activity rates. In comparison to the current study’s findings, the values suggested by Monalisa et al. [34] and the PMD [35] are very low. It is concluded that the current study’s findings are consistent with those of recent similar studies and the Building Code of Pakistan [23].