Probabilistic seismic hazard assessment of Nepal using multiple seismic source models

The potential for devastating earthquakes in the Himalayan orogeny has long been recognized. The 2015 MW7.8 Gorkha, Nepal earthquake has heightened the likelihood that major earthquakes will occur along this orogenic belt in the future. Reliable seismic hazard assessment is a critical element in development of policy for seismic hazard mitigation and risk reduction. In this study, we conduct probabilistic seismic hazard assessment using three different seismogenic source models (smoothed gridded, linear, and areal sources) based on the complicated tectonics of the study area. Two sets of ground motion prediction equations are combined in a standard logic tree by taking into account the epistemic uncertainties in hazard estimation. Long‐term slip rates and paleoseismic records are also incorporated in the linear source model. Peak ground acceleration and spectral acceleration at 0.2 s and 1.0 s for 2% and 10% probabilities of exceedance in 50 years are estimated. The resulting maps show significant spatial variation in seismic hazard levels. The region of the Lesser Himalaya is found to have high seismic hazard potential. Along the Main Himalayan Thrust from east to west beneath the Main Central Thrust, large earthquakes have occurred regularly in history; hazard values in this region are found to be higher than those shown on existing hazard maps. In essence, the combination of long span earthquake catalogs and multiple seismogenic source models gives improved seismic hazard constraints in Nepal.


Introduction
Due to the ongoing collision between the Indian and Eurasian plates since ~65 Ma (Ding L et al., 2005), the Himalayan mountain range has grown into the highest orogenic belt on Earth, experiencing widespread seismicity. Situated in the central segment of the Himalayan orogenic belt, Nepal is one of the most earthquake-prone countries in the world. The classical tectonic divisions of this region from south to north are the Sub-Himalaya, Lesser Himalaya, Higher Himalaya, and Tethyan Himalaya (Yin A, 2006), which are separated by four major E-W trending faults: the Main Frontal Thrust (MFT); the Main Boundary Thrust (MBT); the Main Central Thrust (MCT); and the South Tibetan Detachment System (STDS) (Figure 1). These surface faults have merged with the Main Himalayan Thrust (MHT), the detachment that separates the underthrusting Indian plate from the overriding Himalayan orogeny (Nábělek et al., 2009). The convergence rate between the Indian and Eurasian plates is 30-40 mm/year towards the NNE (Stevens and Avouac, 2015). About half of this rate is accumulating along the MHT (Bettinelli et al., 2006).
The Himalayan orogen has experienced many devastating earthquakes (Table 1), including the 2015 M W 7.8 Gorkha, Nepal earthquake. Geodetic measurements (e.g., Stevens and Avouac, 2015) have identified that the MHT is fully locked at the shallower portion in its interseismic period and the accumulation process of the seismic strain continues until being released by large earthquakes Rajendran et al., 2015;Sapkota et al., 2013). Previous studies have predicted that an earthquake of M W 8.0 at 100 km east to the Kathmandu valley could potentially harm 1.3% of the population and damage 50% of the buildings (Wyss, 2005). The 2015 M W 7.8 Gorkha earthquake caused about 9000 casualties and destroyed over fifty thousands of buildings with huge economic losses (Bilham, 2015). However, the Gorkha Earthquake was considered to be smaller than the great earthquakes anticipated in the central Nepal (Bilham, 2015;Hand and Pulla, 2015) and eventually heightened the likelihood of catastrophic future earthquakes along the Himalayan orogenic belt.
Probabilistic seismic hazard assessment (PSHA) has long been employed in determining appropriate building design code provisions by many earthquake-resistant infrastructure projects (e.g., IBC, 2006). The implementation of seismic design provisions for buildings and other infrastructure is an effective way to reduce seismic risk and consequently to avoid earthquake-related disasters. After the pioneer work of Cornell (1968), PSHA has been widely studied on global, regional, and local scales (e.g., Chaulagain et al., 2015;Nath and Thingbaijam, 2012;Ordaz et al., 2014;Ram and Wang GX, 2013;Sawires et al., 2016;Zhang PZ et al., 1999). Earthquake catalogs (both instrumental and historical) are the primarily basis for PSHA studies. However, the complete part of the instrumental earthquake catalog is very short (e.g., around 55-60 years for this study region) compared to the time required for earthquake generation by tectonic processes (Bilham, 2013). This catalog limitation in PSHA causes low hazard value at some regions where great earhquakes have actu-  (Kijko et al., 2016). This limitation is evident from some cases of big earthqukes, such as the 2008 Sichuan M W 7.9, the 2010 Haiti M W 7.0, and 2011 Japan M W 9.0 earthquakes (Wang ZM et al., 2012). In this study therefore we accumulate all available historical and instrumental earthquake catalogs to assess the probabilistic seismic hazard for this area.
Seismic hazard maps for Nepal have been produced since the 1990's (e.g., Chaulagain et al., 2015;Rahman et al., 2018;Ram and Wang GX, 2013;Zhang ZM et al., 1999). Most of these previous studies estimated the hazard value using only the complete part of the catalogs (about 60 years of data), based on single seismic source models and old attenuation equations. A recent study of the GEM Faulted Earth Project characterized the attribute of the Himalayan Frontal Thrust (MFT) fault in the seismic hazard assessment for the frontal Himalaya (Berryman et al., 2014). In this study, we use a methodology developed by Kijko and Smith (2012) to compute seismicity parameters (e.g., Gutenberg-Richter b value, rate of seismicity). In general, PSHA procedures assume that seismicity parameters are constant over time. But study of earthquake catalogs has shown that these parameters may change over time (Kijko et al., 2016). Thus we calculate the b value and seismicity rate by taking into account the incompleteness of the earthquake catalog since the 12 th century. Moreover, three seismogenic source models (i.e. the smoothed gridded seismicity, the linear source model, and the areal source model) and two sets of ground motion prediction equations (GMPEs) are combined using a standard logic tree structure for capturing model-related epistemic uncertainties. The topographic elevation in the areal source and smoothed gridded seismicity models are also taken into account, while paleoseismic and geodetic records are incorporated into the linear source model. A set of seismic hazard maps for peak ground acceleration (PGA) and spectral acceleration (SA) are produced in this analysis for a better understanding of the spatial variation of seismic hazard.

Earthquake Catalog
Earthquake catalog provides the basic data for delineating seismogenic sources and computing the seismicity parameter, especially the mean seismic activity rate λ, the Gutenberg-Richter b value, and the maximum expected earthquake magnitude M max . We consider the incompleteness of the earthquake catalog which includes both historical earthquakes and the instrumental catalog, with different levels of completeness (Kijko and Smit, 2012;Kijko et al., 2016). Ambraseys and Douglas (2004) and Szeliga et al. (2010) provide information regarding historical major earthquakes for Himalaya and adjacent regions. We also use the histor-ical earthquake catalogs from the National Oceanic and Atmospheric Administration (NOAA), GEM Global Historical Earthquake Archive (GEM-GHEA), and the National Seismological Centre of Nepal (NSC). Instrumental earthquake catalogs between 1906 and 2016 are collected from the U.S. Geological Survey (USGS) and International Seismological Centre (ISC).
The catalogs from different sources are then compiled together and checked manually to remove duplicate events on the basis of location, time, and magnitude. We convert the different types of magnitude (e.g., body wave magnitude, surface wave magnitude, duration magnitude, and local magnitude) contained in the catalog into M W using the different empirical relationships shown in the Table 2. Approximately 2100 events with M W ≥4.0 are recorded in the comprehensive catalog. The seismicity model usually assumes that the occurrence of an earthquake is independent. All foreshocks and aftershocks are therefore excluded from the catalog, deploying a declustering process following the algorithm of Gardner and Knopoff (1974). A total of 900 events (including 13 historical earthquakes with magnitude equal to or greater than 6.0) are included in the final catalog ( Figure 1).
Assessment of different levels of completeness of various subcatalogs for the instrumental earthquake catalog is very important. The reason is that records of large earthquakes are generally complete for longer periods while the catalogs of smaller earthquake are complete for shorter durations. The completeness of each earthquake catalog is governed by such considerations of socioeconomic and historical circumstances, demographic variations, and seismic station coverage. Failure to acknowledge the incompleteness of the catalogs would result in overestimate and underestimate of recurrence rates for small and large earthquakes, respectively. We first assume a threshold magnitude (M 0 ) of M W 4.0 and then estimate completeness by following the visual cumulative method (Tinti and Mulargia, 1985). The modern instrumental catalog is found to be complete since 1964 (Figure 2a). In utilizing the incomplete instrumental earthquake catalog since 1906, the entire catalog is divided into different sub-catalogs based on different levels of completeness. The complete time window for various magnitude intervals (e.g., 4.0-5.0; 5.0-6.0; 6.0-7.0) is assessed using the statistical method of Stepp (1972). The complete time windows are found for different levels of completeness of 4.0, 5.0 and 6.0 to be after 1964, 1925-1964, and 1906-1925, respectively (Figure 2b).

Earthquake Source Models
Proper seismic source zoning is essential in PSHA. Seismogenic sources with different characteristics are generally delineated according to tectonic settings, seismological and geological attrib-  (Kaviris et al., 2008) utes, and fault characteristics. The Himalayan orogen is regarded as a very complex earthquake producing environment. We thus use a model that combines three different seismic sources: areal, linear, and smoothed gridded sources.

Areal Seismic Sources
In general, the area of a seismic source zone is delineated using the basic principle that seismicity within a single source zone is uniform and homogenous. It is therefore assumed that every point within the source has an equal potential to generate a future earthquake (Baker, 2015). Recent studies (Chaulagain et al., 2015;Ram and Wang GX, 2013) have delineated 23 areal seismic source zones ( Figure 1) for the study region based on quantitative analysis of the earthquake distribution, fault information, and seismotectonic settings. In this study, the same (23) seismic areal source zones are utilized. High seismicity clusters are present in western (zones 16, 17 and 18) and east central (zones 9, 10, 11 and 12) Nepal. The frontal edge of the Himalaya along the boundary of southern Nepal shows sparse seismicity.

Linear Seismic Sources
Usually, major earthquakes are associated with active faults, and in seismic hazard assessment these active faults are considered as linear seismic sources. Although, there are numerous methods for earthquake source simulation (e.g., Zhang WB et al., 2006), for the Himalayan-Tibetan region, Styron et al. (2010) identified active faults using remote sensing imagery, aerial photographic interpretation, paleoseismic studies, and field investigations. For our study region, 14 active faults are assigned as linear seismic sources ( Figure 1). Apart from the MFT, there are no well-constrained attributes for the other faults. Therefore, two magnitudearea and magnitude-length scaling empirical relationships of Yen YT and Ma KF (2011) and Wells and Coppersmith (1994) are applied to characterize these faults' attributes. The compilation and evaluation of approximately 72 such kinds of relationships were carried out in the GEM Faulted Earth and Regionalisation Components project (Stirling and Goded, 2012). They found that the relationship of Yen YT and Ma KF (2011) is better in quality because of its updated wide-spread datasets with greater magnitude coverage than that of Wells and Coppersmith (1994) who used older datasets and relatively narrow magnitude coverage.
They recommended the quality score of 1 (best available) and 2 (good) for Yen YT and Ma KF (2011) and Wells and Coppersmith (1994), respectively. Accordingly, we assign the weight for Yen YT and Ma KF (2011) and Wells and Coppersmith (1994) as two-thirds and one-third, respectively. For large magnitude earthquakes, the occurrence rate along these faults is determined using long-term slip rates as: where μ is shear modulus 3.0 × 10 11 dyne/cm 2 , L is rupture length, W is rupture width, is the fault slip rate, and M 0 is the seismic moment. The relationship between the seismic moment and magnitude (M) is that of Hanks and Kanamori (1979): Small scale fault attributes (e.g., dip, dip direction, slip, and slip rate) are not well documented for this study region. We therefore assign to small-scale faults the same values for the above mentioned parameters as those of nearby known faults (MFT). The Main Frontal Thrust (MFT) is one of the well-constrained active faults in this study area (e.g., Ader et al., 2012). The GEM faulted earth project (Berryman et al., 2014) characterized the different attributes of the MFT for all of the three different segments in the study area. In this study, we used the same attributes (e.g., rupture length, width, dip, dip direction, slip, slip rate, and expected maximum magnitude) for the three segments of MFT (i.e. MFT-1, MFT-2, and MFT-3, respectively from west to east).
We deployed two magnitude-frequency distribution models in the linear source model: the Characteristic-earthquake model for the empirically estimated earthquakes with magnitude equal or higher than M W ≥6.5, and the Gutenberg-Richter model for instrumentally recorded earthquakes. The computed linear source parameters (e.g., seismic moment, slip rates, rupture length, rupture) for active faults are illustrated in Table 3 for the Himalayan region.

Smoothed Gridded Sources
As earthquake hypocenters are widely scattered in relation to fault locations, we use a zone-free approach for the spatial smoothing of seismicity (Frankel, 1995;Woo, 1996) to avoid subjectivity in delineating areal seismic sources. This approach is termed "kernel estimation"; it implements spatial smoothing of seismicity and depicts the seismicity's non-uniform spatial distribution, by assuming that seismicity varies from place-to-place. For smoothed gridded seismicity, thousands of 0.1º × 0.1º grid cells are made for the entire study area. At each grid cell, earthquakes of M W ≥4.0 are counted for different magnitude intervals in estimating the seismicity rates for corresponding magnitude intervals. The final seismic activity rate is then obtained by smoothing the calculated seismicity rates using the multivariate kernel probability density function following the equation of Woo (1996) as: where n j is the rate of seismicity in the j-th grid cell, is the smoothed rate of seismicity in the i-th cell, c is the correlation distance accounting for location uncertainties, and Δ ij is the distance between the i-th and j-th cells. The sum is taken over cells j within a distance of 3c of cell i. It is important to note that the resulting smoothed seismicity will be concentrated around the epicenters of the recorded earthquakes for very small c values, while larger c values spread out the seismicity and eventually do not reflect the true spatial variation of seismicity as assumed in the zone-free approach. Therefore, in this study, the c value of 0.5° is assigned, following Frankel (1995).

Maximum Earthquake Magnitude and Focal Depth
The maximum expected earthquake magnitude (M max ) for every seismic source zone is computed on the basis of historical (e.g., 1505 M W 8.2 and 1934 M W 8.2 earthquakes) and recent earthquake data, but also considering great earthquakes that occurred in nearby seismic zones with similar tectonic settings. The estimate of M max is computed following the algorithm of Kijko and Singh (2011). This estimation takes into account the uncertainty in the earthquake magnitude determination. The estimated maximum expected magnitude of an earthquake for each areal source zone, along with its standard deviation, is given in Table 4. The computed M max values vary between 6.5 and 8.6. For smoothed gridded seismicity, the same range of M max is used for the respective grids.
Accurate determination of earthquake location along the mountain belt is usually difficult because of insufficient station coverage. The ISC catalog reveals that earthquakes are distributed throughout the thickness of the Tibetan and Indian crust ranging between 0 and 90 km. Moreover, some recent studies have found that the hypocenters of most of the earthquakes in Nepal are located at depths ranging between 0 and 30 km (Pandey et al., 1999;Priestley et al., 2008). In addition, hypocenter relocation results from a nearby broadband seismic array found that the focal depths for the Gorkha earthquake sequence ranged primarily between 10 km and 20 km (Adhikari et al., 2015;Bai et al., 2016;Letort et al., 2016). Thus, in this analysis, the average focal depth of 15 km is assigned as a constant depth for all the source models. It is essential to note that the topographic elevation of the study area varies significantly from south to north. The elevation of the source zones along the south frontal border is less than 1 km and in the central and northern zones is about 3 km and 5 km, respectively. Seismic intensity attenuates reasonably with respect to distance from the source and thereby the topography has a considerable effect on hazard value. Thus, the topographic elevation is considered for each source zone along with the assigned average focal depth.

Seismicity Model and Parameters
To estimate the earthquake magnitude exceedance rates, a Modified Gutenberg-Richter (MGR)-Poisson model is chosen as a seismicity model for all seismogenetic sources except characteristic earthquakes. Herein, the seismicity is expressed as (Cornell and where λ 0 is the rate of earthquake activity with threshold magnitude M 0 (here M 0 = 4.0), β is a parameter equivalent to the b value (β = 2.303×b), and M max is the maximum expected magnitude for the sources. In computing the λ(M) value, the uncertainty of both M max and β are taken into account.
Conventionally, seismicity parameters are calculated based on the minimum magnitude of completeness (M C ) of the catalog. The Maximum Curvature (MAXC) method of Wiemer and Wyss (2000) is used to evaluate the M C value. The estimated M C value of M W 4.0 is taken as the M 0 because earthquakes smaller than M W 4.0 are not likely to cause any damage to infrastructure, from an engineering point of view.
In this study, the essential seismicity parameters are computed by taking into account the incompleteness of the instrumental catalog and the historical earthquake records. The Gutenberg-Richter b value and its standard error for every source zone are determined following the joint maximum likelihood function of Kijko and Smit (2012) based on the historical and instrumental earthquake catalog. This new approach is the extension of the Aki-Utsu (Aki, 1965;Utsu, 1965) b value estimator which takes into account both incompleteness of the earthquake catalogs and temporal variation of seismicity. It provides straightforward approximations for the standard errors and confidence intervals (Kijko and Smit, 2012). In this method, the whole span of the earthquake catalog is divided into different sub-catalogs and the parameter, β (β = 2.303×b) is calculated as: β where m j i is the sample of n i earthquake magnitudes recorded during the time span of the i-th sub-catalog (i=1, 2, …, s). In addition, maximization of the likelihood function in equation (5) provides the generalized Aki-Utsu -value estimator as: However, in this analysis, the b value and its standard errors are estimated for the whole region, as well as for each areal source zone. For the data of the whole region, the estimated b value is 0.82 and the standard error is 0.12. These values are assigned to zones 4, 5, 6, 8, and 21 because of low seismicity. On the other hand, for rest of the zones, the individual b value and its standard error are assigned for hazard computation (Table 4). In case of smoothed gridded sources, the b value and its standard error of 0.82 and 0.12 are employed to all of the grids.β Once the -value is known, then the average annual rate of seismicity (λ 0 ) at threshold magnitude for each seismic source zone is computed using the joint maximum likelihood function (Kijko and Smit, 2012) as:

Selection of Ground Motion Prediction Equations
The ground motion prediction equations (GMPEs), also known as attenuation relationships, are essential in estimating the ground motion parameters (e.g., PGA, SA). These equations usually give ground motion intensity as a function of many variables such as earthquake magnitude, source-to-site distance, and soil condition. In general, the attenuation relationships are obtained empirically using a statistical regression method on a particular set of strong motion data.
The seismotectonic setting of the Himalayan orogenic belt is very complex. This study area is a tectonically active interplate region (Nath and Thingbaijam, 2012) and the earthquake occurrence pattern of this region is not simple and even rather compound in combination. For example, Nath and Thingbaijam (2012) identified that this study region belongs to the active shallow crustal regime (ACR) with many subduction zone interface (SZI) earthquakes across the Himalaya. Earthquake focal depth and mechanism study (Bai et al., 2017) stated that the occurrence of large earthquakes is fairly associated with subduction interface in this region. Moreover, it was assumed by some previous seismic hazard studies that the Himalaya and its surrounding region belongs to either SIZ or ACR (e.g., Ram and Wang GX, 2013), or both (e.g., Chaulagain et al., 2015). In avoiding the controversy, we consider that the study region behaves as both of these environments. However, the strong ground motion data for the Tibetan-Himalayan region are insufficient and consequently no particular GMPEs have been studied previously.
Usually, this problem is resolved by selecting the GMPEs developed in other regions of similar seismotectonic settings, which could represent sufficiently the ground motion in this region. Therefore, a total of 8 GMPEs are selected in this analysis according to the criteria (e.g., tectonic regime, functional form, regression coefficients, frequency range of the model, journal types, datasets, etc.) of Cotton et al. (2006), and a number of trellis plots (i.e. predicted spectral acceleration (PSA) versus magnitude, PSA versus distance, PSA versus structural periods) as Stewart et al. (2013). Two Next Generation Attenuation (NGA-west2) models: Abrahamson et al. (2014) and Chiou and Youngs (2014) from the global datasets; one model from Europe and Middle East datasets of Ambraseys et al. (2005); and other one from mainly Japanese datasets of Zhao JX et al. (2006) are selected for ACR. The weight assessment and thereby the ranking of different GMPEs are not carried out because of scarcity of strong motion records in the Himalayan-Tibetan region. Due to insufficient observed data, many seismic hazard studies for Asian countries (e.g., Chaulagain et al., 2015;Kolathayar and Sitharam, 2012;Nath and Thingbaijam, 2012;Ornthammarath et al., 2011) used equal weight for all models. In this study, therefore, for each of these four models, equal probabilistic weights (1/4) are assigned. Besides, three models from global datasets covering the large distance, namely BC Hydro 2016 known as Abrahamson et al. (2016), Atkinson and Boore (2003), and Youngs et al., (1997); and one for local dataset of Zhao JX et al. (2006) are selected for SIZ. Equal probabilities (1/4) of weight to each of the models are assigned, as for ACR.
Compatibility among these selected GMPEs is checked when applied in the CRISIS2015 (Ordaz et al., 2015) program. Finally, these two sets of GMPEs are combined together using a logic tree. In the similar way, GMPEs are employed for all three source models.

Logic Tree Structure
Uncertainties are commonly associated with different models used in PSHA. The application of a logic tree allows the capturing of epistemic uncertainties in different input models Sabetta et al., 2005) by deploying alternative models in the hazard estimation. In this study, three types of seismogenic source models and two sets of GMPEs (each consisting of four different GMPEs) are combined using the logic tree. Equal weights (1/3) are assigned to each of the source models, as there is no substantial observed data for the study region in evaluating the ranking of the weight. For the linear source model, equal weight (1/2) is assigned for the Characteristic Earthquake and Gutenberg-Richter seismicity models as there is no reason to prefer one model over the other. Figure 3 displays the standard logic tree framework used for the linear source model in this analysis.

Seismic Hazard Estimation
To estimate the seismic hazards on the basis of computed seismi-city parameters, the seismic hazard module CRISIS2015 (Ordaz et al., 2015) is used. This module has high efficiency in hazard calculation and flexibility in model selection (Danciu et al., 2010). For areal and linear seismic sources, the principle of CRISIS in computing the seismic hazard is based on the standard Cornell approach (Cornell, 1968). The assumption of the Module is that within a source zone, seismicity is evenly distributed by unit area and thereby all the points might be a potential earthquake focus. A spatial integration by subdividing the original sources is performed by the CRISIS program in obtaining accurate hazard value. Once a source is subdivided into sub-sources, the acceleration exceedance rate due to the i-th single source can be computed using the equation: where M 0 and M max are the threshold and maximum magnitudes, respectively; is the probability that the acceleration exceeds a certain value a at the site across distance R ij for an earthquake with magnitude M; and W ij is the weight of each subelement. The above expression also assumes that ∑W ij = 1. The total average acceleration exceedance rate at the site due to the contributions of all the sources, N, within 300 km is calculated as: For smoothed gridded seismicity, once the seismicity parameters and the GMPEs are known for each of the nodes of the grids, the hazard at a given node is calculated based on the effects of the totality of the nodes and their corresponding distances to the site of interest as in equations (8) and (9).

Results
For seismic hazard estimation, the entire study region is divided into small grids with the size of 0.1º×0.1º and the total number of grid cells is approximately 5000. At the center of each grid cell, the PGA and SA values are computed for a referenced bedrock condition. The contribution of all the potential sources to hazard value within the radial distance of 300 km from the center of each grid cell is taken into account. We estimate the hazard value by disaggregating the hypocentral distance into small intervals of 1 km, and the magnitude range (between the minimum and maximum magnitude) into small incremental values of 0.1. Initially, the haz-ard value at each grid cell is computed using all of the eight different attenuation models for all three seismogenic sources separately. Afterward, the final hazard value at each grid cell is obtained by the combination of these individual values, using the logic tree structure. This analysis produces different seismic hazard maps for PGA value and SA value at short (0.2 s) and long (1.0 s) natural periods at 5% critical damping factor for 10% and 2% probabilities of exceedance in 50 years, corresponding to 475 and 2,475 year return periods, respectively. PGA hazard maps for the probabilities of exceedance of 10% and 2% in 50 years are shown in Figure 4a, b. Figure 5a, b illustrates the spatial variations of SA value obtained at 0.2 s and 1.0 s natural periods for 2% probability of exceedance.
The overall spatial distribution of hazard value regarding both PGA and SA is non-uniform. Because of heterogeneous character- istics of the seismotectonic settings at the study region, there are clear lateral variations in the computed seismic hazard maps. The predicted hazard distribution pattern largely follows the background seismicity distribution across the study area. The regions of highest potential hazard are likely distributed along strike of the Indian-Eurasian plate collision zones. However, the hazard estimate for 10% probability of exceedance shows a wide range of values, and the minimum and maximum hazard values are found as low as 0.21g, and as high as 0.64g, respectively (Figure 4a). The zones of highest hazard are identified in the far western and central to northeast central Nepal (Figure 4a) where the rate of seismicity is very high compared to other parts of Nepal (Figure 1). For far western Nepal (i.e. south of Darchula city), we obtain peak hazard values as high as 0.64g. In central Nepal around the Kathmandu valley, peak hazard values range up to 0.59g. The south-ern-most part along the Nepal-India boundary is identified as a low seismic hazard zone. MHT has long been predicted as a potential source for producing great earthquakes in the Himalayan region Rajendran et al., 2015;Sapkota et al., 2013). Large earthquakes in the central portion of the Himalayan orogenic belt may greatly affect the densely-populated cities of Nepal, such as Kathmandu and Pokhara .
Nonetheless, the predicted seismic hazard value is increased significantly with decreasing probability of exceedance. For example, at 2% probability of exceedance in 50 years, the computed hazard level for central and northeast central Nepal is about a factor of 1.8 higher than that of 10% in 50 years; for far western Nepal, the comparable factor is nearly 1.5 (Figure 4b)  as high as 1.02g, respectively. We produce the SA maps (Figure 5a, b) as per the requirement of the modern building design provisions, which are analogous with the PGA maps in terms of hazard distribution pattern. SA value at 0.2 s and 1.0 s natural periods for 2% probability of exceedance in 50 years is found in the order of 0.87-2.38g, and 0.29-1.19g, respectively.
To understand how different source models contribute to the computed seismic hazard values, we estimate the hazard value for each model individually. The resulting hazard maps for these source models show that the levels of maximum and minimum hazard value are quite different from model to model. For instance, the linear source model produces maximum values as high as 0.60g whereas the maximum hazard values for the smoothed gridded and areal source models are 0.71g and 0.54g, respectively. We also find that the effect of sources on hazard value varies spatially. In regions with high hazard, all the source models (areal, smoothed gridded, and linear) produce nearly similar ground motion values ( Figure 6). However, for northern-most, southern-most, and the central west regions where background seismicity is low, different levels of effect are observed among the models. For example, in the southern and central west regions of Nepal (e.g., around the cities of Bara, and Rukum), the linear source model contributes comparatively higher seismic hazards than the other two source models, and consequently augments the combined hazard value levels for these regions (Figure 6). On the other hand, in northern Nepal (near Gorkha city), the contribution from areal and smoothed gridded sources is comparatively higher than that of the linear source model. In the southwest region (near Kailali city), the highest contribution is from the smoothed seismicity source, and the areal source contributes the least hazard value, while the linear source model yields an intermediate effect. Moreover, to infer the attenuation pattern of seismic hazard with respect to various natural periods, the uniform seismic hazard spectra at 10% probability of exceedance in 50 years for different cities is computed (Figure 7). The highest hazard value is obtained at the short natural period of 0.2 s for all cities. With the increase of the natural periods, the hazard value falls rapidly, and at the longer periods of about 1.5 s and above, the values became nearly flattened and exhibit very gentle slopes of attenuation.

Conclusions and Discussions
To produce seismic hazard maps of Nepal, we incorporate the latest datasets, multiple seismogenic source considerations, updated GMPEs, topographic elevation data, well constrained focal depth data, and consideration of inherent uncertainties associated with each of these inputs. An updated uniform declustered earthquake catalog since 12th century is utilized. 23 areal sources, 14 linear sources and hundreds of smoothed gridded sources are combined using a logic tree to account for the epistemic uncertainties. Essential seismicity parameters (e.g., b value and rate of seismicity, λ) are computed using the joint maximum likelihood method (Kijko and Smit, 2012). The incompleteness of the catalog (i.e. different complete sub-catalogs) and the uncertainty in earthquake magnitude determination are considered for seismicity parameters computation. The modified Gutenberg-Richter and Characteristic-earthquake magnitude-frequency models, two sets of GMPEs (four GMPEs for the active shallow crust and four GMPEs for the subduction interface) are utilized. In essence, the combination of past seismicity data (both historical and instrumental), paleoseismic and geodetic data, along with three sesimogenic source models, thus provides more comprehensive and explainable constraints to the computation of seismic hazard values.
The spatial variation pattern of the hazard value shows good agreement with the maps of Zhang PZ et al. (1999), Ram and Wang GX (2013), and Rahman et al. (2018). In addition, the resulting hazard values for the Nepal area are comparable with results for the adjacent regions in India (e.g., Kolathayar and Sitharam, 2012;Nath and Thingbaijam, 2012) and Tibet (Rahman et al., 2018). The rock-site PGA value around the epicenter of the 2015 Gorkha earthquake shown in previous studies (e.g., Chaulagain et al., 2015;Ram and Wang GX, 2013;Zhang PZ et al., 1999) is below 0.4g. The hazard value for the source area of the Gorkha earthquake has been estimated to be low (Chaulagain et al., 2015;Ram and Wang GX, 2013) based on the single source model and the short span earthquake catalog. However, the PGA values obtained in this study are up to 0.45g at the same probability of exceedance (i.e. 10%). We suggest that this difference is probably a result of our consideration of multiple seismogenic sources along with the wide span of our earthquake catalog including historic large earthquakes.
The hazard value is assessed in terms of PGA for 10% and 2% probabilities of exceedance on bedrock level. The spectral acceleration distribution maps at 0.2 s and 1.0 s natural periods are produced for 2% probability of exceedance in 50 years. The region of the Lesser Himalaya from east to west is found to be has high seismic hazard potential. The seismic hazards estimated in this study would help constrain the national building design provisions for the central segment of the Himalayan region. Better seismic hazard estimate in a tectonically active region requires detail paleoseismic, geological, and geodetics studies. Strong ground motion data is essential in order to develop the new attenuation relationships for this tectonically active region.