Rehman et al.
Progress in Earth and Planetary Science(2023) 10:27
https://doi.org/10.1186/s40645-023-00559-x

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Ground motion hazard of the China–
Pakistan Economic Corridor (CPEC) routes
in Pakistan
Qasim Ur Rehman1,2*, Muhammad Waseem3, Waqas Ahmed1, Ihtisham Islam1,4,
Hammad Tariq Janjuhah4*   and George Kontakiotis5

Abstract
Pakistan has seen a burst of infrastructure development recently due to the increased connection between Asia andEast Europe. The China–Pakistan Economic Corridor is a project between China and Pakistan aimed to improve theregional infrastructure that would ultimately enhance the connection between Asia and Eastern Europe. However,
the active tectonics of Pakistan could put this infrastructure at risk if it is not built to the highest hazard prevention
standard. This study reports the ground motion hazard by using the probabilistic seismic hazard assessment approachand the areal seismic source model. The seismic hazard maps of the China–Pakistan Economic Corridor in Pakistan are
derived using the Cornell–McGuire (1968–1976) approach, which takes into account all earthquakes (25AD-2020) that
occurred in Pakistan and nearby regions, the newest ground motion prediction equations, and an updated seismotectonic source model of Pakistan. The final ground motion intensities are attained as peak ground acceleration and
5% damped spectral acceleration at T = 0.2 s and 1.0 s for 475- and 2475-year return periods (estimated for bedrocksite conditions). The results are displayed as color-coded maps that represent the amplitude deviation of groundmotion. From the spatial evaluation of the maps, a peak ground acceleration value of 0.40–0.52 g for the 475-year
return period and a spectral acceleration (0.2 s) value of 1.66–2.13 g for 2475-year return period are mostly observedon the northern and western routes. The central and eastern routes are mostly characterized by a peak ground acceleration value of 0.22–0.24 g for the 475-year return period and a spectral acceleration (0.2 s) value of 0.95–1.13 g due
to diffused seismicity and lower number of faults in this region. The ground motion intensity values obtained in thisstudy can be utilized for the seismic design of all kinds of infrastructure and bridges along the CPEC routes in accordance with the Building Code of Pakistan, the International Building codes, and the load and resistance factor designcodes published by American Association of the State Highway and Transportation Officials.
Keywords Probabilistic seismic hazard, Seismic design, Areal source model, PGA and SA, CPEC route*Correspondence:
Qasim Ur Rehman
qrehman@uoh.edu.pk
Hammad Tariq Janjuhahhammadtariq013@gmail.com
Full list of author information is available at the end of the article© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

1 Introduction
Pakistan has a unique geographical location on theworld map. It has transportation links with its neighbors(China, India, Afghanistan, and Iran), and the Arabian
Ocean connects it to the rest of the world in the south.
The China–Pakistan Economic Corridor (CPEC) projectis the infrastructure program that will link Central Asia
and Europe through Pakistan. The project includes suchas upgrading highways, constructing 125 miles of tunnels, establishing a 2000-mile railway track from Kashgar
to Gwadar, upgrading Gwadar airport, oil and gas pipelines, power grids, and related infrastructure that link
China with Asia, Europe, the Middle East, Africa, andbeyond (Irshad 2015). The highways will act as a channel for connecting billions of people across Asia, Africa,
and Europe. Its main goal is to improve transportation
between China and other regions, boosting trade andeconomic growth (Makhdoom et al. 2018; Tong 2014).
However, the global disaster database specifies that
the African, European, and Asian countries have twice
the comparative average disaster damages. Primarily,
the impact is ten times higher in southern and eastern
Asia than the worldwide average (Chang et al. 2021).
The highways in Pakistan are vulnerable owing to theactive tectonics of the region, topographic relief, geologic structure, human engineering activity, and groundmotion hazard (Fig. 1). The most frequent geologic hazards in the region are slope failures, erosion, avalanches,
earthquakes, and landslides (Rehman et al. 2021). The
landslides adversely affect the people and economy byabolishing the transportation system, buildings, andmajor infrastructure projects. On January 4, 2010, northPakistan experienced a devastating Attabad landslide thatblocked the Hunza River and created a dam upstream.
According to the National Disaster Management Authority report of July 2010, about 16 km of the Karakoram
Highway, i.e., part of the existing CPEC route, was buried, over 115 houses were destroyed, and 20 people were
missing (Karim et al. 2015). Consequently, the region
faced a slowdown in trade between Pakistan and China,
which was about US$ 8.7 billion before the disaster(Cook and Butz 2015). The Pakistan has also experienced
a significant number of earthquakes with varying focal
depths and moment magnitudes (Mw) > 7.0. On October
8, 2005, a massive earthquake of Mw 7.6 hit Pakistan anddestroyed the transportation networks, communication
services, infrastructure, farm fields, and electricity systems. Over 73,000 people died and 3.5 million were left
homeless as a result of this earthquake (Fujiwara et al.
2006). The earthquake resulted in the triggering of various landslides in the vicinity of the epicenter locationnear Muzaffarabad and Balakot. The Hattian Bala landslide having volume 68 × ­106 ­m3 was triggered by the

Page 2 of 31

Kashmir 2005 earthquake and killed more than 1000
people (Basharat et al. 2012). Additional major earthquakes in the region include the 1819 Allah band, 1931
Mach, 1935 Quetta, 1945 Makran, 1997 Harnai, 2001
Bhuj, 2005 Kashmir, 2008 Awaran, 2011 Dalbandin, and2013 Saravan (Bilham 1999; Bilham et al. 2007; Waseem
et al. 2019). This series of earthquakes has resulted innumerous fatalities and infrastructural damage, revealing poorly constructed structures and inadequate seismic
designs (Naseer et al. 2010). These hazards slow down
the economies of developing countries such as Pakistan
through loss of human lives and by destruction of infrastructure (Rehman et al. 2014). The earthquakes with
Mw > 6.5 represent the recurrence of future major earthquakes in the region (Monalisa and Jan 2007). Therefore,
it is necessary to assess the ground motion hazard of thetransportation routes connecting China and Pakistan.
An accurate seismic hazard assessment is a step toward
reducing the impact of future earthquakes. Damages canbe mitigated by advising on earthquake-resistant design,
seismic vulnerability, and regional loss evaluation models. A probabilistic seismic hazard analysis (PSHA) is a
useful quantitative tool for estimating the uncertainty
in the size, location, and ground motion intensity of
future earthquakes in a region and identifying communities at risk from a seismic hazard. The goal of PSHA isto quantitatively estimate, on a probabilistic basis, howmuch shaking (seismic hazard) a site is likely to experience from future earthquakes. The method predicts theground motion intensity at a given location. The probability density function (PDF) is used to model the deviation in ground motion. The rate of earthquakes is thencomputed by taking into account the ground motion
measurements in relation to the threshold magnitude.
The PSHA produces a curve of ground motion amplitudes according to the annual rate. Pakistan’s seismic
hazard assessment has been the subject of collaborative
efforts by various organizations and researchers who
have utilized a wide range of approaches. These include
the research conducted by Ministry of housing and works
for building code of Pakistan (BCP 2007) and by otherresearchers (Mahmood et al. 2020; Monalisa and Jan
2007; Rafi et al. 2012; Sesetyan et al. 2018a; Waseem et al.
2020; Waseem et al. 2019; Waseem et al. 2018; Zaman
et al. 2012; Zhang et al. 1999). The factors that contribute to the variety of adopted approaches in these studies
include: (1) the selection of ground motion prediction
equations (GMPEs) for the region due to the lack of specialized GMPEs for this tectonic regime; (2) the consideration of shallow focus (focal depth up to 50 km) or deep
focus (focal depth greater than 50 km) earthquakes orboth; (3) the choice of a seismotectonic model based on
active fault zones and data (Rahman et al. 2019; Sesetyan

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 3 of 31

Fig. 1 The China–Pakistan Economic Corridor (CPEC) routes overlaid with the major tectonic zones of Pakistan, including the Chaman, Indus
Suture, Main Mantle Thrust, and Karakoram Fault systems. The figure highlights the potential seismic hazards associated with the CPEC routes, as
they traverse through seismically active regions with high rates of tectonic activity

et al. 2018b; Waseem et al. 2020); (4) the availability ofa variety of methodologies, such as Cornell (1968),
McGuire (1976), and the kernel estimate method (Rafiet al. 2012; Waseem et al. 2020, 2019, 2018; Zaman et al.
2012) for seismic hazard assessment. The specialized
GMPEs are not derived from the tectonic regime of Pakistan. Therefore, the present study has utilized the four
random GMPEs from the set of established equations
developed for the similar tectonic regions of the world,
considered both the shallow and deep focus earthquakes,

an updated areal seismotectonic model of Pakistan, andemployed the Cornell–McGuire approach. This approachconsiders a constant seismicity rate that does not change
in the identified seismic sources. Other advantages ofchoosing this standard approach for the present study
include: (1) the availability of a limited-time earthquake
catalog; (2) fault characterization in Pakistan is not very
good; and (3) the ability to effectively assess the maximum magnitude potential of seismic sources.

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

The highway network between China and Pakistan
is linked by several bridges over the Indus River and its
major tributaries. The seismic design of all infrastructure, and particularly the bridges, requires spectral
intensities in compliance with 2012 specifications of
the American Association of State Highway and Transportation Officials (AASHTO 2012) and International
Building Code (IBC) 2021 (Ching and Winkel 2021). The
research evaluates the CPEC routes in Pakistan for all
types of infrastructure in terms of peak ground acceleration (PGA) and spectral acceleration (SA) responses sothat adequate structural engineering mitigation measures may be implemented in time to make this economicroute productive. The paper specifically examines theseismic hazard risk qualitatively in terms of PGA and SA
values at T = 0.2 s (short-period) and 1.0 s (long-period)for 475 years (10% probability of exceedance in 50 years)
and 2475 years (2% probability of exceedance in 50 years)

Page 4 of 31

RPs, respectively. This is accomplished by incorporating
all the required input parameters into the OpenQuake
version 3.10 software engine (Pagani et al. 2023) for hazard computation (Fig. 2).

2 Geologic and tectonic settings of Pakistan
Pakistan lies in a complicated tectonic setting with three
active mountain ranges (Karakoram–Hindukush–Himalaya) intersecting in the north and an active tripe plateboundary (Indian–Arabian–Eurasian) in the south andsouthwest. In the middle section, lies the Sulaiman foldbelt, formed by the oblique collision of Indian and Eurasian plates and comprises several transform faults (Farah
et al. 1979; Farah and DeJong 1984; Jadoon et al. 1994,
1992; Jadoon and Khurshid 1996; Lawrence et al. 1981;
Yeats and Lawrence 1982). The southern section comprises the southwest-oriented Bela–Chaman–Kurram
belt, which connects the Indian plate to the EurasianFig. 2 Map of the China–Pakistan Economic Corridor (CPEC) routes project in Pakistan, from its initiation in 2013 to its projected completion in
2030 (source www.​cpec.​gov.​pk). The CPEC routes connect China’s western region of Xinjiang to Pakistan’s southwestern port of Gwadar, passing
through various regions and cities of Pakistan

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

plate. Along the western boundary of the Indian plate,the Arabian Ocean’s Murray ridge and Owen-fracturezone are generated by the interaction of the three plates(Eurasian–Arabian–Indian). On the south coast of Pakistan, the Makran subduction zone is formed as a result of
the interaction of the Arabian and Eurasian plates next tothe deep port of Gwadar (Malkani and Mahmood 2016).
Figure 1 shows CPEC routes with major tectonic elements and the number of well-known recognized faults
in Pakistan. The Himalayas, Pamir-Hindukush, Karakoram, Sargodha high, Suleiman mountains, Run of Kutch,
and Makran area are the most vulnerable seismic zones
in the world (Waseem et al. 2020). The following methods have been chosen for the hazard computation of
CPEC routes.

3 Methods
This study is carried out by taking into account the mainshock events delineated by coordinates 58°–83° N and
22°–40° E and by setting up a 15-km fishnet grid in ArcGIS software (Fan et al. 2022) on either side of the CPEC
routes, which run from Khunjerab to Gwadar in Pakistan. Furthermore, the updated seismic source model ofPakistan developed by Waseem et al. 2020 is used. The
Gardner and Knopoff (1974) algorithm (Gardner andKnopoff 1974) is used to de-cluster the homogenizedearthquake catalogue (25AD-2020) in Zmap (ver. 6.0)
software by using the MATLAB R2013a platform. The
Open Quake engine, version 3.10, is used to compute theseismic hazard. The source model, GMPEs, and groundmotion intensity parameters are used as input parameters. The source model includes the vertices of eachseismic source zone, as well as a description of the tectonic regime, a-value and b-value, maximum and lowest threshold magnitude, and upper and lower seismic
focal depth encountered in each seismic source zone. The
choice of mean (50%) and mean plus one standard deviation (84th percentile) has due significance in the determination of the focal depths of respective seismic sourcezones. However, the 84th percentile is opted due to thereason; ground motion hazard is abnormally higher due
to the outdated choice of using the mean (50%). The
four updated ground motion prediction equations are
assigned to distinct tectonic regimes (i.e., active shallow crust, subduction interface, and intraslab), and equal
weights are assigned to each GMPE on the logic tree.
The traditional PSHA areal source initially computes
hazard curves for various ground motion and spectral
acceleration intensities at designated areas of interest.
The hazard maps can then be interpolated at 10% and 2%
probability of exceedance in 50 years. The Open Quake
engine uses Eq. (1) to determine the rate of earthquakes.

Page 5 of 31

R(b) =

νi

P[B > b/m, r]fm(M) · fR/M (r, m)dr dm(1)whereas R(b)—the annual rate of earthquakes fabricating
earthquake amplitude “B” greater than “b”, vi is the annual
rate of earthquakes in a seismic source greater than theminimum threshold magnitude.
P[B > b/m, r] is the probability that an earthquake
at a particular location, at distance “r” with magnitude“m” produces an amplitude “B” [e.g., peak ground velocity (PGV), peak ground displacement (PGD), PGA, etc.]
greater than “b”; fm(M) · fR/M are the probability density
functions for earthquake distance and magnitude.
The annual rate of earthquakes given by the equation
can be strained over all the seismic sources by summation. The earthquake occurrence process is Poissonian,and based on this assumption, Eq. (2) estimates the probability that ground motion "B" will exceed "b".
P[B > b/m, r] = 1 − eR(b)(2)

The methods include the collection and compilation ofearthquake data catalogs; data homogenization to a single magnitude scale; declustering, completeness analysis;
geological development of seismic source models; source
characterization; determination of Guttenberg–Richter(G–R) parameters of seismic source zones; selection of
GMPEs; and epistemic uncertainty of the computations
by incorporating a logic tree.
3.1 Data processing for the earthquake catalogue

The first step in predicting the probability of seismic hazard in a zone of interest is to compile an accurate, up-todate, and homogenized earthquake catalogue. A reliable
catalog records all earthquakes that occur in the region
and must be homogenized on a single magnitude scale.
It is used to compute seismic zone activity rates and to
identify the seismic source model. The composite andhomogenized strong motion data compiled by Khanet al. (2018) for Pakistan from 25AD-2016 are used inthis study (Khan et al. 2018). The earthquake catalogue
was further updated for the years 2017–2020 by including data from several online sources, such as the International Seismological Centre, the United States Geological
Survey, and local networks demarcated by coordinates
58°–83° N and 22°–40° E. The data from these sources
are recorded on a variety of magnitude scales, including surface–wave magnitude (Ms), local magnitude (ML),
body-wave magnitude (mb), and moment magnitude(Mw). The data are then homogenized to the moment
magnitude scale using Zara et al. (2014) relationships,
which were validated by Khan et al. (2018) for this region.

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

For the conversion, the bilinear Eqs. (3), (4), (5), and (6)are employed. For the conversion of Ms to Mw, bilinear
Eqs. (3) and (4) are used,

Mw = 0.58 × Ms + 2.46
For

(3)

3.5 ≤ Ms ≤ 6.0

Mw = 0.94 × Ms + 0.36

(4)Page 6 of 31

dividing line in the model to separate shallow and deepfocus earthquakes. The catalogue includes 5506 shallow
and 1131 deep-focus earthquake events. Figure 3 shows
the deep-focus earthquakes (focal depth > 50 km) thatare mainly concentrated in the northern and lower western CPEC routes, and shallow-focus earthquakes (focal
depth < 50 km) are mainly concentrated in the northern,
central, and western CPEC routes.
3.3 Completeness analysis

For

6.1 < Ms ≤ 8.2

For the conversion of Mb to Mw, bilinear [Eq. (5)] is used,

Mw = 0.93 × Mb + 0.45

(5)

For 4.0 ≤ Mb ≤ 6.2
For conversion of ML to Mw [Eq. (6)], (Zare et al. 2014) isused,

Mw = 1.01 × ML − 0.05
For

(6)

4.0 ≤ ML ≤ 8.3

The different magnitude scale data are converted to themoment magnitude scale for the reason that the newest
GMPEs utilize the moment magnitude as the analytical variable. The complete earthquake catalog included
33,210 events with moment magnitudes greater than 4.0
from 25 AD to 2020.
3.2 Declustering

To avoid the overestimation of hazard analysis, thedeclustering technique is used, which separates dependent events (foreshocks and aftershocks) from independent events (the main shock). The Poissonian assumption
is used to calculate the probability of a hazard recurrence, which takes into account the separate events that
occur arbitrarily in time and space. The dependent events
are spatially and temporally dependent on the independent events. The Gardner and Knopoff (1974) approachis used in the ZMAP 6.0 package (Wiemer 2001) on the
MATLAB R2013a platform to remove the dependent
events (Wiemer 2001). The approach was developed for
the Southern California region, but it is also applicableto other parts of the world with similar tectonic settings.
This algorithm has been used in several studies in Pakistan, e.g., Building code of Pakistan (2007), Zaman et al.
(2012), Waseem et al. (2018), and Mahmood et al. (2020).
Using Gardner and Knopoff ’s (1974) technique, inthe declustered earthquake catalogue, a total of 6637
independent events were found that corresponded to
the respective defined seismic source zones. To delineate the source model, a focal depth of 50 km is used as a

The recurrence relationships used to calculate the annual
activity rate of earthquakes are based on the completeness periods of earthquake magnitudes. The Tinti and
Mulgaria (1985) procedure is used to perform the completeness analysis of the de-clustered, homogenizedearthquake catalogue. To perform the completeness analysis, all independent events were grouped into 0.5 unitof magnitude bins (e.g., 4.00–4.50, 4.51–5.00, 5.01–5.5,
5.51–6.00, 6.01–6.50, 6.51–7.00, 7.01–7.50, and 7.51–
8.00). Every magnitude range is subjected to a completeness check, and the completeness years are determinedso that the data in the completeness periods can be considered to compute the activity rate in each source zone.
The year 1902 was a completeness period of earthquakes
for the 7.51–8.00 magnitude range, and this period isfrom 1902 to 2020. Because the data prior to 1902 wereinsufficient and would have led to an underestimation ofactivity rates, only the earthquakes that occurred during
the completeness periods were used. The total numberof events was then plotted against the years on a graph,and the best fit line was drawn across the data points todetermine the completeness period (Table 1). Figure 4
shows the completeness periods of earthquake moment
magnitudes.
3.4 Geological development of seismic source model

The classical Cornell (1968) approach requires a seismic
source model as an input parameter for the computation of hazard analysis. As shown in Fig. 5, a modified
areal seismic source model developed by Waseem et al.
(2020) for Pakistan is used in this study. An aerial source
model signifies the spatial distribution of seismicity patterns in areas with unknown fault locations (Sesetyan
et al. 2018b). In Pakistan, faults are not very well constrained in terms of slip rates, dip angles, rupture, geometry, and other parameters. Therefore, in the absence ofwell-marked faults in the region, area source zones havebeen used in the current study. Faults do, in fact, show
the variation of real seismic hazard, but due to the limitations on fault data, they were not used. Area sources
are defined based on the fault trends and by using local
geologic information. The used areal model describes
the probable seismic source zones based on tectonics,

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 7 of 31

Fig. 3 The earthquake catalogue utilized for probabilistic seismic hazard analysis (PSHA) computation of China–Pakistan Economic Corridor (CPEC)
routes in Pakistan, a earthquake events with focal depth up to 50 km b earthquake events with focal depth greater than 50 km

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Table 1 The completeness periods of the seismic catalogue,
which indicate the moment magnitude ranges of 0.5 units
above the minimum threshold moment magnitude of 4.0, andthe duration of time for which the catalogue is considered to becomplete
Moment magnitude (Mw) ranges

Completeness period

4.00–4.50

1999–2020

4.51–5.00

1980–2020

5.01–5.50

1978–2020

5.51–6.00

1977–2020

6.01–6.50

1969–2020

6.51–7.00

1962–2020

7.01–7.50

1906–2020

7.51–8.00

1902–2020

8.01–8.50

1902–2020

active fault zones, and seismicity patterns. The modified
source model includes the 34 active shallow crust seismic
source zones, four deep subduction intra-slab seismic
source zones, and one deep subduction interface seismic source zone. Active shallow seismic source zones areassociated with seismicity up to a focal depth of 50 km,
whereas deep subduction source zones are associatedwith seismicity greater than a focal depth of 50 km. As
shown in Fig. 6a, b, the shallow and deep seismic source
models with CPEC routes are established separately witha borderline focal depth of 50 km. In the modified seismotectonic model, deep seismic source zones overlappedwith shallow active source zones. Each seismic sourcezone (SSZ) boundary and shape must be defined by theknowledge decree. The seismically high and low zones
can be differentiated based on well-defined active faults.
Whereas the region lacks marked active faults, the seismicity pattern of each SSZ and geologic expert opinion is
employed to demarcate each deep SSZ. In this model, thegeological information of the region is prioritized overseismicity information to distinguish adjacent SSZ. The
Karakoram ranges, for example, were separated from the
Hindukush ranges, while the Himalayan region is separated from the Kohistan Island arc sequence. The seismic source model presented in this work also considers
the known active faults in Pakistan and adjacent regions(Waseem et al. 2020). However, the seismic source zones
were identified by taking into account the region’s tectonic features, and a summary of noticeable tectonic features is provided in the following section (Fig. 1).
The geology and tectonics of North Pakistan are characterized by the coexistence of three terranes: (1) the
Karakoram block (SSZ 2, 3); (2) the Kohistan Island Arc(SSZ 5); and (3) the Indian Plate (SSZ 9, 10, 13) (Searle
et al. 1999). The Karakoram block is separated fromPage 8 of 31

Kohistan Island Arc to the south by the Main Karakoram
Thrust. In the north and south, the Kohistan Island Arc is
confined to the Main Karakoram Thrust and Main Mantle Thrust. The Indian Plate lies south of the Main Mantle
Thrust. The KB contains sedimentary rocks (Paleozoic–
Mesozoic) in its northern wedge, Karakoram batholith rocks (Cretaceous-Miocene) in its central wedge,
and a metamorphic complex (late Paleozoic-Cenozoic)
in its southern wedge. In this part, the rocks are highly
deformed due to the pre- and post-Indian-Eurasian collision with granitic intrusions of the Eocene-Miocene age(Zanchi and Gaetani 2011). The Karakoram block has
overall low seismicity, as evidenced by a Mw 6.0 Gilgit
earthquake, but the rocks at high elevations have residual
strength due to active tectonics and several collisional
episodes of the India–Asia collision (Rehman et al. 2021).
Even the low moment magnitude value in the region cancause catastrophic hazards to infrastructure. The KB has
several faults, including but not limited to the Karakoram
fault system, North Pamir fault, South Pamir fault, Tirch
Mir fault, and Reshun fault. The area has high induced
seismicity in the north and northwest regions due to thepresence of the Pamir-Karakoram (SSZ 2, 29) and Hindukush ranges (SSZ 1,32,33) (Waseem et al. 2018). The
KIA has metasedimentary and volcanic rocks (Cretaceous–Paleocene) in the north, along with Main Karakoram Thrust: the Kohistan batholith (late Cretaceous) and
Lahor granite intrusion in the center, and mafic–ultramafic complexes (Early Cretaceous) in the south, along
the Indus suture zone (Jagoutz and Schmidt 2012; Petterson 2010). The Pamir and Hindukush mountains stretch
along this arc and have experienced shallow-intermediate
focus earthquakes (Mukhopadhyay and Dasgupta 2015).
The Nanga Parbat syntaxial (SSZ 6) bend in the Kohistan
Island Arc is rising at a rate of 6 mm/yr., indicating active
tectonics in the region (Treloar et al. 1991). The Ladakharc is the eastern extension of Kohistan Island Arc. Basedon the cross-cut connection of Jutal dykes in the collisional regime, the Kohistan Island Arc is believed to havecollided with Karakoram block around 75 Ma in the Late
Cretaceous (Searle et al. 1999). Around 50–55 Ma ago,
the Indian plate collided with the Kohistan Ladakh arc(Bignold and Treloar 2003; Petterson and Windley 1985;
Rehman et al. 2011).
The Himalayas and Hazara Kashmir syntaxis (HKS)(SSZ 10) were formed due to the collision of the
Indian and the Eurasian plate around 50–55 Ma ago(Ding et al. 2016). In Pakistan, the Indian plate Himalayan sequence south of the Indus suture zone may be
divided into three tectonostratigraphic divisions: (1)the Greater Himalaya; (2) the Lesser Himalaya; and(3) the Sub-Himalaya (Qasim et al. 2018). The Greater

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 9 of 31

Completeness Analysis [Mw 4.00-6.00]
1800
1600

Cumulative no. of events

1400
1200
1000
800

1999

600
400

1980
1978
1977

200
0
1500

1600

1700

1800

1900

2000

2100

Years
4.00<Mw<4.50

4.51<Mw<5.00

5.01<Mw<5.50

5.51<Mw<6.00

Completeness Analysis [Mw 6.01-8.50]
200
180

Cumulative no. of events

160
140
120
100

1969

80
60

1962

40

1906

20
0

25

525

1902
1025

1525

2025

2525

Years
6.01<Mw<6.50

6.51<Mw<7.00

7.01<Mw<7.50

7.51<Mw<8.50

Fig. 4 Completeness periods obtained through regression analysis using the Tinti and Mulargia (1985) method for CPEC routes in Pakistan byconsidering the moment magnitude range of 0.5 units above the minimum threshold moment magnitude of 4.0

Himalayas are comprehended to the north of Main
Central Thrust (MCT). The MCT’s position in Pakistan
is disputed; nevertheless, Chaudhary and Ghazanfar(1990) suggested the Batal fault that passes through the
Kaghan valley is equivalent to the MCT. In the northand south, the Lesser Himalayas are surrounded by the
MCT and main boundary thrust (MBT), respectively(Chaudhry and Ghazanfar 1990). The region has highseismicity due to the presence of anticlinorium HKS
and several active faults, including the sinistral Jhelum fault, Nathiagali thrust (NT), and Panjal thrust(PT). In the north and south, the Sub-Himalayas (SSZ15) are limited by the Himalayan frontal thrust (HFT)
and the MBT. The major active faults include the dextral Kalabagh fault (KB), the Riasi Thrust (RT), and theSalt Range Thrust (SRT). The zone primarily consists

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 10 of 31

Fig. 5 The seismic source model that is used for the hazard analysis of the China–Pakistan Economic Corridor (CPEC) routes in Pakistan (Modified
after Waseem et al. 2020). The model provides a representation of the potential interface subduction zone, intraslab subduction zone and shallow
active tectonic zone along the CPEC routes in Pakistan

of Cenozoic rocks with limited exposure to Cambrianlayers.
The Indian plate platform, which represents the stable
continental shelf, is separated into two tectonic zones:(1) the Indian shield zone (Nagar Parkar ridge) (SSZ 24)
and (2) the Indus platform tectonic zone (SSZ 16). The
Indian–Eurasian collision produced reverse faults asa result of the compression (Chaudhry and Ghazanfar1990). The presence of Jacobabad and Kandkot highs, theSukhar rift, and the Nagar Parkar ridge have resulted inhigh seismicity in the region. Major earthquakes include
893 Shahbandar (Mw 8.0), 980 Mansura (Mw 6.0), 1668
Samawani (Mw 6.5), 1819 Allah bund (Mw 7.7) and 2001
Bhuj (Mw 7.9). The Indus trough and plains, on the other
hand, have diffused seismicity compared to the Indianshield zone, possibly because of less active faults.

The late Cretaceous–Neogene India–Eurasia collision deformed the western boundary of the Indian
plate, which is marked by the Bela-ophiolite mélangezone (SSZ 23) (Zaigham and Mallick 2000). The Chaman Ornach-Nal transform fault zone (SSZ 20), whichis composed of the Bela ophiolite belt, separates the
Indian plate in the east from the Makran Chagai trenchof the Eurasian plate in the west (Lawrence et al. 1981).
The Kirthar range is surrounded to the east by the Jacobabad Khairpur high (SSZ 22) and the Indus platform.
Pakistan’s central and southern regions are divided into
four major tectonic provinces: (1) the Bela Chaman Kurram fault zone and ophiolite belt/Axial belt (SSZ 17, 23);(2) the Indian plate platform (SSZ 16); (3) the Chagai
Makran trench arc system (SSZ 26, 27, 28); and (4) the
Indian Ocean (SSZ 25, 26).

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 11 of 31

Fig. 6 The individual shallow active and deep interface/intraslab seismic source models used for CPEC (China–Pakistan Economic Corridor) routes
in Pakistan. a The shallow active source model takes into account shallow earthquake events encountered within the uppermost 50 km of the
Earth’s crust b the deep interface/intraslab seismic source model considers deep earthquake events encountered at depths greater than 50 km

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

The ophiolite axial belt is divided into the Suleiman range, Sibi Syntaxis, Kirthar range, Bela-ophiolitemélange zone, and Ornach-Nal fault zone. The Suleiman Range (SSZ 19) is a fold and thrust belt formed by
the Indian plate rotating clockwise and colliding with
the Eurasian plate. Because of active tectonics and thepresence of the Kingri, Waziristan-Loralai, and Muslim
Bagh faults, the area has high seismicity. In the region,intermediate and shallow-focus earthquakes are common. Waseem et al. (2019) have found this range to havethe highest seismic hazard potential. Sibi syntaxis (SSZ
21) represents the foreland syntaxial bend with dextral transpression and basement faults. Kachhi (1909)
and Sharig (1931) recorded high induced seismicity inearthquakes of Mw 7.1 and 6.8, respectively. The Kirthar
range (SSZ 23) is a north–south-oriented external foreland thrust fold belt contained by passive roof duplex
transpression. The Kirthar fault separates the Kirthar
and Khude fold belts. The zone was hit by the 1931 Mach
earthquake (Mw 7.3). The Bela ophiolite mélange (SSZ23), a suture zone, represents the internal thrust and foldbelt with sinistral transpression. The zone is bounded
to the east by faults including the Diwani, Mor, and Pab
faults, and to the west by the sinistral strike-slip OrnachNal fault zone. This fault is further interpolated andmarks the Makran subduction zone (SSZ 26) boundary in
the Indian Ocean. The zone contains the faults, including Chaman, Ghaziabad, and Ornach-Nal faults. Basedon InSAR data, the slip rate of the Ghazaband fault is
16.3 + 2.3 mm/yr (Fattahi and Amelung 2016). The active
sinistral Chaman fault (SSZ 20), with a length of morethan 850 km, forms a transpressional border between theEurasian and Indian plates in the northwest of Pakistan
and extends into Afghanistan, joining the dextral Herat
fault (SSZ 30,31) and the Pamir fault system (MKT in
Pakistan) (Lawrence et al. 1981). Based on InSAR data,
the slip rate of the Chaman fault is 8 + 3.1 mm/yr (Fattahi and Amelung 2016). The induced high seismicity inthe zone is described by the 1935 Quetta earthquake (Mw
7.6). The Chagai–Makran trench arc system is subdivided
into two tectonic blocks, namely (i) the Ras Koh-Chagaiarc (SSZ 27) and (ii) the Makran flysch basin (SSZ 26).
The Chagai and Mashkhel in the Sistan region represent
the tectonic setting of the Island Arc and Forearc, withdeep subduction megathrusts in the trench. The area
contains the Dalbandin and Ras Koh faults, in addition toseveral unidentified faults. The major earthquakes in thisblock include the 2011 Dalbandin earthquake (Mw 7.2)
and the 1983 Iran earthquakes (Mw 7.0).
The Makran flysch basin is divided into two tectonic
zones: (i) Central Makran and (ii) the Makran coastal

Page 12 of 31

trough. The Central Makran zone has intermediate subduction megathrusts represented by the Panjgur andHoshab Faults, which induced the 1929 earthquake (Mw6.5). The faults such as Ormara, Makran Coast faults,
and Makran subduction thrust faults all mark shallow subduction megathrusts in the Makran Coastal
Trough zone (SSZ 26). The induced high seismicity ofthis active tectonic zone is characterized by the 1945
Markran earthquake (tsunami) (Mw 8.2) and the 2013
Awaran earthquake (Mw 7.7). The Darvaz fault zone (thePamirs) is considered as SSZ 34. Deep seismicity (focal
depth > 50 km) is concentrated mostly in Pakistan’s
Himalaya–Hindukush–Karakoram mountainous area,
demarcated by zones (SSZ 35,36,37,38) and the Makran
subduction zone (SSZ 39) (Fig. 6b).
3.5 Source characterization

The seismic source zones are characterized by identifying maximum magnitude potential, the earthquake count,and the G–R parameters (Gutenberg and Richter 1956)
in each source zone. The maximum magnitude potentialof each SSZ is computed as: (1) by finding the maximum
magnitude value encountered in each source zone withthe addition of 0.5 units to handle the probable estimation
errors or (2) by taking the fault-based magnitude potential of the source zones. Waseem at el. (2020) have shownthat both methods produce magnitudes in a similar series.
As a result, the maximum historic earthquakes approachis used to calculate the maximum magnitude potential ofseismic source zones with an addition of 0.5 units to thehighest observed moment magnitude value of each sourcezone. The key advantage of using this approach is that it
does not abnormally reduce the hazard level.
3.6 Gutenberg–Richter (G–R) parameters of seismic source
zones

After completeness analysis, a declustered earthquake
catalogue and earthquake count of each seismic sourcezone are utilized to calculate the G–R parameters of each
SSZ following the G–R recurrence law (Gutenberg andRichter 1956). It is attained for all defined source zones
by the execution of least square regression analysis on
the logarithm of the cumulative number of earthquakes
and their respective moment magnitude values in eachsource zone. The standard G–R recurrence law accommodates an infinite number of earthquake magnitudes(Kramer 1996). As a result, the truncated G–R recurrence law is employed to predict the anticipated activity
rate in each source zone. Equation (7) gives the truncatedlaw formula:Rehman et al. Progress in Earth and Planetary Sciencemin(2023) 10:27

exp−β· Mw −Mw
 = v0 ·max
min
1 − exp−β·(Mw −Mw )(7)

whereas v0 = exp(α − β · Mwinf ) λ = activity rate, Mwmax
and Mwmin—the maximum and minimum moment magnitude bounds of the source. The G–R parameters are

Page 13 of 31a, b, α, and β. However, α = 2.303*a and β = 2.303*b. The
a-value suggests the cumulative seismicity of each zone.
The b-value is a ratio of large to small magnitudes, anda value around 1 indicates a seismically active zone.
A lower b-value indicates a higher potential for major
future earthquakes and vice versa.

Table 2 The table presents the characterization of seismic source zones based on the geological, tectonic, and seismological data,
including the location, area, maximum magnitude, seismicity rate, and recurrence intervals of earthquakes
Seismic sourcezone (SSZ)

No. ofevents (N)a valueb value

λ value

α value

β value

vo

Mmin

M1

Mmax

1

266

4.401

0.800

15.85310.1371.843

15.860

4.0

7.3

7.8

2

77

3.907

0.794

5.388

8.999

1.828

5.389

4.0

7.4

7.9

3

53

4.503

0.953

4.869

10.370

2.196

4.869

4.0

6.4

6.9

4

47

2.822

0.631

1.973

6.500

1.455

1.973

4.0

7.5

8.0

5

37

2.529

0.5814

1.600

5.825

1.338

1.599

4.0

7.4

7.9

6

13

2.270

0.621

0.611

5.228

1.430

0.611

4.0

6.4

6.9

7

51

3.974

0.850

3.720

9.152

1.959

3.720

4.0

6.7

7.2

8

37

2.999

0.685

1.808

6.907

1.578

1.808

4.0

6.9

7.4

9

165

4.585

0.884

11.117

10.5602.037

11.122

4.0

8.5

9.0

10

33

3.041

0.733

1.285

7.004

1.688

1.285

4.0

7.5

8.0

11

47

3.557

0.767

3.070

8.193

1.768

3.070

4.0

7.7

8.2

12

166

4.361

0.830

10.95210.044

1.912

10.956

4.0

7.1

7.6

13

10

2.953

0.814

0.495

6.802

1.876

0.495

4.0

5.7

6.2

14

86

4.283

0.865

6.618

9.864

1.993

6.619

4.0

7.5

8.0

15

61

4.309

0.912

4.544

9.923

2.102

4.544

4.0

6.6

7.1

16

144

4.567

0.896

9.621

10.519

2.063

9.624

4.0

7.1

7.6

17

101

5.078

1.035

8.672

11.695

2.383

8.674

4.0

7.0

7.5

18

71

4.330

0.894

5.645

9.973

2.060

5.647

4.0

7.7

8.2

19

74

5.721

1.196

8.646

13.175

2.754

8.649

4.0

6.8

7.3

20

95

3.921

0.789

5.826

9.031

1.817

5.828

4.0

7.7

8.2

21

73

4.210

0.882

4.786

9.697

2.032

4.787

4.0

7.5

8.0

22

88

3.065

0.611

4.18

7.059

1.407

4.181

4.0

6.4

6.9

23

157

5.455

1.085

13.026

12.5632.498

13.031

4.0

8.2

8.7

24

46

3.231

0.727

2.091

7.441

1.675

2.091

4.0

7.8

8.3

25

33

3.978

0.854

3.65

9.162

1.966

3.650

4.0

6.1

6.6

26

204

2.191

0.668

0.33

5.047

1.539

0.329

4.0

8.2

8.7

27

47

3.973

0.846

3.886

9.151

1.948

3.886

4.0

6.6

7.1

28

45

2.969

0.668

1.978

6.839

1.539

1.978

4.0

6.8

7.3

29

252

4.594

0.835

17.931

10.5801.923

17.940

4.0

7.0

7.5

30

25

3.378

0.806

1.427

7.781

1.856

1.427

4.0

7.5

8.0

31

45

3.860

0.852

2.834

8.890

1.962

2.834

4.0

7.0

7.5

32

542

5.667

0.998

47.19513.051

2.299

47.228

4.0

7.7

8.2

33

137

4.486

0.856

11.442

10.331

1.973

11.446

4.0

7.6

8.1

34

249

4.924

0.915

18.256

11.3402.108

18.265

4.0

7.3

7.8

35

391

4.114

0.703

19.911

9.475

1.620

19.922

4.0

7.9

8.4

36

81

4.546

0.957

5.214

10.470

2.204

5.215

4.0

7.5

8.0

37

36

2.699

0.620

1.65

6.217

1.429

1.649

4.0

7.4

7.9

38

38

3.302

0.747

2.062

7.605

1.720

2.061

4.0

7.7

8.2

39

36

2.820

0.631

1.972

6.494

1.453

1.972

4.0

7.5

8.0

Mmax = M1 (maximum observed magnitude in each source zone) + 0.5 units

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 14 of 31

SSZ 1

100

N = 10 4.4017- 0.8004Mw

10
1
0.1
0.01
0.001

4

5

6

Mw

7

8

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

G-R model for shallow seismic source zones of Highways network of CPEC
SSZ2

100

N = 10 3.9078- 0.7941Mw

10
1
0.1
0.01
0.001

4

5

1
0.1
0.01
0.001

4

5

Mw

6

7

8

SSZ5

100
10

N = 10 2.5296- 0.5814Mw

1
0.1
0.01
0.001

Cum. no. of earthquakes/year

N = 10 4.503- 0.9539Mw

10

4

5

Mw

6

6

7

8

SSZ4

100

7

8

Earthquake data

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

SSZ3

Mw

100
10

N = 10 2.8227- 0.6319Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

SSZ6

100

N = 10 2.2701- 0.621Mw

10
1
0.1
0.01
0.001

4

5

Mw

6

7

8

Fitting

Fig. 7 The Gutenberg–Richter model plots, illustrating the expected earthquake magnitude for seismic source zones (1–39) in Pakistan. The modelis based on historical seismic data and can be used to estimate the probability of an earthquake of a certain magnitude occurring in a specificsource zoneRehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 15 of 31

SSZ 8

10

N = 10 2.9992- 0.6855Mw

1
0.1
0.01

Cum. no. of earthquakes/year

4

5

Mw

6

7

8

SSZ9

100

N = 10 4.5856- 0.8849Mw

10
1
0.1
0.01
0.001

4

5

Mw

6

7

8

10

N = 10 3.5579- 0.7677Mw

1
0.1
0.01
4

5

Mw

6

7

8

Earthquake data
Fig. 7 continued100
10

N = 10 2.9992- 0.6855Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

SSZ10

100
10

N = 10 3.0414- 0.7331Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

SSZ12

SSZ11

100

0.001

Cum. no. of earthquakes/year

Com. no. of earthquakes/year

0.001

Cum. no. of earthquakes/year

SSZ 7

100

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

G-R model for Shallow Seismic Source Zones of CPEC100

N = 10 4.3615- 0.8305Mw

10
1
0.1
0.01
0.001

4

Fitting

5

Mw

6

7

8

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 16 of 31

10

N = 10 2.9536- 0.8147Mw

1
0.1
0.01
0.001

4

5

6
Mw

7

8

N = 10 4.309- 0.9129Mw

1
0.1
0.01
4

5

Mw

6

7

8

10

N = 10 5.0785- 1.0351Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

Earthquake data
Fig. 7 continued

1
0.1
0.01
0.001

SSZ17

100

N = 10 4.2835- 0.8657Mw

10

Cum. no. of earthquakes/year

10

SSZ14

100

SSZ15

100

0.001

Cum. no. of earthquakes/year

SSZ13

100

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

G-R model for shallow seismic source zones of Highways network of CPEC

4

5

6
Mw

7

8

SSZ16

100

N = 10 4.5676- 0.8961Mw

10
1
0.1
0.01
0.001

4

5

Mw

6

7

8

SSZ18

100

N = 10 4.3305- 0.8947Mw

10
1
0.1
0.01
0.001

4

5

Fitting

Mw

6

7

8

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 17 of 31

Cum. no. of earthquakes/year

SSZ19

100

N = 10 5.7212- 1.1961Mw

10
1
0.1
0.01
0.001

4

5

6

Mw

7

8

Cum. no. of earthquakes/year

SSZ21

100

N = 10 4.2108- 0.8827Mw

10
1
0.1
0.01
4

5

Mw

6

7

8

SSZ23

100

N = 10 5.4552- 1.0851Mw

10
1
0.1
0.01
0.001

4

5

6

7

8

Cum. no. of earthquakes/year

0.001

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

Cum. no. of earthquakes/year

G-R model for shallow seismic source zones of Highways network of CPEC
SSZ20

10

N = 10 3.9218- 0.7891Mw

1
0.1
0.01
0.001

4

Fig. 7 continued

Mw

6

7

8

SSZ22

100
10

N = 10 3.0652- 0.611Mw

1
0.1
0.01
0.001

4

5

6
Mw

7

8

SSZ24

100
10

N = 10 3.2311- 0.7277Mw

1
0.1
0.01
0.001

4

5

Mw

Earthquake data

5

Fitting

6
Mw

7

8

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 18 of 31

G-R model for shallow seismic source zones of Highways network of CPEC

SSZ26

10

Cumu. no. of earthquakes/year

100

N = 10 3.7987- 0.8541Mw

1
0.1
0.01

Cumu. no. of earthquakes/year

4

5

Mw

6

7

8

N = 10 3.9739- 0.8461Mw

1
0.1
0.01
4

5

10

Mw

6

7

8

0.01
0.001

4

6

N = 10 4.5944- 0.8352Mw

0.1
0.01
4

5

Mw

6

Mw

8

SSZ28

100
10

N = 10 2.9699- 0.6684Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

SSZ30

1

7

8

Earthquake data
Fig. 7 continued

0.1

SSZ29

100

0.001

1

Cumu. no. of earthquakes/year10

N = 10 3.8685- 0.7247Mw

10

SSZ27

100

0.001

100

Cumu. no. of earthquakes/year

0.001

Cumu. no. of earthquakes/year

Cumu. no. of earthquakes/year

SSZ25

100

N = 10 3.3788- 0.8061Mw

10
1
0.1
0.01
0.001

4

5

Fitting

Mw

6

7

8

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 19 of 31

G-R model for shallow seismic source zones of Highways network of CPEC
SSZ32

100
10

N = 10 3.8604- 0.852Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

100

Cumu. no. of earthquakes/year

Cumu. no. of earthquakes/year

SSZ31

N = 10 5.6671- 0.9983Mw

10
1
0.1
0.01
0.001

4

5

Mw

1
0.1
0.01
4

5

6

7

8

Mw

SSZ35

100
10

N = 10 4.1143- 0.7038Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

Earthquake data
Fig. 7 continued

Cumu. no. of earthquakes/year

N = 10 4.4861- 0.8569Mw

Cumu. no. of earthquakes/year

Cumu. no. of earthquakes/year
Cumu. no. of earthquakes/year

100

0.001

7

8

SSZ34

SSZ33
10

6

100
10

N = 10 4.9242- 0.9157Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

SSZ36

100
10

N = 10 4.5464- 0.9573Mw

1
0.1
0.01
0.001

4

Fitting

5

6
Mw

7

8

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 20 of 31

Cumu. no. of earthquakes/year

SSZ38

SSZ37

100
10

N = 10 2.6998- 0.6206Mw

1
0.1
0.01
0.001

4

5

6
Mw

7

8

Cumu. no. of earthquakes/year

Cumu. no. of earthquakes/year

G-R model for shallow and deep seismic source zones of Highways network of CPEC

100
10

N = 10 3.3026- 0.7471Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

SSZ39

100
10

N = 10 2.8202- 0.6313Mw

1
0.1
0.01
0.001

4

5

Mw

6

7

8

Earthquake data

Fitting

Fig. 7 continued

The threshold magnitude, Mw 4.0, was regarded as theminimum engineering magnitude that might cause damage
and is thus important in earthquake engineering. Table 2
shows the calculated G–R parameters for each seismic
source zone. Figure 7 shows the Gutenberg–Richter model
for seismic source zones (1–39) for the CPEC routes.
3.7 Ground motion prediction equations

The reliable ground motion prediction equations for aparticular location can be used to generate an accurate
estimate of ground motion and spectral acceleration
intensity values based on different rates of exceedance. Such an estimate is grounded on strong-motion
data in addition to the geology and tectonics of thearea. Ground motion prediction equations have notbeen developed in regions such as Pakistan; therefore,
selecting suitable GMPEs based on the region’s geology and tectonics is challenging. In Pakistan, there arethree broad tectonic settings: (1) active shallow crustal
regions such as the Himalayas, Sulaiman, Kirther, andChaman regions; (2) the Makran subduction interface
zone in the southwestern part of Pakistan; and (3) the
subduction intraslab zone in the Pamir, Hindukush,
Karakoram, and Himalayas (Waseem et al. 2020). The
current study’s criterion for selecting globally recognized GMPEs is based on similar tectonic settings.
In this study, the ground motion prediction equations developed for similar tectonic regions have been
selected randomly from the set of equations. Theseequations are assigned equal weights in the logic tree
and more than one equation is assigned to different tectonic settings of Pakistan. Strong motion data are also
scarce due to a lack of different data sharing agencies
in Pakistan responsible for earthquake data acquisition.
Only PGA values for unknown recording locations areavailable without the metadata. These PGA values areRehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 21 of 31

established for shallow earthquakes in northern Pakistan. Likewise, the data from the subduction interface
zone along the Makran Subduction Zone in Pakistan are
also not available. The only study about the evaluationof global ground motion prediction equations is available by Waseem et al. (2022) that was published after
the submission of the current study. Therefore, for thisstudy, random selection from the set of well-establishedground motion prediction equations from literature isconsidered. The ground motion data used in the analysis is divided into shallow and deep earthquakes basedon a 50-km focal depth dividing line. For the subduction
interface/intraslab regions, earthquake events with focal
depths greater than 50 km are considered; however,for the active shallow crust regions, earthquake eventswith focal depths up to 50 km are considered. The four
GMPEs are selected at random, and each tectonic zone
is assigned two GMPEs. The GMPEs of Boore et al.
(2014) and Akkar et al. (2014) are considered for shallow active areas (SSZ 1 to 34, except SSZ 26), whereas
Youngs et al. (1997) and Zhao et al. (2006) are considered for the subduction interface/intraslab (SSZ 26,
35–39).
The GMPE of Boore et al. (2014), developed using nextgeneration attenuation models, predicts the PGA for shallow crustal earthquakes and has a moderate functional
form. The required site parameters are Vs30 (average shear
wave velocity from surface up to a depth of 30 m), depth (inmeters) to 1 km/s shear wave velocity layer and magnitude.
The secondary variables are depth to the top of the ruptureand hypocentral depth. The function form of Boore et al.
(2014) for aleatory variability is presented in Eqs. (8, 9, 10,
11, 12):(8)σ = ∅2 + τ 2 τ1
τ = τ1 + (τ2 + τ1 )(Mw − 4.5)
τ
2

 �
�

 ∅�Mw , rjb ���ln V /V∅ = ∅ Mw , rjb − �∅V ln( (V2 /Vs,30) )
2
1

�
 �
∅ Mw , rjb − �∅VMw ≤ 4.5
4.5 < Mw < 5.5
Mw ≥ 5.5(9)
Vs,30 ≥ V2
V1 ≤ Vs,30 ≤ V2
Vs,30 ≤ V1

(10)

∅(Mw ) �
 rjb ≤ R1


� �

r
jb
ln
�
�
R1 �  R < r ≤ R∅ MW , rJb = ∅(Mw ) + �∅R  �
1
2
jb
R

2
ln

/R1



∅(Mw ) + �∅R

(11)


 ∅1
∅(Mw ) = ∅1 + (∅2 − ∅1 )(Mw − 4.5)
∅
2

Mw ≤ 4.5
4.5 < Mw < 5.5
Mw ≥ 5.5(12)where ∅ is intra-event, τ is inter-event variability,
R1 = 110, R2 = 270, ∅R = 0.10, ∅V = 0.07, V1 = 225,
V2 = 300, ∅1 = 0.695, ∅2 = 0.495, τ1 = 0.398 and τ2 = 0.348.
The GPME of Akkar et al. (2014) forecasts the pseudoresponse continuum for 5% damping up to 4.0 s and isdeveloped for active shallow crust regions. The equations are derived by using a reference data set of groundmotions in Europe (RESORCE) and the Middle East
region, which is comprised of 1041 records from 221
earthquakes. The required distance parameter is ‘Joyner–
Boore’ distance, the rupture parameter is magnitude, andthe site parameter is Vs30. The functional form of Akkar
et al. (2014) is presented in Eqs. (13,14,15) (for M ≤ 6.75)
is:
ln (Y ) = ln (YREF ) + ln[S Vs,30 , PGAREF(13)

where ln(YREF) is the reference ground motion model:� a1 + a2 (Mw − c1 ) + a3 (8.5 − Mw )2 + [a4 + a5 (Mw − c1 )] ln R2 + a2 + a8 FN + a9 FR Mw ≤ 6.75
6�
(14)
ln (YREF ) = a + a (M − c ) + a (8.5 − M )2 + [a + a (M − c )] ln R2 + a2 + a F + a F M > 6.75
1
7
w
1
3
w
4
5
w
1
8 N
9 R w
6



And ln[S9 Vs,30 , PGAREF is the nonlinear site amplification function:�n Vs,30/
�
VREF
�n 
�
b1 ln Vs,30 VREF + b2 ln 
��
� �
V
s,30
ln S Vs,30 , PGAREF =(PGAREF +C)
/
VREF

�
�


 b ln min (Vs,30 ,VCON )







�

1

�

VREF



PGAREF +C�Vs,30 ≤ VREF
Vs,30 > VREF

(15)

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

In these equations, R is the source-site distance measure; Y is the ground motion intensity measure in g; FN
and FR are the dummy variables used to represent thedifferent fault classes; VREF and VCON are the reference
Vs,30 in the nonlinear site model and the limiting Vs,30
after which the site amplification is a constant of 750 m/s
and 1000 m/s, respectively. The values agreed for othercoefficients are, a8 =  − 0.191, a9 = 0.0937, a5 = 0.2529,
a6 = 7.5, a7 =  − 0.5096, a2 = 0.0029 for all distance metrics. a1 = 1.85329, a3 =  − 0.02807, a4 =  − 1.23452,
c1 = 6.75, VREF = 750 ms , VCON = 1000 ms , b1 =  − 0.41997,
b2 =  − 0.28846, c = 2.5, n = 3.2.
The GMPEs of Youngs et al. (1997) and Zhao et al.
(2006) were developed for subduction intraslab andinterface regions such as Japan, Cascadia, and Alaska.
The required site parameter for GMPE Youngs et al.
(1997) is Vs30, rupture parameters are magnitude andfocal depth. The function of Youngs et al. (1997) is presented in Eq. (16) is:Cln PGA =C1∗ + C2 M + C3∗ ln rrup + e

C4∗ − C 2∗ M
3

+ C5 Zss + C8 Zt + C9 H

(16).
whereasIn this equation, PGA is measured in g, C1 = 0.2418,
C2 = 1.414, C3 =  − 2.552, C4 = ln(1.7818), C8 = 0.3846,
C9 = 0.00607, and σ = 1.45 − 0.1M . The other coefficients are not needed for prediction in rock. The tectonic
type of the equation is interface and intraslab with a focal
depth, H, between 10 and 229 km. In this equation, theeffect of depth and tectonic type is significant.
The GMPE of Zhao et al. (2006) implements the equations for the subduction interface. The required rupture parameters are the magnitude and focal depth.
C1∗ = C1 + C3 C4 − C3∗ C4∗ , C3∗ = C3 + C6 Zss , C4∗ = C4 + C7 Zss

Page 22 of 31

The function form of Zhao et al. (2006) is presented in
Eq. (17): loge y =aMw + bx − loge (r) + e(h − hc )δh(17)+ FR + SI + SS + SSL loge (x) + Ck
wherer = x + C exp (dMw )where y is in cm/s2, δh = 1 when h ≥ hc and 0 otherwise, a = 1.101, b =  − 0.00564, c = 0.0055, d = 1.080,
e = 0.014112,
SI = 0.000,
SS = 2.607,
SSL =  − 0.528,
CH = 0.293, C1 = 1.111, C2 = 1.344, C3 = 1.355, C4 = 1.420,
h = focal depth, hc = 15 km (best depth effect for shallow events), σ = 0.604 (intra-event) and τ = 0.398 (interevent). The equation classifies events into three source
types; (1) Crustal, (2) Interface SI, (3) Slab SS and SSL. FR isabout faulting mechanisms. Hard rock NEHRP site class
A, Vs,30 > 1100 m/s, CH is used. The equation provides areasonable distribution with respect to magnitude anddepth.
All the GMPEs employed are well recognized and areused to compute seismic hazard analysis. Figure 8 shows
the logic tree approach used in the study.
3.8 Epistemic uncertainties

The epistemic uncertainty in the computation of seismic
hazard analysis is concerned with the lack of knowledge
about earthquakes and ground motions. These uncertainties are well expressed using a logic tree, and the weightassigned to the respective logic tree branch illustrates theconfidence level of the analyst. The logic tree branches
are mainly associated with the approach adopted for thehazard analysis, source model, the selection of GMPEs,
and source characterization. Figure 8 shows the logic

Fig. 8 The logic tree approach, which assigns equal weights to each GMPE. This approach implies that each equation has an equal likelihood ofproviding accurate predictions of ground motion

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 23 of 31

Fig. 9 Seismic hazard maps of CPEC routes in terms of PGA for 475-year return period. a The province-wise maps show the SA values for eachprovince along the CPEC routes, b the district-wise maps provide a more granular view of the seismic hazard by showing the SA values forindividual districts

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 24 of 31

Fig. 10 Seismic hazard maps of CPEC routes in terms of PGA for 2475-year return period. a The province-wise maps show the SA values for eachprovince along the CPEC routes b the district-wise maps provide a more granular view of the seismic hazard by showing the SA values for individual
districts

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 25 of 31

Fig. 11 Seismic hazard maps of CPEC routes in terms of SA at T = 0.2 s for 2475-year return period. a The province-wise maps show the SA values
for each province along the CPEC routes b the district-wise maps provide a more granular view of the seismic hazard by showing the SA values forindividual districts

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

tree approach used in the study. The logic tree approachis only developed for attenuation equations in the current
study. Each attenuation equation is assigned equal weight(0.5) to preserve equal uncertainty in the assessment.

4 Results and discussion
The province- and district-wise seismic hazard maps for
the CPEC routes in terms of PGA for 475- and 2475-year
RPs are shown in Figs. 9 and 10. For engineering bedrocksite conditions, these maps are in agreement (i.e., Vs30
of 760 m/s). The color-coded maps signify the amplitude deviation of ground motion parameters. According to the overall spatial evaluation of seismic hazard
maps, the northern and western CPEC routes are prone
to the highest seismic hazard levels, i.e., the attained
PGA values are in the range of 0.40–0.52 g for 475-year
RP. This is due to the proximity of multiple active faults
in the Hindukush, Himalayan, Karakoram, and Sulaiman
mountains, whereas the central and eastern CPEC routes
have comparatively lower seismic hazard levels for 475year RP (Fig. 9a). The western CPEC route, which passes
through the districts of Ziarat, Quetta, Sibi and is next
to the Sulaiman fold belt and the Bella Chaman Kurram
fault zone, is likewise seismically hazardous, with PGA
values ranging from 0.40 to 0.52 g for 475-year RP. Intermediate and shallow focus earthquakes are frequent inthis region due to the presence of faults such as Kingri,
Waziristan-Loralai, Muslim Bagh, Chaman, and Ghaziabad. The northern CPEC route passes through Ghizer,
Diamer, Kohistan, Swat, Shangla, Mansehra, Haripur,
Islamabad, Attock, Swabi, and Peshawar districts, whichis adjacent to the Pamir–Karakoram Mountains, Hindukush–Himalaya ranges, Kohistan Island arc, and suturezones. This route is characterized by high seismicity, with
PGA values up to 0.52 g for 475-year RP (Fig. 9b). The
region is tectonically active due to the Indian–Eurasianplate collision and has numerous regional faults such as
Karakoram, North Pamir, South Pamir, Tirch Mir, Reshun, MMT, MCT, MBT, sinistral Jhelum, Nathiagali,
Panjal thrust, Khairabad, and Himalayan frontal thrust.
The CPEC route passing through the Mianwali, Zhob,
Kachhi, Kalat, Khuzdar, Panjgur, Kech, and Gwadar
districts bordering the Chagai arc and Makran subduction zone is recognized as the second most hazardous
route, with PGA values ranging between 0.32–0.40 g for
475-year RP (Fig. 9b). The well-known Dalbandin and
Ras-Koh faults, as well as a significant number of unidentified faults, are located in the Chagai and Mashkhel
areas. The 2011 Dalbandin earthquake (Mw 7.2) and the1983 Iran earthquake (Mw 7.0) indicated moderate seismicity in the area. In central Makran region, Panjgur andHoshab faults have generated several intermediate subduction megathrusts. The Makran coastal trough zone isPage 26 of 31

characterized by shallow subduction megathrusts such as
the Ormara fault, Makran coast fault, and Makran subduction thrusts (Fig. 1). The zone has high seismicity,
as evidenced by the 1945 Makran earthquake (Mw 8.2)
and the 2013 Awaran earthquake (Mw 7.7). The central
CPEC route is characterized by lower hazard values, i.e.,
PGA values of 0.22–0.24 g, except in a section of the D.G.
Khan district where the ground intensity is approaching
0.40 g for 475-year RP (Fig. 9b). The eastern CPEC route
passing through the Indian shield region has been recognized as the less hazardous route for 475-year RP, which
may be due to scattered seismicity and a limited numberof faults in the region. Figure 10 presents the provinceand district-wise color-coded maps for the 2475-year RPwith PGA values ranging between 0.43 and 0.93 g.
The spectral acceleration values (SA) computed at
T = 0.2 s (short-period) and T = 1.0 s (long-period) for
a 2% risk of exceedance in 50 years (2475-year RP) provided the hazard definition in line with the AASHTO(2012) (Specifications 2012) and IBC (2021) (Ching
and Winkel 2021) criteria. Figures 11 and 12 present
the province- and district-wise color-coded SA maps at
T = 0.2 s and 1.0 s for 2475-year RP. These maps are inagreement with engineering bedrock site conditions(i.e., Vs30 of 760 m/s). The spectral acceleration maps
reveal that the hazard is not uniformly spread over CPEC
routes in Pakistan, with increased hazards observed in
the northern and western CPEC routes (Fig. 11a). The
presence of faults such as Tirch Mir, Reshun, MBT, Kalabagh, Khairabad, east–west oriented MMT, and MKT in
the region results in the highest ground intensity values
encountered in the northern CPEC route (Fig. 11b). The
maximum ground intensities encountered in the western CPEC route passing through Ziarat, Quetta, Sibi,
Kechhi, and Khuzdar districts are owing to the presenceof active Chaman, Ghaziabad, and Chaman faults (Fig. 2).
For the 2475-year RP, the SA values (0.2 s) and SA (1.0 s)
range between 0.95–2.13 g and 0.24–0.64 g, respectively(Figs. 11, 12).
The current study has utilized an updated earthquake
catalogue (25AD-2020), reorganized completeness periods, introduced four updated GMPEs, computed activity
rates in each source zone, incorporated the 84th percentile of focal depths in each source zone, established a
15-km fishnet grid on either side of the CPEC routes andalso used the latest OpenQuake version 3.10 software
engine for hazard computation, which makes it exclusive to other studies carried out for the territory of Pakistan, e.g., Building Code of Pakistan (2007), Zaman et al.
(2012), and Waseem et al. (2020). Ground motion values
are calculated by taking into account both shallow anddeep-focus earthquakes. Instead of considering the average focal depth (50%) value within each source zone, theRehman et al. Progress in Earth and Planetary Science(2023) 10:27

Page 27 of 31

Fig. 12 Seismic hazard maps of CPEC routes in terms of SA at T = 1.0 s for 2475-year return period. a The province-wise maps show the SA values
for each province along the CPEC routes, while b the district-wise maps provide a more granular view of the seismic hazard by showing the SA
values for individual districts

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27

hypo depth of earthquakes within a source zone is determined by using the 84th percentile of the focal depths(the 50% + one standard deviation of focal depths). Taking the average focal depth (50%) value in a source zone
has the effect of abnormally increasing the hazard levelwithin each source zone.
The CPEC routes connect Khunjerab in North Pakistan to Gwadar in South Pakistan and extend throughout
the country via their western, central, and eastern routes(Fig. 2). Therefore, a comparison of earlier studies for the
territory of Pakistan can be furnished to the CPEC routes
in Pakistan based on obtained PGA values. The PGA
agreeing to a 10% probability of exceedance in 50 years(475-year RP) is a standard parameter for hazard computation, and hence, a comparison can be established basedon this parameter. The maximum PGA value obtained inthis study is 0.52 g for the period of 475-year RP, whichis consistent with the values obtained by Sesetyan et al.
(2018a) and higher than the values obtained by Waseem
et al. (2020). They used shear wave velocity of 760 m/s
as a site condition, EZ-FRISK software version 7.62, andthree sigma of ground motion prediction equations tocapture the variation in ground motions. Sesetyan et al.
(2018a) have used a shear wave velocity of 750 m/s as asite condition and three standard deviations of groundmotion prediction equations. The current study has used
760 m/s shear wave velocity and three standard deviations of equations that are comparable with Waseem
et al. (2020) and Sesetyan et al. (2018a) in terms of theseparameters. The BCP (2007) study used 760 m/s shear
wave velocity and one standard deviation of groundmotion prediction equation. Complete information is notavailable for all the mentioned studies; therefore, it may
be safely assumed that they have followed the variationin ground motion and rock site conditions to evaluate
the ground motion hazard. This study is also comparable to the studies of the BCP (2007), and the Global Seismic Hazard Assessment Programme (GSHAP) by Zhang
et al. (1999) due to similarity in the hazard computation
approach and source models. However, BCP has classified Pakistan into five seismic zones, with the highest
attainable PGA value marked as greater than 0.32 g andno upper bound value is provided for the 475-year RP.
The equivalent PGA values obtained in this study varyfrom 0.22–0.52 g for the CPEC routes in Pakistan, utilizing an aerial source for 475 year RP (Fig. 9). The spatial
distribution of these PGA values obtained for the CPEC
routes in Pakistan can be equally helpful in Pakistan’s
revised seismic zonation map. According to Zhang et al.
(1999), Rafi et al. (2012), Sesetyan et al. (2018a), andWaseem et al. (2020), the attained PGA values range
between 0.04–0.50 g, 0.004–0.22 g, 0.02–0.50 g, and
0.15–0.40 g, for 475-year RP, respectively. The results ofPage 28 of 31

Zhang et al. (1999) and Sesetyan et al. (2018a) agree withthe upper bound of attained PGA values in the current
study. However, the results of Waseem et al. (2020) areslightly lower than the upper bound but in agreementwith the lower bound of attained PGA values for the territory of Pakistan.
The ground motion intensity values obtained in thisstudy (i.e., PGA and SA) may be used for all types ofinfrastructure, particularly the spectral intensity values
for the seismic design of the bridges along with the CPEC
routes under AASHTO (2012) and IBC (2021) guidelines.
The PGA and SA maps can be a way forward to adopt
structural mitigation measures for all kinds of infrastructure being built around and across the CPEC routes inPakistan to make these routes secure and productive.

5 Conclusions
The result of the probabilistic seismic hazard assessment
offers rationalized seismic hazard maps for the CPEC
routes in Pakistan. This research leads to the following
conclusions:
1. The highest seismic hazard is apparent in the northern and western CPEC routes in Pakistan due toactive tectonic features in the vicinity of these sections. The industrial zones, transport infrastructure,and tunnels built along these routes are vulnerable,
even to small ground motion intensities.
2. The western CPEC route in Pakistan, which passes
through Ziarat, Quetta, and Sibi districts adjacent to
the Sulaiman fold belt and the Bella Chaman Kurram
fault zone, is the most seismically hazardous route,with attained PGA values ranging between 0.40 and
0.52 g for 475-year RP and SA (0.2 s) values ranging
between 1.66 and 2.13 g for 2475-year RP. The routes
with highest PGA values are more prone to seismic
risk and vulnerability. This highway section contains
many active faults in the neighborhood (i.e., the Chaman fault, Ghaziabad fault, Muslim Bagh fault, andWaziristan-Loralai fault). The planned CPEC infrastructure projects like hospitals, educational institutes, and roads build along this route are more vulnerable to seismic hazard.
3. The continuation of the western CPEC runs adjacent to the Makran subduction zone and Chagai arc,
considerably passing through Zhob, Kachhi, Kalat,
Khuzdar, Panjgur, Kech, and Gwadar districts arecharacterized as the second most hazardous route,i.e., PGA value ranges between 0.32 and 0.40 g for
475-year RP and SA (0.2 s) values ranges between
1.39 and 1.66 g for 2475-year RP. The planned CPEC
projects like the Gwadar airport and energy plants

Rehman et al. Progress in Earth and Planetary Science(2023) 10:27(coal, wind, solar) can be designed by following theattained accelerations values for earthquake-resistant
structures.
4. The central and eastern CPEC routes are generally
characterized by lower hazard values, i.e., PGA values
range between 0.22–0.24 g for 475-year RP and SA(0.2 s) values range between 0.95 and 1.13 g for 2475year RP, except in a highway section of D.G. Khan
district where the ground acceleration approaches
0.40 g for 475-year RP may be due to the presenceof Sulaiman fold and thrust belt (Zindapir anticlinalzone).
5. The peak ground acceleration (PGA) values achieved
at 10% and 2% probability of exceedance in 50 years
and spectral acceleration values (SA) obtained at
T = 0.2 s and T = 1.0 s for 2% probability of exceedance in 50 years provides hazard definition in accordance with the AASHTO LRFD (2012) and IBC (2021)
recommendations and can be used for seismic designof any type of infrastructure built along the CPEC
routes.
Acknowledgements
This work has been carried out as a part of PhD. dissertation at the National
Centre of Excellence in Geology, University of Peshawar, Pakistan. The firstauthor highly acknowledges the online earthquake data sources from wheredata is collected, alike International Seismological Centre, United States
Geological Survey, the local networks including Water and Power Development Authority Pakistan, Pakistan Meterological Department and the National
Geophysical Data Center. Additionally, the comments and recommendationsof the editor and anonymous reviewers are thankfully acknowledged as having performed a constructive review, improving the quality of the article.
Author contributions
Conceptualization, Q.R., M.W., W.A.; methodology, Q.R., M.W., W.A.; software,
Q.R.; validation, Q.R., M.W., W.A., I.I, H.T.J., G.K.; formal analysis, Q.R.; investigation,
Q.R., M.W., W.A., I.I, H.T.J., G.K.; resources, M.W., W.A.; data curation, Q.R., M.W.,
W.A.; writing—original draft preparation, W.A.; writing—review and editing,Q.R., M.W., W.A.; I.I, H.T.J., G.K.; visualization, Q.R., M.W., W.A, I.I, H.T.J, G.K.; supervision, M.W., W.A. All authors read and approved the final manuscript.
Funding
This research received no external funding.
Availability of data and materials
The data used in this work is available on request.

Declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The author declares no competing interests.
Author details
1
National Centre of Excellence in Geology, University of Peshawar, Peshawar 25120, Pakistan. 2 Department of Earth Sciences, University of Haripur,
Haripur 22610, Pakistan. 3 Department of Civil Engineering, University

Page 29 of 31

of Engineering and Technology, Peshawar 25120, Pakistan. 4 Departmentof Geology, Shaheed Benazir Bhutto University Sheringal, Upper Dir 18000,
Khyber Pakhtunkhwa, Pakistan. 5 Department of Historical Geology‑Paleontology, Faculty of Geology and Geoenvironment, School of Earth Sciences,
National and Kapodistrian University of Athens, Panepistimiopolis, Zografou,
15784 Athens, Greece.
Received: 24 May 2022 Accepted: 5 May 2023

References
AASHTO-LRFD (2012) Bridge design specifications, 6th ed. American Association of State Highway and Transportation Officials, Washington, DC. ISBN:
978-1-56051-523-4
Akkar S, Sandıkkaya MA, Bommer JJ (2014) Empirical ground-motion models
for point-and extended-source crustal earthquake scenarios in Europe
and the Middle East. Bull Earthq Eng 12(1):359–387. https://​doi.​org/​10.​
1007/​s10518-​013-​9461-4
Basharat M, Rohn J, Ehret D, Baig MS (2012) Lithological and structural control
of Hattian Bala rock avalanche triggered by the Kashmir earthquake 2005,
sub-Himalayas, northern Pakistan. J Earth Sci 23(2):213–224. https://​doi.​
org/​10.​1007/​s12583-​012-​0248-3
BCP. Building Code of Pakistan with Seismic Provisions (2007) Minist. Hous.
Work. Government of Pakistan Islamabad
Bignold SM, Treloar PJ (2003) Northward subduction of the Indian Platebeneath the Kohistan island arc, Pakistan Himalaya: new evidence fromisotopic data. J Geol Soc Lond 160(3):377–84. https://​doi.​org/​10.​1144/​
0016-​764902-​068
Bilham R (1999) Slip parameters for the Rann of Kachchh, India, 16 June 1819,
earthquake, quantified from contemporary accounts. Geol Soc Lond
Spec Publ 146(1):295–319. https://​doi.​org/​10.​1144/​GSL.​SP.​1999.​146.​01.​18
Bilham R, Lodi S, Hough S, Bukhary S, Khan AM, Rafeeqi SFA (2007) Seismic
hazard in Karachi, Pakistan: uncertain past, uncertain future. Seismol Res
Lett Seismol Soc Am 78(6):601–613. https://​doi.​org/​10.​1785/​gssrl.​78.6.​601
Boore DM, Stewart JP, Seyhan E, Atkinson GM (2014) NGA-West2 equations forpredicting PGA, PGV, and 5% damped PSA for shallow crustal earthquakes. Earthq Spectra Earthq Eng Res Inst 30(3):1057–1085. https://​doi.​
org/​10.​1193/​07011​3EQS1​84M
Chang M, Cui P, Dou X, Su F (2021) Quantitative risk assessment of landslides
over the China–Pakistan economic corridor. Int J Disaster Risk Reduct
63:102441. https://​doi.​org/​10.​1016/j.​ijdrr.​2021.​102441
Chaudhry MN, Ghazanfar M (1990) Position of the Main Central Thrust in thetectonic framework of Western Himalaya. Tectonophysics 174(3–4):321–
329. https://​doi.​org/​10.​1016/​0040-​1951(90)​90329-7
Ching FDK, Winkel SR (2021) Building codes illustrated: a guide to understanding the 2021 international building code. Wiley. ISBN: 978-1-119-77240-8
Cook N, Butz D (2015) The dialectical constitution of mobility and immobility:
recovering from the Attabad Landslide disaster, Gojal, Gilgit-Baltistan,
Pakistan. Contemp South Asia 23(4):388–408. https://​doi.​org/​10.​1080/​
09584​935.​2015.​10909​50
Ding L, Qasim M, Jadoon IAK, Khan MA, Xu Q, Cai F, Wang H, Baral U, Yue Y(2016) The India–Asia collision in north Pakistan: insight from the U–Pb
detrital zircon provenance of Cenozoic foreland basin. Earth Planet Sci
Lett 455:49–61. https://​doi.​org/​10.​1016/j.​epsl.​2016.​09.​003
Fan F, Xiao C, Feng Z, Chen Y (2022) Land-planning management based on
multiple ecosystem services and simulation in tropical forests. J Environ
Manag 323:116216. https://​doi.​org/​10.​1016/j.​jenvm​en.​2022.​116216
Farah A, DeJong KA (1984) Geodynamics of Pakistan. Geological Survey ofPakistan
Farah A, Lawrence RD, DeJong KA (1979) An overview of the tectonics ofPakistan. Mar Geol Oceanogr Arab Sea Coast Pak. Van Nostrand Reinhold
Company, New York, pp 161–176
Fattahi H, Amelung F (2016) InSAR observations of strain accumulation andfault creep along the Chaman Fault system, Pakistan and Afghanistan.
Geophys Res Lett 43(16):8399–8406. https://​doi.​org/​10.​1002/​2016G​L0701​
21
Fujiwara S, Tobita M, Sato HP, Ozawa S, Une H, Koarai M, Nakai H, Fujiwara M,
Yarai H, Nishimura T, Hayashi F (2006) Satellite data gives snapshot of theRehman et al. Progress in Earth and Planetary Science(2023) 10:27

2005 Pakistan earthquake. Eos Trans Am Geophys Union 87(7):73–77.
https://​doi.​org/​10.​1029/​2006E​O0700​01
Gardner JK, Knopoff L (1974) Is the sequence of earthquakes in Southern
California, with aftershocks removed, Poissonian? Bull Seismol Soc Am
Citeseer 64(5):1363–1367
Gutenberg B, Richter CF (1956) Earthquake magnitude, intensity, energy, andacceleration: (second paper). Bull Seismol Soc Am 46(2):105–145. https://​doi.​org/​10.​1785/​BSSA0​46002​0105
Irshad MS (2015) One belt and one road: Dose China–Pakistan economic corridor benefit for Pakistan’s economy? J Econ Sustain Dev 6(24):66
Jadoon IAK, Khurshid A (1996) Gravity and tectonic model across the Sulaiman
fold belt and the Chaman fault zone in western Pakistan and eastern
Afghanistan. Tectonophysics 254(1–2):89–109. https://​doi.​org/​10.​1016/​
0040-​1951(95)​00078-X
Jadoon IAK, Lawrence RD, Lillie RJ (1992) Balanced and retrodeformed
geological cross-section from the frontal Sulaiman Lobe, Pakistan: duplex
development in thick strata along the western margin of the Indian Plate.
Thrust Tectonics 66:343–356
Jadoon IAK, Lawrence RD, Lillie RJ (1994) Seismic data, geometry, evolution,and shortening in the active Sulaiman fold-and-thrust belt of Pakistan,southwest of the Himalayas. Am Assoc Pet Geol Bull 8(5):758–774
Jagoutz O, Schmidt MW (2012) The formation and bulk composition of modern juvenile continental crust: the Kohistan arc. Chem Geol 298:79–96.
https://​doi.​org/​10.​1016/j.​chemg​eo.​2011.​10.​022
Karim D, Karim I, Daveel S, Khan A (2015) Damming hunza river by massive
Attabad landslide, story of a risk management initiative from Hunza,
Pakistan. J Mt Area Res 5:1–9. https://​doi.​org/​10.​53874/​jmar.​v5i0.​76
Khan S, Waseem M, Khan MA, Ahmed W (2018) Updated earthquake catalogue for seismic hazard analysis in Pakistan. J Seismol 22(4):841–861.
https://​doi.​org/​10.​1007/​s10950-​018-​9736-y
Kramer SL (1996) Geotechnical earthquake engineering. Pearson Education,
India
Lawrence RD, Yeats RS, Khan SH, Farah A, DeJong KA (1981) Thrust and strike
slip fault interaction along the Chaman transform zone, Pakistan. Geol
Soc Lond Spec Publ 9(1):363–70. https://​doi.​org/​10.​1144/​GSL.​SP.​1981.​
009.​01
Mahmood K, Ahmad N, Khan U, Iqbal Q (2020) Seismic hazard maps ofPeshawar District for various return periods. Nat Hazards Earth Syst Sci
20:1639–1661. https://​doi.​org/​10.​5194/​nhess-​20-​1639-​2020
Makhdoom AS, Shah AB, Sami K (2018) Pakistan on the roadway to socio-economic development: a comprehensive study of China Pakistan Economic
Corridor (CPEC). Gov Res J Polit Sci 6(6):66
Malkani MS, Mahmood Z (2016) Fluorite from Loralai-Mekhtar and Celestite
from Barkhan, Dera Bugti, Kohlu, Loralai and Musakhel districts (Sulaiman
Foldbelt) and Karkh area of Khuzdar district (Kirthar Range): a glimpse on
Tectonic and Sedimentary Mineral Resources of Indus Basin Pak. Geol.
Surv. Pakistan. Inf Release 980:1–16
Monalisa KAA, Jan MQ (2007) Seismic hazard assessment of the NW Himalayanfold-and-thrust belt, Pakistan, using probabilistic approach. J Earthq Eng
11(2):257–301. https://​doi.​org/​10.​1080/​13632​46060​10312​43
Mukhopadhyay B, Dasgupta S (2015) Seismic hazard assessment of Kashmirand Kangra valley region, Western Himalaya, India. Geomat Nat Hazards
Risk 6(2):149–183. https://​doi.​org/​10.​1080/​19475​705.​2013.​832405
Naseer A, Khan AN, Hussain Z, Ali Q (2010) Observed seismic behavior of buildings in northern Pakistan during the 2005 Kashmir earthquake. Earthq
Spectra Earthq Eng Res Inst 26(2):425–449. https://​doi.​org/​10.​1193/1.​
33831​19
Pagani M, Silva V, Rao A, Simionato M, Johnson K (2023) OpenQuake engine
manualPetterson MG (2010) A review of the geology and tectonics of the Kohistan
island arc, north Pakistan. Geol Soc London Spec Publ 338(1):287–
327.https://​doi.​org/​10.​1144/​SP338.​14
Petterson MG, Windley BF (1985) RbSr dating of the Kohistan arc-batholith in
the Trans-Himalaya of north Pakistan, and tectonic implications. Earth
Planet Sci Lett 74(1):45–57. https://​doi.​org/​10.​1016/​0012-​821X(85)​90165-7
Qasim M, Ding L, Khan MA, Jadoon IA, Haneef M, Baral U, Cai F, Wang H, Yue
Y (2018) Tectonic implications of detrital zircon ages from lesser Himalayan Mesozoic-Cenozoic strata, Pakistan. Geochem Geophys Geosyst
19(5):1636–1659. https://​doi.​org/​10.​1002/​2017G​C0068​95

Page 30 of 31

Rafi Z, Lindholm C, Bungum H, Laghari A, Ahmed N (2012) Probabilistic seismic
hazard of Pakistan, Azad-Jammu and Kashmir. Nat Hazards 61(3):1317–
1354. https://​doi.​org/​10.​1007/​s11069-​011-​9984-4
Rahman MM, Bai L, Khan NG, Li G (2019) Probabilistic seismic hazard assessment for Himalayan–tibetan region from historical and instrumental
earthquake catalogs. Earthqu Multi-Hazards Around Pac Rim 2:161–181.
https://​doi.​org/​10.​1007/​s00024-​017-​1659-y
Rehman G, Ahmad S, Khan SD, Ali F, Ali TH, Khan SF (2014) Threat of glacial lake
outburst flood to Tehsil Gupis from Khukush Lake, District Ghizer, Gilgit
Baltistan, Pakistan. Nat Hazards 70(2):1589–1602. https://​doi.​org/​10.​1007/​
s11069-​013-​0893-6
Rehman HU, Seno T, Yamamoto H, Khan T (2011) Timing of collision of theKohistan–Ladakh Arc with India and Asia: debate. Isl Arc 20(3):308–328.
https://​doi.​org/​10.​1111/j.​1440-​1738.​2011.​00774.x
Rehman QU, Ahmed W, Waseem M, Khan S, Farid A, Shah SHA (2021) Geophysical investigations of a potential landslide area in Mayoon, Hunza
District, Gilgit-Baltistan, Pakistan. Rud Zb Bull 36(3):66
Searle MP, Khan MA, Fraser JE, Gough SJ, Jan MQ (1999) The tectonic evolutionof the Kohistan-Karakoram collision belt along the Karakoram Highway
transect, north Pakistan. Tectonics 18(6):929–949. https://​doi.​org/​10.​1029/​
1999T​C9000​42
Sesetyan K, Danciu L, Demircioglu MB, Giardini D, Erdik M, Akkar S, Gülen
L, Zare M, Adamia S, Ansari A, Arakelyan A (2018a) The 2014 seismic
hazard model of the Middle East: overview and results. Bull Earthq Eng
16:3535–3566. https://​doi.​org/​10.​1007/​s10518-​018-​0346-4
Sesetyan K, Demircioglu MB, Duman TY, Can T, Tekin S, Azak TE, Fercan ÖZ(2018b) A probabilistic seismic hazard assessment for the Turkish territory—part I: the area source model. Bull Earthq Eng 16(8):3367–3397.
https://​doi.​org/​10.​1007/​s10518-​016-​0005-6
Specifications LBD (2012) American Association of State Highway and Transportation Officials (AASHTO). Washington, DC, USA
Tong L (2014) CPEC industrial zones and China–Pakistan capacity cooperation.
Strateg Stud JSTOR 34:174–84
Treloar PJ, Potts GJ, Wheeler J, Rex DC (1991) Structural evolution and asymmetric uplift of the Nanga Parbat syntaxis, Pakistan. Himalya Geol Rundschau 80(2):411–428. https://​doi.​org/​10.​1007/​BF018​29374
Waseem M, Lai CG, Spacone E (2018) Seismic hazard assessment of northern Pakistan. Nat Hazards 90(2):563–600. https://​doi.​org/​10.​1007/​
s11069-​017-​3058-1
Waseem M, Khan MA, Khan S (2019) Seismic sources for southern Pakistan andseismic hazard assessment of Karachi. Nat Hazards 99(1):511–536. https://​doi.​org/​10.​1007/​s11069-​019-​03755-5
Waseem M, Khan S, Asif Khan M (2020) Probabilistic seismic hazard assessmentof pakistan territory using an areal source model. Pure Appl Geophys
177:3577–3597. https://​doi.​org/​10.​1007/​s00024-​020-​02455-7
Waseem M, Rehman ZU, Sabetta F, Ahmad I, Ahmad M, Sabri MMS (2022)
Evaluation of the predictive performance of regional and global groundmotion predictive equations for shallow active regions in Pakistan. Sustainability MDPI 14(13):8152. https://​doi.​org/​10.​3390/​su141​38152
Wiemer SA (2001) Software package to analyze seismicity: ZMAP. Seismol Res
Lett Seismol Soc Am 72(3):373–82. https://​doi.​org/​10.​1785/​gssrl.​72.3.​373
Yeats RS, Lawrence RD (1982) Tectonics of the Himalayan thrust belt in northern Pakistan. masc
Youngs RR, Chiou SJ, Silva WJ, Humphrey JR (1997) Strong ground motion
attenuation relationships for subduction zone earthquakes. Seismol Res
Lett Seismol Soc Am 68(1):58–73. https://​doi.​org/​10.​1785/​gssrl.​68.1.​58
Zaigham NA, Mallick KA (2000) Bela ophiolite zone of southern Pakistan:
tectonic setting and associated mineral deposits. Geol Soc Am Bull
112(3):478–489. https://​doi.​org/​10.​1130/​0016-​7606(2000)​112<​478:​
BOZOS​P>2.​0.​CO;2
Zaman S, Ornthammarath T, Warnitchai P (2012) Probabilistic seismic hazard
maps for Pakistan. In: Proceedings of the 15th world conference earthquake engineering WCEE, pp 1–10
Zanchi A, Gaetani M (2011) The geology of the Karakoram range, Pakistan:
the new 1: 100,000 geological map of Central-Western Karakoram. Ital
J Geosci Società Geologica Italiana 130(2):161–262. https://​doi.​org/​10.​
3301/​IJG.​2011.​09
Zare M, Amini H, Yazdi P, Sesetyan K, Demircioglu MB, Kalafat D, Erdik M,
Giardini D, Khan MA, Tsereteli N (2014) Recent developments of theRehman et al. Progress in Earth and Planetary Science(2023) 10:27

Middle East catalog. J Seismol 18(4):749–772. https://​doi.​org/​10.​1007/​
s10950-​014-​9444-1
Zhang P, Yang Z, Gupta HK, Bhatia SC, Shedlock KM (1999) Global seismic hazard assessment program (GSHAP) in continental Asia. http://​hdl.​handle.​
net/​2122/​1387
Zhao JX, Zhang J, Asano A, Ohno Y, Oouchi T, Takahashi T, Ogawa H, Irikura
K, Thio HK, Somerville PG, Fukushima Y (2006) Attenuation relationsof strong ground motion in Japan using site classification based on
predominant period. Bull Seismol Soc Am Seismol Soc Am 96(3):898–913.
https://​doi.​org/​10.​1785/​01200​50122

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Page 31 of 31