Seismic hazard assessment of the Shillong Plateau, India

Abstract The Shillong Plateau is an earthquake-prone region in the northeastern India. Based on regional seismotectonic studies, we present here a deterministic seismic hazard assessment (DSHA) and maps of peak horizontal accelerations (PHA) for three largely populated districts – the East Khasi hills, the Ri-Bhoi, and the West Garo hills – within the Shillong Plateau. The hazard analysis methodology is based on the analysis of 72 earthquake sources (active faults) located within 500 km seismotectonic region around the plateau. Using an average sample log-likelihood approach, suitable ground motion prediction equations (GMPEs) are identified. As a variation in hypocentral distances can affect the ranks (or weights) of selected GMPEs, DSHA is performed separately for the three selected districts. Analyses show that the northern part of the East Khasi hills, eastern part of Ri-Bhoi district and the West Garo hills districts exhibit the highest PHA value of 0.46 g at site class A (hard rocks). In addition, response spectra for the Shillong city, Nongpoh, and Tura indicate that the maximum spectral acceleration reaches 0.67 g, 0.77 g, and 0.64 g at 0.1 s, respectively. These assessments indicate that the Barapani, Oldham, and Dauki faults influence significantly the seismic hazard in the studied region.


Introduction
The northeastern region of India (NE-India) is located between two plate boundaries: the Indian-Burmese plate boundary from the east and the Indian-Eurasian plate boundary from the north. The lithospheric deformation along these plate boundaries is responsible for the origin of faults within NE-India. The Shillong Plateau (SP) formed due to fault movements in the southwestern portion of NE-India (Bilham and England 2001) spreads over approximately 25,000 km 2 and comprises of the Indian state of Meghalaya. Meghalaya is divided into eleven districts (Figure 1), among which the East Khasi hills district (where the capital city, Shillong is located) and the West Garo hills district are densely populated. Stress release along the faults led to several major to great earthquakes (EQs) in the past as shown in Figure 2. This makes the SP an active seismic region.
The 1869 Cachar EQ (marked as 13 in Figure 2) caused significant damages to public and governmental buildings in the Shillong city and Nongpoh, the headquarters of the Ri-Bhoi district. During the 1897 Assam EQ (marked as 18 in Figure 2),  many buildings were significantly damaged or destroyed in the Shillong city, and in Cherrapunji located in the East Khasi hills district. About 500 to 600 people were reported dead in Cherrapunji due to landslides caused by the EQ. Landslides had been triggered by the 1897 Shillong EQ in Nongpoh and in Tura located in the West Garo hills district. In Tura, apart from the landslide, the material from the failed slope got accumulated, blocking the water way and forming lakes; masonry houses on the other hand also suffered major damages (Oldham 1882;Bilham 2008; Baro and Kumar 2015). The 1930 Dhubri EQ caused shaking up to intensity VII (the Rossi-Forel scale) in Cherrapunji and the Shillong city and up to intensity VIII in Tura leading to severe damages of buildings and infrastructure. To reduce damages from future EQs, it is essential to quantify the seismic hazard potential of the SP and to transfer the knowledge to engineers and disaster management authorities to undertake actions. A deterministic (scenario-based) seismic hazard analysis (DSHA) is presented in this paper for three selected districts within SP -East Khasi hills, Ri-Bhoi, and West Garo hillswhere the most populated citiesthe Shillong city, Cherrapunji, Nongpoh, and Turaare located. DSHA follows the following steps: (i) to identify and characterize EQ sources (or tectonic faults) located around the site of interest; (ii) to determine the maximum possible EQ magnitude on each source; (iii) to select ground motion prediction equations (GMPEs), which mimic the observed ground acceleration in the best way; and (iv) to determine peak horizontal acceleration (PHA) for selected sites as well as spectral accelerations for different frequency (or period) of seismic waves. We present these research steps of the DSHA in the following sections.

Seismotectonics of the Shillong Plateau
Numerous faults located within and around the SP have been identified, and a seismotectonic map of the SP has been developed by Baro and Kumar (2017). This map ( Figure 2) covers major part of northeast India, portion of the eastern Himalayas, plains of Bengal and Bangladesh within a circle of 500 km radius with its center located at the geographical point (25.63 N, 91.25 E). The entire seismotectonic region is divided into four seismic source zones, which differ from each other in terms of rupture characteristics, tectonic features, thickness of overburden, geology, and the rate of plate movement: (1) the Shillong Plateau-Assam Valley Zone (SPAVZ), (2) the Indo-Burma Ranges Zone (IBRZ), (3) the Bengal Basin Zone (BBZ), and (4) the Eastern Himalayas Zone (EHZ). Baro and Kumar (2017) identified 72 tectonic faults within the seismotectonic region ( Figure 3; Table 1). The prominent tectonic features within each of the above four seismic source zones are summarized below, and illustrated in Figures 2 and 3.

Shillong Plateau -Assam Valley zone (SP-AVZ)
The SP located within the SP-AVZ is surrounded by numerous tectonic faults including the Dauki fault, the Dapsi thrust (marked as 11 in Figure 3), the Oldham (or Brahmaputra) fault, and the Kopili fault. The Dauki Fault is an EW-trending south dipping normal fault located on the southern boundary of the SP (Baro and Kumar 2015). The 90 $ 100 km long Dapsi thrust is located to the west of the Dauki fault and trends NW-SE (Kayal 2008). In addition, the Dhubri fault, which was the source of the 1930 Dhubri EQ, lies to the west of the SP. According to Bilham and England (2001), the Oldham fault, situated at the northern edge of the SP was responsible for the plateau's uplift 60 million years ago and was also the source for the 1897 Assam EQ. However, Rajendran et al. (2004) argues that the Brahmaputra fault located to the north of the SP, rather than the Oldham fault, is the source of the Assam EQ. The 300-400 km long and 50 km wide, northwest dipping Kopili fault is located to the northeast of the SP within the Assam Valley, and it was responsible for the 1869 Cachar EQ, the 1943 Assam EQ, and the 2009 Assam EQ (Kayal et al. 2010). Recent paleoseismic investigations by Kumar et al. (2016) revealed several EQs (MW !7.0) generated by the Kopili fault during the past 900 years.  Table 1.) SP is shown as the shaded area (modified after Baro and Kumar 2017). Indo-Burma ranges zone (IBRZ) The IBRZ, situated to the east and southeast parts of the SP-AVZ comprises of the Indo-Burmese mountains and the Tripura folds. A number of faults, such as the Kaladan fault, the Naga fault, the Churachandpur-Mao fault (CMF; also known as the Eastern boundary thrust fault), and the Kabaw fault, are associated with the Indo-Burmese plate boundary. According to Sikder and Alam (2003), the Kaladan fault extends from the Arakan coast in the south (exhibiting thrust faulting) to Nagaland in the north (exhibiting strike-slip faulting). The Kaladan fault extends toward north as the 400 km long Naga thrust (Wang et al. 2014). The CMF lays to the east of the Naga thrust in the Imphal valley and exhibits presently strike slip faulting (Kundu and Gahalaut 2013). The 170 km long CMF is capable to produce Mw 7.6 EQ (Wang et al. 2014), and the 280 km long strike-slip Kabaw fault located to the east of the CMF can produce even higher magnitude EQs in the future.

Bengal Basin zone (BBZ)
The BBZ is located to the west of the IBRZ and is composed of the deep alluvial deposits carried by the rivers Ganges and Brahmaputra during the Paleogene-Neogene times (Mohanty and Walling 2008). Compared to the SP-AVZ and IBRZ, this zone is seismically less active. Kayal (2008) attributes this trait to the deep alluvial deposits as well as the far-off distance of the BBZ from any of the plate boundaries. However, intraplate EQs had occurred within the BBZ in the past including the 1918 Srimangal EQ (mb7.6) on the NE-SW trending Sylhet fault. The 1935 Pabna EQ (MW6.2) is another intraplate EQ occurred in the 500 km by 25 km Eocene Hinge Zone of the BBZ (Figure 3) located to the south of Kolkata (Kayal 2008;Nath et al. 2014). Other major faults within this zone are the Garhmayna-Khandaghosh fault (GKF) and the Jangipur-Gaibandha fault (JGF) shown in Figure 3.

Eastern Himalayas zone (EHZ)
The EHZ is situated to the north of the SP-AVZ and comprises of 2500 km long Himalayan mountain range spreading from west to east (Kayal 2008). Several tectonic Note that names of some of the faults are unknown in the literature, and hence marked in the Table as Fxxx (see Figure 3 for the location of the faults). M obs is the maximum observed EQ magnitude; M p is the maximum possible EQ magnitude; and h is the focal depth of the observed EQ.
faults in this zone, such as the Indus Suture Thrust (IST) and the Main Central Thrust (MCT), run through the entire length of the Himalayas. Kayal (2008) divided the northeastern part of the Himalayas into the Trans-Axial Himalaya, the Central Himalaya, the Lesser Himalaya, and the Outer Himalaya. The Outer Himalaya separates from the Lesser Himalaya by the MBT (Ni and Barazangi 1984). The MCT, located to the north of the MBT, lies between the Lesser Himalayas and the Central Himalaya. The IST lies to the north of the MCT. No major EQs are reported in the recent times in the northeastern segment of IST even though it has potential to produce an EQ of Mw 8.5 (Srivastava et al. 2015).

Earthquake catalogue and the maximum possible earthquake magnitude
In this study, we use the EQ catalogue compiled by Baro and Kumar (2017) from three different sources: the United States Geological Survey, the India Meteorological Department, and the National Disaster Management Authority. The catalogue after removal of foreshocks and aftershocks contains 2359 EQs occurred between the years of 1411 and 2015 with the magnitude greater than 2.5. The EQ catalogue was divided into four sub-catalogues corresponding to the four seismic source zones discussed in the previous section. The completeness of each sub-catalogue was checked with respect to both magnitude and time, and the magnitude of completeness of each subcatalogue was determined to be 4.0 (Baro and Kumar 2017). Ground shaking due to an EQ is influenced by its magnitude, and hence, an assessment of the maximum possible EQ magnitude (Mp) is of importance in seismic hazard analysis. We estimate Mp for each of the 72 seismic sources. There are several statistical methods for Mp assessments at a fault using EQ catalogues (e.g. Kijko and Sellevoll 1989;Kijko and Graham 1998). Unfortunately, the number of EQs associated with individual faults does not allow for a reliable assessment of the Mp using statistical methods. Moreover, Mp cannot be solely determined from EQ catalogues (e.g. Holschneider et al. 2011), and thus an additional information should be used for this aim. This information includes the fault plane geometry and the maximum observed displacement on the fault; long-term geological and paleo-seismological assessments to reduce the errors in calculations of Mp; and others. In addition, numerical modeling of regional fault dynamics is a powerful alternative to overcome the shortcomings associated with the short time duration (about 100þ years) of recorded EQ catalogues (e.g. Ismail-Zadeh et al. 2015, Ismail-Zadeh et al. 2017

Ground motion prediction equations
An important step while performing seismic hazard analysis is to assess the level of ground shaking at the site of interest in terms of ground motion parameters. As regional EQ ground motion measurements are scarce and influenced by soil conditions, these ground motion parameters are obtained using appropriate GMPEs. Several region specific GMPEs have been developed so far to estimate the ground motion parameters (or seismic hazard level) at a site. Nath et al. (2005), Das et al. (2006), Raghukanth and Iyengar (2007), Sharma et al. (2009), Nath et al. (2009, Baruah et al. (2009), Gupta (2010), NDMA (2010) and Anbazhagan et al. (2013) developed GMPEs for the Himalayas and NE-India considering recorded or synthetic EQ data as well as their combination. The existence of several GMPEs could lead to the uncertainty associated with the degree of suitability of the GMPE for the study area. This uncertainty could be addressed by applying the logic tree approach, wherein appropriate weights or ranks suitable for the study area are added to the GMPEs ). Nath and Thingbaijam (2011) pointed out that many earlier studies paid a little attention to the selection of proper of GMPEs, and they adopted the efficacy test method by Scherbaum et al. (2009) for the selection and ranking of the most suitable GMPEs applicable to different regions of India.
Following Nath and Thingbaijam (2011), initially we have selected the GMPEs proposed by Hwang and Huo (1997), Toro (2002), Atkinson and Boore (2006), Raghukanth and Iyengar (2007) and Nath et al. (2009), and checked whether these GMPE can be applied within the seismotectonic region of the radius 500 km, and in the magnitude range of 4.0 to 8.7. Except for the GMPE by Toro (2002), the remaining GMPEs are either developed for shorter distance ranges or for lesser magnitude ranges. Therefore, for the present seismic hazard analysis we use the following GMPE (Toro 2002): where y is the spectral acceleration (in g); M w is the moment magnitude; R JB is the Joyner-Boore distance, that is, the shortest distance from a site to the surface projection of the rupture surface; e e is epistemic uncertainty; e a is aleatory uncertainty; and A 1 to A 6 are the constants (presented in Table 2 by Toro et al. 1997). R JB can be transformed into the hypocentral distance using the focal depth of EQs (listed in Table 1, column 6). Also following Nath and Thingbaijam (2011), we have selected the GMPE by Kanno et al. (2006) suitable for the Himalayas and the Indo-Myanmar regions: where r is the hypocentral distance in km; B 1 , B 2 , and B 3 are the regression coefficients; and e 2 is the error between observed and predicted values. Values of these coefficients for different oscillation periods are taken from Kanno et al. (2006). In addition, we use two other GMPEs developed by NDMA (2010) and by Anbazhagan et al. (2013) for the Himalayas. GMPE proposed by NDMA (2010) is represented as: where g is the acceleration due to gravity; C 1 to C 8 are the coefficients specific for northeast India; f 0 ¼ maxðln r = 100 À Á ; 0Þ, and e 3 is the standard error. Values of the coefficients for different periods are taken from NDMA (2010). Anbazhagan et al. (2013) developed a GMPE for the Himalayan region considering both recorded and synthetically generated ground motion data for 14 significant EQs, occurred at different segments of the Himalayan belt. Synthetic ground motions were developed for the past EQs with no ground motion records, which also included the 1897 Assam EQ. The GMPE developed by Anbazhagan et al. (2013) is as follows: where b is a decay parameter; D 1 , D 2 , and D 3 are regression constants; and r represents the standard error. Values of the coefficients for different periods are taken from Anbazhagan et al. (2013). Table 2 lists the GMPEs used in the present study with the magnitude range and distance range up to which each of the selected GMPE is valid. GMPEs considered in this study have been developed on the basis of limited ground motion data, and for this reason, the GMPEs have a standard error term. A logic tree approach for GMPEs provides a possibility to select more than one GMPE in seismic hazard analysis so that the standard error in the proposed seismic hazard value can be minimized. However, a relative weightage of selected GMPEs will influence the standard error term, and hence, a suitable weight should be assigned to each selected GMPE before using the logic tree ). Scherbaum et al. (2009) proposed to assess the appropriateness of each selected GMPE in terms of the Kullback-Leibler divergence D, which is the measure of how one probability distribution (e.g. for the 'true' ground motion) diverges from an expected probability distribution (e.g. for the ground motion predicted by a GMPE model) as: where f is the probability distribution resembling 'true' observations, which are actually unknown; g is the probability distribution representing a selected GMPE model, E f is the expected value from a given model. The Kullback-Leibler divergence becomes smaller with the probability distribution g approaching the probability  (2006) Equation (2) 5.5-8.2 r 450 NDMA (2010) Equation (3) 4.0-8.5 r 500 Anbazhagan et al. (2013) Equation (4) 5.3-8.7 r 300 distribution f. Scherbaum et al. (2009) suggested to calculate the second term in Equation (5) using the average sample log likelihood (LLH): where x 1 , x 2 , … , x N are samples of ground motion values obtained using a GMPE model g x i ð Þ . Following Scherbaum et al. (2009), a normal distribution with the probability density function (PDF) over range of x i (i ¼ 1, … , N) determined as: where m is the mean (expectation) of the normal distribution, r is the standard deviation, and r 2 is the variance. For this PDF, the value of LLH can be estimated as: Using the LLH value for each GMPE model, the data support index (DSI) can be calculated to evaluate the percentage by which the weight of a model is increased or decreased by the data used. The weight (w j ) of GMPE model g j (j ¼ 1, … , M) can then be calculated as: and the DSI of model g j as: where w unif is the uniform weight w unif ¼ 1=M, and M ¼ 4 in this study. To generate a set fx i g of the ground motion parameter (PHA) using a selected GMPE, we follow the following procedure. For each of three cities located in the districts of East-Khasi Hill, Ri-Bhoi, and West Garothe Shillong city (25 34 0 N, 91 53 0 E), Nongpoh (25 51 0 N, 91 49 0 E), and Tura (25 30 0 N, 90 12 0 E), (1) a minimal distance (we refer to it here as an 'epicentral distance') between the city and the nearest fault within each of four seismic source zones is calculated; (2) the maximal epicentral distance is assumed to be the radius of the seismotectonic region (i.e. 500 km) for each of four source zones; and (3) x i values (PHA values) are calculated from relevant GMPEs based on the range of epicentral distances (from minimal to maximal distances) and on the highest possible EQ magnitude (M p ) within each seismic source zone. From Table 1, the highest M p is 8.7, 8.0, 8.3, and 8.2 for SP-AVZ, IBRZ, BBZ, and EHZ, respectively. For example, for the Shillong city, the minimal epicentral distance from the IBRZ, BBZ, and EHZ source zones are 120 km, 60 km, and 160 km, respectively. Since the Shillong city is located within SPAVZ, the minimal epicentral distance is 30 km (the distance from the Shillong city to the Barapani fault). Hence, the range of epicentral distances for the Shillong city is 30 to 500 km, 120 to 500 km, 60 to 500 km, and 160 to 500 km from the SP-AVZ, IBRZ, BBZ, and EHZ seismic source zones, respectively. Similarly, the set of values fx i g is generated in the case of Nongpoh and Tura; the minimum epicentral distance within SP-AVZ is calculated from Nongpoh and Tura with respect to the Dapsi fault and the Barapani fault, respectively.
Based on the x i values so obtained, the LLH value for the GMPE model g j as well as w j , and DSI k are calculated using Equations (8)-(10). Only those GMPEs, which give positive DSI values, are identified and ranked in the order of the highest to the lowest DSI values. New w i values are then calculated for the GMPE models having positive DSI to obtain the weights for relevant GMPEs, which are used for seismic hazard assessments within each of four seismic source zones. Tables 3(a-c) present the LLH values, DSI values, ranks and weights (w i ) of different GMPEs obtained for each seismic source zone for the (a) East Khasi hills, (b) Ri-Bhoi and (c) West Garo hills districts, respectively.
For the three districts within the source zones of IBRZ, BBZ and EHZ, the GMPEs developed by Toro (2002), NDMA (2010) and Anbazhagan et al. (2013) are found to be suitable. Within SP-AVZ, the GMPEs developed by Toro (2002) and NDMA (2010) are suitable for the East Khasi hills and Ri-Bhoi districts, and the GMPEs by Toro (2002), NDMA (2010), and Anbazhagan et al. (2013) are suitable for the West Garo hills district. Note that the GMPEs by Toro (2002), NDMA (2010), and Anbazhagan et al. (2013) were developed for the site conditions with the shear wave velocity (V s ) of 1.8 km s À1 , 3.6 km s À1 , and 3.6 km s À1 , respectively. According to the NEHRP site classification scheme (BSSC 2004), the sites with V s > 1.5 km s À1 are considered to belong to site class A (a hard rock).
The ground motion predicted by the GMPE can be compared to that recorded during EQs to select the suitable prediction equations, when both the recorded and predicted data are assessed at the site class A. Here, we compare PHA recorded during two EQs at station Nongstoin (see Figure 1 for its location in the SP, at site class A according to PESMOS) with those predicted at the site of the station by three GMPEs (see Table 4). Also, Table 4 lists the data related to other recording stations, but located at site class B for comparison with the recording at station Nongstoin. For the Mw5.9 EQ occurred at a hypocentral distance of 408 km from the Nongstoin station, the recorded PHA is 8.88 cm s À2 , and predicted PHAs for the GMPEs developed by Toro (2002), NDMA (2010), and Anbazhagan et al. (2013) are 7.63, 17.07, and 8.77 cm s À2 , respectively. Meanwhile, for the Mw6.2 EQ occurred at a hypocentral distance of 200 km from the same station, the recoded PHA is 19.43 cm s À2 , and the PHAs predicted by the same three GMPE are 8.59, 18.39, and 8.28 cm s À2 , respectively. The PHA values obtained from Toro (2002) and Anbazhagan et al. (2013) are closely matching to the recorded maximum acceleration value for the first EQ, and the PHA values obtained from NDMA (2010) to the recorded maximum acceleration value for the second EQ. A difference between recorded and predicted values can be attributed to GMPEs due to limitations and uncertainties of the equations. Meanwhile, the comparison provides a confidence that GMPEs, which give positive DSI, predicts PHA values closer to recorded data, and thus, they can be used in seismic hazard assessments.

DSHA
Seismic hazard can be defined as a seismic phenomenon that 'may cause loss of life, injury or other health impacts, property damage, loss of livelihoods and services, social and economic disruption, or environmental damage' (UNISDR 2009). In seismological and EQ engineering community, seismic hazard is assessed using strong ground motion parameters, such as, peak ground acceleration and/or seismic intensity. This assessment is based on the knowledge of seismic wave excitation at the source, seismic wave propagation, its attenuation, and site effect in the region under consideration. Meanwhile, a comprehensive seismic hazard assessment should combine the results of seismological, geomorphological, geological, and tectonic investigations and modeling (e.g. Ismail-Zadeh 2014). At present, two basic methods of seismic hazard assessment are in use: probabilistic and deterministic. The probabilistic hazard analysis determines the rates of exceedance of certain levels of ground motion estimated over a specified period of time, and considers uncertainties in EQ source, magnitude, path, and site conditions (e.g. Cornell 1968). Meanwhile there is some criticism of the probabilistic assessment, and it is related to the validity of a point source model, ground motion uncertainties in the mathematical formulation of the method, and incapability to correctly model the dependencies between large numbers of uncertain random parameters (Panza et al. 2011). DSHA is based on specified EQ scenarios. Namely, for a given EQ, the attenuation of seismic energy with distance is assessed to determine the level of ground motion at a particular site using the available knowledge on EQ sources and wave propagation processes. Although the frequency of the ground motion is usually not addressed in DSHA, the method is considered to be robust for an assessment of seismic hazard and useful for a decision-making (Babayev et al. 2010).
The origins of DSHA could be traced back to nuclear power industry applications. It is in association with the fact that ground shaking due to a scenario-based EQ in the vicinity of critical infrastructure (e.g. nuclear power plants, dams, etc.) is to be Table 4. Comparison between recorded PHA and PHA predicted by GMPEs by Toro (2002), NDMA (2010) and Anbazhagan et al. (2013)  seriously considered. DSHA performance includes an identification of seismic sources located around a site of interest, which could influence the seismic hazard scenario at the site. In DSHA, an EQ is assumed to occur at the point of the fault closest to the site of hazard assessment. Thus, the shortest distance from the site to the source is calculated. Based on the information of past EQs, which had occurred around the site, M P is determined. Considering the identified sources, the shortest distance, and M P , the ground motions, which could occur at the site, are assessed. Several studied on DSHAs were performed across India. For example, Desai and Choudhury (2014) assessed seismic hazards for Greater Mumbai, Kumar et al. (2013) for the Lucknow city, and Parvez et al. (2003) for the entire country and the adjoining areas. Here, we use the scenario-based seismic hazard assessment for the SP.

Peak horizontal acceleration (PHA)
The M p values obtained in this work are employed in the selected GMPEs to perform DSHA for the study area. To obtain an unbiased measure of minimum hypocentral distances between the faults and study area, the SP is divided into 58 Â 21 cells along the longitude and latitude respectively. The size of each cell is 0.05 Â 0.05 . The minimum hypocentral distances between the 72 faults and the top left corner of each cell are determined. These minimum hypocentral distances and M p values are then employed in the selected GMPEs. The PHA values are estimated at site class A from the selected GMPEs at the top left corners of each cell. The PHA values are calculated taking into account the weights assigned to each of the GMPEs from the LLH calculations as discussed earlier. Because the seismotectonic region consists of four seismic source zones, four sets of PHA values are obtained in this study at each cell based on weights of each GMPE for that seismic source zone. The highest PHA value among the four is assigned to each cell to produce a PHA contour map. Figure 4 presents the results of the DSHA for three districts of (a) East Khasi hills, (b) Ri-Bhoi, and (c) West Garo hills. The maximum PHA values at Shillong city, Nongpoh, and Tura are 0.36 g, 0.46 g, and 0.33 g respectively.
The PHA values obtained on the basis of DSHA are juxtaposed with PHA contours maps developed by the NDMA (2010) based on probabilistic seismic hazard analysis (PSHA). It should be mentioned that while DSHA assesses the worst scenario and presents the values of the ground motion related to the worst scenario, PSHA allows determining the rates of exceedance of the specified ground motion within a desired time. Therefore, a direct comparison between the DSHA and PSHA results is meaningless.
The NDMA (2010) contours over the SP had PHA values of 0.40 g to be exceeded in 50 years at 0.02 probability, 0.45 g at 0.01 probability, and 0.55 g at 0.005 probability at the oscillation period 0s. Note that the contours developed by NDMA (2010) are for the entire SP at different percentages of probability of exceedance and not for individual cities as attempted in this study. The PHA values estimated in this study are further compared to the hazard maps for Manipur (an Indian state east to Meghalaya) developed by Pallav et al. (2012) at 0.02 probability of exceedance in 50 years. Pallav et al. (2012) found that the PHA values range from 0.4 to 1.1g. The highest PHA value of 0.46 g within the East Khasi hills district is obtained on the northern side of the district as shown in Figure 4(a). Similarly, the highest PHA value of 0.46 g is obtained for the eastern part of the Ri-Bhoi district [Figure 4(b)]. These high PHA values are due to the presence of the Barapani fault at that location. However, at other locations within the two districts, the Oldham fault is found to influence significantly the PHA values. Similarly, for the West Garo hills district the Oldham fault is responsible for the high PHA value toward the north. However, toward south of the district, the Dauki fault influences the PHA values. The Barapani fault, the Oldham fault, and the Dauki fault are all located within the SP-AVZ. This shows that the tectonic and geological setting of the SP-AVZ is such, that this zone has the potential to cause the highest seismic hazard among the four seismic source zones. This is in agreement with the results obtained in Baro and Kumar (2017), where a similar conclusion was derived based on the assessments of the b-value, the return period, and the probability of exceedance.

Response spectra
Further attempts are done to develop response spectra for the Shillong city, Nongpoh, and Tura. Response spectra are calculated at top left corners of the cells at the minimum hypocentral distance from each of the above 72 faults and for the corresponding M p values. The hypocentral distance and M p are then used in selected GMPEs with their appropriate weights. The coefficients in GMPEs (Equations (1)-(4)) are the function of the period of oscillations. This way, DSHA is performed for different periods, and the response spectra are developed for above three locations. Developed response spectra at site class A for the Shillong city, Nongpoh, and Tura are shown in Figure 5. The spectral acceleration (SA) at Nongpoh at 0 s, 0.1 s and 0.2 s are 0.46 g, 0.77 g and 0.55 g respectively, which are the highest values among the three sites. An increase in the spectral acceleration values can be observed after 0.2 s at Tura.
The response spectra developed in this study is in agreement with those developed earlier (Eurocode 8 2005; NDMA 2010) (see Figure 6). For example, the Eurocode 8 (2005) are matching well up to 0.1s for all the above three locations. However, the  predicted SA at 0s obtained for Nongpoh is higher than those the estimated on the basis of Eurocode 8 (2005). For the Shillong city, the response spectra from NDMA (2010) give the SA of 0.45 g and 0.80 g at 0 s and 0.1 s, respectively, which are higher than the values obtained in this study, while those for the periods greater than 0.1 s matches well the results of the present study [ Figure 6(a)]. A similar trend can be observed for Nongpoh [ Figure 6(b)], where the higher values of 0.59 g and 0.93 g are attained at 0 s and 0.1 s, respectively, and a closer match with the present studies for the periods greater than 0.2 s. For Tura [ Figure 6(c)], the obtained SA values of 0.35 g and 0.63 g at 0 s and 0.1 s, respectively, are similar to the present study, while not well matching the results of the present study at the periods greater than 0.2 s.

Conclusion
DSHA is performed for the three districts -East Khasi hills, Ri-Bhoi, and West Garo hillswithin the Shillong Plateau. We have obtained that the PHA values range from 0.27 g to 0.46 g within these districts. For the East Khasi hills and the West Garo hills districts, the northern parts show higher PHA values of 0.46 g. For the Ri-Bhoi district, higher PHA value of 0.46 g is observed in its eastern part. In terms of engineering parameters, the higher PHA values, the higher seismic hazards. We consider that the higher PHA values in the eastern part of the Ri-Bhoi district and the northern part of the East Khasi hills district are associated with the Barapani fault. Meanwhile, the Oldham fault influences the PHA values for the remaining areas within the two districts. The higher PHA value in the northern part of the West Garo hills district is also due to the Oldham fault. Therefore, we propose that the Oldham fault is a potential source of high seismic hazard in the SP.
The DSHA maps for the East Khasi hills, Ri-Bhoi, and West Garo hills districts have been compared to the PHA contours maps developed by NDMA (2010) and Pallav et al. (2012), and the comparison shows a reasonable agreement. The response spectra developed for the Shillong city, Nongpoh, and Tura have also been compared to the response spectra developed using Eurocode 8 (2005), and the GMPE given by NDMA (2010). Except of some deviations in the response spectra associated with the use of multiple GMPEs and LLH weights in this work, the comparison shows a good match with our results. Findings from the present study provide detailed seismic hazard assessment for the three densely populated districts of the SP.