1 Introduction

Seismic hazard studies on a regional scale are carried out using either a deterministic seismic hazard analysis (DSHA) or probabilistic seismic hazard analysis (PSHA) approach. DSHA considers a worst-case scenario for a particular seismic source, referred as the Maximum Credible Earthquake (MCE), to estimate the level of ground shaking (Reiter 1990; Kramer 1996; Bommer 2002). On the other hand, PSHA considers the likelihood of various earthquakes from multiple potential seismic sources, each having a range of uncertainty in source characteristics (e.g., rupture length, distance, fault dip, maximum magnitude, slip rate) (Abrahamson 2000; Bommer 2002). Uncertainty is treated explicitly, and the annual probability of exceeding a specified ground motion is computed.

Deterministic earthquake ground motion simulation is a common technique for estimating seismic-induced ground shaking (Wang 2015). The technique uses numerical methods and models that incorporate the physics of an earthquake source and the propagation of seismic waves (Taborda and Roten 2015). Deterministic scenarios are useful for simulating worst-case events that could affect an area (McGuire 2001). An example is the Great Southern California ShakeOut, in which an M 7.8 scenario earthquake, rupturing the southern segment of the San Andreas Fault, was simulated (Bielak et al. 2010). Other examples are the simulation of scenario earthquakes along the Central Marmara Fault (CMF) and North Boundary Fault (NBF) in Turkey (Pulido et al. 2004), and the simulation of an Mw 6 scenario earthquake in the Grenoble Valley (France) by Stupazzini et al. (2009).

For the region exposed to the 2005 Kashmir earthquake in northern Pakistan, only PSHA approaches have been used till date. MonaLisa et al. (2008) carried out a PSHA for the seismically active NW Himalayan Fold and Thrust Belt in Pakistan. On a larger scale, PMD and NORSAR (2007), Zaman et al. (2012) and Sultan (2015) conducted a probabilistic seismic hazard assessment covering entire Pakistan. However, in PSHA the near-surface effects of an earthquake, such as topographic amplification, are often not included which significantly influence the spatial distribution of seismic response. Depending on the surface geometry, topography could lead to scattering and/or focusing of propagating seismic waves (Lee et al., 2008, 2009a). Previous studies show that topography amplifies seismic-induced ground shaking at mountain ridges, while it de-amplifies in valleys (Hartzell et al. 1994; Spudich et al. 1996; Bouchon et al. 1996; Assimaki et al. 2005; Lee et al. 2009a, b; Hough et al. 2010; Kumagai et al. 2011; Takemura et al. 2015; Khan et al. 2017, 2020). The directions of incident seismic waves also have an impact on the amplification of ground shaking (Pedersen et al. 1994; Ashford and Sitar 1997; Ashford et al. 1997; Shafique et al. 2008; Meunier et al. 2008; Qi et al. 2010; Xu et al. 2013). Topography-induced scattering of seismic energy and its effect on ground shaking (de-)amplification has been studied with numerical simulations (Lee et al. 2008; Takemura et al. 2015).

The Mw 7.6 2005 Kashmir earthquake raised concerns for another major earthquake in the region (Durrani et al. 2005). The earthquake occurred on the Muzaffarabad Fault (MF) (Hussain et al. 2009), but the nearby Murree Fault (MBT) and Punjal Fault (PF) are also potential sources for a major earthquake (MonaLisa et al. 2009). Wyss (2006) and Bilham and Hough (2006) predict large magnitude future earthquakes for the region with a high number of expected casualties. In this paper, we use spectral element method (SEM), for the first time to estimate ground motion for the larger Muzaffarabad region, the results of which are applicable both on regional and local scales. SEM has been widely used in simulating 3D seismic response (e.g., Komatitsch and Vilotte 1998; Lee et al. 2008; Stupazzini et al. 2009; Chaljub et al. 2015). This approach allows simulating seismic waves by incorporating source, path and site effects (e.g., Atkinson and Macias 2009; Stupazzini et al. 2009; Bielak et al. 2010). Moreover, it also allows incorporation of topographic effects, which has not been addressed in the previous studies.

The approach adopted in this study is a mix of DSHA and PSHA, where it considers the worst-case scenario in the form of the 2005 Kashmir earthquake as source and makes use of multiple potential sources (by changing its location). Moment tensor (MT) solution of the 2005 Kashmir earthquake (Dziewonski et al. 1981; Ekström et al. 2012) is used as a source for a total of 25 scenario earthquakes using DSHA, which includes the actual 2005 Kashmir earthquake. The results from these scenarios are combined in a probabilistic seismic hazard map for the region showing the probable ground displacement due to topography and its variability.

2 Study area

The city of Muzaffarabad lies in an active tectonic zone and is surrounded by several major active faults (Fig. 1). This region is prone to large magnitude earthquakes (MonaLisa et al. 2009) as it is situated on the collision zone where the Indian plate subducts underneath the Eurasian plate at a rate of 31 mm/year (Bettinelli et al. 2006). Apart from the 2005 Kashmir earthquake along the Muzaffarabad Fault (MF), several other earthquakes greater than Mw 5 have been recorded, as well as two (Mw 6 and Mw 6.4) aftershocks of the 2005 Kashmir earthquake (Tahirkheli 2010). Studies indicate that the energy stored in the Himalayas gives a high probability of earthquakes of Mw >8.0 in the future (Durrani et al. 2005; Stevens and Avouac 2016). The rugged topography of the study area, with close to 4 km elevation difference, makes it prone to topographic (de-)amplification (Lee et al. 2009b; Hough et al. 2010; Shafique and van der Meijde 2015).

Fig. 1
figure 1

Major faults in the region surrounding the city of Muzaffarabad, overlaid on an ASTER GDEM elevation model (resampled to 270 m spatial resolution). The moment tensor (MT) solution of the 2005 Kashmir earthquake is indicated by the red beachball icon. To test the sensitivity of the region to earthquakes along the regional faults, the same MT is placed at different locations along the regional faults. The locations of these hypothetical sources are indicated by red dots and labeled with source number for reference

Full size image

3 Materials and methods

We use SEM for simulating 3D seismic wave propagation to estimate ground shaking. SEM was originally developed by Patera (1984) for computational fluid dynamics. It was adapted for 3D seismic wave propagation by Komatitsch and Vilotte (1998) and Komatitsch and Tromp (1999). SEM incorporates a 3D physical model, considering the velocity, density and free surface topography for a full elastic waveform simulation including all possible waves, based on source characteristics of an earthquake (Chaljub et al. 2007). We used SPECFEM3D Cartesian (Computational Infrastructure for Geodynamics 2016) software as SEM simulator in this study.

Before a SEM simulation is run, a 3D mesh model is developed to represent the physical characteristics of area under investigation (Casarotti et al. 2008). The mesh is composed of hexahedra elements that are isomorphous to a cube (Komatitsch et al. 2002). The mesh elements have material and structural properties that define how it reacts to applied conditions such as an earthquake. We used Cubit v.13.0 software, developed by (Sandia National Laboratories 2011) in conjunction with GeoCubit (Computational Infrastructure for Geodynamics 2016) for meshing. GeoCubit is developed at the Istituto Nazionale di Geofisica e Vulcanologia (INGV) in Italy. This software automates mesh generation in three steps. The first step is the construction of a 3D geometry. GeoCubit creates points from a DEM, combines these points to make lines and eventually creates a surface by combining the lines. After creating the topographic surface, a 3D geometrical model is created. The second step is meshing the geometry with hexahedral elements. The last step involves defining boundary conditions, such as free and absorbing surfaces, and mesh export for simulation with SEM.

The surface topography of the geometric model is based on Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global DEM data. It is retrieved from the Global Data Explorer, courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, http://gdex.cr.usgs.gov/gdex/. The geometric model extends to a depth of 40 km. Introduction of two tripling layers to the mesh increases the spatial resolution of the mesh from 2430 m at the bottom depth of 40 km to 270 m at the topographic surface. There are no seismic waves speed models available for this area: Tomographic velocity models are too coarse; the best resolution is 1 degree (Johnson and Vincent 2002). Seismic lines are not available for the region, and nearby seismic lines (Bhukta and Tewari 2007) cannot provide sufficient detail on the Muzaffarabad region. We therefore assign constant wave speeds (Vp = 2800 m/s, Vs = 1500 m/s) and density (ρ = 2300 km/m3) throughout the model, which is representative of upper crustal conditions (Johnson and Vincent 2002; Bhukta and Tewari 2007). The mesh can resolve frequencies up to approximately 5.5 Hz (Khan et al. 2017). The topographic surface was set as free surface; the remaining 5 plane surfaces are set as absorbing surfaces. In order to estimate topography-related (de-)amplification, we also construct a mesh model with same properties but without topography (i.e., having a plain top surface). Amplification maps are based on the ratio between maps with and without topography, respectively, of modeled maximum displacement for each model point.

The seismic point source is the MT solution of the 2005 Kashmir earthquake. In order to observe the impact of possible future earthquakes, the same MT is placed at 24 different locations at approximately 6-km intervals along major faults (MF, MBT, PF) in the study area (Fig. 1). A total number of 25 simulations were performed to give a fair representation for every potential epicenter point. More simulations could be performed by reducing the interval; however, this will get computationally more expensive. We believe that with an interval of 6 km, each scenario earthquake would have only 3 km error in location on each side on the causative fault. We consider this to be acceptable for earthquakes of this magnitude (with rupture lengths of tens of kilometers). All simulations are performed in SPECFEM3D with a delta source time function of half duration ``0'' at a duration of 35 s. Since SEM is efficient in simulating low frequency ground displacement and has limited capability in simulation of high frequency accelerations (Dhanya et al. 2016), our analysis utilizes the shakemaps obtained from these simulations in the form of peak ground displacement (PGD) (Fig. 2). A mean of PGDs (Fig. 3a) obtained from all 25 simulations is calculated together with the corresponding standard deviations (Fig. 3b). The mean and standard deviation are subsequently used to create a seismic hazard map (Fig. 3c).

Fig. 2
figure 2

Peak ground displacement (PGD) obtained from 25 simulations with changing location of the moment tensor (MT). Simulation 0 has the original source location of the 2005 Kashmir earthquake. The source is represented by a yellow circle and number corresponding to the source locations shown in Fig. 1. It can be observed that in each image the spatial pattern of the PGD amplitudes is influenced by the topography and the changing source location

Full size image
Fig. 3
figure 3

a Mean peak ground displacement from all 25 simulations, showing high average values on ridges (red tones) and low average values in valleys (blue tones). b The standard deviation of peak ground displacement, indicating higher variation on ridges compared to valleys. c Topography-based seismic hazard map for the study area derived from mean PGD (a) and its standard deviation (b). The red area indicates a relatively high hazard, yellow indicates medium, and green indicates a relatively low hazard area. The saturation represents the standard deviation, and whitish tones therefore represent uncertainty (more white is more uncertainty, and vice versa). d Digital elevation model of the study area

Full size image

The seismic hazard map is shown as a hue, saturation and value (HSV) color composite. The mean PGD was put as hue with a color ramp of green to red, and standard deviation (rescaled to a range of 0–1 and inverted) as saturation. The value was uniformly given a value of 1. This results in a map (Fig. 3c) that has areas with low to high mean PGD in green to red tones, respectively, and with low standard deviation areas in darker tones and high standard deviation areas in lighter (white) tones.

To obtain amplification maps, the PGD shakemaps obtained from plain surface models are subtracted from the shakemaps of models with topography. Amplification maps are presented for each location along different faults (Fig. 4). The average of all amplification maps was used with its corresponding standard deviation to show the variability of PGD amplitudes with respect to topography and changing source location (Fig. 5), indicating regions that might experience strong topographic effects.

Fig. 4
figure 4

Difference of peak ground displacements (∆PGDs) obtained from 25 simulations. Simulation 0 is the original 2005 Kashmir earthquake. Each source is represented by a yellow circle and a number corresponding to the source locations shown in Fig. 1. The pattern of amplification (positives values) and reduction (negative values) follows the topography and is influenced by the changing location of the MT

Full size image
Fig. 5
figure 5

a is the digital elevation model of the area; b is the mean peak ground displacement difference (∆PGDs) from all 25 simulations, showing an overall amplification on ridges (red) and de-amplification (blue) in valleys. c Standard deviation of ∆PGDs, indicating a higher variation in (de-)amplification on ridges when compared to valleys

Full size image

4 Results

In total, 25 simulations are carried out (Fig. 2), in which the MT solution of the 2005 Kashmir earthquake moves from its actual position to imaginary positions along major faults in the region (Fig. 1). In general, the PGD amplitude decreases with increasing distance from the MT solution (Fig. 2), except for places where topographic (de-)amplification occurs. It can be observed that, by comparing Fig. 2 with Fig. 1, in each model output the spatial pattern of the PGD amplitudes is influenced by the local topography. From Fig. 2, it can be seen that not only the topography itself, but also the location of the source relative to topographic features has an influence for the observed variation in amplitude: The direction of incoming waves with respect to the orientation of the slopes has an impact on the occurrence of possible (de)amplification. A ridge with high values of PGDs gets lower, or vice versa, when the location of source is changed.

The mean PGD of all 25 simulations (Fig. 3a) shows a radial pattern with high values (approximately 320 cm) in the central part of the study area, while reaching lower values (approximately 100 cm) at the edges. A relation between the topography (Fig. 3d) and the mean PGD (Fig. 3a) can be observed visually. In general, higher mean PGD values can be observed for elevated areas (ridges) compared to low-lying surrounding areas (valleys).

Depending on the direction of wave propagation, the exposure to incoming waves and distance from the MT locations, the PGD values vary between the 25 simulations. This is visible in the standard deviation of the 25 PGD simulations, which shows a variability between approximately 25 cm and 200 cm.

The central part of the study area has overall higher values. The deepest parts of valleys have relatively low values, but where the topography is higher, also the standard deviations increases. The largest standard deviations are located along the Muzaffarabad fault in the elevated areas to the northeast (Fig. 3b).

A seismic hazard map based on the mean PGD (Fig. 3a) and corresponding standard deviation (Fig. 3b) is derived in Fig. 3c. This HSV map shows the mean PGD in green to red hues: Areas with a high mean PGD are considered to be high hazardous areas, and vice versa. The standard deviation of the PGD is shown in saturation: Areas with low standard deviation appear as (intense) saturated tones and areas with high standard deviation appear as (whitish) unsaturated tones. This is, for example, clearly visible in the region N of Muzaffarabad which has high average PGD with high standard deviation, compared to an area east of Muzaffarabad that is more prominently visible with a high PGD and low standard deviation.

Although topographic effects are visible in Figs. 2 and 3, it is difficult to quantify the impact of topography on the PGD. To explore the topographic effect, the topographic amplification effect is obtained from models with and without topography (Fig. 4). The resulting pattern of amplification (positives values) and de-amplification (negative values) appears to follow the topography and is also influenced by the changing location of the MT. The (de-)amplification pattern is also found to have dependency on the angle of incident of waves with the slopes (Fig. 3). Depending on the location of the MT, the same slope can either be amplified or de-amplified. For example, the eastern region has a de-amplified response for the original MT location (Fig. 4-0), but gets an amplified response when the location of MT is changed (e.g., Fig. 4-8, 13, 23). On average, the topography increases the modeled PGD with 125 cm.

When all simulations shown in Fig. 4 are averaged (Fig. 5b), a correlation with topography (Fig. 5a) becomes visible: Ridges are subjected to an overall amplified response while valleys show a de-amplified response. The standard deviation (Fig. 5c) shows that the difference in amplification depends on topography and ranges from 5 cm in valleys and flat areas to 125 cm on some of the ridges. When slopes surrounding Muzaffarabad are facing away from a simulated source, they show an amplified seismic response. Similarly, when slopes are facing toward a MT location, they are in the shadow zone of seismic energy and show a de-amplified seismic response. For example, the slopes facing W-SW around N 34°30′ −E 73°30′ have an amplified response when the source is located near the faults in NE of the area (Fig. 4-12, 13, 14, 15, 16). The standard deviation (Fig. 5c) also shows higher values on same slopes, indicating high variability in amplification due to changing location of the earthquake source.

5 Discussion

Seismic hazard analysis is carried out by estimating ground motion for hypothetical earthquakes in the area of Muzaffarabad, Pakistan, with the MT solution of the 2005 Kashmir earthquake. This source is moved to 24 different locations along major faults in the region, at approximately 6 km intervals. The rough terrain of the study area has an impact on seismic response in the area, as observed in Fig. 2. During the 2005 Kashmir earthquake (Fig. 2-0), ridges experienced higher amplitudes of PGD compared to valleys (Khan et al. 2020). The area encompasses multiple major faults that may cause an earthquake in the future, and such an event could occur at any location along these faults. Our results predict higher amplitudes on ridges and lower amplitudes in valleys. It also shows that areas closer to the earthquake source will shake more compared to farther area, unless changed by topographic (de-)amplification. The central part of the areas has relatively high PGD values due to the presence of most MT solutions nearby. The values become lower toward the boundaries of the area because fewer MT solutions influence the result. An exception, however, is a topography-related amplification effect in the northeastern part of the area, which also has relatively more MT solutions nearby. The high mountains in the north are facing away from most of the MT locations and therefore do not trap energy. As a result, these mountains show relatively little amplification and have a lower standard deviation than the ridges surrounding Muzaffarabad City. A large amount of variation is expected on the ridge tops while comparatively low variation is expected from the valleys, as evident from standard deviation of PGDs in Fig. 3b. Most likely reason for this is that slopes in one direction have low response because of being in shadow zone, but for the opposite direction shows maximum amplification. It is why the deepest part of valleys have relatively low standard deviation values compared to elevated areas. The pattern of amplitudes is thus a combined effect of the topography and source location.

We also investigate how much the area is subjected to topographic (de-)amplification, compared to the hypothetical case when there would be no topographic variation (i.e., plain surface). The amount of (de-) amplification, beside the slope, geometry and height of topographic features, also depends on its location with respect to the seismic source (MT). It was found that the PGD of models with topography can have an approximate factor 1.7 increase in magnitude, when compared to plain surface models. The spatial extent of areas with a high PGD may get reduced by topography controlled destructive interference or a seismic shadow zone. Areas near the MT should have high amplitudes of PGDs, but due to topography these amplitudes can be de-amplified (e.g., the bluish areas in Fig. 4). The ridges, on the other hand, have an amplified seismic response (e.g., the reddish regions in Fig. 4). This amplified response can be observed even when a ridge is far from the MT. In principle, the amplification factor is distance independent and will only vary due to directional variations. The mean amplification in terms of difference (∆PGD) ranges shows the same phenomena of amplification (reddish) on the ridges and de-amplification (bluish) in the valleys (Fig. 5b). Like the standard deviation for PGD, standard deviation for ∆PGD also shows high variation on the ridges, with low variation in the valleys, due to changing MT location (Fig. 5c). Key to this analysis is that ridges show higher absolute values of PGD and higher amplification values. But there is a strong dependency on the location of the earthquake and the incoming direction of the seismic energy with respect to the orientation of the topography. This results in increased standard deviations for the average PDG and ∆PGDs for ridges in comparison with values for valleys or flat areas. This makes predicting the effect of earthquakes on slopes much more uncertain but also makes very clear that topography can lead to very strongly enhanced seismic impact.

Our seismic hazard map (Fig. 3c) integrates the mean PGD and standard deviation values in a single map. This map summarizes the topographic and potential source location effect on seismic-induced ground shaking in the study area. Its scale is relative, meaning that reddish tones indicate more hazardous areas and green tones indicate safer (but not per definition safe) areas. In our analysis, we defined that an area with high PGD and low standard deviation is more hazardous than an area with high PGD but with high standard deviation. Despite the fact that such a location might experience a higher absolute PGD in any of the 25 simulations, the overall effect on the area is lower since it also experiences much lower amplitude events. The validity of this reasoning can be seen in the areas directly east and north of Muzaffarabad. The area east has similar maximum amplitudes as the area to the north, but much lower standard deviation. There are therefore many more events that will lead to high values in this region compared to the area to the north (despite having similar or even larger maximum amplitudes). One could argue that an area with high PGD and high standard deviation is more hazardous compared to high PGD but with low standard deviation, but based on statistics the chance of a destructive event is larger for the area to the east than to the north. This way we also add information on the certainty that an area might be hit by a high impact event, compared to the average PGD map alone.

We have applied in our modeling always the same CMT solution, suggesting similar faulting mechanism and orientation on all faults. Although the faulting mechanism in the area is predominately thrust, the orientation is clearly more variable. This might have an impact on the energy radiation pattern of the each of the modeled earthquakes. We have considered 4 possible scenarios and evaluated the impact of each choice. The first option was to use the moment tensors of earthquakes on each respective fault. However, these events are relatively small in amplitude and would therefore not represent the worst-case scenario required for the hazard map. The second possibility could be the use of a synthetic moment tensor. However, this would then not be realistic since it will not be clear if that pattern of energy radiation is representative for the area. A third possibility would be the use of same moment tensor with its strike orientation parallel to the respective fault in the area (rotated from the Kashmir earthquake CMT). However, such modification would still be unrealistic since then the radiation pattern will also be rotated and there is no way to check if that is more realistic than any other energy radiation pattern. The fourth, and last, option is the use of the 2005 Kashmir earthquake moment tensor and use it as is along the faults in the region. However, in such a case we would have the same radiation pattern for all earthquakes, which is likely to be unrealistic. We have finally opted to go for scenario 4. The main purpose of this study was seismic hazard analysis by developing a methodology that adopts a deterministic approach for producing a probabilistic seismic hazard map. For this purpose, we want to use a worst-case scenario earthquake. It should represent the maximum recent historic magnitude for the region with a source mechanism that is known to occur. We accept that the actual radiation pattern of the energy could be somewhat different for each of the modeled scenarios. Since for other fault the radiation pattern would be an unknown anyway, we believe that using the CMT of the 2005 Kashmir earthquake is a realistic simplification to obtain a worst-case scenario seismic hazard map for the region.

In the Third United Nations (UN) World Conference on Disaster Risk Reduction, held in Sendai, Japan, from 14 to 18 March 2015, the UN member states adopted the “Sendai Framework for Disaster Risk Reduction 2015–2030.” The United Nations Office for Disaster Risk Reduction UNDRR has been tasked to support the implementation, and the follow-up activities and review of the Framework (UNDRR 2015). The Sendai Framework sets four priorities areas for action: understanding disaster risk, strengthening disaster risk governance to manage disaster risk, investing in disaster risk reduction for resilience and enhancing disaster preparedness for effective response, and to “Build Back Better” in recovery, rehabilitation and reconstruction (UNDRR 2015). The seismic hazard map and the methodology adopted in this study could be directly used in understanding disaster risk by identifying vulnerable areas and level of exposure. In areas of high topographic relief, there is an increased probability of co-seismic landsliding. The methodology adopted and results obtained from the scenario earthquake modeling and the corresponding derived seismic hazard map, jointly with other causative factors like anthropogenic factors and geology, can be used in identifying potential landslides zones. The information then can be used for risk reduction, preparedness, and taking mitigation measure to avoid and/or minimize the potential economic, social and environmental consequences of earthquakes that may occur in the future. Similarly, it can also be helpful for engineers and town planners in designing and planning infrastructure (such as bridges, buildings, town planning etc.) in the area so that it may be able to sustain in case of any future seismic event.

6 Conclusion

This study simulated the seismic response of possible future earthquakes in the area surrounding the city of Muzaffarabad (Pakistan). Totally, 25 simulations were carried out, in which the MT solution of the 2005 Kashmir earthquake was moved from its original and actual position to imaginary positions along major faults in the region. Generally, the PGD amplitude decreases with increasing distance from the MT solution, but is superimposed by a clearly visible topographic amplification effect. Not only the topography itself, but also the relative location of the source with respect to a slope has an influence on the observed variation in amplitude. The mean PGD of all 25 simulations varies between approximately 100 cm at the borders of the study area and approximately 320 cm for the elevated parts in the center of the study area. Depending on the direction of propagation, the exposure to incoming waves, and the average distance from a MT location, the PGD values vary between simulations. This becomes visible in the standard deviation of the PGD simulations, which shows a large variability from approximately 25 cm up to almost 200 cm.

We presented a seismic hazard map, based on an integration of the mean PGD and its standard deviation, summarizing the topographic and potential source location effect on seismic-induced ground shaking in the study area. It provides a classification from hazardous to safe in relative terms (so not per definition safe in absolute sense). Combining the hazard of multiple future scenarios, this map indicates potential seismic hazard and is a guide in earthquake preparedness.