Seismic risk assessment for the North Eastern Region of India by integrating seismic hazard and social vulnerability

1

ABSTRACT

The present study aims at conducting a comprehensive seismic risk assessment for the North Eastern Region of India at regional and sub-regional levels by integrating probabilistic seismic hazard and social vulnerability assessments. Bedrock-level peak ground acceleration varied from 0.14 to 0.69g for the return period of 475 years. Using PCA, the social vulnerability index (SVI) was generated considering district-level socioeconomic indicators. Built environment quality, illiteracy, access to amenities, dependent population, and employment opportunities contributed to high SVI. Most vulnerable districts were concentrated in the Brahmaputra floodplains, Tripura fold belt, and Imphal valley. At the regional level, significant parts of Assam, Meghalaya, Arunachal Pradesh, and Tripura lie in moderate to very high-risk zones. At the sub-regional level, Nagaland accounts for the highest proportion of areas in high to very high-risk zones. The findings will aid site-specific resilient infrastructure design, disaster risk reduction, and effective resource allocation for the risk-prone areas.

KEYWORDS:

• PSHA
• vulnerability index
• cluster analysis
• disaster risk reduction
• geospatial analysis

1. Introduction

The North Eastern Region (NER) of India, located in the Eastern Himalayas, is one of the most seismo-tectonically active regions in the world (Das et al., 2012) and is placed in seismic zone V (highest-risk zone; except Sikkim) as per seismic zonation map of India (IS-1893, 2016). Geographically, the region has a convergent plate boundary between the Indian and Eurasian plates in the north resulting in the formation of the Himalayan fold mountains and active subduction under the NNE-SSE trending Indo-Burmese ranges along the eastern boundary of the country, which are responsible for region’s complex geology and tectonic setting (Agrawal et al., 2022; Baro et al., 2018). Along with these features, the NER also consists of features such as Brahmaputra Valley Plain, Bengal Basin, Mishmi Hills, Shillong-Mikir Plateau, and Tripura Fold Belt (Verma & Kumar, 1987).

The influence of hazards in a specified time and space may vary in terms of environmental, physical, and socioeconomic losses (Agrawal et al., 2021). This region has suffered extensive loss of lives and damage to property due to several significant earthquakes, including two great events of moment magnitude (Mw) >8.0, namely, the 1897 Shillong earthquake and the 1950 Assam earthquake (Bilham & England, 2001; Mrinalinee Devi & Bora, 2016; Thingbaijam et al., 2008). Recently, on 28th of April 2021, an Mw 6.0 earthquake occurred near Dhekiajuli in Assam, India, leading to ground cracking and the collapse of several houses. The consequences of a seismic event are not only restricted to structural collapse but can also cause soil liquefaction, sand boiling, and landslides (Dixitet al., 2012a; Dixit et al., 2016; Gupta et al., 2021; Mase, 2020). The topographic features may further amplify the ground motion in a seismic event (Agrawal & Dixit, 2022b). The region has seen massive infrastructure development and unplanned fast urbanization during the previous few decades. According to the country’s most recent Census (Census of India, 2011), NER has 14 towns and urban agglomerations, each having a population of more than 100,000 people. Compared to the time of occurrence of previous significant earthquakes, the increased population density has enhanced the region’s susceptibility to seismic hazards. Therefore, to minimize potential damages or loss of lives and for effective disaster risk management, reduction, and mitigation strategies against any future earthquake, it is essential to quantify the seismic risk by assessing the seismic hazard and vulnerability at the regional and sub-regional levels, which will aid the concerned authorities in undertaking appropriate actions. Mathematically, seismic risk can be defined as the product of hazard and vulnerability (D’Amato et al., 2022; Frigerio et al., 2016).

The seismic hazard assessment (SHA) deals with the quantitative estimation of ground motion at a particular site in a specific time interval, either deterministically or probabilistically. It can be expressed as peak ground acceleration (PGA), spectral acceleration, or other ground motion parameters (Dixitet al., 2012b; Kramer, 1996; Mase & Sugianto, 2021). The probabilistic seismic hazard assessment (PSHA) is based on the total probability theorem. It considers various uncertainties associated with earthquake magnitude size, its spatial and temporal occurrence, and the predictive relationship lacking in the deterministic approach (Cornell, 1968; Mase, 2022; McGuire, 2008). In the past, several researchers like Bhatia et al. (1999), Nath & Thingbaijam (2012), NDMA (2010), Parvez et al. (2003), and Sitharam & Kolathayar (2013) have conducted seismic hazard studies at the national scale using probabilistic and deterministic approaches. SHA studies are also available at the regional scale, such as S. Das et al. (2006), Das et al. (2016), Ghosh & Chakraborty (2017), Sharma & Malik (2006), and Ghione et al. (2021) for the NER, Sitharam & Sil (2014) for the state of Tripura and Baro et al. (2018, 2020) for the Shillong Plateau, Meghalaya, and Bahuguna & Sil (2020) for the state of Assam along with several other studies contributed different aspects of seismic hazard in the region (Dixit et al. 2016; Raghukanth & Dash 2010; Raghukanth et al. 2011,2009). However, these studies have only concentrated on the hazard component and have ignored the vulnerability aspect, an equally important issue.

Vulnerability is the potential for damage or loss of societal assets. It can be described as the characteristics and situation of people or assets in a region that make it susceptible to negative consequences of a hazard (Cutter et al., 2008, 2003). It could be a physical, socioeconomic, environmental vulnerability, or a combination. The physical vulnerability deals with the building stock in the vicinity and its susceptibility to hazards. In contrast, socioeconomic vulnerability (SV) refers to a community’s current socioeconomic conditions and demographic features that can affect the community’s ability to prepare for, respond to, manage and cope with the effects of hazards (Ogie & Pradhan, 2019). The SV is generally measured in terms of the social vulnerability index (SVI), generated based on various indicators such as population composition, employment opportunities and occupation, literacy, quality and condition of the built environment, and other basic facilities that influence the resilience and coping capacity of a community. Several risk assessment studies have been conducted considering either socioeconomic component (Armaş, 2008; Cerchiello et al., 2018; Martins et al., 2012) or physical vulnerability alone (Del Gaudio et al., 2017; Fuentes et al., 2021; Laguardia et al., 2022).

Natural hazard risk assessment considering socioeconomic vulnerability involves multivariable selection and analysis using different approaches. Gupta & Dixit (2022a) have performed a flood risk assessment for the Assam region using the MCDA-AHP approach. Lee et al. (2019) integrated three novel and efficient models involving SWARA, radial basis function (RBF), and Teaching–Learning-Based optimization (TLBO) methods to assess the seismic vulnerability of the Tehran region, the capital of Iran. Banica et al. (2017) used Principal Component Analysis (PCA) for indicators’ processing and correlation analysis. With the help of the AHP, a seismic vulnerability index was generated for Iasi, one of the major earthquake-prone cities in Romania. SV studies for natural hazards have been performed for different hazard-prone regions by several researchers like Ge et al. (2013) for the Yangtze River delta, China, Siagian et al. (2014) for Indonesia, and Gautam (2017) for the Nepal region.

Similarly, Derakhshan et al. (2020), Frigerio et al. (2016), Zhang et al. (2017) and Mase (2020) have performed SV studies for seismic hazards in different parts of the world. In the Indian scenario, including NER, seismic hazard-vulnerability studies primarily focused on the physical or built environment (Baruah et al., 2020; Dutta et al., 2021; Joshi et al., 2019). Furthermore, only a few studies have performed a seismic risk assessment, considering SV or coping capacity in the Indian scenario (Agrawal et al., 2021; Jena et al., 2021). However, no comprehensive study evaluating seismic hazards and SV for seismic risk assessment are available for Eastern Himalayas.

The main goal of the present study is to generate a seismic risk map of NER and analyze risk scenarios at regional and sub-regional levels by performing PSHA in terms of PGA at the bedrock level and SV assessment using PCA. The present study is the first such attempt to assess the seismic risk for the NER of India (except Sikkim) in Eastern Himalayas, considering SV at regional and sub-regional levels. In the subsequent sections, a detailed discussion regarding the study area, methodology, and results are made.

2. Topography, geology, and seismotectonic setting of the study area

The study area lies between 21°56ʹN to 29°27ʹN latitude and 89°41ʹE to 97°24ʹE longitude, consisting of the Seven-Sister States, namely Arunachal Pradesh, Assam, Manipur, Meghalaya, Mizoram, Nagaland, and Tripura, of the NER excluding Sikkim (Figure 1). Sikkim does not have a common boundary with any of these seven states. The study area shares its boundary with neighboring countries like Bhutan, Tibet (China) in the north, Myanmar in the east, Bangladesh in the southwest, and the neighboring state, West Bengal, in the west. The study area covers ~253,353 km² of the area, and more than 45 million people (Census of India, 2011) reside there.

Figure 1. Demographic profile of the study area.

The population density ranges from 17 person/km² (for Arunachal Pradesh) to 398 person/km² (for Assam) in the region, and a large chunk of the population resides in the alluvial-plain regions (the Brahmaputra, Tripura piedmont, and Imphal valley plains). The region has a predominantly rural population (except Mizoram), as shown by the histogram in Figure 1. Agriculture and allied activities contribute the most to the economy (Gupta & Dixit, 2022b). In recent years, the region is also witnessing an increase in tourism activities. Therefore, there is a shift of housing types from traditional to modern-multi story buildings, increasing the region’s vulnerability.

Physiographically, the entire study area can be classified into the Eastern Himalayas, Brahmaputra Valley Plains, Shillong-Mikir Plateau, Barak Valley, Arakan-Yoma, and Naga Hills (Agrawal & Dixit, 2022b; Das et al., 2012; Gupta & Dixit, 2022a; Verma & Kumar, 1987). The study area is characterized by varying geological units ranging from Precambrian rocks in the upper reaches and Shillong Plateau region to recent Quaternary sediment deposits of 3 to 5 km thickness in Brahmaputra Valley Plains (Nandy, 2001). The seismotectonic map of the study area, as shown in Figure 2a, was prepared from the Geological Survey of India (GSI, 2000) and available literature (Das et al., 2012; NDMA, 2010; Thingbaijam et al., 2008).

Figure 2. a) Seismotectonic map of the study area. b) Seismic source zonation of the study area showing the spatial distribution of declusterd earthquake.

In the northern portion, from south to north, the major faults are the Main Frontal Thrust (MFT), Main Boundary Thrust (MBT), and Main Central Thrust (MCT), which are north dipping listric thrusts. MFT, also known as Himalayan Frontal Thrust (HFT), is the southernmost Himalayan thrust, separating the Brahmaputra basin from the Outer Sub-Himalayas (Figure 2a). MBT separates the Outer Sub-Himalayas from the Lesser Himalayas, which is bounded by MCT in the north (Gupta et al., 2022; Kayal, 2008).

The northeast portion of the study area, the Eastern Syntax zone, is the tri-junction point of the Indian, Eurasian, and Burmese plates and has a high crustal deformation zone due to plate tectonics, which was the place of the 1950 Assam earthquake. The major seismotectonic features in this region are the NW-SE trending Mishmi Thrust, Lohit Thrust, and the right-lateral Po Chu Fault (Figure 2a). The NE-SW trending Naga Thrust and Disang Thrust are significant features in the central portion of the study area, forming the southeastern boundary of the Brahmaputra valley plains in the region (Figure 2a). The Shillong Plateau results from ‘pop-up tectonics’ and is confined by Oldham Fault in the north and the E-W trending Dauki Fault in the south (Kayal et al., 2006). The Oldham Fault is a south-dipping fault that was the source of the 1897 Assam earthquake (England & Bilham, 2015), whereas the north-dipping Dauki Fault was the source of the 1923 Meghalaya earthquake of magnitude 7.1 (Baro et al., 2020). The Dhubri Fault, the source of the 1930 Dhubri earthquake of Mw 7.1, is another prominent tectonic feature in the region, forming the western boundary of the Plateau (Kayal, 2008). The NW-SE trending Kopili Fault, located east of the Oldham Fault, separates the Mikir Hills from the Shillong Plateau (Agrawal & Dixit, 2022a). This fault generated 1869 Cachar and 1943 Assam earthquakes of magnitudes greater than 7.0 (Dasgupta, 2011).

Along with these features, Eocene Hinge Zone, Sylhet, Churachandpur-Mao Fault, and Shan Sagaing Faults are other critical features in the southwest and south of the study area responsible for several smaller to moderate size earthquakes (Thingbaijam et al., 2008). The study area is divided into five distinct source zones in terms of fault properties, seismic source, geology, and plate tectonics for seismic source zonation (Das et al., 2012; Thingbaijam et al., 2008), shown in Figure 2b.

3. Methodology

3.1. Seismic hazard assessment

In the present study, PSHA was performed at the regional and sub-regional levels of the study area. PSHA involves seismic source zonation, earthquake recurrence relation, and the ground motion attenuation equations to produce seismic hazard curves regarding ground motion level and an associated annual frequency of being exceeded. PSHA begins with preparing the earthquake catalog and homogenization to a common scale. The earthquake catalog was prepared considering a seismo-tectonic region of a 500 km radius, centered on the coordinate point of 26.16°N and 93.28°E, shown in Figure 2b (Anbazhagan et al., 2015). The earthquake data was compiled from international and national databases such as USGS, International Seismological Centre (ISC), NCS India, and Bhukosh-GSI (Baro & Kumar, 2017; Gupta et al., 2021; Thingbaijam et al., 2008). A total number of 8959 earthquake events, with magnitude, hypocentral depth, and time and location of occurrence, were collected from 1760 to 2021 (261 years) from these databases. The obtained earthquake catalog comprised events in different magnitude scales, such as body-wave magnitude (m_b), surface-wave magnitude (M_S), local magnitude (M_L) and moment magnitude scale (Mw), and converted to a common scale using orthogonal regression relation (Table 1, Appendix A). The homogenized earthquake catalog was further processed to eliminate the duplicate and dependent events from the dataset by performing declustering of the catalog by adopting the classical window technique of Gardner & Knopoff (1974) using an open-source software ZMAP (v 7.0) (Wiemer 2001; Anbazhagan et al., 2019; Sitharam & Sil, 2014). After declustering, the dataset contains 6599 events within the considered seismotectonic region, among which 4837 events are larger than Mw 3.5 (Figure 2b). The completeness of the homogenized-declustered earthquake dataset with respect to time was also checked using the statistical analysis proposed by Stepp (1972) discussed in detail in Appendix A. Based on the completeness study, the earthquake catalog of the recent 70 years was considered to develop the Gutenberg-Richter (G-R) recurrence relation for each source zone (Appendix A). The largest possible earthquake, M_max, was evaluated based on the conventional incremental value method, IVM (Anbazhagan et al., 2019; Gupta, 2002), and the procedure suggested by Kijko & Sellevoll (1989). The seismic activity rate of an individual fault was evaluated by following an approach similar to Raghukanth & Iyengar (2006) and NDMA (2010). The detailed discussion and methodology of homogenization, declustering, completeness study, maximum magnitude calculation, and fault-level deaggregation are mentioned in Appendix A.

Table 1. Relation between different magnitude scales.

Scale	Data points	Relation		R^2
mb to Mw	349	Mw = (1.008 ± 0.018) mb – (0.095 ± 0.086)	2.4 ≤ mb≤ 6.9	0.896
ML to Mw	277	Mw = (0.919 ± 0.022) ML + (0.286 ± 0.085)	2.5 ≤ ML ≤7.0	0.864
MS to Mw	130	Mw = (0.715 ± 0.031) MS + (1.796 ± 0.151)	3.0 ≤ MS ≤ 7.2	0.803

3.1.1. Selection of ground motion prediction equation (GMPE)

Knowing site-specific attenuation relation is significant for evaluating ground motion parameters (Mase et al., 2018). The ground-motion prediction equations (GMPEs) are region-specific and should preferably be developed. However, sometimes due to the limited quality and quantity of available strong ground motion data, it is not always possible to develop own ground motion prediction model. In such scenarios, various researchers have shown in the past that one can use the previously developed model for the same region or other region based on similar seismotectonic characteristics (Anbazhagan et al., 2015; Bahuguna & Sil, 2020; Das et al., 2016; Ornthammarath et al., 2011; Pallav et al., 2012; Tanapalungkorn et al., 2020).

In the present study, six GMPEs suitable for the Himalayan region were adopted (Table 4) based on the available literature (Anbazhagan et al., 2019; Ramkrishnan et al., 2021) and validated through the recorded strong motion data for the study area obtained from COSMOS strong-motion center (www.strongmotioncenter.org) (Sitharam & Sil, 2014). For this purpose, calculated PGA vs. hypocentral distances was obtained for different combinations of magnitude and focal depth using the selected GMPEs (Table 2), as shown in Figure 3a-3e. The observed PGA values from strong ground motion records for the same magnitude – focal depth (Mw, h) combination are also plotted in Figure 3a-3e. The graph shows that for all Mw-h combinations, ANBA2013 and NATH2012 predicted PGA values less than the observed ones. ATBO2003 overestimates the PGA for events with smaller focal depths (Figure 3b, ,) and underestimates for larger depths (Figure 3a and 3e). BAHU2020 generally underestimates the PGA values, making it the lower bound in a few cases except for Mw 7.3, h 90 combinations (Figure 3).

Figure 3. Comparison of GMPEs with the observed PGA values for a different combination of moment magnitude and hypocentral distance (km): (a) Mw 5.0, h 43, (b) Mw 5.9, h 15, (c) Mw 6.0, h 17, (d) Mw 6.0, h 34, (e) Mw 7.3, h 90.

Table 2. List of selected GMPEs.

Sl. No.	GMPE	Abbreviation	Remark
1	Jain et al., 2000 non-subduction zone: $ln(PGA)=-3.443+0.706M-0.8028ln(R)$ with SE = 0.44 subduction zone: $ln(PGA)=-0.332+0.00233R+0.59ln(R)$ with SE = 0.59 Where PGA in g, R is the shortest source-to-site distance, and SE is the standard error	JAIN2000	For Central-Himalayan region
2	Atkinson & Boore, 2003$logY=c_1+c_{2}M+c_{3}h+c_{4}R-g\cdot logR+s_l(c_5{S}_C+c_6{S}_D+c_7{S}_E)$Where Y in cm/s^2,$R=\sqrt{D_{fault}^2+{\mathrm{\Delta }}^2}$, $\mathrm{\Delta }=0.00724\times {10}^{(0.507M)}$, c1 = 2.991, c2 = 0.03525, c3 = 0.00759, c4 = −0.00206, g = 10^(1.2−0.18M), σ1 = 0.20 (intra-event) and σ2 = 0.11 (inter-event) for interface events (h < 50 km) and c1 = −0.04713, c2 = 0.6909, c3 = 0.01130, c4 = −0.00202, g = 10^(0.301−0.01M), σ1 = 0.23 and σ2 = 0.14 for in-slab events (h > 50 km) and (c5,c6,c7) = (0.19, 0.24, 0.29) for all events. $s_l$ is frequency-dependent constant. SC,SD, and SE are equal to zero for site class B (NEHRP), VS,30 > 760 m/s.	ATBO2003
3	Nath et al., 2012$ln(P)=9.143+0.247M-0.014(10-M{)}^3-2.67ln(r_{rup}+32.9458e^{(0.0663M)})$where P = PGA (g), $r_{rup}$= fault-rupture distance (km). Standard Deviation = 0.330.	NATH200	For the Shillong region of NE India
4	Anbazhagan et al., 2013$log(y)=-1.283+0.544M+b\cdot log(X+e^{(0.381M)})+\sigma$where Y = Spectral acceleration (SA) (g), $X=\sqrt{R^2+h^2}$, where R = epicentral distance (km), h = focal depth (km), b is decay parameter (−1.792), and σ = 0.283 (for 0 s period).	ANBA2013	For Himalayan region
5	Ramkrishnan et al.,2021$logy=-2.135+0.437M-1.099log(X+e^{(-0.080M)})\pm 0.549$where y = PGA (g), X = hypocentral distance (km), and SE = ±0.549	RAMK2021	For Central-Himalayan region
6	Bahuguna & Sil, 2020$ln(PGA)=6.680+1.134M-0.001R-0.7098lnR$where PGA in g, R = hypocentral distance (km) (ranges from 100–400 km), and standard deviation = 0.53	BAHU2020	For Assam region

Table 3. List of common social vulnerability indicators and their variables.

Sl. No.	Indicator	Variables
1	Population composition	P01	Population density
2		P02	Rural population (%)
3		P03	Female population (%)
4		P04	Population belongs to the socially backward class (%)
6	Age	Age01	Age less than 07 (%)
7		Age02	Age group of 07 to 60 (%)
8		Age03	Above the age of 60 (%)
9	Literacy	L01	Effective literacy rate
10		L02	Illiterate (%)
11		L03	Illiterate female (%)
12	Employment status and occupation	EO01	Population belongs to MW^1 class (%)
13		EO02	Female population belongs to MW^1 class (%)
14		EO03	Population belongs to the OMW^2 class (%)
15		EO04	Female population belongs to the OMW^2 class (%)
16		EO05	Population belongs to MrW^3 class (%)
17		EO06	Female population belongs to MrW^3 class (%)
18		EO07	Population belongs to the OMrW^4 class (%)
19		EO08	Female population belongs to the OMrW^4 class (%)
20		EO09	Non-permanent employment (%)
21		EO10	Female population with non-permanent employment (%)
22		EO11	Non-working population (%)
23		EO12	Non-working female population (%)
24	Building material	BM01	With brick or stone roof (%)
25		BM02	With kutcha roof (%)
26		BM03	With kutcha wall (%)
27		BM04	With kutcha floor (%)
28	House condition	HC01	Residential houses in dilapidated condition (%)
29		HC02	Residential cum other houses in dilapidated condition (%)
30	Family size	HH01	Houses with 4-5 households (%)
31		HH02	Houses with 6 or more households (%)
32	Amenities	A01	Houses with no electricity and have dependence upon kerosene or other oil as a source of light (%)
33		A02	Houses with no water source within or near the premises of the house (%)

¹MW: Main Workers; Workers who worked for more than six months in the reference period

OW: Other Workers; Workers other than cultivators, agricultural laborers, or household workers. e.g. Government servants, municipal employees, teachers, bankers, trade, commerce, etc.

²OMW: Other Main Workers; Main workers falling under OW

³MrW: Marginal Workers; Workers who worked for less than six months

⁴OMrW: Other Marginal Workers; Marginal workers falling under OW

Table 4. Selected variables based on PCA.

Factor	Extracted variables	Eigenvalue	Variance explained (%)	Weightage factor (wi)
1	HC01, A01, L02, P02, L03, Age01, HC02, BM03, EO12, HH02, EO11, BM02, and A02	19.528	37.87	0.42
2	BM01, EO03, EO07, HH01, Age03, and P04	2.080	27.35	0.30
3	EO02, EO06, EO05, EO10, and EO01	1.026	25.31	0.28

The PGA values obtained using the JAIN2000 model are well within the range of observed values for focal depths more than 45 km, and the RAMK2021 model was found to be suitable for focal depths less than 45 km (Figure 3). As a result, the GMPEs RAMK2021 (Ramkrishnan et al., 2021) for seismic source zones with an average focal depth of less than 45 km and JAIN2000 (Jain et al., 2000) for seismic source zones with an average focal depth ≥45 km were employed in this study.

3.1.2. PSHA for NER of India

To evaluate seismic hazards at the bedrock level, using a probabilistic approach, the entire study area was divided into grids of the size of 0.2° × 0.2°. Each grid center was considered a site of interest at which the seismic hazard in terms of PGA was evaluated by considering all the active seismic sources within a radius of 500 km (Anbazhagan et al., 2019).

The procedure followed for PSHA assumes that an event within a seismic source follows a stationary Poisson process (Kramer, 1996). The probability of ground motion parameter, Y, exceeding a specified level, y, in a specified period T, at a given site is expressed as

(1) $P(Y>y)=1-exp({\mu }_{y}T)$(1)

Where µ_y is the mean annual rate of exceedance as detailed below

(2) ${\mu }_y=\sum _{i=1}^n{N}_i(m_{min})\underset{m}{\int }\underset{r}{\int }P(Y>y|m,r)p_{R|M}(r|m)p_M(m)drdm$(2)

Where n is the number of faults, N_i(m_min) is the annual frequency of events on an i^th fault having m ≥ m_min, p_M(m) is the probability density function (PDF) corresponding to the magnitude, p_R|M (r|m) is the conditional PDF corresponding to hypocentral distance (r), and P(Y > y| m, r) is the probability of exceedance of ground motion parameter, Y, over y, for given magnitude m occurring at a distance r from the site.

Thus, µ_y incorporates the temporal, spatial, and magnitude uncertainty of a future event and ground motion uncertainty produced by the seismic event at the site. Equation 2(2) ${\mu }_y=\sum _{i=1}^n{N}_i(m_{min})\underset{m}{\int }\underset{r}{\int }P(Y>y|m,r)p_{R|M}(r|m)p_M(m)drdm$(2) shows the summation of individual contributions of 21 faults (i = 1, 2, 3, …, 21) for the SHA at each site to obtain the annual exceedance of PGA. All the calculations mentioned above were performed using MATLAB.

To produce the hazard curve considering all the sources, µ_y of a particular site of interest was summed up and plotted against the target PGA level, y. A typical hazard curve obtained for Shillong City is shown in Figure 4. Using the hazard curves, the PGA value for 10% probability of exceedance in 50 and 100 years of return period and 2% probability of exceedance in 50 years of return period are obtained. Correspondingly, thematic maps were produced using ArcGIS.

Figure 4. Seismic hazard curves (Shillong city).

3.2. Social vulnerability assessment

3.2.1. Dataset and its processing

SV depends on the indicators like population, age, gender, literacy, employment status, stock of built structures and their occupancy, etc. (Cutter et al., 2003; Depietri, 2020; Fatemi et al., 2017; Siagian et al., 2014; Wood et al., 2010). For the present study, data regarding the vulnerability indicators, comprising 54 variables, were collected at the district level for all 82 districts from India’s 15^th housing and population census (Census of India, 2011) based on a detailed and critical literature review (Agrawal et al., 2021; Armaș & Gavriș, 2013; Cutter et al., 2008, 2003; Derakhshan et al., 2020; Izquierdo-Horna & Kahhat, 2020; Schmidtlein et al., 2011). These selected variables influence SV and positively or negatively impact the same (Martins et al., 2012). These selected factors influence people’s resilience, coping capacity, and recovery ability after a regional hazard (Frigerio et al., 2016). The district-level SVI was calculated based on the methodology proposed by (Cutter et al., 2003), which considers social, economic, and demographic indicators for accessing the vulnerability of a place (Agrawal et al., 2021; Cutter et al., 2008).

The collected dataset was first normalized to transform the values into a standard scale (Dintwa et al., 2019). Multi-collinearity analysis was then performed to remove the variables with very high correlation, and a subset of 33 variables was retained and utilized further for the analysis (Table 3).

3.2.2. Factor extraction using PCA and generation of SVI

The principal component analysis (PCA) is a prevalent factor analysis and data reduction technique that reduces a large set of variables to a smaller set of orthogonal components (principal components), explaining most of the variance among the applied dataset (Jolliffe, 2002). In factor analysis, the sample adequacy of the applied dataset is checked using KMO (Kaiser-Meyer-Olkin) and Bartlett’s test. The KMO statistics values range from 0 to 1, and values close to 1 suggest an adequate dataset to proceed with factor analysis or PCA. For PCA, the KMO value should be at least 0.6; for PCA results to be reliable, the value should be greater than 0.8 (Chakraborty et al., 2020). Bartlett’s sphericity test was also employed to confirm the factor selection, which tests the null hypothesis for the correlation matrix. The output was in the form of the p-value of chi-squared statistics. The small p-value (<0.05) suggests that the pairwise correlation matrix is not an identity matrix, and the dataset is suitable for PCA (Sharma, 1996). In the present study, a KMO value of 0.897 (>0.8) and a p-value of 0.000 (<0.05) in Bartlett’s test were obtained, which indicates sufficient data adequacy for statistical analysis.

The dataset was then entered into the PCA using IBM SPSS 26.0 software. The varimax rotation with Kaiser Normalization and eigenvalue rules, i.e., factors with eigenvalue >1.0, was utilized for factor selection. The number of factors extracted was further confirmed by Cattell’s scree plot method (Cattell, 1966), as shown in Figure 5. Three factors comprising 24 variables together explaining ~ 90.5% variance of the employed dataset were extracted (Table 3 and 4).

Figure 5. Scree plot.

After the factors were extracted, the factor scores obtained from PCA corresponding to the factor were aggregated to generate the composite SVI score at the district level using Equation 3(3) $SVI=\sum _{i=1}^{n_f}{w}_i\times \mathrm{F}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{i}$(3) (Ge et al., 2013).

(3) $SVI=\sum _{i=1}^{n_f}{w}_i\times \mathrm{F}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{i}$(3)

Where n_f is the number of factors and

$w_i=\frac{varianceexplainedbyfactori(\mathrm{\%})}{totalvarianceexplainedbyallthefactors(\mathrm{\%})}$

The obtained SVI scores were mapped in the GIS environment to produce the SVI map for the study area. The SVI scores are normalized to a scale of 0 to 100 using Equation 4(4) $SV{I}_j^N=\frac{\left(SV{I}_j-SV{I}_{min}\right)}{\left(SV{I}_{max}-SV{I}_{min}\right)}\times 100$(4) for ease of interpretation (Chakraborty et al., 2020) and classified into five vulnerability classes using the quantile classification scheme. In the present study, the comparative assessment of SV at a sub-regional scale was also performed by generating the SVI for each state individually.

(4) $SV{I}_j^N=\frac{\left(SV{I}_j-SV{I}_{min}\right)}{\left(SV{I}_{max}-SV{I}_{min}\right)}\times 100$(4)

3.3. Spatial cluster analysis

Spatial autocorrelation was performed to analyze the autocorrelation of the dataset throughout the study area, and Global Moran’s I value, which ranges from −1 to 1, was obtained (Karuppusamy et al., 2021). A spatial statistical tool for hotspot analysis (Getis-Ord Gi*) in ArcGIS was employed to identify the spatial clusters within a specific area (Brandt et al., 2020). This tool examined the rank of each district within the context of the rank of neighboring districts. The hotspots were located based on the values of statistically significant z for 99, 95, and 90% confidence levels. Based on the calculated Z-scores and p-values, the districts of the study area were classified and grouped as spatially significant hot or cold spots. Typically, the hotspots exhibit higher z scores and lower p scores (Al-Dogom et al., 2018). In this analysis, a zone of indifference was selected for the spatial relationship conceptualization, and a threshold distance of 71,542 m was used. False discovery rate (FDR) correction was applied to identify spatial clusters at the local level better.

3.4. Seismic risk assessment

For the seismic risk assessment, the derived hazard values in PGA for the 475-year return period were considered and normalized to generate the seismic hazard index (SHI) and classified into five hazard classes. The classified SHI map was then multiplied with SVI scores to quantify the risk in the region. The values of 1 to 5 were assigned to SVI and SHI classes such that 1 represents the very low class and 5 represents the very high class. The generated risk matrix is shown in Figure 6, where values range from 1 (lowest) to 25 (highest) after the multiplication. The risk matrix values are classified into five groups, viz. very low (1–2), low (3–6), moderate (7–10), high (11–16), and very high (17–25) following Derakhshan et al. (2020) and Frigerio et al. (2016). Following the same, a seismic risk zonation map for the study area was prepared.

Figure 6. Risk matrix.

4. Result and discussion

In the present study, the results of PSHA are presented in terms of PGA at bedrock level, which was obtained from the hazard curve, that is, PGA vs. mean annual rate of exceedance. The SVI was generated by applying the FA and PCA as the factor extraction method. The results of SHA and SVI were then integrated to prepare the risk zonation map for the study area.

4.1. Seismic hazard assessment

The probabilistic approach of SHA was adopted, and the results are presented as bedrock-level PGA (Figure 7). For this purpose, seismicity parameters (a and b) and magnitude of completeness (M_C) were estimated. Using G-R law (Gutenberg & Richter, 1944), the estimated b value ranges from 0.77 (zone 2) to 1.02 (zone 5). The value ranges from 3.67 to 4.90, as shown in Figure A.3 and Table A.1. The obtained M_C value of the considered dataset varies from 3.4 (zone 3) to 3.7 (zone 5) for different seismic source zones (Table A.1). The estimated seismic parameters and M_C values are compared with the previous studies for the same study area in Table 5. The observed lower bound of the b value is found to be similar to that of Thingbijam et al., (2008) and NDMA (2010), while the upper bound is comparable to values reported by Sharma & Malik (2006) and Das et al. (2016) as listed in Table 5. Lower b values are observed for the seismic source zone in the northeastern, northern, and eastern portions of the considered seismotectonic region (in increasing order) because these regions experience many higher magnitude earthquakes. On the other hand, a high value is primarily observed for seismic source zone 5, where the smaller size seismic events are more frequent (Figures 2b and Figure A.3).

Figure 7. Spatial distribution of PGA at the bedrock level for a return period of (a) 475 years, (b) 950 years, and (c) 2475 years.

Table 5. Comparison of the estimated seismic parameters with previously reported values.

Parameter	Das et al., 2016	NDMA, 2010	Thingbaizam et al., 2008	Sharma & Malik, 2006	Bahuguna & Sil, 2020	Present study
a	1.68–5.76	-	4.47–8.59	-	0.15–4.52	3.67–4.89
b	0.43–1.07	0.66 ± 0.03–0.80 ± 0.02	0.61–1.36	0.42–1.04	0.18–0.9	0.77–1.02
Mc	3.6–4.6	-	4.0–4.8	-	-	3.4–3.7

The PGA values at bedrock level (in g) for 10% probability of exceedance in 50 and 100 years and 2% probability of exceedance in 50 years were obtained from the derived seismic hazard curve at each site of interest. The obtained PGA values range from 0.14 to 0.69 g and 0.17 to 0.86 g for a 10% probability of exceedance in 50 years (return period of 475 and 950 years), respectively (Figure 7a and 7b) and from 0.22 to 0.93 g for 2% probability of exceedance in 50 years (Figure 7c). In all cases, higher PGA values are observed in the NER’s northern, northeastern, and western portions. In contrast, on moving south-southeast wards, the PGA values decrease gradually (Figure 7), which can be attributed to the influence of MFT, MCT, Dauki Fault, Mishmi, and Lohit Thrust in the region’s northern, central, and western portions.

The calculated PGA values for some important cities in the region are also compared with those from previous studies, as shown in Table 6. The calculated PGA values in the present study are relatively lower than those reported by Nath & Thingbaijam (2012) and are on the higher side compared to other studies (Table 6). However, in the case of Aizawl, Imphal, and Kohima, the PGA values for low probability of exceedance (return period = 475 years) are comparable with that of (NDMA, 2010) as mentioned in Table 6. In contrast, at Aizawl and Imphal, the calculated PGA values are less than that reported by Sharma & Malik (2006) and higher than that by Sil et al. (2013). At Guwahati and Shillong, the calculated PGA values are in a higher range than that reported by Ghione et al. (2021) and comparable with Bahuguna & Sil (2020) for a higher return period. Such variations can be attributed to the selection of different seismic source zones and ground motion attenuation models, which can be considered a limitation of the PSHA method. The predicted seismic hazard maps in the present study can act as an input in sustainable infrastructure planning, urban planning, and risk mitigation and disaster preparedness plans.

Table 6. Comparison of estimated PGA values obtained for important cities with previously reported values.

Sl. No.	City	PGA (g)
		Present Study^‡	NDMA, 2010 ^†	Nath & Thingbaijam, 2012 ^‡	Other studies
1	Guwahati	0.60–0.85	0.23–0.40	0.66–1.40	0.35 (Ghione et al., 2021) ^⁎; 0.46–0.92 (Bahuguna & Sil, 2020) ^†
2	Agartala	0.32–0.47	0.12–0.20	0.25–0.60	0.22 (Das et al., 2016) ^⁎; 0.11–0.20 (Sil et al., 2013) ^‡
3	Aizawl	0.20–0.32	0.22–0.45	0.45–1.20	0.3 (Sharma & Malik, 2006) ^⁎; 0.1–0.17 (Sil et al., 2013) ^‡
4	Imphal	0.28–0.42	0.30–0.55	0.70–1.40	0.4 (Sharma & Malik, 2006) ^⁎; 0.18–0.8 (Pallav et al., 2012) ^‡
5	Shillong	0.61–0.83	0.25–0.45	0.72–1.30	0.55 (Ghione et al., 2021) ^⁎
6	Itanagar	0.55–0.88	0.28–0.45	0.70–1.20	0.44 (Sharma & Malik, 2006) *; 0.18 (Das et al., 2016) *
7	Kohima	0.25–0.40	0.25–0.55	0.60–1.30	0.5 (Sharma and Malik, 2006)*; 0.15 (Das et al., 2016) *

For return period of * 475 years; ‡ 475–2475 years; † 475–4950 years

4.2. Social vulnerability assessment

Based on PCA, using Kaiser Criteria of factor retention, three components (factors) with eigenvalues greater than 1.0 were extracted that together explain ~90.5% of the total variance (Table 4). The selection of these three components was also confirmed by the scree plot, shown in Figure 5 (Sharma, 1996). Factor 1 comprising of 13 variables, reflecting housing condition (HC01 & HC02), lack of basic amenities (A01 & A02), literacy (L02 & L03), population composition and age (P02 & Age 01), quality of build environment (BM02& BM03), employment opportunities available (EO11 & EO12), and the number of households (HH02) as listed in Table 3 and 4. It imitates the relative living condition and socioeconomic status at the district level, explaining the highest variance (37.87% of the dataset) among all three factors. The variables such as HC01, A01, L02, and P02 have the highest loading.

The spatial distribution of SV in terms of Factor 1 is shown in Figure 8a. The districts of Assam and Tripura generally show a high to very high class of vulnerability corresponding to this factor. These districts have a high percentage of houses in dilapidated condition, with a lack of amenities reflective of poor living conditions, which may impact the community’s ability to cope with the impact of the disaster and the reason for high to the very high class of vulnerability. The study also finds that high illiteracy and dependent population further contribute to increased vulnerability in these districts, which is consistent with other SV studies (Dintwa et al., 2019; Frigerio et al., 2016)

Figure 8. Spatial distribution of social factors (a) Factor 1, (b) Factor 2, (c) Factor 3 and (d) SVI.

Factor 2 comprises six variables, governed by building material, aged population, and period of employment (Table 4). It explains 27.35% of the total variance in the dataset, the second-highest after Factor 1. The poor built quality, the duration, and the nature of employment opportunities available in a region are vital features governing the resilience capacity of the people of the region to the negative consequences of hazardous events. The spatial distribution of vulnerability in Factor 2, as shown in Figure 8b, reveals that these variables are relatively more dominant in districts of Tripura and Assam (mainly in upper Assam and Barak valley). The houses with poor built quality are more likely to get damaged in case of disaster and thus positively impacting the vulnerability.

Similarly, the aged population is dependent upon the younger generation. It may have mobility constraints, impacting the ability to quickly move away from the region impacted by the hazard in case of disaster and affecting the region’s vulnerability. All these factors enhance a region’s vulnerability and decrease the society’s resilience and coping capacity in a hazardous event.

Factor 3 represents the employment and percentage of the female population involved in agriculture and other related activities (Figure 8c; Table 4). The study area is relatively less urbanized and has agriculture and allied activities as a significant source of their income. Agriculture and small-scale household industries are low-paying jobs, and the female population of this region is found to be mainly involved in it. After the disaster, the non-permanent marginal workers having relatively low-paying jobs are more likely to lose their jobs due to disruption in daily activities and businesses (Morrow, 1999). The spatial distribution of Factor 3 shows that the districts of Assam, Meghalaya, Mizoram, and Manipur, having a high percentage of the population dependent on agriculture and small-scale industries, are under a highly vulnerable class.

The spatial distribution of SVI scores at the district level is also presented as a GIS-based choropleth map (Figure 8d). The presented SVI map will help to advance the understanding of socioeconomic vulnerability. The index scores were generated by combining the obtained component scores after the PCA using Equation 3(3) $SVI=\sum _{i=1}^{n_f}{w}_i\times \mathrm{F}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{i}$(3) and normalized using Equation 4(4) $SV{I}_j^N=\frac{\left(SV{I}_j-SV{I}_{min}\right)}{\left(SV{I}_{max}-SV{I}_{min}\right)}\times 100$(4) before entering into the GIS environment and classified using the quantile classification scheme in ArcGIS (Figure 8d). The adopted classification scheme and their areal coverage are listed in Table 7. Figure 8 reveals that the region’s factors and composite SVI scores are not evenly distributed. The most vulnerable districts are concentrated in the Brahmaputra flood plains, Tripura fold belt, and Imphal valley. Similar to the SV studies performed by Martins et al.(2012), Siagian et al. (2014), and Chakraborty et al. (2020) for their corresponding study areas, the finding of the present study also suggests that the characteristics like the quality of the built environment, illiteracy, access to amenities, dependent population, and employment opportunities are the major contributor to the relative level of SV of the districts of NER.

Table 7. Classification scheme of SHI, SVI, and Risk zone at a regional scale and corresponding areal distribution.

Value Range	Class	Area (%)
SHI
0.14–0.25	Very Low	18.31
0.26–0.36	Low	14.82
0.37–0.48	Moderate	7.75
0.49–0.59	High	20.48
0.60–0.70	Very High	38.64
SVI
0–0.013	Very Low	27.69
0.013–0.056	Low	20.50
0.056–0.165	Moderate	17.96
0.165–0.399	High	14.56
0.399–1.000	Very High	19.29
Risk Zone
01–02	Very Low	9.36
03–06	Low	40.40
07–10	Moderate	16.12
11–16	High	17.62
17–25	Very High	16.50

Overall ~19% area of the study area falls under the very high class of SVI, concentrated in Assam only (covering 15 out of 27 districtsas per Census of India, 2011), whereas 21 districts (12 of Assam, 2 of Manipur, and 4 of Tripura), covering ~15% area of the study area are identified with high class of SVI (Figure 8d). The states of Assam and Tripura reside the highest share of the population (~69.5% and ~8%, respectively) in NER (Census of India, 2011). Along with the relatively high population density, an average of ~60% population is non-working, a largely dependent population and ~85% of houses are made of weak building materials and poor quality, with a large share of the rural population in districts of these states (Census of India, 2011). Such districts also lack basic amenities like the availability of drinking water and electricity. These factors are responsible for the concentration of high to a very high class of SV in these states compared to the rest of NER. Considering the seismic risk, the proposed SVI can be used to locate socioeconomically most vulnerable communities using the maps and prepare an effective seismic risk reduction and mitigation policy for the region.

4.3. Cluster analysis

The spatial cluster analysis in Figure 9 represents hot spot and cold spot analysis. It is observed that hot spots with a 99% confidence level are found in central Assam and the lower part of Barak Valley, and hot spots with a 95% confidence level are observed in upper Barak valley and some parts of lower, northern, and upper Assam. A significant part of the lower Assam and northern part of Mizoram are located in the hot spot with a 90% confidence level. 17.14% of the total study area emerges as a hot spot with an average SVI score of 0.329. The cold spots with 95% confidence are mainly predominant in Nagaland, and cold spots with 90% confidence interval are observed in the northern part of Arunachal Pradesh. It is found that the spatial variation of cold spots spread over about 8.58% of the study area, with an average SVI score of 0.177. With an average SVI score of 0.208, the remaining areas are considered non-significant in terms of the hotspot and cold spot analysis, which are the northwestern and southern parts of the study area.

Figure 9. Spatial cluster analysis of SVI.

4.4. Seismic risk assessment considering SV at the regional level

Considering the SV, the seismic risk index at the regional level was obtained by integrating the SVI with the SHI and classified into five risk categories using the developed risk matrix (Table 7; Figures 6 and 10). The SHI map considering a 10% probability of exceedance in 50 years, was prepared. The spatial distribution of SHI shows that ~59% area of NER falls under the high to very high class (Figure 11), and these tracts are primarily concentrated in Arunachal Pradesh, Assam, and Meghalaya (Figure 10a).

Figure 10. (a) Seismic hazard class; (b) Risk zonation.

Figure 11. Comparison of areal distribution (%) of different classes of SV, SH, and Risk indexes in NER.

The presented risk map shows that SHI significantly influences the spatial distribution of seismic risk along with SVI. The upper reaches of the Himalayas in the northern portion of the study area (districts of Arunachal Pradesh) possess high to very high SHI and very low to moderate SVI and seismic risk index. The western portion of the study area, covering parts of Assam and Meghalaya, has moderate to very high SVI and SHI, and the seismic risk index ranges from high to very high (Figures 8d, 10aand 10b). For the central, southeast, and southern portion of the study area (central Assam, Manipur, Mizoram, Nagaland & Tripura), the risk index ranges from high to very low class, decreasing while moving towards the southeast, similar to the distribution of SHI in the region (Figure 10) even though for the same area, the SVI varies from very high (in the central part) to very low with a random distribution of moderate to low classes (Figure 8d). More than 50% of the total area is associated with either moderate, high, or very high-risk classes (Figure 11).

4.5. Distribution of seismic risk considering SV at sub-regional level

In the present study, the seismic risk was also conducted at the sub-regional level for all seven sister states of NER. State-wise seismic hazard classification was performed, and SVI mapping was done for each state (sub-regional level). The final risk zonation for each state was obtained by integrating an individual state’s SH class and SVI.

The SHI map considering a 10% probability of exceedance in 50 years was prepared for the seven states and classified into five classes very low (Class 1), low (Class 2), moderate (Class 3), high (Class 4), and very high (Class 5), shown in Figure 12. The PGA value for the Assam region ranges from 0.292 to 0.671 g, and 32.89% of the total area of Assam lies in a very high seismic hazard class (Figure 13). High PGA values are found in the western Assam region. From the seismotectonic map of the study area, it can be observed that lower Assam is situated in the proximity of significant faults like Dhubri, Oldham, Kulsi, Dhansiri-Kopili Fault, and thrusts like MBT, MFT, etc. (Bahuguna & Sil, 2020). Very low seismic hazard class is observed for the eastern part of Central Assam and some northeastern Assam regions. For Arunachal Pradesh, the PGA value ranges from 0.334 to 0.692 g, and approximately 51.15% of the total area of Arunachal Pradesh falls under the very high seismic hazard class (Figure 12). Among 16 districts of Arunachal Pradesh, four districts lie in the very high seismic hazard class. The presence of major thrusts like MCT, MFT, Lohit, and Mishmi Thrust increases the seismicity of these districts and the entire region of Arunachal Pradesh (S. Das et al., 2006). Meghalaya region experienced the 1897 Great Assam earthquake, responsible for the rise of the northern edge of the Shillong plateau (England & Bilham, 2015). In the present study, the PGA value of the Meghalaya region lies from 0.500 to 0.682 g, and the region of Garo and Khasi Hills falls, located in the western and southwestern part, under a very high seismic hazard zone covering about 65.14% of the total area of Meghalaya (Figures 12 and 13). Major tectonic features like Dauki, Dhubri Fault, Dapsi Thrust, Barapani shear zone, Umrangso, and NW-SE trending Kopili Faults contribute to its high seismicity (Baro et al., 2020).

Figure 12. Seismic hazard classification map at sub-regional level.

Figure 13. Comparison of areal distribution (%) of different seismic hazard classes at the sub-regional level.

For Mizoram, it is found that 64.47% of the total area of Mizoram lies in the very low to a low seismic hazard class; only 7.74% of the area is in the very high seismic hazard zone, and the range of PGA value varies from 0.144 to 0.3 g. The entire central and southern portion of Mizoram shows very low to a low seismic hazard class, and only a small portion of northern Mizoram falls under a very high seismic hazard zone. It lies near the Mat Fault and Tuipui fault (Sitharam & Sil, 2014). The lowest and highest PGA values observed for Nagaland are 0  and 0.308 g, respectively, and about 59.85% of the total area of Nagaland falls under a moderate to very high seismic hazard zone (Figure 13). The Naga and Disang Thrust are major seismotectonic features in this region, and comparatively high PGA values are observed in the northern part of Nagaland. Approximately 40% of the total area of Tripura falls under a very high to high seismic hazard class with PGA values ranging from 0.222 to 0.351 g, mainly concentrated in the northern part of the Tripura region showing the influence of the Sylhet fault and Arakan Trench (Sil et al., 2013).

Similarly, for Manipur, the overall PGA value remains low compared to other states, varying from 0.172 to 0.308 g. Some parts of western Manipur fall under a very high seismic hazard zone. It is found that about 44.93% of the total area is characterized by moderate to a very high seismic hazard class and 55.07% of the total area of Manipur is in a very low to low seismic hazard zone, mainly in the eastern part (Figures 12 and 13). Many active faults like the Dauki and Kopili Fault surrounding the Manipur region contribute to the region’s seismicity (Pallav et al., 2012).

The present study calculated SVI for the seven states and classified it into five classes (Figure 14). For the Assam region, it is observed that 46.71% of the total area falls under the high to very high SVI class, and it constitutes central, southern, northeastern, and two districts of western Assam. In Arunachal Pradesh, more than 75% of the total area of Arunachal Pradesh lies under the very low to moderate SVI class (Figure 15), and high to very high SVI is observed for the southeast part of the state. Most of Meghalaya, approximately 67.97% of the total area of Meghalaya, is classified in the very low to moderate SVI class. The West Garo Hills in the west and Jaintia Hills in the east of Meghalaya state are located under a high to very high SVI class (Figure 14). The northern and southwestern part of Mizoram state falls under the high to very high SVI class, constituting approximately 47.17% of the total area of the Mizoram region. The SVI classification shows that western and eastern Mizoram is under a very low to low SVI class. The northern and eastern parts of Nagaland, comprising five districts, showed high SVI values, covering 45.96% of the total area of Nagaland. Very high SVI values can be observed for the southwestern region of Nagaland, and it is observed that 79.06% of the total area of Nagaland state can be classified into moderate to very high SVI classes (Figures 14 and 15). Out of four districts of Tripura high to very high SVI range is observed for two districts covering 48.28% of the total Tripura state area. More than 60% of the Manipur area is characterized by moderate to very high SVI values covering eight districts, and only 38.79% of the total area falls under the very low to low SVI class. From the present study, it is observed that the area that comes under high to very high SVI class is more in states like Mizoram, Nagaland, and Tripura as compared with the rest of the states of NER (Figures 14 and 15). Higher SVI values are mainly associated with high population density, low employment opportunities, lower literacy rates, poor housing conditions, etc. (Fatemi et al., 2017).

Figure 14. Distribution of SVI classes at the sub-regional level.

Figure 15. Comparison of areal distribution (%) of different SVI classes at the sub-regional level.

A risk zonation map for each of the seven states was prepared by integrating the state-wise seismic hazard class and SVI map and classified into five risk zone (Figure 16). High to very high-risk class is observed for western, central, northern, and southern parts of Assam and comprises 35.53% of the total area of Assam (Figure 17). The northeast of Assam and the southern part of Central Assam can be classified as a very low to moderate risk class (Figure 16). For Arunachal Pradesh, four classes of risk zone are identified, and 65.60% of the state’s total area covers moderate to very high-risk classes (Figure 17), mainly concentrated in the western, southern, and southeastern parts of Arunachal Pradesh. A significant area of Meghalaya state, 85.25% of the total area of Meghalaya, can be classified under moderate to very high-risk zone. It mainly includes the state’s central, southeastern, and western regions (Figures 16 and 17).

Figure 16. Risk zonation map at the sub-regional level.

Figure 17. Comparison of areal distribution (%) of different risk classes at the sub-regional level.

According to the risk zonation map of Mizoram, the eastern and southern part of the state shows a very low to low-risk class, whereas only a small portion of the northern region falls under the very high-risk zone. About 60.50% of the total area of Mizoram comes under the very low to low-risk class, and 13.04% falls under high to very high-risk zonation (Figure 17). The northern and southwestern parts of Nagaland state constitute a high to very high seismic risk zone, covering about 51.47% of the state’s total area (Figures 16 and 17). The western and entire southeastern part of Nagaland shows very low to moderate seismic risk covering 48.53% of the total area of Nagaland. However, among all the seven states, the highest proportion of area, 51.47%, falls under Nagaland’s high to very high-risk class. For Tripura, a major part of the state is considered under the very low to low-risk class, including the two districts (Figure 16). It is found that 37.50% of the total area falls under the high to very high-risk class and the very low to moderate risk class covers 62.49% area (Figure 17). Four risk zones are obtained for Manipur state ranging from low to very high, and 81.94% of the total area is characterized by a low to moderate risk zone (Figure 17). The low-risk zone areas are identified in the western, northeastern, and southeastern parts of Manipur, and a small area of the very high-risk zone is present in the extreme western part of the state. Risk zonation depends on the spatial variation of the study area’s seismic hazard class and SVI; for example, in Arunachal Pradesh, the minimal area has a very low seismic hazard class, but more than 40% of the area has low to very low SVI (Figures 13, 15 and 17). As a result, approximately 34% of the area falls under the low-risk class, presented in Figures 16 and 17 (Frigerio et al., 2016). The detailed classification of SH, SVI, and risk zone is given in Table 8.

Table 8. Classification scheme of SHI, SVI, and Risk zone and corresponding areal distribution.

		SHI		SVI		Risk Zone
State	Class	Value Range	Area (%)	Value Range	Area (%)	Value Range	Area (%)
Arunachal Pradesh	Very Low	0.334–0.405	1.61	0–0.196	27.23	1.0–2.0	0.00
	Low	0.406–0.477	3.08	0.196–0.407	15.42	3.0–6.0	34.40
	Moderate	0.478–0.549	17.59	0.407–0.701	34.32	7.0–10.0	16.52
	High	0.550–0.621	26.58	0.701–0.949	10.34	11.0–16.0	33.97
	Very High	0.622–0.692	51.15	0.949–1	12.69	17.0–25.0	15.11
Assam	Very Low	0.292–0.368	14.00	0–0.258	16.57	1.0–2.0	5.09
	Low	0.369–0.444	16.45	0.258–0.364	24.73	3.0–6.0	35.85
	Moderate	0.445–0.520	18.01	0.364–0.521	11.97	7.0–10.0	23.52
	High	0.521–0.596	18.65	0.521–0.709	34.51	11.0–16.0	22.22
	Very High	0.597–0.671	32.89	0.709–1	12.21	17.0–25.0	13.31
Manipur	Very Low	0.172–0.198	26.96	0–0.027	18.80	1.0–2.0	0.00
	Low	0.199–0.224	28.11	0.027–0.121	19.96	3.0–6.0	60.48
	Moderate	0.225–0.250	21.48	0.121–0.141	21.14	7.0–10.0	21.47
	High	0.251–0.277	13.08	0.141–0.200	29.89	11.0–16.0	17.34
	Very High	0.278–0.308	10.37	0.200–1	10.22	17.0–25.0	0.71
Meghalaya	Very low	0.500–0.536	3.67	0–0.001	8.22	1.0–2.0	0.00
	Low	0.537–0.573	9.68	0.001–0.545	23.36	3.0–6.0	17.95
	Moderate	0.574–0.609	9.50	0.545–0.788	36.39	7.0–10.0	32.23
	High	0.610–0.645	12.01	0.788–0.886	16.90	11.0–16.0	32.19
	Very High	0.646–0.682	65.14	0.886–1	15.13	17.0–25.0	17.62
Mizoram	Very Low	0.144–0.175	36.95	0–0.133	6.43	1.0–2.0	9.86
	Low	0.176–0.207	27.52	0.133–0.490	37.05	3.0–6.0	50.64
	Moderate	0.208–0.238	15.82	0.490–0.639	9.35	7.0–10.0	26.09
	High	0.239–0.270	11.96	0.639–0.854	30.83	11.0–16.0	7.90
	Very High	0.271–0.300	7.75	0.854–1	16.34	17.0–25.0	5.50
Nagaland	Very Low	0.203–0.224	19.03	0–0.027	11.51	1.0–2.0	13.18
	Low	0.225–0.245	21.12	0.027–0.074	9.48	3.0–6.0	18.51
	Moderate	0.246–0.266	23.11	0.074–0.086	9.80	7.0–10.0	16.83
	High	0.267–0.287	23.18	0.086–0.407	45.97	11.0–16.0	27.64
	Very High	0.288–0.308	13.56	0.407–1	23.25	17.0–25.0	23.83
Tripura	Very Low	0.222–0.248	13.22	0–0.309	21.20	1.0–2.0	26.64
	Low	0.249–0.273	23.88	0.309–0.310	12.98	3.0–6.0	26.84
	Moderate	0.274–0.299	23.24	0.310–0.311	17.53	7.0–10.0	9.02
	High	0.300–0.325	22.73	0.311–0.944	28.62	11.0–16.0	30.73
	Very High	0.326–0.351	16.93	0.944–1	19.67	17.0–25.0	6.78

The significance of the spatial distribution of seismic hazard class, SVI, and risk zonation mapping at the sub-regional level is that it gives a clear understanding of the areas that require immediate and most attention in terms of mitigation and management of hazards at the district level. From the risk zonation at a regional level, it appears that the entire state falls under the high-risk zone for Meghalaya, which shows the uniform spatial distribution of seismic hazard class and SVI. However, this is not the case in general, and it became challenging to prioritize and differentiate the area district-wise for making disaster management policies. So, at sub-regional risk zonation mapping of Meghalaya shows that out of seven districts as per Census of India (2011), one district is under very high-risk class, one is under high, one is for low, one district shows moderate to high, one low to high-risk class, one low to moderate risk and one district shows low to very high-risk zone (Figures 10a and 16). Similarly, the other states of NER show more variation at the sub-regional level than the regional level of a seismic hazard class, SVI, and risk zonation mapping. From the sub-regional mapping, areas of critical importance and the extent of the affected area can be obtained so that the allocation of resources can be done effectively.

5. Conclusion

The present study conducted a seismic risk assessment considering PSHA and SV at regional and sub-regional levels for NER. The steps of PSHA include the preparation of an updated earthquake catalog and homogenization, delineating the study area into seismic source zones, and estimating PGA values using existing GMPEs. Updated seismic hazard maps for the return periods of 475, 950, and 2475 years were generated. The seismic hazard map for the return period of 475 years is classified into five classes to generate seismic hazard classes for regional and sub-regional levels. It is found that 59.11% of the total area of the NER comes under the high to very high seismic hazard class. The states of Assam, Arunachal Pradesh, and Meghalaya exhibit relatively higher PGA values. At the sub-regional level, it is observed that 77% of the total area of Meghalaya comes under high to very high seismic hazard zones. The spatial variation of seismic hazard class correlates well with the surrounding tectonic framework of the study area at regional and sub-regional levels. Socioeconomic indicators were considered for the assessment of SV, and based on PCA, three factors comprising 24 socioeconomic variables were extracted and used to generate an SVI map. The SVI map was further classified into five SVI classes, showing that 66.14% of the total study area lies under the very low to moderate SVI class, whereas 33.86% of the area contributes to the high to very high SVI class. Assam and Tripura show high to very high SVI values, and very low to low SVI classes were observed for Arunachal Pradesh and Nagaland. The spatial cluster analysis shows high social vulnerability patterns in Assam.

From the obtained seismic risk zonation map at the regional level, about 50.23% of the total study area shows moderate to very high-risk zone consisting of mainly Assam, Meghalaya, and some parts of Arunachal Pradesh and Tripura. At the sub-regional level, about 85.25% of the total area of Meghalaya can be grouped under moderate to very high-seismic risk zone, including the central, southeastern, and western regions of the state. Among all the seven states, Nagaland shows the highest proportion of area,that is 51.47%, falling under the high to very high-seismic risk class.

The produced risk maps are easy to interpret and essential for urban planners and decision-makers for the identification of the most seismically vulnerable areas . This study can also aid in the design of site-specific resilient infrastructure. It can assist the state authorities and local bodies in preparation for disaster risk mitigation planning and effective resource allocation in the areas of critical importance.

Data Availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request.

CRediT authorship contribution statement

Conceptualization: JD; Methodology: NA, LG, JD; Formal analysis and investigation: NA, LG; Validation: NA, LG, JD; Visualization: NA, LG; Writing - original draft preparation: NA, LG; Writing - review and editing: JD, SKD; Resources: JD; Supervision: JD.

Disclosure Statement

The Coalition for Disaster Resilient Infrastructure (CDRI) reviewed the anonymised abstract of the article, but had no role in the peer review process nor the final editorial decision.

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

The Article Publishing Charge (APC) for this article is funded by the Coalition for Disaster Resilient Infrastructure (CDRI). The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Appendix A

A.1 Homogenization and declustering of earthquake catalogue

The obtained earthquake catalog comprised events in different magnitude scales, like body-wave magnitude (m_b), surface-wave magnitude (M_S), local magnitude (M_L), and moment magnitude scale (Mw). Since all magnitude scales except Mw saturates for higher magnitude, the catalogue was homogenized to a standard magnitude scale (i.e., Mw) by developing correlations between different magnitude scales and Mw using an orthogonal regression approach (Wason et al., 2012). There are 349, 277, and 130 data points for the orthogonal regression between m_b – Mw, M_L – Mw, and M_S – Mw, respectively (Figure 18).

Figure 18. Conversion from (a) m_b to Mw, (b) M_L to Mw, and (c) M_S to Mw.

For PSHA, it is essential to have a mutually exclusive set of events. Therefore, the homogenized earthquake catalog was further processed to eliminate the duplicate and dependent events from the dataset. For this purpose, the declustering of the earthquake catalog was performed by adopting the classical window technique of Gardner & Knopoff (1974) using an open-source software ZMAP (v 7.0) (Wiemer, 2001; Anbazhagan et al., 2019; Sitharam & Sil, 2014). After declustering, the dataset contains 6599 events within the considered seismotectonic region, among which 4837 events are of magnitude larger than Mw 3.5.

A.2 Completeness of catalogue

The completeness of the data in terms of magnitude and time must be examined before estimating seismic parameters for PSHA. The magnitude of completeness (M_C) is the lowest magnitude above which the catalogue, in a selected space-time window, is considered to be complete (Rydelek & Sacks, 1989; Wiemer & Wyss, 2000). In the present study, the M_C values of each considered seismic zone are obtained through (Wiemer & Wyss, 2000) maximum curvature method (MAXC). Various scholars worldwide have adopted a similar procedure (Anbazhagan et al., 2019; Baro & Kumar, 2017; Borah et al., 2021; Gupta et al., 2021; Woessner & Wiemer, 2005). The open-source software ZMAP (v 7.0) by Wiemer (2001) was used for the calculation. M_C values range from 3.4 (zone 3) to 3.7 (zone 5) for the considered study area, as listed in Table A.1.

Further, the completeness of the homogenized-declustered earthquake dataset with respect to time was also checked using the statistical analysis proposed by Stepp (1972), as shown in Figure A.2. To do so, the number of earthquakes per decade is classified into six magnitude ranges, namely, (i) 3≤ Mw < 4; (ii) 4≤ Mw < 5; (iii) 5≤ Mw < 6; (iv) 6≤ Mw < 7; (v) 7≤ Mw< 8; and (vi) Mw ≥ 8. In this method, the earthquakes were assumed to follow the Poisson distribution model, and the variance of the sample mean for a definite unit time interval is obtained using the Equation A.1.

(A.1) ${\sigma }_{\lambda }^2=\left(\lambda /n\right)$(A.1)

Where λ is an unbiased mean rate for each definite unit time interval of the sample, and n is the number of unit time intervals equal to sample length (T) in years (Stepp, 1972). Then, the standard deviation (σ_λ) can be obtained by Equation A.2

(A.2) ${\sigma }_{\lambda }=\sqrt{\lambda /T}$(A.2)

As per Stepp (1972), σ_λ will be equal to 1/√T for each time interval if λ is constant. Therefore, σ_λ for each magnitude class is plotted corresponding to each time interval (Figure A.2). The dataset of a given magnitude interval is recognized as complete to the extent that plotted data points are parallel to line 1/√T for that duration. According to Figure A.2, the duration of completeness is 60, 70, 100, 120, 150, and 260 years for Mw 3.0–4.0, 4.0–5.0, 5.0–6.0, 6.0–7.0, 7.0–8.0, and ≥8.0, respectively.

Figure 19. Completeness of the earthquake catalog with time.

A.3 Estimation of seismic parameters

The seismicity of a region can be described by seismic parameters a and b, which correlate with the rate of occurrence of an event of a particular size. The distribution of event sizes in a given period is best described by a most widely accepted G-R recurrence law (Kramer, 1996), as given in Equation A.3.

(A.3) $log(N)=a-b(\mathrm{M}\mathrm{w})$(A.3)

Where N represents the number of cumulative events per year greater than an event of a given magnitude, a and b are constants of regression, known as seismic parameters.

Based on the completeness study, the earthquake catalogue of the recent 70 years was considered to evaluate the recurrence relation for each source zone (Figure A.3a-e). The total number of earthquakes above the magnitude of completeness are 566, 431, 722, 2247, and 871 in the source zones 1, 2, 3, 4, and 5, respectively. The obtained values of seismic parameters a and b are summarized in Table A1.

Table A1. Seismic parameters and M_C values.

Source zones	MC (using MAXC)	Seismic parameters		R^2
		a	b
1	3.60	4.06 ± 0.36	0.84 ± 0.07	0.96
2	3.50	3.67 ± 0.36	0.77 ± 0.07	0.95
3	3.40	4.23 ± 0.33	0.85 ± 0.06	0.97
4	3.50	4.70 ± 0.19	0.86 ± 0.03	0.98
5	3.70	4.90 ± 0.37	1.02 ± 0.07	0.97

Figure 20. Gutenberg-Richter relation for each source zone (a) Zone 1, (b) Zone 2, (c) Zone 3, (d) Zone 4, and (e) Zone 5.

A.4 Estimation of maximum magnitude (M_max)

The largest possible earthquake, M_max, that a seismic source can produce is an essential input parameter for PSHA as ground motion intensity is influenced by earthquake magnitude. In this study, M_max was evaluated based on the conventional incremental value method, IVM (Anbazhagan et al., 2019; Gupta, 2002), and the procedure suggested by Kijko and Sellevoll (1989), typically referred to as KS89. It is based on the doubly truncated G–R relation given by Equation A.4 and can be applied only when the region’s seismic parameter ‘b’ is known (Kijko, 2004).

(A.4) $M_{max}=m_{max}^{obs}+\mathrm{\Delta };\text{break}Where\mathrm{\Delta }=\frac{E_1(n_2)-E_1(n_1)}{\beta exp(-n_2)}+m_{min}exp(-n)$(A.4)

M_max is the calculated maximum magnitude, $m_{max}^{obs}$ is the observed maximum magnitude associated with each fault, n is the number of events above M_C in the region, and $m_{min}$ denotes the minimum magnitude.

Based on the M_C value, m_min in the present study was taken as 3.5. $n_1=\frac{n}{1-exp(-\beta (m_{max}-m_{min})}$, $n_2=n_{1}exp[-\beta (m_{max}-m_{min})]$, and E₁(n_i) is an exponential integral function, which can be approximated as $E_1(n_i)=\frac{n_i^2+a_1{n}_i+a_2}{n_i(n_i^2+b_1{n}_i+b_2)}exp(-n_i)$, where a₁ = 2.334733, a₂ = 0.250621, b₁ = 3.330657 and b₂ = 1.681534 (Abramowitz and Stegun 1970).

On the other hand, in IVM, which is relatively simple and applied by many researchers, M_max is obtained by adding a constant value of 0.5 to $m_{max}^{obs}$ value of each seismic source (Anbazhagan et al., 2019; Bahuguna & Sil, 2020; Gupta, 2002). Values of M_max calculated by both methods are given in Table A2.

Table A2. Seismic source zone characterization of the study area.

Source zones	Fault name	Fault ID	αi	χi	Observed Mmax	Calculated Mmax
						IVM	KS89
1	Main Central Thrust (MCT)	1	0.49	0.59	6.8	7.3	6.86
	Main Boundary Thrust (MBT)	2	0.51	0.41	6.8	7.3	6.86
2	Siang Fault (SiF)	3	0.18	0.16	6.6	7.1	6.63
	Lohit Thrust (LT)	4	0.20	0.40	6.3	6.8	6.32
	Mishmi Thrust (MT)	5	0.62	0.43	7	7.5	7.06
3	Shan-Sagaing Fault (SSF)	6	1.00	1.00	7.6	8.1	7.96
4	Eastern boundary thrust and Kabaw Fault (EBT & KBF)	7	0.42	0.71	7.3	7.8	7.49
	Chaurachandpur-Mao Fault (CMF)	8	0.09	0.15	7.2	7.7	7.35
	Naga Thrust (NT)	9	0.25	0.06	7.3	7.8	7.49
	Disang Thrust (DT)	10	0.24	0.08	7	7.5	7.10
5	Tista Fault (TF)	11	0.15	0.02	4.5	5	4.50
	Dhubri Fault (DF)	12	0.09	0.07	6.7	7.2	6.87
	Dhansiri-Kopili Fault (DKF)	13	0.08	0.11	5.7	6.2	5.72
	Samin Fault (SF)	14	0.05	0.13	7.1	7.6	7.64
	Oldham Fault (OF)	15	0.07	0.13	8.1	8.6	8.55
	Dhudhnoi Fault (DhF)	16	0.06	0.11	8.1	8.6	8.54
	Barapani Shear Zone (BSZ)	17	0.03	0.06	5.2	5.7	5.20
	Dauki Fault (DaF)	18	0.18	0.13	7.1	7.6	7.64
	Madhupur Blind Fault (MBF)	19	0.04	0.03	7.1	7.6	7.64
	Sylhet Fault (ShF)	20	0.09	0.11	6.5	7	6.60
	Arakan Trench (AT)	21	0.14	0.09	6.8	7.3	7.01

A.5 Fault deaggregation

Estimating the seismic activity rate of an individual fault is a significant step in conducting PSHA of any region. For this purpose, an approach similar to Raghukanth and Iyengar (2006) and NDMA (2010) was adopted in this study. A conservation property was heuristically used to develop recurrence relations for different seismic regions of the study area. The number of earthquakes per year with Mw≥ m_min, i.e., N(m_min = 3.5) in a region, was calculated from the G–R relation of that region using Equation A.3. Since all these events are associated with the faults within the seismo-tectonic region, it should be equal to the sum of the number of earthquakes occurring on individual faults, i.e. $=\sum _{i=1}^n{N}_i(m_{min})$, where N_i(m_min) is the annual frequency of events of Mw ≥ m_min on the i^th fault in the region, (i = 1,2,3, …,n). N_i(m_min) of any fault depends on the fault length and past seismic activity of that fault. The evaluation of N_i(m_min) involves two basic assumptions: (i) longer faults will have a higher capacity to rupture into smaller segments, and (ii) shorter faults may be more active in producing relatively smaller events. Correspondingly, the recurrence relation for i^th source is obtained using the following Equation A.5

(A.5) $N_i(m_{min})=0.5({\alpha }_i+{\chi }_i)\cdot N(m_{min})$(A.5)

Where ${\alpha }_i=L_i/\sum L_i$ is the weighing factor for length of i^th fault (L_i),${\chi }_i$ is another weighting factor defined as the ratio of the number of EQs associated with i^th fault to the total number of EQs in the region. In this study, 21 active seismic sources were identified and the detail of each fault in terms of the number of EQs associated with it, its length, weighing factors, and the evaluated M_max, are given in Table A.2. The b-value of each fault is considered equal to that of the associated seismic zone, and the annual rate of events of magnitude ≥ Mw (λ_m) is obtained using Equation A.6.

(A.6) $N_i(m)=N_i(m_{min})\left[1-\frac{1-e^{\{-\beta (m-m_{min})\}}}{1-e^{\{-\beta (m_{max}-m_{min})\}}}\right];\text{break}\mathrm{f}\mathrm{o}\mathrm{r}{m}_{min}\le m\le m_{max}$(A.6)

Where, m_min is the minimum threshold magnitude, m_max is the maximum potential magnitude of the fault i, and β = 2.303b. The individual fault-level recurrence relations are shown in Figure 21a-214.

Figure 21. Fault level recurrence relation for source (a) Zone 1, (b) Zone 2, (c) Zone 3, (d) Zone 4, and (e) Zone.