Seismic hazard assessment of the Shillong Plateau using a probabilistic approach

Seismotectonic processes of the Shillong Plateau (SP) have been influenced by the Himalayan orogeny, the India-Burmese subduction, and the Bengal Basin evolution leading to high seismic activity in the region. With a goal of assessing seismic hazards in the SP and providing a scientific information to engineering and disaster risk management communities, a probabilistic seismic hazard analysis is employed to determine hazard in highly-populated districts of SP and particularly in Shillong, Nongpoh, and Tura cities, located within the districts. This analysis is based on the use of historical and instrumentally recorded regional earthquakes since 1411 and deals with uncertainties related to earthquake magnitudes, rupture locations, and the frequency of ground motion exceedance. Individual hazard curves indicate that the Barapani fault possesses the highest frequency of seismic hazard for Shillong city and Nongpoh, and the Eocene hinge zone and Dauki faults are responsible for the highest frequency of seismic hazard at Tura. The results of the hazard assessment together with those obtained earlier using a scenario-based approach demonstrate that although the Oldham fault located near Tura can produce a great, but rare earthquake, few other nearby faults are capable of producing smaller magnitude events with a higher probability of occurrence. ARTICLE HISTORY Received 29 May 2019 Accepted 2 October 2020


Introduction
The northeastern (NE) India is tectonically active region intertwined with numerous faults, which have been the source of several strong to great damaging earthquakes. The regional seismicity has encouraged scientists to understand the seismic hazard potential of the region. As earthquakes are associated with uncertainties related to their location, size and ground motion exceedance, probabilistic seismic hazard (PSH) analysis could assist in addressing these uncertainties and provide an information on the exceedance of a prescribed value of ground shaking due to an earthquake within a certain duration of time and with a defined probability of the exceedance. Sharma and Malik (2006), Anbazhagan et al. (2009), NDMA (2010), Pallav et al. (2012), and Sil et al. (2013) used PSH approach to assess seismic hazard for either the entire region or its different parts. In addition, Parvez et al. (2017) performed neo-deterministic seismic hazard assessment of India by computing synthetic seismograms with input data set consisting of structural models, seismogenic zones, focal mechanisms and earthquake catalogues.  performed deterministic (scenario-based) seismic hazard (DSH) assessment for the Shillong Plateau (SP) located in the southwestern part of NE-India and encompassing the state of Meghalaya (Figure 1(a)). They found that the Barapani, the Oldham, and the Dauki faults were the most-seismic hazard causing sources with estimated peak horizontal acceleration (PHA) of up to 0.47 g within the SP, which was higher than that recommended by the Indian code of seismic design of structure (IS 1893(IS 2016. This DSH analysis considered the worst scenarios for each seismic source based on its maximum possible earthquake magnitude, and did not account for the likelihood of the scenarios to occur, e.g., during the design life of constructed building and other infrastructure in the region. The present work uses the probabilistic approach to seismic hazard analysis of the SP to overcome the limitation associated with DSH analysis. The principal goal of this study is to determine seismic hazards in three populated districts of the SP, providing a useful and usable scientific information to engineering applications and disaster risk management communities. Therefore, the objectives of the study are; (i) to assess earthquake hazards in the region using PSH analysis, (ii) to analyze the results obtained by probabilistic (in this study) and scenario-based assessments , and (iii) to provide a broader knowledge on seismic hazards in the region accounting for various scenarios and uncertainties in the assessments. This would permit decision-makers to consider seismic safety strategies based on the scientific knowledge on ground shaking and potential losses due to regional earthquakes, on the cost-benefit analysis, and planning for mitigation strategy.

Seismotectonics
The SP is surrounded by faults, which have generated several major (7.0 Mw 7.9) and one great (Mw ! 8.0) earthquakes for the recorded period of time. The SP is considered to have been formed as a result of tectonic movements between the Indian, Eurasian and Burma plates, approximately 60 My ago (Johnson and Alam 1991). The movements of the tectonic plates led to devastating earthquakes, and three human settlements, such as Shillong city, Nongpoh, and Tura located within the SP had suffered severe damages during these earthquakes (Oldham 1882;Bilham 2008; Baro and Kumar 2017a).
The northeastern region of the Indian plate converges with the Eurasian plate at the rate of 5 cm y À1 , and the SP inside the main Indian plate rotates in a clockwise direction at the rate of 1.15 per My (Vernant et al. 2014) leading to strain-stress localization. Seismic sources in the region within the 500-km radius from the center of the SP (Figure 1(b)) were identified by Baro and Kumar (2017b). We note that only known seismic sources, which generated seismicity from 1411 until 2015 are included in the analysis. Based on seismicity and tectonic patterns, regional geology, rupture characteristics, and the rate of the regional movement, this seismotectonic region was split into four seismic source zones: the Shillong Plateau -Assam Valley Zone (SP-AVZ), the Indo-Burma Ranges Zone (IBRZ), the Bengal Basin Zone (BBZ), and the Eastern Himalaya Zone (EHZ) (Figure 1(b,c); Baro and Kumar 2017b). We discuss here briefly the tectonic setting of the seismic source zones following .
The IBRZ is located to east of the SP-AVZ and comprises of the north-south trending faults in the vicinity of the tectonic boundary between the Indian and Burma plates. The faults in this zone originated as a result of the collision and subduction process between the two plates. According to Wang et al. (2014), the Kabaw fault is a strike-slip fault having a capacity to generate earthquake of Mw8.4, and the Churachandpur-Mao Fault (CMF) located west of the Kabaw fault could generate an earthquake of Mw7.6. To the west of the CMF, the Naga thrust fault is located, which continues as the Kaladan fault towards the south. The Kaladan fault shows a strikeslip mechanism northward and a thrust mechanism southward (Maurin and Rangin 2009).
Southwestward, the IBRZ changes from mountainous terrain to the plains of the Bengal Basin, filled-up with sediments dating around 60 to 2.6 Mya (Mohanty and Walling 2008). Two most prominent tectonic sources within this seismic source zone are the Sylhet fault and the Eocene Hinge Zone. Nandy (2001) suggested that the Sylhet fault was likely to be the source for the 1918 Srimangal earthquake (m b 7.6). According to Nath et al. (2014) and Kayal (2008), the Eocene Hinge Zone located at a depth of 4.5 km beneath Kolkata City, might have been the source for the 1935 Pabna earthquake (Mw6.2). Other major faults that exist within this zone are the Garhmayna-Khanda Ghosh fault, Jangipur-Gaibandha fault, Pingla fault and Sainthia-Bahmani fault.
The EHZ is located to the north of the SP-AVZ and comprises of the northeastern part of the Himalayan mountain range. The Himalayas originated approximately 60 Mya when the Indian plate started subducting under the Eurasian plate (Johnson and Alam 1991). The subduction led to the formation of thrust faults within the zone. Some of these faults are the Indus Suture Thrust (IST), the Main Central Thrust (MCT), and the Main Boundary Thrust (MBT). Some studies such as Mukhopadhyay (2011) suggest that several past earthquakes could be traced back to the Main Himalayan Thrust, located between the underthrusting Indian plate and the overlying Himalayas. However another school of thought exists which suggests that the MCT is presently in a dormant state (Ni and Barazangi 1984;Kayal 2014). Nevertheless, it is unequivocal that the northeastern section of the MCT and MBT have remained dormant, leading to the speculation of the existence of the Assam Gap (Khattri 1987), which is also called as the Northeast Seismic Gap (Mittal et al. 2012). It is also speculated that an earthquake of Mw8.5 might occur in this region of the seismic gap (Srivastava et al. 2015).
The diversity among the four seismic source zones provides a difficulty to estimate the future seismic hazard potential of the SP (Baro and Kumar 2017b). For example, faults situated in the BBZ and the IBRZ produced major earthquakes in the past; meanwhile, strong earthquakes occurred in the IBRZ, but was not yet observed in the BBZ (Figure 1(b)). Hence, to account seismic contributions from the both seismic source zones, the probability of occurrence of different size (magnitude) earthquakes should be considered. The faults located within the studied seismic source zones are of varying lengths (Figure 1(c), following nomenclature of faults as , and different size earthquakes can occur at different locations within the same fault. Thus, it is essential to consider the probability of rupture at any point within a fault, i.e., at different hypocentral distances from the site. The presence of major earthquakes within these zones also raises the question of probability of occurrence of such earthquakes during the lifetime of a structure. The above mentioned uncertainties associated with determining the seismic hazard potential of the SP could be addressed by PSH assessment.

Hazard analysis
We use the regional declustered earthquake catalogue (Baro and Kumar 2017b), where foreshocks and aftershocks of main earthquake events have been removed using the static window method of Reasenberg (1985). This catalogue containing 2359 earthquakes for the period from 1411 to 2015 and is split into four sub-catalogues corresponding to four seismic source zones introduced in Sect. 2. Baro and Kumar (2017b) estimated the mean annual activity rate, aand b-values of the Gutenberg-Richter (GR) relationship, and the maximum earthquake magnitude for each seismic source zone using a relevant sub-catalogue and the Kijko et al. (2016) method. In addition,  estimated the maximum possible magnitude (M p ) at each of 72 faults within the seismotectonic region (four seismic zones; see Table A1 in Appendix).
A PSH anlysis considers uncertainties associated with earthquake magnitude, hypocentral distance, and ground motion exceedance at a given site due to a particular fault rupture (Reiter 1989). Hereinafter, we follow the seismic hazard assessment approach described in EM-11101110 (1999). This approach allows for assessing the frequency of ground motion exceedance at a given site during a certain time interval. Using the Poisson distribution model, the probability of exceedance P E (z) of a ground motion level z (or a spectral acceleration) at a site of interest for an exposure time period "t" is given the following equation: where v z is the annual frequency of ground motion exceedance at the site of interest, which can be determined considering the uncertainties as where K is the total number of seismic sources in the seismotectonic region (K ¼ 72 in this study); k k ðm i Þ is the frequency of earthquake magnitude m i to occur on each source k (i ¼ 7 values of magnitudes considered in this study); PðR ¼ r j jm i Þ is the probability that during an earthquake of magnitude m i the rupture occurs on source k at the hypocentral distance r j (j ¼ 14 sub-faults on each source fault considered in this study) from the site of interest; and PðZ>zjm i r j Þ is the probability that an earthquake of magnitude m i occurring at the hypocentral distance r j will cause ground motion Z exceeding a known threshold value of z. Below, we present how the three uncertainties, associated with magnitude, distance and ground motion exceedance, can be assessed.

Frequency of different magnitudes (m i ) of EQs on each faultk k ðm i Þ
For each seismic source zone, the frequency of occurrence of earthquakes (N) of magnitude mi and greater can be presented in the truncated exponential form of the GR (frequency-magnitude) relationship (EM-1110 1999) where M is the magnitude; b is the b-value in the GR relationship; m 0 is the lowest earthquake magnitude of interest to the calculation (m 0 ¼ 4 in this study); and aðm 0 Þ is the frequency of earthquake occurrences with magnitude m 0 and greater, that is, the ratio of the total number of earthquakes of magnitude m 0 ! 4 in the catalogue to the duration of the catalogue's completeness, which is 80 years. For SP-AVZ, IBR, BBZ, and EHZ, the frequency of occurrence aðm 0 Þ is estimated to be 3.28, 17.01, 2.375, and 6.81, respectively.
To assess the frequency of occurrence of earthquakes on individual fault within each seismic zone (N k ), we use the following relationship proposed by Iyenger and Ghosh (2004): where r k ¼ length of fault k summation of all fault lengths within the seismic sourse zone , and d k ¼ total number of EQs on fault k summation of the total number of EQs in the seismic sourse zone : Table A1 lists 72 faults along with their characteristics, such as their starting and end point coordinates, maximum observed and possible earthquake magnitudes, the number of earthquakes on each faults (including the total number of earthquakes per seismic source zone), the fault's length (including the total length of faults per seismic source zone), and parameters r k and d k for each fault.
The interval between magnitudes m 0 and M p is divided into I-1 parts (in this study, I ¼ 7). For each point of the division (i ¼ 2, … , I À 1), the values ofN k ðM>m i þ Dm=2Þ and N k ðM>m i ÀDm=2Þ are estimated considering Dm ¼ ðM p À m 0 Þ=ð2nÞ and noting that N k ðM>m 1 ÀDm=2Þ ¼ N k ðM>m 1 Þ and N k ðM>m I þ Dm=2Þ ¼ N k ðM>m I Þ: At each of the points and for all 72 faults, the frequency of occurrence k k ðm i Þ of earthquake magnitude m i on fault k can be then calculated as A typical calculation of the frequency of occurrence k k ðm i Þ for the Oldham fault is presented in Table A2 (in Appendix), and a plot of the k k ðm i Þ versus m i is presented in Figure 2. As expected, the frequency of earthquake occurrences decreases with the increase in magnitudes and reaches lowest value of 0.000019 for magnitude 8.7.

Uncertainty in hypocentral distance
The conditional probability distribution function of the hypocentral distance R given that the magnitude M ¼ m i for a rupture segment uniformly distributed along a fault can be estimated from the following equations (Der Kiureghian and Ang 1977; Iyenger and Ghosh 2004): where r is the shortest distance (in km) between the site and a rupture segment; L is the length of the fault (in km); X(m i ) is the rupture length (in km) corresponding to earthquake magnitude m i ; and notations d, D, L 0 , and r 0 are explained graphically in Figure 3. X m i ð Þ is calculated corresponding to m i for each of the 72 faults, using the following equation (Iyenger and Ghosh 2004): Here, the minimum of two terms in Eq. (9) defines the length of the rupture. To calculate the conditional probability distribution function (Eqs. (6) -(8)) with respect to distance between the rupture segment and the site, each of the 72 faults is divided into 14 sub-faults. For each sub-fault, (i) the distance between the site and the location of the sub-fault center þ(the sub-fault length)/2 on the fault, and (ii) the distance between the site and the location of the sub-fault center -(the sub-fault length)/2 on the fault is considered. The difference between the conditional probabilities at the two locations is taken as the probability of rupture to occur at distance from the sub-fault center. Each of the three districts (east Khasi Hill, West Garo Hill, and Ri-Bhoi) is then divided into grids of 0.05 x 0.05 and the distance between the top left corner of each grid to each sub-fault's center is used to calculate the probability of rupture. Typical calculations considering the Oldham fault with m i ¼ 4.0 at the Shillong city center is summarized in Table A3 (in  Appendix), and the results of the calculations are presented in Figure 4. The cumulative probability of rupture varies between the hypocentral distances of 77.87 km to 170.54 km as 0.038 and 1.0, respectively.

Uncertainty in ground motion exceedance
Uncertainties are considered in PSH analysis by estimating the probability of exceedance of the ground motion parameter above a threshold value for a given m i and r j . According to EM-1110EM- (1999, the conditional probability distribution for exceedance of the ground motion PðZ > zjm i r j Þ is estimated using the following log-normal distribution whereFðUÞ is the cumulative distribution function of a truncated normal distribution; U ¼ ð ln Z À E½ln zÞ=S½ln z; Z is the ground motion of interest to this study; z is the ground motion calculated from the selected ground motion prediction equation (GMPE) based on m i and r j ; E represents the mean level of ground motion given by the selected GMPE; and S represents the standard error in the ground motion prediction by the selected GMPE (presented below). The functionFðUÞ in Eq. (10) can be represented as where F(U) is the standard cumulative normal distribution. In this way, for various combinations of m i and r j , the probability that ground motions assessed by GMPEs exceed the initially defined range of Z is calculated using Eq. (10). For the present work, the values of Z vary from 0.025 g to 0.80 g in increments of 0.025 g considering all 72 faults. Once the value of PðZ > zjm i r j Þ is calculated, using formulas (5)-(10) the overall frequency of exceedance of a ground motion to occur can be determined.
In this study, we employ the following GMPEs : where, z is the spectral acceleration (in g), M is the moment magnitude, R M ¼ r þ 0:089 exp ð0:6MÞ, r is the hypocentral distance, e a is aleatory uncertainty, e e is epistemic uncertainty, and A 1 to A 6 are the constants (presented in table 3 by Toro et al. 1997 where g is the acceleration due to gravity, f 0 ¼ maxð ln ðr=100Þ, 0Þ, e is the standard error; and B 1 to B 8 are the constants specified for NE-India (NDMA 2010); and where b d is a decay parameter, C 1 , C 2 and C 3 are regression coefficients, and r e is the standard error .
The termSð ln zÞ in Eq. (10) is estimated asSð ln zÞ ¼ S 1 þ S 2 þ S 3 , where S i is the product of the weight of the GMPE and the value of the standard error term in Toro et al. (2002) for S 1 , in NDMA (2010) for S 2 , and in Anbazhagan et al. (2013) for S 3 .
The weights of GMPEs for each seismic source zone are listed in Table A4 (in Appendix).
GMPEs assist in an assessment of the level of the ground motion at a given site due to an earthquake. In this study, the ground motion parameters are calculated at the sites corresponding to shear wave velocity of 1500 m s À1 or to class A (hard rock) as per the NEHRP site classification scheme (BSSC 2004). Because measured data on the ground motion are limited, a standard error term is introduced in most GMPEs to account for variation in ground motions. When several GMPEs are employed in a hazard assessment, a logic tree approach helps to minimize the standard error by using weighting (ranking) coefficients prescribed to each GMPE (Delavaud et al. 2009). Note that as the highest maximum possible magnitudes differ within each source zone, the ranks of the same GMPEs vary with the zone. Moreover, the GMPE's weights can be influenced by the difference in the hypocentral distances from the selected sites -Shillong, Nongpoh, and Turato each of the seismic source zones .

Hazard curves
A hazard curve presents the frequency of exceedance of a specific level of ground motions defined over a range of the motions. As discussed in Sect. 3, the values of the frequency of exceedance versus corresponding ground motion due to a fault rupture is calculated at each site using Eqs. (5)-(10). Repeating this procedure for all faults considered, individual hazard curves for the faults can be determined at each site. Figure 5 illustrates the hazard curves for various contributing faults at Shillong city (a), Nongpoh (b), and Tura (c). All hazard curves show that the stronger a ground motion (PHA) is, the lower is the frequency of occurrence of the that motion in a specified time period. The hazard curves are presented not only for the sources close to the sites considered, but also for those sources capable producing significant ground motion at the sites. The Barapani fault possesses the highest frequency of exceedance for any possible ground motion at Shillong city followed by the Dauki fault, the Kopili fault, the Oldham fault, the Eocene Hinge Zone and other faults ( Figure 5(a)). Also the Barapani fault is the highest hazard-causing fault for Nongpoh ( Figure 5(b)). The Dauki and Dhubri faults show the highest frequency of exceedance until the peak ground acceleration is lower than 0.225 g, and hence, both the faults generate the highest seismic hazard at Tura (Figure 5(c)). At greater PHA values, the Eocene Hinge Zone shows a higher frequency of exceedance.
An individual fault-based hazard curve accounts for ground motion at the site of interest from all possible earthquake magnitudes to occur at all possible location on that fault. Hence, a PSH assessment provides an aggregated assessment of all possible scenarios of ground motions with estimated frequency of exceedance of the motions. A deaggregation procedure allows for identifying the earthquake scenarios at a source fault (for instance, magnitude-distance pairs), which have contributed to the hazard curve at the site of interest. While all the contributions are accounted to get a hazard curve, the deaggregation procedure determines relative contribution to the seismic hazard from the range of magnitude (M), distance (R), and the ground motion scatter or the random error, which is defined as the number of standard deviations (often referred to as 'epsilon'), by which the ground motion deviates from its median value. The highest contributing combination can help in understanding the selected return period or the probability of ground motions being exceeded (McGuire 1995;Bazzurro and Cornell 1999;Kumar et al. 2013;Sokolov and Ismail-Zadeh 2015). Meanwhile, it should be noted that if small magnitude earthquakes occur much more frequently than bigger magnitude events near the site of interest, the number of small events may be large enough to cause many exceedances of selected ground motion level due to highlevel positive deviations from the mean values of predicted ground motion. Thus, the relative contribution to hazard from small events may be comparable with the contribution from much larger, but less frequent events (Sokolov and Ismail-Zadeh 2015).
As the number of standard deviations has already taken into consideration during the event-by-event calculation of seismic hazard, as described in Sect. 3, we perform here deaggregation in terms of the magnitude and the source-to-site (hypocentral) distance considering the full range of error. The deaggregation procedure for the case of the 10% probability of exceedance in 50 years (similarly for the case of the 2% probability of exceedance) assumes that the value of seismic hazard (suppose z1) is obtained for the fault causing the highest hazard (which can be seen in Figure  5(a-c)). Depending upon the considered magnitude and hypocentral distance bins (in Figure 6) for this fault (sections 3. 1.1 and 3.1.2), the frequency of exceedance is determined for each combination of the magnitude and hypocentral distance. If the contributions from all magnitude-distance combinations are added, it will give a total seismic hazard of z1, which will have the 10% probability of exceedance in 50 years. These frequencies of exceedance in each magnitude-hypocentral distance bin are expressed in terms of percentage contributions to the seismic hazard value. As a result of the deaggregation, a 3-D histogram is plotted with the magnitude, hypocentral distance, and percentage contribution to the seismic hazard as its axes. For this work, the deaggregation histograms are developed using the procedure described above for the worst seismic hazard causing faults for Shillong, Nongpoh, and Tura (as per Figure 5(a-c)), and the results are presented in Figure 6(a-d). The Barapani fault provides the highest contribution of 13.7% to seismic hazard at Shillong city, if an earthquake of Mw6.1 occurs at 24.67 km hypocentral distance (Figure 6(a)). At Nongpoh, the highest contribution from the Barapani fault is 13.9% due to an earthquake of Mw6.1, which might occur at 10.18 km hypocentral distance ( Figure 6(b)). At Tura, the highest contributions of 41.05% from the Dauki fault is possible due to an earthquake of Mw7.6 occurring at 69.76 km hypocentral distance (Figure 6(c)). Considering the Dhubri fault (another high hazard-causing fault for Tura), a contribution of 26.3% is possible at Tura due to an event of Mw7.6 occur at 80.41 km hypocentral distance (Figure 6(d)). The deaggregation procedure helps to identify the size and the location of a rupture area within a hazard-causing fault, which will produce a highest PHA among all possible seismic scenarios.
A summation of hazard curves for a site, from all faults considered, gives the total hazard curve indicating the overall frequency of exceedance of a particular ground motion. At the period of seismic wave oscillation 0 s, the frequency of exceedance is 0.00174 for the spectral acceleration 0.2 g, which gives a return period of 575 years for Shillong city (Figure 7(a)). Thus, the annual probability of exceedance of 0.2 g at Shillong city for the next 50 years is 8.33% according to Eq. (1). At oscillation periods of 0.4 s, 1 s and 2 s, the frequency of exceedance are estimated to be 5.69 Â 10 À4 , 1.97 Â 10 À4 and 9.28 Â 10 À7 , respectively. The probability of exceedance in 50 years at oscillation periods 0.4 s, 1 s, and 2 s are 2.93%, 0.98% and 0.0046%, respectively. Similarly, the frequency of exceedance of the spectral accelerations can be assessed for Nongpoh (Figure 7(b)) and Tura (Figure 7(c)).

Hazard assessment maps
Probabilistic hazard maps can be plotted by presenting seismic hazard values (e.g., PHA) at selected grid points in the study area, for a specified probability of exceedance and for a desired period. In this section, the seismic hazard maps are developed for the districts of the East Khasi hills, Ri-Bhoi, and the West Garo hills because of highly populated cities located in the districts, although the same technique can be applied to study seismic hazards of the entire Meghalaya state. PHA values are calculated for each of the three districts using 247, 273, and 323 grid points, respectively (laid out in a grid size of 0.05 x 0.05 ) and the numerical code by Kumar et al. (2013). Figure 8 shows the seismic hazard maps for the selected districts at 2% probability of ground motion exceedance in 50 years. The highest PHA value of 0.54 g and 0.48 g reaches in the northern part Khasi hills and the central Ri-Bhoi (Figure 8(a,b)). This indicates that these parts of the districts are seismically more hazardous than the remaining parts. The southern part of the West Garo exhibits a PHA value of 0.22 g (Figure 8(c)). PHA values at 2% probability of exceedance in 50 years (with a return period of 2475 years) in Shillong city, Nongpoh, and Tura are 0.28 g, 0.52 g, and 0.20 g, respectively. Figure 9 shows the hazard maps for the same three districts at 10% probability of ground motion exceedance in 50 years (with a return period of 475 years). Similarly, the northern East Khasi hills, the eastern Ri-Bhoi, and the southern West Garo hills exhibit highest seismic hazard, and the maximum PHA values at Shillong city, Nongpoh and Tura are estimated to be 0.19 g, 0.32 g, and 0.11 g, respectively. Bhatia et al. (1999), Das et al. (2006), and NDMA (2010) developed PSH assessment maps for the SP for 10% probability of exceedance in 50 years, and showed that the PHA values over the SP vary between 0.25 g and 0.30 g, and hence in good agreement with our assessment results. However, Bhatia et al. (1999) and Das et al. (2006) did not present seismic hazard values in the study area. Similarly, NDMA (2010) assessment considered the NE-India as one seismic source zone.

Spectral acceleration
While PSH maps provide probable seismic scenarios at specific period of oscillation, uniform hazard spectra (UHS) present seismic hazards at various oscillation periods for a specific frequency of exceedance ( Figure 10) and provide a design response spectrum for structural analysis. The spectral acceleration (SA) values at Shillong city at period 0 s are 0.28 g for 2% and 0.19 g for 10% probabilities of exceedance, and the peak SA values for 2% and 10% probabilities of exceedance are attained at period 0.1 s as 0.50 g and 0.33 g, respectively (Figure 10(a)). The SA values for Nongpoh at period 0 s are 0.52 g and 0.32 g for 2% and 10% probabilities of exceedance, respectively (Figure 10(b)), and for Tura are 0.20 g and 0.11 g at the same probabilities of exceedance (Figure 10(c)).
In addition to UHS obtained here for Shillong, Nongpoh and Tura, response spectra presented in IS 1893 (2016) is also plotted in Figure 10 (see also Table 1). IS 1893 (2016) presents the SA values for (i) a maximum credible earthquake (MCE), which is considered to be the most severe earthquake to occur in a region with the ground motion that has 2% probability of exceedance in 50 years, and (ii) for a design basis   GEOMATICS, NATURAL HAZARDS AND RISK earthquake (DBE), which can be expected to occur during the construction design life with the ground motion that has 10% probability of exceedance in 50 years. Figure 10 shows that the SA values based on IS 1893 (2016) are higher than those derived in this work, and the reason for that is site conditions considered in this work and in IS 1893 (2016). The latter defines hard soil/rock as a site condition having the standard penetration test N-value above 30, which belongs to site class B/C based on the NEHRP site classification scheme (BSSC 2004). On the other hand, UHS values derived in the present study corresponds to the NEHRP site class A condition. Thus, the site conditions for the present findings and those given in IS 1893 (2016) correspond to two different site conditions, and hence a direct comparison of the obtained results for site class A with the building code IS 1893 (2016) data is inappropriate as the values of the spectral acceleration will be much higher in the case of hard soil/rock sites compared to hard rock sites. The findings obtained in this work can be used for building design only after a refinement of the results based on an assessment of the local soil's contribution to the ground shaking. In case of a building's location on rocks related to site class A, the present findings can be used directly.
Further, we compare the spectral acceleration obtained in this work with those obtained earlier by NDMA (2010) and  for engineering bedrock conditions corresponding to site class A and three locations -Shillong City, Nongpoh, and Tura (see Table 1). Although the present study and study by NDMA (2010) use a probabilistic SHA, there are some differences in their assessment methodologies. For example, NDMA (2010) carried out PSHA for India, and the Shillong Plateau was a part of one of thirty-two source zones. NDMA (2010) developed PGA contour maps for the entire country, where only average values within a comparatively small area, like Shillong Plateau, can be extracted. As a result, the same spectral accelerations are listed in Table 1 for the three locations. Another source of the discrepancy in the results can come from the fact that a source zone is defined differently in both cases, e.g., in the case of the present study, the seismic activity, slip rate, and tectonic setting are considered. Moreover, more than one GMPE are employed in our analysis of ground motions.
Despite the differences in the methodologies, our results fit rather well the previous estimates (NDMA 2010), especially when the average values of the spectral accelerations obtained in this work for the three locations (a sum of the spectral accelerations determined for a specific period for Shillong City, Nongpoh, and Tura divided by three) are compared to the values obtained in NDMA (2010) (see the numbers in parentheses next to the NDMA's numbers in Table 1). The present assessment refines the NDMA (2010)'s results at the regional/local scales highlighting the importance of regional seismic hazard assessments for engineering aims. The values of spectral accelerations based on DSHA (Table 1;  tend to be closer to those values obtained by PSHA in the case of 2% probability of exceedance in 50 years than in the case of 10% probability, because average return periods are normally longer for large earthquakes.

Discussion: scenario-based versus probabilistic seismic hazard assessments
Compared to deterministic seismic hazard maps, probabilistic maps do not present the ground shaking at site, but instead present a level of ground shaking, which can be exceeded with a certain probability, within a certain period. The PSH maps provide a low bound of seismic hazard useful for engineering purposes (Ismail-Zadeh 2018). Keeping this in mind, we compare here results of the PSH assessment performed in this work with the results of deterministic (scenario-based) seismic hazard (DSH) assessment of the SP obtained by . We note that both DSH and PSH analyses are based on the same earthquake catalogue, same b-values and maximum possible magnitude values, and the same GMPEs. However, the pattern of seismic hazard for the West Garo hills district significantly differs in the case of DSH and PSH assessments. Namely, the highest PHA values for the scenario-based hazard assessment were calculated for the north-eastern part of the district, and the values decrease to its south-western part (Figure 11(a)). Meanwhile, the pattern of seismic hazard values based on 10% (Figure 11(b)) and 2% (Figure 11(c)) probability of ground motion exceedance in 50 year shows the highest ground motion to be associated with the southern part of the district. The Oldham fault, located to the north of the SP and generated an earthquake of Mw8.0 (1897 Shillong earthquake), is very likely to influence the pattern of seismic hazard in the district. In the DSH assessment, its contribution is significant because of the possible maximum earthquake magnitude, the fault had produced. Meanwhile, the great earthquakes at the fault have lower frequency of reoccurrence, and hence the contribution of the Oldham fault to the probabilistic assessment of seismic hazard is smaller compared to other faults, such as the Eocene hinge zone, Dauki fault and Dhubri fault (Figure 5(c)).
To show the influence of the Oldham fault to the results of the DSH assessment, the same DSH calculations (as in  has been performed by removing the Oldham fault from the analysis. In this case, the DSH and PSH assessment results for the West Garo hills district show a similar trend of the highest seismic hazard towards the south with gradually decreasing seismic hazard towards the north (Figure 11(d)). Hence, we conclude that the effect of the fault generating great earthquakes is prominent in DSH assessments. However, the probability of rupture of the Oldham fault to generate a great earthquake is smaller (but not negligible) during the design life of a structure.
Consideration of the maximum possible earthquake magnitude at a fault leads to the worst-case scenario in seismic hazard quantification, and this might generate difficulties in decision making related to building construction/retrofitting and disaster risk management. Besides an economic component associated with construction/retrofitting costs, the determination of the maximum possible earthquake magnitude has significant uncertainties (e.g., Kijko 2004;Kijko and Singh 2011). There is no widely accepted procedure for the determination of the maximum earthquake magnitude in the region. But even though the maximum possible magnitude could be determined exactly, a recurrence time of events with this magnitude, or slightly below it, is unknown. It is because of a lack of observations on large earthquakes occurred on natural faults. Note that some attempts were made to estimate statistically the recurrence times of large earthquakes on individual faults (e.g., Z€ oller et al. 2007).
The term of "recurrence time" is perceived differently in scientific and non-scientific communities dealing with seismic hazards and risks. Earthquake recurrence time is referred to an average time between two subsequent big events (of the same magnitude) normally evaluated based on several known occurrences of big events. Therefore, a recurrence time should be used cautiously in hazard assessment, because a great earthquake may happen in a region within a much shorter (compared to a recurrence time determined) period since the last big event. For example, Chile experienced two great earthquakes of Mw9.5 in 1960 (38.24 S;73.05 W) and of Mw8.8 in 2010 (36.29 S;73.24 W) with an epicentral distance of 217 km between two events and within time interval 50-years although such events are extremely rare. The major 2003 (Mw7.6) and 2013 (Mw7.8) Scotia Sea earthquakes occurred at the same fault segment (Vall ee and Satriano 2014). The time difference between the 1812 Wrightwood (estimated Mw7.5) and 1857 Fort Tejon (estimated Mw7.9) earthquakes in California was much shorter than the estimated average recurrence time of about 330 years (Jacoby et al. 1988). The 193486.59 E) and the 2015 Nepal Mw7.8 earthquake (28.23 N;84.73 E) occurred at a spatial distance of about 238 km and temporal distance of 81 years, although based on a comprehensive analysis of geological, paleo-seismological, geophysical and geodetic observations since Late Holocene, the return times of great earthquakes in eastern Nepal were estimated to range between 750 ± 140 and 870 ± 350 years (Bollinger et al. 2014).
Similarly, the recurrence time of great earthquakes in the Shillong Plateau region cannot be known with the accuracy, to communicate to the earthquake engineering community or disaster risk management authorities. However, the information about ground shaking due to a great earthquake in the region, such as the 1897 Assam earthquake, is essential for design/retrofitting of important building constructions in the region. Meanwhile, a determination of all possible seismic scenarios expected during the design life of infrastructure is vitally important for a seismic safety assessment. Such hazard assessments should include various uncertainties related to determination of magnitudes of the events, their recurrence times, ground motion due to site conditions etc. PSH analysis allows to overcome some of the issues related to uncertainties, although without scenarios of large events, the analysis will be not comprehensive (e.g., Sokolov and Ismail-Zadeh 2015).
When speaking about seismic safety as the goal of hazard assessments, an important issue is an analysis of physical and social vulnerability, exposure, and hence risk assessment (as a convolution of hazard, vulnerability and exposure) allowing to elaborate strategic countermeasure plans for the disaster risk mitigation. Seismic risk assessments facilitate a proper choice of safety measures, ranging from building codes and insurance to establishment of rescue-and-relief resources (Ismail-Zadeh 2018). Most of the practical problems require estimating risk for a specific territory (e.g., populated cities or districts) and within this territory separately for the objects such as lifelines, sites of vulnerable constructions, etc. Multi-scale hazard and risk assessments are helpful for decision-makers at local, municipal, and national levels. Therefore, the results of the PSH analysis presented here together with the results on scenario-based hazard analysis  for three most-populated districts of the Meghalaya state is of important for decision-making.

Conclusion
In this study, the PSH analysis has been performed for the districts of the East Khasi hills, Ri-Bhoi, and the West Garo hills, located within the SP. Hazard curves are developed at the sites of Shillong city, Nongpoh, and Tura located within the East Khasi hills, Ri-Bhoi and the West Garo hills districts, respectively. Hazard curves show that the Barapani fault possesses the highest frequency of exceedance at Shillong city and Nongpoh. At Tura, both Dauki and Dhubri faults are found to be the highest hazard causing faults.
The northern part of the East Khasi hills district, eastern part of Ri-Bhoi district and southern part of the West Garo hills district show high PHA value. This suggests that these parts of the districts are prone to higher seismic hazards. Within these three districts, the PHA values range from 0.22 (0.11) g to 0.46 (0.33) g at 2% (10%) probability of exceedance in 50 years. A comparison of our results with those developed by Bhatia et al. (1999), Das et al. (2006), and NDMA (2010) show a reasonable agreement. Also, uniform hazard spectra for Shillong city, Nongpoh and Tura obtained in this study are found comparable with those in NDMA (2010). A comparison between the results of DSH and PSH assessments for the same region performed by  and in this work shows that a contribution of the Oldham fault to seismic hazard of the region is prominent in the case of DSH assessments, and much smaller in case of PSH assessments because of the recurrence time of large events at the Oldham fault.   (1) Sub-fault no.; (2) hypocentral distance with respect to sub-fault center (r); (3) X(m i ) (km); (4) L o (km); (5) r o (km); (6) hypocentral distance with respect to subfaults length / 2 þ subfault center location; (7) P(B < r); (8) hypocentral distance with respect to sub-faults length / 2sub-fault center location; (9) P(C > r); (10) probability of rupture at sub-fault center.