Earthquake vulnerability in the Himalaya by integrated multi-criteria decision models

Himalayan mountains are one of the most seismo-tectonically active zones on the surface of the earth. Recurring moderate and high magnitude earthquakes are not uncommon in this region. This paper maps the earthquake vulnerability in the region using integrated multi-criteria decision models. Factors which may influence vulnerability in a region can be categorized in social, geotechnical, structural, and physical parameters. We have used the analytical hierarchy process (AHP) approach to determine the weights of various parameters, which were further used to develop earthquake vulnerability maps for the study area using the VIseKriterijumska Optimizacija I Kompromisno Resenje (VIKOR), and Grey Relational Analysis (GRA) methods. There is good correlation between the vulnerability estimated through AHP-VIKOR and AHP-GRA methods. Our analysis indicates that more than 12% area may be under high to very-high vulnerability, whereas more than 44% population, and about 43% buildings are highly vulnerable to hazards due to earthquakes. The results would be useful for various hazard mitigation and infrastructure planning agencies working in the region.


Introduction
The Himalayas are one of the world's most seismo-tectonically active region on the earth. Continuous collision of the Indian and the Eurasian plates result in accumulation of the strain energy (Dewey et al. 1989), which is frequently released in the form of earthquakes. The extreme destructive power of earthquakes endangers life and properties (Xiwei 2010). The available historical database (Rahman et al. 2017) indicates that many catastrophic earthquakes have occurred in this region (Table 1). These destructive earthquakes cause substantial casualties and financial loss in the affected neighboring countries. One of the most recent and devastating earthquakes, i.e., the Gorkha earthquake (M ~ 7.8), 1 3 occurred in the central Himalaya on 25th April 2015. This earthquake caused approximately 9000 casualties, about 22,000 injured, destroyed thousands of buildings and caused substantial economic losses (Bilham 2015).
Although there have been significant developments in understanding earthquakes and related processes in the recent past, it is nearly impossible to predict earthquakes in space and time with any practical use. Lack of implementing the Government's standard infrastructure development codes often results in poor quality infrastructure. The rapid urbanization with improper management and poor constructions, besides high population density in cities, poses a significant risk due to earthquakes, leaving the modern cities more Table 1 Major earthquakes (historical and recent) in the Himalayan seismo-tectonic zone ( taken from USGS; NOAA; Bhatia et al. 1999;Zhang et al. 1999;Ambraseys and Douglas 2004;Szeliga et al. 2010) . Date  vulnerable and exposed to such severe natural disasters (Asadi et al. 2019). Identifying vulnerable areas, developing good quality infrastructure, and early warning systems may reduce earthquakes' primary and secondary impacts on society. Due to better planning, strict implementation of building codes, and the rule of law, loss of life and casualties are relatively lower. In contrast, financial losses are usually higher in developed countries, whereas the situation is reversed in developing countries (Ebert and Kerle 2008). Furthermore, the underprivileged sections of society may be more vulnerable because of the poor quality of their houses, which may not withstand even moderate earthquakes (Schilderman 2004). The women population may be even more vulnerable than men because of cultural practices and social norms (Ruddock 2007). With rapid population growth and urbanization, it is essential to study urban and rural areas' structural and social vulnerability. The development of good health care facilities with a better communication network may help in minimizing the casualties in pre-and post-earthquake situations (Merciu et al. 2018).
One of the main reasons for the lack of availability of earthquake vulnerability maps may be the limited availability of datasets (Rashed and Weeks, 2003). However, a variety of data is now available in the public domain freely, which allows us to evaluate earthquake vulnerability based on various parameters. Recent studies have used MCDM methods such as Analytic Network Process (ANP), and Analytical Hierarchy Process (AHP) approach, and the geographical information system (GIS) (Alizadeh et al. 2018a) to study earthquake risk and vulnerability. An integrated approach using Artificial Neural Network (ANN) and AHP was used to develop the earthquake risk assessment map for the Banda Aceh (Jena et al. 2019). The VIseKriterijumska Optimizacija I Kompromisno Resenje (VIKOR) method was used to estimate seismic vulnerability (Jena et al. 2020), whereas the AHP and the Technique for Order Preference by Similarity to Ideal Solution (TOPSIS) were used to develop earthquake hazard and risk maps of Küçükçekmece, Istanbul (Nyimbili et al. 2018). AHP and GIS were used to assess the seismic vulnerability of the residential houses in urban areas of Tabriz city (Alizadeh et al. 2018a) and for the school buildings in Tehran (Panahi et al. 2014). Chakraborty and Joshi (2016) used AHP to map a multi-disaster scenario, whereas the AHP and GIS were used to map the socio-economic and structural vulnerability and risk index for Kolkata (Nath et al. 2015). Duzgun et al. (2011) used the simple additive weighting (SAW) method to assess the earthquake vulnerability for Eskisehir, Turkey. In contrast, a combination of physical, social, and systemic parameters was used to evaluate earthquake vulnerability in Victoria, British Columbia, using AHP (Walker et al. 2014). Karaman and Erden (2014) created high-resolution earthquake hazard maps using an AHP and GIS for Istanbul. For the Himalayan region, several studies have been done to estimate the landslide and flood susceptibility mapping; however, not much has been done to map the seismic vulnerability except the earthquake hazard model for Sikkim Himalaya done through fuzzy logic AHP (Pal et al. 2008).
Accurate mapping of earthquake vulnerability is always challenging, particularly on a large scale. Considering the importance of the risks related to impacts of the earthquake on neighboring densely populated areas, we apply an integrated MCDM model to assess the earthquake vulnerability in the Himalayan seismotectonic zone. The historical seismicity, geological features, lithology of the area, slope, population density, quality of the infrastructure, location of health facilities, and availability of good communication networks are few of the main factors often considered for mapping the earthquake risk and vulnerability. These primary factors can be broadly classified into social, geotechnical, structural, and physical parameters, making this study suitable for applying multi-criteria MCDM models. The MCDM methods depend on data, which have their associated uncertainties. Therefore, we used the VIKOR and GRA methods to optimize these uncertainties and provide relatively robust results. To prepare an earthquake vulnerability map for the study region, two MCDM models have been integrated, i.e., the analytical hierarchy process (AHP), VIseKriterijumska Optimizacija I Kompromisno Resenje (VIKOR) technique (AHP-VIKOR), and the analytical hierarchy process (AHP), Grey Relational Analysis (GRA) technique (AHP-GRA). The AHP method evaluates the priority and criterion ranking and relies on knowledge-based available literature and the experts' opinion. The VIKOR technique helps to determine the compromise ranking and comprises the solution acquired through the initial weights. The Grey Relational Analysis (GRA), integrated with AHP, evaluates the complicated uncertainty among the multi-criteria in a given system and optimizes it. We estimated the weights of the variables using AHP, whereas VIKOR and GRA provide the ranking of the criteria/alternative. The methods have been integrated with GIS for better visualization of the results.
In this study, we have attempted to address questions of societal benefit by applying the data-driven integrated MCDM models. We have proposed a novel method by coupling two approaches discussed earlier, previously not used to address such problems. Furthermore, we evaluate our new approach of assessing earthquake vulnerability with other popular MCDM models. In summary, the main objective of the paper is to assess the earthquake vulnerability for the Himalayan region by applying a novel integrated MCDM method for large-scale regions by applying a number of publicly available datasets that directly or indirectly impact the assessment of the earthquake vulnerability.

Study area
The study area covers the Himalayan mountains extending from northern to the north-eastern border of India, is probably one of the most seismo-tectonically active features. The study area is spread in several countries such as India, China, Bangladesh, Myanmar, Pakistan, and Bhutan, with about 976661km 2 , inhabited by more than 90 million people living in the vicinity (Fig. 1). Therefore, it is essential to map the earthquake vulnerability to mitigate future potential events. The region's prominent tectonic feature is the Himalayan Frontal Arch (HFA), stretching in the NW to SE direction and dominated by several earthquakes that occur in this region. Other prominent tectonic features include the Main Central Thrust (MCT), Indus Tsangpo Suture (ITS), and the Main Boundary Thrust (MBT), which extend throughout the Himalayan region. Active collision and subduction of Indian and Eurasian plates have resulted in considerable accumulation of strain energy, which is released in the form of moderate to large earthquakes from time to time (Dewey et al. 1989). The region falls under seismic zones IV and V in India's seismic hazard zonation map due to high seismic activities in this region (BIS 2002). The expected peak ground acceleration (PGA) for the Himalayan region is between 0.10 and 0.40 g, with a 10% probability of exceedance in 50 years .
Most of the earthquakes in this region are shallow earthquakes with focal depths less than 30 km; however, a few deep earthquakes have also been reported, indicating breaking of the subducting Indian plate in the mantle. Furthermore, the stress released along the Himalayan region subduction zone is not homogeneous in space and time (Bilham 2015), making this practically impossible to predict the earthquakes with any practical use. Numerous earthquakes of magnitude varying from moderate to large occurred in the last 200 years. The recent historical earthquakes, such as the 1934 Bihar-Nepal earthquake 1 3 (Mw ~ 8.4) and the 1950 Assam-Tibet (Mw ~ 8.6) earthquake, along with the 2015 Gorkha earthquake (Mw ~ 7.8), proved highly devastating for the region. The Gorkha earthquake of 2015, which ruptured a segment of the Main Himalayan Thrust fault (MHT) (Rupakhety 2018), resulted in 8790 deaths and 22,300 injured persons (NPC 2015). These earthquakes pose a significant threat to the densely populated northern India. Therefore, it is important to understand the associated hazard and risks in the region. In the following sections, we discuss the data, methodology, and results obtained in this study.

Data and methodology
This study uses several datasets available in the public domain including geological data, historical and recent seismicity, population, infrastructure, social, and healthcare facilities data for estimating earthquake vulnerability in the study region. The datasets have been compiled from various sources available in the public domain. Since the study area is distributed in several countries, we have collected data from multiple international agencies and homogenized them. The social data, population, and literacy rate datasets have been collected from the Office of the Registrar General & Census Commissioner, India (https:// censu sindia. gov. in/), the Bangladesh Bureau of Statistics (http:// www. bbs. gov. bd/), the China National Bureau of Statistics (http:// www. stats. gov. cn/ engli sh/), the National Statistics Bureau, the Royal Government of Bhutan (http:// www. nsb. gov. bt/ main/ main. php), the Myanmar Population and Housing Census (http:// themi mu. info/ census-data), the Central Bureau of Statistics Nepal (https:// cbs. gov. np/), and the Pakistan Bureau of Statistics (http:// www. pbs. gov. pk/). The structural and physical characteristics datasets were obtained from USGS and OpenStreetMap (https:// www. opens treet map. org/), whereas the availability of road and railways transport networks were collected from DIVA-GIS (https:// www. diva-gis. org/). The earthquake magnitude and epicenter data were obtained from NOAA (https:// www. noaa. gov/), whereas the dataset for active faults is taken from the global earthquake model (https:// www. globa lquak emodel. org/). The SRTM elevation cen. uni-hambu rg. de/) were also used. The datasets were last updated on 27th November 2020. Kernel density and Euclidean distance were estimated for the layers and the distances. The PGA was evaluated using the equation developed by Panjamani et al. 2016 for the Himalayan region. The methodology is discussed in detail in the following sections.

Multi-criteria decision-making model
Multiple criteria decision-making (MCDM) is a mathematical tool for supporting the subjective assessment of decision-makers (Zavadskas and Turskis 2011). This method helps examine the complexity of the problem and judge several alternatives based on the appropriate selection of alternatives (Malczewski and Liu 2014). MCDM also helps decisionmakers store, modify, analyze, and visualize data. On the other hand, a geographic information system (GIS) is a platform where the decision-maker can explore and realize the alternatives' desirability. An expert's opinion is required to understand the relative importance of various parameters, which sometimes leads to several uncertainties (Jankowski and Nyerges 2001;Meng and Malczewski 2015). We have applied MCDM models in an integrated manner, i.e., the AHP-VIKOR and AHP-GRA have been used to map the earthquake vulnerability in the entire Himalayan seismicity zone. The AHP is used to evaluate the priority and criterion ranking. The VIKOR and GRA help allocate rank to a social, geotechnical, structural, and physical vulnerability and produce the final earthquake vulnerability map.

AHP method
Analytic Hierarchy Process (AHP) (Saaty 1980) is a commonly used, powerful, and simple MCDM that includes objective and subjective factors that help decision-makers manage many complex problems. In summary, the method requires a pairwise comparison matrix created with the criterion scores. Comparison between each alternative is performed with some specific criterion (x 1 , x 2 …x n ). The scores scaled between 1 and 9 are assigned for equally important and extremely important criteria (Saaty 1980), based on expert's opinion and/or on the basis of available literature. The criteria weights (w i ….w n ) are then evaluated using the normalized matrix with the condition that The weighted sum vector is formulated as W i = ∑x ij w j, where x ij indicates the ith class rank for the jth layer, and W j is the weight for the jth layer. Then, the consistency index (CI) and consistency ratio (CR) is given by where n is the number of criteria being used, λmax is the principal eigenvalue, and RI is a randomness indicator that depends on the dimension of the comparison matrix. Boulos (2003) estimated the RI's value for the matrices with a dimension between 1 and 15. For a consistent matrix, CR's acceptable value must be less than equal to 0.10. CR > 0.10 indicates that the inputs obtained from the expert's opinion and/or based on available literature are unreliable and need to be reviewed. After developing the comparison matrix, the priorities, consistency ratio, and rank can be estimated. In this study, assessment of social, CI = max − n ∕ (n − 1); and CR = CI∕RI physical, geotechnical, and structural vulnerabilities have been calculated using the technique discussed above for the study area.

VIKOR method
The Vise Kriterijumska Optimizacija I Kompromisno Resenje (VIKOR) method (Opricovic 1998;Mardani et al. 2016) is often used to solve MCDM problems based on the aggregating LP-metric function. This approach is considered a powerful tool to address the MCDM problems, especially when the decision-maker is unaware of the system design at the foundation level. In case of contradictory criteria, this method chooses between sets of available alternatives (Chiu et al. 2013) and properly optimizes the multi-criteria complex systems (Opricovic and Tzeng 2004). Firstly, the best ( X * i ) and worst ( X i ) values for each criterion from the normalized decision matrix, j = 1, 2….m is computed as , where X i and X * i are the minimum and maximum value of X ij . The value of the matrix is X ij , and X i indicates the alternative of i th criteria. Secondly, the S j and R j, which signify the utility measures and regret the alternative (X j ) measures, are computed as: where W i is a criterion weight that expresses their importance relatively, in this case, i = 1, 2…n. After that, values of Q j , j = 1, 2…j is estimated where S * =min (S j ), S = max (S j ); and R * = min (R j ), R= max (R j ). S* and S indicate the minimum and maximum measures of utility, whereas R* and R represent the minimum and maximum measures of regret, respectively. ν is the strategical maximum group utility weight, while 1-ν is measured as an individual regret weight. Typically, the ν value is considered 0.5 (Liu et al. 2013), but if ν > 0.5, the Q j index rank can be determined that incline toward the majority agreement. In the case of ν < 0.5, the Q j value indicates the majority's negative attitude.

Grey relational analysis (GRA) method
The Grey Relational Analysis (GRA) model was developed to efficiently address MCDM problems with inaccurate and inadequate information (Deng 1982(Deng , 1988Hamzaçebi and Pekkaya 2011;Wu and Peng 2016). Normalization of data is required to avoid the unit difference and make the dataset unbiased (Haq et al. 2008). In this case, the reference sequence is X 0j, and the original data Y ij is attributed j of alternative i. Translation of Y ij into the comparability sequence X ij is possible as follows to maximize and minimize the response, Δ ij, the deviation sequence of the reference and the comparability sequence are computed as Δ ij =|X 0j -X ij |, and the grey relational coefficient is calculated as where Δ min and Δ max indicate the minimum and maximum values of the deviation sequence, ζ is the distinguishing or identification coefficient, and the value of ζ belongs to 0 to 1; usually, taken as 0.5. The purpose ζ is to expand or compress the range of the grey relational coefficient. Finally, the grey relational grade is computed as Here, W j is the attribute weight j and ∑ n j=1 W j = 1 and Γ (X 0 , X i ) is the grey relational grade (GRG). GRG signifies the correlation level between the reference and the comparability sequence and represents overall quality characteristics (Haq et al. 2008).
We adopted an integrated MCDM approach to the problem, i.e., AHP-VIKOR and AHP-GRA, using 26 parameters to estimate four vulnerability maps. The selected parameter was categorized into four types of vulnerability, i.e., social, geotechnical, structural, and physical vulnerability. The social, geotechnical, structural, and physical vulnerability maps were generated using AHP through experts' opinions and knowledge. The VIKOR and GRA approaches were used to rank the four vulnerability layers and develop the earthquake vulnerability maps for the region (Tables 6 and 7). The weighted sum technique has been used for all the layers with rank and weight values to develop the final vulnerability maps. The weighted sum approach is formulated as V = ∑ n i=1 L(AHP). Where V is the final vulnerability layer, L(AHP) is the layers acquired by the AHP approach, and i is the number of layers (Jena et al. 2020). The classification of each layer is done as suitable conditions (2 classes), medium domain (1 class), and unsuitable conditions (2 classes) to standardize the criteria. Based on this vulnerability, layers are scaled from 1-5. The scaled maps were categorized into five classes, i.e., very low, low, moderate, high, and very highly suitable. The final vulnerability maps were compared to evaluate the area, population, and the number of buildings under the vulnerable zones and their statistical relationships.
We used Pearson correlation to evaluate the performance of the applied integrated MCDM models. The Pearson Correlation has been applied in bivariate data, which produced a correlation coefficient (r) that measures the strength and direction of linear relationships between pairs of continuous variables which are normally distributed. The Pearson's correlation coefficient (Swinscow 1997) is defined as: where xi and yi are the values of x and y for the ith individual and x and y are the mean values of x and y. The value of the correlation coefficient varies from -1 to + 1. The positive and negative value shows a positive or negative linear correlation, respectively.

Results and discussions
The Himalayan mountain chain covers a large area spread in several countries. Estimating earthquake vulnerability is a multi-criteria problem and substantially includes several parameters that are likely to affect the region in case of an earthquake, such as population density, quality of buildings and infrastructure, distance from the seismo-tectonically active geological features, historical and recent seismic activities, availability of the health care facility, communication, and transport network, etc. Therefore, the multi-criteria approach is better suited for the evaluation of such problems. In the following sections, we discuss our results obtained for the region using the techniques discussed above.

Social vulnerability
The earthquake hazard is more multi-faceted and riskier in urban areas because of rapid urbanization without proper planning and management (Rahman et al. 2015). Unplanned growth in cities results in the inadequate distribution of health and infrastructure facilities besides environmental issues and societal problems (Hassanzadeh et al. 2013). This directly impacts society's most vulnerable groups, i.e., women, children, the physically disabled, and the people below the poverty line. People are constantly moving towards urban areas because of poor livelihood in rural areas, pressuring the local resources. People are shifting to temporary shelters that are ill-maintained and enhance the risk of more casualties due to the secondary effects of earthquakes, such as the collapse of buildings and infrastructures. However, an increase in the literacy rate enhances awareness among people about such disasters and their impacts on society (Rygel et al. 2006) and affects earthquakes vulnerability (Martins et al. 2012). Information about populated sites, religious places, parks, and visiting places should be timely updated to minimize the risk factor (Wisner et al. 2003;Alizadeh et al. 2018b). Parameters included in estimating social vulnerability are shown in Fig. 2, and the classes of the parameters with AHP calculation are shown in Table 2. As the Himalayan region is not uniformly inhibited, the estimated social vulnerability map appears to be nonuniform (Fig. 3), as expected. The risk is higher in high population density areas of the central Himalayan region and in the westernmost part of Himalaya, where few major cities are located. Kathmandu valley and its surrounding region with relatively high population density and literacy rate and hosting several famous and religious places appear to be under very high vulnerability zone, which is also true for cities such as Islamabad, Rawalpindi, etc. With a high population and low literacy rate, certain areas are relatively more vulnerable socially, whereas areas with low population density, such as Bhutan, are under relatively very low to low social vulnerability. With a chain of high mountains and harsh climatic conditions, the Tibetan region is relatively sparsely populated. Therefore, these areas fall under very low to low social vulnerability as per our analysis. The area within northern India with moderate to high population density appears very low to moderately vulnerable considering the parameters used in our study. On the other hand, Bangladesh and Myanmar encompass a small fraction of the Himalayan region extending into their territory. Bangladesh, one of the densely populated countries in the region, appears to be under low to moderate social vulnerability. An average population density of 83 per km 2 and a literacy rate of about 90% make Myanmar relatively low socially vulnerable as per our study.
Overall ~ 55.46% of the study area falls under very low social vulnerability, 32.86% low, 7.52% moderate, 3.27% high, and the rest is under very high vulnerability. The social vulnerability as a function of the population comprises 15.81% of the population under very low social vulnerability, 37.42% low, 14.71% moderate, 16.05% high, and 16.01% (Table 8). Therefore, our concern should be moderate to a very high degree of social vulnerability.

Physical vulnerability
The time is taken to respond to an emergency lead to the concept of "Golden Hour". The inaccessibility of transport networks, communication blockage, destruction, or lack of healthcare facilities and services increases the mortality rate. Thus, distance to major healthcare services and hospitals, lack of road, rail, and air transport networks, lack of community services such as fire service, disaster management teams, etc. plays a vital role in mitigating and imparting proper rescue during an earthquake and related secondary impacts (Hosseini et al. 2009;Manshoori 2011;Karimzadeh et al. 2014;Alizadeh et al. 2018b). Usually, underground gas pipes are destroyed due to the ground motion during an earthquake, leading to increased casualty rates in the nearby areas. Therefore, we consider proximity to the gas station as a secondary parameter in our study of physical vulnerability.
The various parameters and their respective rankings estimated using AHP are shown in Table 3. Physical vulnerability (Fig. 3) estimates using parameters discussed earlier indicate that 32.14% of the study area lies under very low physical vulnerability, 23.97% low, 18.25% moderate, 16.50% high, and 9.14% very high vulnerability. The geotechnical vulnerability in terms of the population indicates that ~ 77.35% of the total population is under very low vulnerability, 16.05% low, 5.71% moderate, and a negligible population falls under high to very high vulnerability (Table 8). Thus, it is evident that more than 90% of the total population is under very low to low physical vulnerability. Similarly, most populated major cities in the Himalayan region are marked under very low physical vulnerability (Fig. 3), except that a few are under moderate to high vulnerability. The Tibetian Fig. 2 Flowchart illustrating method adopted to compute various vulnerability maps for the study area area with hilly terrain makes road and rail networks construction challenging and makes the area high to very high physical vulnerability. Therefore, we can suggest that the countries in the region should try their best to develop adequate medical facilities and good communication and transport networks essential for disaster mitigation.

Geotechnical vulnerability
The seismicity data indicates the possible earthquake zones in a particular region (Soe et al. 2009). Thus, distance from active fault plays a significant role, and closeness to the fault enhances risk due to seismic hazard (Alizadeh et al. 2018a;Soe et al. 2009). Peak ground acceleration is an important parameter that provides ground shaking information related to earthquake magnitude, source to site distance, and lithology (Soe et al. 2009). Lithology, topography, cracks, and fractures control probable earthquake locations (Hosseini et al. 2009;Alizadeh et al. 2018b). Well-documented active faults and complex tectonics make the Himalayan region one of the most active zones for the occurrence of catastrophic earthquakes and making this region geotechnical more vulnerable (Fig. 4). Estimating the geotechnical vulnerability has been done using parameters such as proximity to the fault and earthquake epicenters, earthquake density, peak ground acceleration, lithology, and elevation. The classes of the factors contributing to the geotechnical vulnerability with AHP calculation are shown in Table 4. The estimated geotechnical vulnerability indicates that the presence of prominent tectonic features such as the Main Central Thrust (MCT), Indus Tsangpo Suture (ITS), and Main Boundary Thrust (MBT) makes the region moderate to highly vulnerable. The geotechnical vulnerability also reveals that Kathmandu valley and Islamabad may be under a very high threat of geotechnical vulnerability. As discussed earlier, Nepal recently faced one of the destructive earthquakes and having several active faults, making the region moderate to very highly vulnerable. In India, most cities are under low to moderate vulnerability except Darjeeling, which is under high threat because of nearby active faults. Based on our analysis of the geotechnical vulnerability, we find that 11.20% of the total area falls under very low vulnerability, 27.71% low, 31.69% moderate, 21.84% high, and 7.55% under very high vulnerability. 3.71% of the total population lies under a very low geotechnical vulnerability, 22.37% low, 29.96% moderate, 28.50% high, and 15.46% very high vulnerability (Table 8). About 44% population lies under high to a very high degree of vulnerability, distributed in more than 29% of the area. Thus, a considerable percentage of people in the region are under threat of high to very high geotechnical vulnerability.

Structural vulnerability
Proper planning and management of land use and infrastructure are required to minimize earthquake vulnerability (Aghataher et al. 2008). The damage during an earthquake will be less for well-constructed sites, stadiums, museums, etc. On the other hand, if the construction quality is poor, then the probable damage and loss of life may be even more during an earthquake if that occurs during holidays when more people are likely to visit those places (Alizadeh et al. 2018a, b). As mentioned earlier in the social vulnerability, people are usually moved to the urban areas for better livelihood. Due to lack of resources, they prefer to construct non-engineered houses that are most vulnerable to disaster. Similarly, vulnerabilities at hotels, transport terminals, dams are multi-faceted (Brown et al. 2017). Dams are not inherently safe against seismic events (Wieland 2016). Lack of maintenance of old structures makes them risky and prone to structural failure during earthquakes (Orchiston, 2012;Roark et al. 2000). The structural vulnerability ware computed, ignoring variations in parameters during day and night (Martins et al. 2012). The weight of the parameters included in the structural vulnerability with the classes is shown in Table 5. It is observed that most of the metropolitan cities in the study area fall under moderate to a very high degree of structural vulnerability (Fig. 4). As discussed earlier, due to harsh climatic conditions in the Himalayan region, people move to urban areas for better earnings and usually building non-engineered structures to reside, making the cities relatively more Red zones are highly vulnerable, while the green zones are relatively low vulnerable based on the geotechnical criteria structurally vulnerable. Several cities in Nepal and Pakistan which are relatively populated may be under high structural vulnerability. Whereas, due to sparse anthropogenic activities and the development of infrastructures, the Tibetian region may be under a very low structural vulnerability. In the Indian Himalayan region, many cities may be under moderate to high structural threat due to the unplanned development and illmaintenance infrastructures (Tables 6, 7).
The result also reveals that out of the total area, ~ 25.47% area falls under very low vulnerability, 30.98% low, 31.95% moderate, 7.46% high and 4.14% very high structural vulnerability area. On the other hand, with respect to population, 1.08% of the population in the study area is under very low vulnerability, 13.19% low, 48.55% moderate, 23.37% high and 13.82% of the population are under very high risk due to structural vulnerability (Table 8). Interestingly, the very low and low structural vulnerability comprises more than 55% of the area. In contrast, the number is slightly less than 15% in terms of population, and the reason is primarily because of harsh climatic and geographical distribution in the Himalayan region, as discussed above. Therefore, most of the region has a significantly lower population density with low structural vulnerability.

Comparison of the AHP-VIKOR and AHP-GRA vulnerability maps
The final vulnerability maps (Fig. 5) of the region were developed using two integrated MCDM methods of AHP-VIKOR and AHP-GRA. A comparison table has been presented, including the percentage of vulnerable areas, population, and buildings estimated by analyzing AHP-VIKOR and AHP-GRA (Table 9). The final vulnerability maps show almost similar patterns of vulnerability levels for the region. The percentage difference in both maps helps us to analyze the performance of the applied models. For better understanding, we adopted the criteria that if the percentage difference is less than 2%, then both results are 'closely matching,' 2%-5% match was classified as 'moderately match,' and greater than 5% match was classified as 'poorly matching' area (Table 8). The result indicates that out of the total study area with a low, high, and a very high degree of vulnerability, there is a good match in models derived by the two methods, whereas moderate vulnerability areas show moderate match and a relatively poor match very low vulnerability areas. In terms of the vulnerable population, moderately match is seen for low, moderate, and a very high degree of vulnerability, whereas very low and high vulnerability areas, are poorly matching. In contrast, AHP-VIKOR and AHP-GRA for the vulnerable buildings show poorly matching in a very low, low, high, and very high level of vulnerability and closely matched for moderate vulnerability. An estimated percentage of vulnerable areas, populations, and buildings have been evaluated (Table 9) by aggregating the result of AHP-VIKOR and AHP-GRA. The Pearson's correlation coefficient was also computed to evaluate the performance of the used integrated MCDM models. The result shows that the correlation coefficient (r) for the total vulnerable area is 0.9766 (Fig. 6a), which indicates a robust correlation; for the total number of vulnerable buildings, the correlation is about 0.8226 (Fig. 6c). In contrast, a moderate correlation (0.7362) was observed for the total population under threat (Fig. 6b). Therefore, our analysis indicates that earthquake vulnerability maps estimated using AHP-VIKOR and AHP-GRA are well correlated, and their performance is similar to a great extent.

Sensitivity analysis
To understand the impacts of the criteria used in the integrated models, we conducted the sensitivity analysis for the integrated MCDM techniques. For the VIKOR coupled AHP models, sensitivity analysis for influence levels of the v value in the VIKOR method can be done to determine which alternatives are affected significantly by maximum group utility and minimum individual regret (Table. 10). Similarly, for the AHP-GRA, we have done sensitivity analysis for influence levels of the ζ (Distinguishing Coefficient) value in the GRA method to understand the value of the used criteria. The result suggests that in the AHP-VIKOR model, the social characteristics and physical characteristics were not affected by the value of v. The social aspects have a higher Qj value than other characteristics involved in the earthquake vulnerability given the maximum group utility and minimum individual regret. In contrast, physical characteristics show low value in terms of maximum group utility and minimum individual regret. The ranking of geotechnical characteristics and structural characteristics were not affected by the value of v. Similarly, for AHP-GRA, the ranking of the four criteria was not affected by the value of ζ. Increment or decrement in the value of ζ did not impact the grey relational grade value of social characteristics. In contrast, the value of grey relational grade for other aspects techniques. Red zones are highly vulnerable, while the green zones are relatively low vulnerable to earthquake hazards keeps increasing with the increment of ζ, which satisfies the fact that no matter what the value ζ occupies, the ranking is always the same (Kuo et al. 2007). Therefore, both the models validate our ranking procedure.

Conclusions
Identifying vulnerability in larger areas susceptible to severe damage due to earthquakes is one of the most important factors for seismic hazard mitigation. This study estimated earthquake vulnerability for the Himalayan region using integrated MCDM techniques, i.e., AHP-GRA and AHP-VIKOR. Twenty-six parameters were used to estimate social, geotechnical, structural, and physical vulnerability classes.
Having an advance estimate of earthquake vulnerability may help the hazard mitigation and planning agencies identify the most affected areas during an earthquake, and can be utilized to plan in advance in case of any destructive earthquake that might affect the region.