A geospatial analysis of multi-hazard risk in Dharan, Nepal

Natural hazard risk assessment generally focuses on a single hazard type, such as earthquakes, landslides, or floods. This emphasis tends to consider physical processes in isolation. However, most locations are simultaneously at risk to multiple, interacting hazards that generate cascading effects or synergies. Although scholars have proposed a multi-hazard risk framework based on probabilities, the quality and quantity of data required for such an approach are often unavailable in developing countries. Using geospatial and socioeconomic data, this study represents a first step in assessing multi-hazard risk in the city of Dharan, Nepal. Three hazards—landslides, floods, and earthquakes—were considered for an integrated hazard assessment using statistical methods and the Analytic Hierarchy Process (AHP). We employed a Social Vulnerability Index (SoVI) to create a vulnerability map of the study area, which was then combined with a multi-hazard hazard map to produce a total risk map. Our results indicate that eastern Dharan along the Seuti River and southwestern Dharan on the left bank of the Sardu River are at high risk to multiple hazards. Central Dharan and the hills in the western portion of the city are categorized as low risk areas. Data limitations, such as availability and spatial resolution, did not allow for dynamic modeling; however, our results identified the spatial extent of low to high risk areas, which can inform future disaster planning. For example, the methodology and results of this study could assist in the development of disaster risk reduction programs and policies. ARTICLE HISTORY Received 8 March 2019 Accepted 20 December 2019


Introduction
Human and economic losses to natural hazards have escalated in recent decades (Guha-Sapir et al. 2004;Bouwer 2011). Disaster risk reduction (DRR) programs have addressed such losses by identifying and assessing risk in specific geographic locations CONTACT Sanam K. Aksha sanam7@vt.edu (Van Westen et al. 2014;Garcia-Aristizabal et al. 2015). The primary objective of risk assessment is to quantify the probability of human and economic losses due to the occurrence of a 'hazardous event', defined as a process that causes loss or harm to human and economic systems (Wohl 2000). Since hazards occur at the intersection of socio-ecological systems, reducing hazard loss necessitates a holistic assessment of risk that considers integrated nuances of the physical, built, and social environments of a place. Hence, risk assessment, and ultimately vulnerability reduction, are processes that require multidisciplinary knowledge of various coupled physical and social processes to calculate the cumulative level of risk posed by hazards. Studies of natural hazard risk frequently emphasize the impacts of an individual hazard, such as landslides (Althuwaynee et al. 2014;Pellicani et al. 2017), floods (Kazakis et al. 2015;Kabenge et al. 2017), earthquakes (Theilen-Willige 2010; Dhar et al. 2017), drought (Lehner et al. 2006), sea level rise (Hinkel et al. 2014), tropical cyclones (Hoque et al. 2018), or wildfires (Adab et al. 2013). However, in reality, many locations are prone to multiple natural hazards that can occur simultaneously or manifest in a set of cascading events (Kappes et al. 2012). Mitigation and planning for such complex outcomes requires a consideration of combined hazard risk. However, efforts to assess multi-hazard risk are impeded by a multitude of barriers, including lack of a common definition for a multi-hazard risk (epistemological issues) (Marzocchi et al. 2012); developing a common approach for integrating different hazards (methodological issues) (Tate et al. 2010); availability of intensive data (data scarcity issues) (Gallina et al. 2016); confinement of risk assessments within disciplinary boundaries (research domain issues) (Kappes et al. 2012;Barrantes 2018); and the inherent uniqueness of each place in generating geographically-specific hazards conditions and outcomes (place-based issues) (Johnson et al. 2016). As a result, few studies have examined integrating multiple hazards in risk and impact assessments (Van Westen et al. 2014;Barrantes 2018).
Despite these challenges, multi-hazard approaches offer many benefits. First the integration of natural hazard risk and social vulnerability offers a more realistic assessment of potential impact because social, economic, and cultural elements are considered simultaneously along with physical geography (Marzocchi et al. 2012;Johnson et al. 2016). Additionally, modeling the spatio-temporal overlap of hazards while considering the possibility of cumulative synergistic and cascading effects improves upon the modeling of single hazard types, which tends to underestimate or overestimate total risk (Kappes et al. 2012;Van Westen et al. 2014). Moreover, the inclusion of socioeconomic factors, such as income, education, ethnicity, and elderly populations, provides much needed insights on levels of human capacity and resilience before a disaster happens (Greiving et al. 2006;Tate et al. 2010). Since the predisaster context serves either to accentuate or attenuate the impacts of a disaster, a multi-hazard approach that considers the spatial, demographic, and physical contexts and their various linkages and feedbacks could play a significant role in mitigating human and economic losses.
A spatial approach must be employed to effectively model multi-hazard risk. A spatial approach in risk assessment facilitates DRR by providing critical information on hazard source areas, possible impact zones, and the geographic distribution of (vulnerable) populations and infrastructure located in hazardous areas (Greiving et al. 2006). Spatial approaches can also help to identify optimal locations for disaster mitigation infrastructure and assist in evacuation, response, resource allocation, and evidence-based policymaking. Thus, the use of remotely sensed data, geographic information systems (GIS) software, and publicly available social, demographic, and economic data have potential to support modeling efforts in terms of increased spatial resolution, computing capacity, rigor in quantitative techniques, and sharing of data for public benefit (Wohl and Oguchi 2004;Bishop et al. 2012;Hoque et al. 2018) This study introduces a model for spatial multi-hazard risk assessment applied to Dharan, Nepal-a location for which spatial data availability is limited in terms of quality, quantity, and access. The purpose is to use the relevant, publicly available geospatial data to assist local decision and policy makers in the efficient deployment of resources for DRR by developing a procedural model for generating a composite risk map. The overall study objectives are to: (1) produce individual hazard assessments for the rapidly growing city of Dharan, Nepal (i.e., for earthquakes, floods, and landslides); (2) calculate social vulnerability for the city of Dharan; and (3) combine results from objectives 1 and 2 into a comprehensive multi-hazard risk assessment for the city of Dharan.
Multi-hazard risk assessment methods can be broadly categorized into two approaches: (a) assessing individual hazards in a particular geographic location independently (Greiving et al. 2006;Gr€ unthal et al. 2006;Carpignano et al. 2009;Schmidt et al. 2011); and (b) evaluating possible interactions and/or cascade effects among the different possible hazardous events (Nadim and Liu 2013;Garcia-Aristizabal et al. 2015;Zhang and Zhang 2017). Most studies acknowledge that hazard interactions (i.e., triggering and/or cascade effects) are important considerations in hazard risk assessment. However, very few studies have explicitly considered cascade effects and interactions among hazards (Kappes et al. 2012;Marzocchi et al. 2012;Gill and Malamud 2014;Van Westen et al. 2014;Liu et al. 2015). Such studies require substantial amounts of input data for dynamic modeling, and sometimes the complexity of hazard 'chains' that can be foreseen discourage analysts from considering such interactions and triggering effects in a comprehensive multi-hazard analysis (Komendantova et al. 2014;Van Westen et al. 2014;Garcia-Aristizabal et al. 2015;Liu et al. 2015). Since no historical hazard occurrence data and limited physical and social data are available for the study area of Dharan, this study adopts the first approach to identify the spatial distribution of multi-hazards risk in the study area.
In the hazards and geography literature, risk (R) has been understood as the combination of hazard exposure (H) and societal vulnerability (V). This relationship can be expressed in a 'pseudo-equation' of R ¼ H Â V (Wisner et al. 2004). Thus, there is no risk if a hazard and vulnerable population do not interact in a particular location. Vulnerability is defined as the susceptibility of people and communities exposed with their social, economic, and cultural resources to cope with and recover from a hazardous event (Hewitt 1997;Mileti 1999;Cutter et al. 2003). The predisposition of social inequalities, livelihoods, degree of social protection, and community and built environment attributes help to determine the level of impact of any disaster regardless of source or physical processes (Comfort et al. 1999;Wisner et al. 2004).
However, the majority of risk assessment research either focuses on a single hazard type or warrants little attention to vulnerability of society to natural hazards. To address the gap, this study combines an integrated spatial assessment of hazards, based on the approach of Greiving et al. (2006), with the hazards of place model developed by Cutter et al. (2000) (Figure 1). Greiving et al. (2006) consider all spatially relevant hazards that produce total risk in a particular location, and Cutter et al. (2003) Social Vulnerability Index method considers various social, economic, and demographic indicators that influence vulnerability of a society. Together, this hybrid model conceptualizes risk as the joint product of (1) spatially relevant hazards in a particular place, and (2) level of vulnerability present in the social systems in a particular place. The combination of both models provides a holistic assessment risk.

Study area
Nepal is susceptible to a multitude of natural hazards, ranging from frequent, regularly occurring hazards (e.g., floods, landslides, and avalanches), to less frequent but higher magnitude hazards such as earthquakes. Floods, for example occur regularly during heavy precipitation events and the annual monsoon season (June-September), whereas large-scale earthquakes occur periodically. The 2015 Gorkha earthquake claimed more than 9,000 lives and destroyed tens of thousands of houses. The combination of these serial and sporadic hazards makes apparent the high level of risk and low level of disaster preparedness that characterize the country of Nepal (Aksha et al. 2018). The city of Dharan (192 km 2 area and located at 26 51 0 N, 87 13 0 E) is situated in the Sunsari District of eastern Nepal, approximately 600 km southeast of the capital Kathmandu ( Figure 2). Dharan is one of three major urban centers in eastern Nepal, with the latest 2011 census reporting the population at 137,705 (CBS, 2012). Dharan is situated at the foothills of the Siwalik range and is characterized by the presence of very young sedimentary rocks such as mudstones, shale, sandstone, and conglomerates. Furthermore, the Main Boundary Thrust (MBT) runs along the north side of Dharan, placing the entire city at seismic risk. The MBT is an active thrust running east-west along the Himalayas that is capable of initiating major earthquakes at any time (Upreti 2001).
Dharan is rapidly growing in population and areal extent, with many settlements expanding on fan deposits of the Sardu and Seuti Rivers that flank the city ( Figure  2). The two rivers frequently flood during the annual monsoon season. A large agglomeration of squatter settlements has also emerged along the banks of Seuti and Sardu Rivers. This highly vulnerable area contains approximately 6,500 households scattered across several administrative wards. Risk to dry landslides is severe in Dharan due to riverbank cutting by the Sardu and Seuti Rivers. Serial flooding and wet landslides during the monsoon season complicate this risk by continuously deteriorating agriculturally productive land and posing constant, increasing risk to settlements that have sprung up as part of urban sprawl processes (Dharan Municipality 2014). Efforts to reduce vulnerability are further exacerbated by weak institutional memory of past disaster events, scarcity of financial and capital resources, and limitations in quantitative, geospatial, and socioeconomic data and their applications.

Data and sources
For our hazard, vulnerability, and risk assessment of Dharan, we collected remotely sensed imagery, hundreds of landslides geolocations, and building footprint data with  respect to floodplain locations, topographic data, environmental variables, triggering factors, and social data (Table 1). These data were obtained from publicly available sources such as the most recent 2011 census (CBS, 2012), socio-economic data available at the Dharan sub-metropolitan office, and several field investigations (acquiring spatial administrative data, verification of landslides, interviews with key informants). Data were managed in a GIS environment and analyzed using several statistical and geospatial methods. High resolution WorldView-3 imagery (1.24 m resolution) dated 3 February 2015 was obtained from the DigitalGlobe Foundation to produce a land use land cover map. Monthly rainfall data for the 2016 calendar year were obtained from the Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks (PERSIANN) available at the Center for Hydrometeorology and Remote Sensing, University of California (http://chrsdata.eng.uci.edu). River flow data were not available for the Seuti and Sardu Rivers because no gauging stations are installed. Based on visual interpretations of WorldView-3 images and Google Earth, an inventory of active landslides was prepared. Landslide locations were field-verified in summer 2016. The methodological framework depicted in Figure 1 guided our overall analysis. Based on previous studies, a spatial database of 12 physical variables (aspect, plan curvature, distance from faults, distance from streams, elevation, flow accumulation, land use, lithology, normalized difference vegetation index (NDVI), rainfall intensity, slope, and topographic wetness index (TWI)) ( Table 1) was prepared considering their influence on multi-hazard risk. Each spatial layer was transformed into a grid spatial database at a pixel size of 30 Â 30 m. All raster layers were referenced using World Geodetic System (WGS)-1984 and the Universal Transverse Mercator (UTM) zone 45 North. Each hazard (landslide, flood, earthquake) was modeled individually based on causative factors, and the models were subsequently integrated to generate an integrated hazard map of the study area. Finally, social, demographic, and economic data (as documented in Table 1) were used to assess social vulnerability at the ward level, which was overlaid on the integrated hazard map to produce a composite risk map of Dharan.

Hazard assessment
Three types of natural hazards (i.e., landslides, floods, and earthquakes) were included in this study based on hazard risk exposure. First, each individual hazard was mapped separately based on selected conditioning factors (Table 1). Subsequently, the three layers were overlaid to prepare an integrated hazard map. These processes are detailed below.

Landslide hazard assessment
Landslide hazard assessment helps to understand the basic characteristics of slopes and their susceptibility to failure. Hillslope failure results from interacting factors that weaken shear strength and/or increase shear stress. Influencing factors are often topographic of origin, such as slope angle, aspect, elevation, and surface curvature (Devkota et al. 2013;Bathrellos et al. 2017). Furthermore, local (hillslope) hydrology, soil moisture, and climatic triggering events are key factors. Rivers and streams are also important factors in triggering landslides because they both govern and can change slope cutting movements across space, time, and event (Van Westen et al. 2008). Land use and vegetation cover influence slope failure by affecting slope stability, hydrological regimes (via the obstruction and interception of precipitation), changes in runoff speed, and perturbations in groundwater flow (Vanacker et al. 2003;Van Westen et al. 2008). Underlying lithology and geology (e.g., location of faults) are also important to consider (Youssef et al. 2016).
We employed a general linear model (GLM) to evaluate landslide hazards in Dharan. GLMs have been employed by other researchers to assess landslide hazards (e.g., Mousavi et al. 2011;Devkota et al. 2013;Youssef et al. 2016). For purposes of using landslide presence/absence as a binomial predictor in a logistic regression model, 71 non-landslide points were randomly generated in ArcGIS 10.3 (Mousavi et al. 2011;Devkota et al. 2013). These points were randomly divided into training (70%) and validating (30%) datasets for model fitting and validation. Based on visual interpretation of WorldView-3 satellite images, the 71 landslide locations were geolocated and 63 were field-verified in July 2016 ( Figure 3). Due to inclement weather and steep, unstable terrain, eight landslide locations could not be field-verified but were assessed using Google Earth imagery.
Based on a review of the literature, we selected ten relevant variables for potential inclusion in the landslide model. These variables included aspect, plan curvature, distance from faults, distance from streams, elevation, land use, lithology, normalized difference vegetation index (NDVI), slope, and topographic wetness index (TWI) ( Table 1). Slope angle, slope aspect, plan curvature, and TWI were derived from a DEM using Spatial Analyst tools in ArcGIS 10.3. We also measured distance from the fault-line and two rivers (Sardu and Seuti) using the Euclidean Distance tool in ArcGIS. High-resolution satellite images were classified into four classes (forest, urban area, water body, and agriculture/open/barren) to obtain a land use map of the study area (Table 1). Aspect, land use, and lithology were defined as categorical variables, and the remaining variables were defined as nominal variables. Values for each landslide and non-landslide point were extracted and imported into R version 3.3.3 (R Core Team 2017) to perform a binomial logistic regression using the GLM function in R Studio. We assessed multicollinearity using a cut-off value of 0.8. As a result of this assessment, two variables (distance from faults and slope) were removed prior to modelling.

Flood hazard assessment
Flood hazard assessments are generally based on meteorological, hydrological, and geomorphological data and typically use hydraulic modeling to predict flood depth based on rainfall-runoff relationships (Vojtek and Vojtekov a 2016). However, many developing nations, including Nepal, lack adequate sets of such data (Skilodimou et al. 2019). The use of hydraulic modeling is limited in these contexts and instead a GIS-based flood hazard analysis can be used to estimate flood hazards (Rawat et al. 2017;Shafapour Tehrany et al. 2017).
Distance from drainage networks is crucial as flood inundation results from the overflow of drainage channels. Areas closer to drainage channels have relatively higher risk of inundation than areas farther away (Butler et al. 2006). Water flows from higher elevations and slopes to lower elevations and slopes. Lower slopes decrease the speed of runoff and are more quickly inundated than higher slopes (Kabenge et al. 2017). High values of flow accumulation indicate concentrated flows and consequently higher flood hazard zones (Kazakis et al. 2015). The geology of an area in terms of soil type and parent rock formation may amplify or extenuate the magnitude of flood events. Soil type determines the infiltration and water holding capacity of an area, thus affecting flood susceptibility (Tehrany et al. 2013). Land use and land cover play roles in interception of precipitation, infiltration, evapotranspiration, and underground water holding capacity of an area. Natural vegetation cover such as forests and grassland reduces the speed and amount of runoff more than built-up areas (Yalcin and Akyurek 2004). Higher rainfall intensity can result in quickly reaching infiltration capacity; thus, the likelihood of flood increases in an area where the amount of rainfall increases (Ouma and Tateishi 2014).
To generate the GIS-based FHI, we extracted elevation, flow accumulation, and slope from a DEM. Geological information was obtained from engineering geological maps of the study area. Euclidean distance was used to calculate distance from drainage networks and a modified Fournier index was applied to calculate rainfall intensity from rainfall measurements (De Luis et al. 2010). The relative importance of each variable was determined based on the Analytic Hierarchy Process (AHP), detailed in section 3.2.4.

Earthquake hazard assessment
Earthquake hazard assessment is the process of evaluating the design motion of the ground during earthquakes, since most of the damages are caused by intense ground shaking or structural collapse. In general, two approaches are used, scenario and probabilistic, which require historical data, geology, tectonics, fault activity, and seismic source models (Giardini 1999;Chaulagain et al. 2015). Probability of exceedance of some particular value of ground motion parameter, most commonly peak ground acceleration (PGA), is a widely adopted method in regional and national earthquake hazard assessment (Papadopoulou-Vrynioti et al. 2013;Chousianitis et al. 2016). At the local level, microzonation is a common tool which primarily includes the effect of only sedimentary layers (Allen and Wald 2009). While impacts are significantly influenced by topographic amplification based on local geological conditions, there are other factors beyond sedimentary layers (Hough et al. 2010).
Advancements in satellite technologies and computing capacity have also allowed earthquake hazard assessment using geomorphometric variables based on DEMs (Geiß and Taubenb€ ock 2013). Several studies have analyzed geomorphic/topographic features in earthquake-prone areas, primarily including five criteria-variables: distance from faults, elevation, flow accumulation, lithology, and slope (Wald and Allen 2007;Allen and Wald 2009;Theilen-Willige 2010;Dhar et al. 2017). We adopted the same variables since they are relevant to our study area as well (Table 1). Earthquake damage varies mainly due to local lithological and hydrogeological conditions. The presence of tectonic structures such as faults and folds can influence seismic hazards and their secondary effects such as land subsidence, liquefaction, and building collapse (Dhar et al. 2017). Higher slope angles and elevations contribute to mass movements more than lower slopes and elevations. Higher flow accumulation is characterized by young unconsolidated deposits and high groundwater tables, which increase the chance of liquefaction during earthquakes (Theilen-Willige 2010).

Weighting criteria-variables using AHP
To integrate various variables in a spatial decision-making process and determine their relative importance, weighting and ranking of the variables (Table 1) for each individual hazard type (landslide, flood, and earthquake) is required (Bathrellos et al. 2012). For the landslides hazard assessment, a statistical method, described above, was adopted (Althuwaynee et al. 2014;Youssef et al. 2016). However, to determine the weight of each variable for flood and earthquake hazards, AHP method was used (Saaty 2008). AHP is a widely used a multi-criteria, decision-making tool that serves a basis for integrating and determining the importance of criteria-variables (Althuwaynee et al. 2014;Kazakis et al. 2015;Kabenge et al. 2017;Hoque et al. 2018). We prepared a pairwise comparison matrix for flood and earthquake hazards. Three experts (local officer from Dharan municipality, a professor from Tribhuvan University in Nepal, and the first author) have filled in the matrix. The relative importance of criteria was assessed from 1 to 9 implying low to high importance, respectively (Saaty 2008) (Table 2). A consistency ratio (CR) was used to justify the consistency of comparisons in the pairwise comparison matrix, and comparisons were considered consistent if CR was less than or equal to 0.1 (Hoque et al. 2018).

Integrated hazard assessment
Susceptibility to individual hazards was first identified based on the weightings of relevant variables. Subsequently, we overlaid the individual hazard maps using the Weighted Overlay Function tool in ArcGIS 10.3 to prepare an integrated hazard map of the study area. Dharan experiences floods and landslides every year, and there was a damaging earthquake in 1985. Since the study area is at constant risk to all three hazards, we assumed that each had the same relative importance and employed equal weights when preparing the integrated hazard map (Collins et al. 2009). Each hazard maps were classified and then ranked 1-5, using the Jenks Natural Break classification method provided in ArcGIS (1 ¼ low, 5 ¼ high) (Devkota et al. 2013;Hoque et al. 2018).

Vulnerability assessment
Vulnerability is primarily controlled by social, economic, and demographic factors that coalesce to influence the capacity of individuals and communities to mitigate, cope with, and reduce disaster risk (Hewitt 1997;Comfort et al. 1999;Juran and Trivedi 2015). In general, a population with relatively greater socioeconomic status, access to resources, and built environment attributes is less vulnerable and typically performs better during and after disaster events. Based on these concepts, we calculated and mapped social, demographic, and economic data (Table 1) to assess social vulnerability at the ward level of Dharan (total of 25 wards).
We adapted the Social Vulnerability Index (SoVI) method of Cutter et al. (2003) to quantify, map, and assess social vulnerability to hazards in the study area. Social, economic, and demographic variables were obtained from the most recent Nepal census as well as from the city office of Dharan. Among the variables, 18 total such as age, built environment, education, ethnicity, family structure, gender, level of employment, occupation, renters, socioeconomic status, special needs, and urban/rural were included (Table 1). Data were aggregated at the ward level and then normalized and standardized for further processing. Principal component analysis (PCA), using the Dimension Reduction tool in SPSS version 22.0, was used to reduce the 18 variables into a smaller number of more meaningful components (Aksha et al. 2019). Varimax rotation and Keiser criterion were employed to identify components with eigenvalues Moderate importance One variable is slightly favored over another 5 Strong importance One variable is strongly favored over another 7 Very strong importance One variable is very strongly favored over another 9 Extreme importance The evidence favors one variable over another and is of the highest possible order of validity 2, 4,6,8 Intermediate values between the two adjacent judgements When compromise is needed higher than 1. Based on cardinality, all components were summed to obtain SoVI scores of the study area. The values were then brought into an ArcGIS environment to create a vulnerability map.

Overall risk assessment
The total risk of Dharan to three natural hazards (landslide, flood, and earthquake) was based on spatial intersections of hazard risk and social vulnerability in the study area ( Figure 2). The integrated hazard map and social vulnerability layers were multiplied in ArcGIS to obtain a multi-hazard risk map. The Jenks natural breaks method was applied to obtain five classes of risk ranging from low to high (Devkota et al. 2013;Hoque et al. 2018).

Individual hazard assessments
The binomial logistic regression for prediction of landslide presence and absence revealed that distance to streams, elevation, lithology (Siwalik), and lithology (Midland) to be significant predictors (p < 0.05) of landslide presence or absence in Dharan; only variables with statistical significance at the 95% confidence level were included in the risk assessment map. Based on the coefficients of each significant predictor, a landslide hazard susceptibility map was prepared ( Figure 4A). The results indicate that 45% of the study area is located in high hazard zones for landslides, meaning that these areas exhibit relatively high susceptibility to landslide occurrence. The forested areas in the south and stable hills in the north are delineated as regions that exhibit relatively low landslide risk. The results are significantly influenced by the lithology (Midland) variable, which has the highest coefficient values among others in the logistic regression equation (Table 3). Due to such a heavy influence of lithology (Midland), the central part of Dharan is classified as a high landslide hazard zone ( Figure 4A). The relative importance of seven indicators of the Flood Hazard Index (FHI), obtained via Analytic Hierarchy Process (AHP) ( Table 4), were used as inputs for the weighted overlay function in ArcGIS. The thematic layers were then overlaid to prepare a final flood hazard map ( Figure 4B). The map shows that roughly 47% of the study area lies in low flood hazard zones. The majority of high flood hazard zones are near the Seuti and Sardu Rivers. The hills in east, north, and west, as well as comparatively elevated regions show low risk to flood hazards.
An earthquake hazard map ( Figure 4C) was prepared for the study area using five different thematic maps: distance from faults, elevation, flow accumulation, lithology, and slope. AHP determined the weightage of the variables and those values were used while overlaying layers in GIS (Table 5). The map indicates that high hills in the north have relatively lower seismic risk, while central Dharan is categorized as a high risk earthquake hazard zone. Since the Main Boundary Thrust (MBT) passes just north of Dharan-and since a few local faults are also present in the study area-a majority of the human settlements in Dharan lie in high earthquake hazard areas. Forested areas in south are classified as a medium hazard zone, likely because their low elevation allows for soil liquefaction due to high groundwater table movement (Theilen-Willige 2010).
The integrated hazard map-prepared by combining the individual landslide, flood, and earthquake hazard maps-is presented in Figure 4D. The map suggests that approximately 2% of the study area lies in a high risk area while 7% is classified as low risk, 26% as low to medium risk, 30% as medium to high risk, and 35% as medium risk. The areas adjacent to rivers are predominantly under high risk to hazards, whereas forested lands in the north and south are categorized as low risk areas.

Vulnerability assessment
PCA reduced 18 variables (Social Data section in Table 1) into four components that explained 86.62% of the variance of the data. The four components are Socioeconomic Status, Education and Built Environment, Ethnicity, and Disability. Socioeconomic Status explained 45.12% of the variance; Education and Built Environment explained 19.74%; Ethnicity explained 10.98%; and Disability explained 10.78%. Based on cardinality, the components were summed to calculate SoVI scores for the 25 wards of Dharan that were assessed in this study, which were then mapped ( Figure 5).
Of the 25 wards, SoVI analyses categorized three as highly vulnerable, six as medium-high, seven as medium, six as low-medium, and three wards as low ( Figure  5). The social vulnerability map reveals that wards 20, 21, and 22 are the most vulnerable, while wards 1, 4, and 5-which are all located in the core of Dharan-are classified as the least vulnerable.

Total risk assessment of dharan
The final risk map ( Figure 6) portrays the greatest risk areas as eastern Dharan along the Seuti River, as well as southwest Dharan along the left bank of the Sardu River. The southern part of the study area is classified as low risk. Southern Dharan is relatively lower in elevation, flatter, and forested. The slope along the Sardu River and its headwater is also categorized as low risk. Central Dharan, the oldest part of the city, is classified as low risk.

Discussion
The overall objective of the study was to assess multi-hazard risk in Dharan, Nepal, using socioeconomic data and geospatial techniques. Based on the literature and methodology used in similar studies (e.g., Greiving et al. 2006;Collins et al. 2009;Van Westen et al. 2014), we generated individual landslide, flood, and earthquake hazard maps. Next, we weighted and combined these maps to delineate low to high hazard risk areas of Dharan based on the combined influence of the three hazards. Finally, we incorporated vulnerability into a final composite risk map, using a SoVI map generated using the SoVI method. Landslide presence in the study area was most predicted by lithology (Siwalik) ( Table 3 and Figure 4A) (3.05, p < 0.01). Siwalik group is characterized by the presence of loose materials comprised of grained sandstones, mudstones, siltstones, and shales (Upreti 2001;Chamlagain 2009). Such lithology renders the area more susceptible to slope failures on stream banks, the formation of steep cliffs, differential erosion resulting in scraped ridges, and frequent landslides on steep slopes.  The flood hazard assessment ( Figure 4B) indicates that the higher flood risk areas are distributed along the banks of Sardu and Seuti Rivers. Both rivers are characterized by active channel shifting due to loss of high sediment loads during the monsoon season (Stoffel et al. 2016) combined with human encroachment in the form of buildings along the riverbanks as well as commercial sand mining and stone quarrying (Dixit 2003;Sudmeier-Rieux et al. 2012). In the winter, during low flow of the rivers, commercial exploitation of the Sardu and Seuti riverbeds deepens their channels, which later fill during monsoon season. One example of simultaneous riverbank encroachment and expansion of human settlements towards the rivers is clearly visible in Figure 7.
Earthquake hazard assessment determined that most of Dharan (16 of 25 wards) is under medium to high, and high risk ( Figure 4C). However, in reality, this measure merely determines relative risk among political geographies that all lie in a high earthquake prone region (i.e., the entire country of Nepal is considered to be at high seismic risk (Chamlagain 2009)). Thus, all wards are at high risk, but some are riskier than others.
From a multi-hazard point of view, almost all densely populated areas of Dharan are classified as medium, medium to high, or high hazard risk zones ( Figure 4D). This is likely due to the compound nature of hazard risk in the study area (i.e., most areas are susceptible to multiple hazards that also have the propensity to cascade). For example, while any given location may exhibit low risk to one natural hazard (e.g., floods in northeastern parts of the study area in Figure 4B), it may simultaneously be susceptible to other natural hazards or hazard risk (see same area in Figure 4D).
Social vulnerability analyses revealed that central Dharan is less vulnerable than other parts of the city ( Figure 5). Central Dharan is the oldest part of the city and has better access to piped water, sewerage infrastructure, and electricity (Dharan Municipality 2014). In fact, this part of the city has no dependence on firewood as a fuel source, which helps to make the environment clean and reduce health problems, especially related to indoor air pollution. Central Dharan also boasts greater educational attainment (which increases awareness on hazards, risk, and vulnerability) and greater access to critical resources such as preparedness materials, emergency shelters, and other services provided by government and non-governmental institutions. Interestingly, the western part of the study area is classified as low risk even though the area is characterized by relatively poorer built environment attributes (e.g., fewer electricity connections, less sewerage infrastructure, and more low quality external walls) and a greater dependence on agriculture for livelihood and firewood as fuel source (Dharan Municipality 2014). Areas classified as medium and medium to high vulnerability are the fastest growing regions of the study area, and they are expanding towards the banks of Sardu and Seuti Rivers. Figure 7A and B provide context for such expansion along the banks of Seuti River (southeastern part of the study area) between 2004 and 2018.
The distribution of risk to these hazards in the study area ( Figure 6) is partly explained by the steep mountainous terrain, rapid unplanned urbanization, and encroachment on the Sardu and Seuti Rivers. Rapidly and poorly planned urbanization in combination with high river bed exploitation and river bank invasion represent significant ongoing pressures (Sudmeier-Rieux et al. 2012;Liu et al. 2016). As a result, natural hazards are more complex and uncertain due to dynamic interactions among increasing populations and land use development, existing ecologies and ecosystem services, and environmental degradation for housing, livelihood, and commercial purposes.
In general, results confirm the literature the distribution of risk is not solely attributed to the presence of hazardous physical agents (high slope gradients, active riverbanks, presence of fault lines) but result of interaction with social environment and human vulnerability (Hewitt 1997;Mileti 1999;Cutter et al. 2000;Wisner et al. 2004). Central Dharan being less risky than other areas of the city also aligns with the literature. Central Dharan is relatively farther from the rivers, contains fewer slums, and acts as the city's economic hub. Interestingly, the hills in western Dharan, which have experienced numerous landslides, are classified as a low risk zone. While landslide risk is high in the western hills ( Figure 4A), low composite risk is explained by a stable geology in the upper region, low density of human settlements, relatively high levels of forest-cover, and relatively low levels of social vulnerability ( Figure 5). However, while the western hills boast these attributes that collectively assign it as a low risk zone, active landslides and expanding squatter settlement on the banks of the Sardu River could render this as a riskier area in the near future. Prospective considerations such as these should be incorporated to guide future development and risk management in Dharan.
The composite risk map ( Figure 6) should not be considered a final product on the status of multi-hazard risk and vulnerability in Dharan. Rather, our aim was to begin a conversation on how to assess better the dynamic nature of risk and vulnerability in Nepal, which requires a prima facie acknowledgement that simultaneous risk to multiple natural hazards is necessarily overlaid upon a canvas of diverse social, economic, political, and built environment factors. This argument especially holds true in the urban core of Dharan, which exhibits exposure to multiple natural hazards throughout the calendar year and a rapidly expanding population, but a social, economic, and built environment profile that helps to attenuate these risks. However, since physical and social phenomena are dynamic, especially in Dharan, continuous measurements of variables against these baselines are required in order to inform future projections of risk assessment in the study area. In the meantime, it is evident, after overlaying the administrative map, that the most recent and newly developing settlements tend to be classified as medium to high risk. This is particularly important in a country such as Nepal that has exhibits high levels of rural-to-urban migration. Dharan is a hotspot for migrants within Nepal, which means that these relatively vulnerable areas will likely become more populated, more developed, and therefore more risky and vulnerable in the future.
We set out to determine how geospatial technologies can inform multi-hazard risk in a data-poor environment, using the city of Dharan, Nepal as a case study. Our experience here echoes the work of others (Gr€ unthal et al. 2006;Schmidt et al. 2011;Kappes et al. 2012;Marzocchi et al. 2012;Hernandez 2014;Van Westen et al. 2014) in that one of the major challenges of modeling multi-hazard risk was the requirement of extensive amounts of data (Bathrellos et al. 2017;Skilodimou et al. 2019). However, such comprehensive data are often absent in the developing world (Barrantes 2018). Even if the data exist, they are often low quality, patchy in their spatial extent, and difficult to access due to government and institutional barriers. Our model partially addresses this gap by using publicly available geospatial and social data gathered through fieldwork. It is always a great methodological challenge to extract useful information from a large set of datasets for environmental decision making (Runfola et al. 2017). With the help of geospatial techniques and a decision support system, our model is able to represent hazard processes in a spatial platform and conduct analyses to compare them. Although our work evaluates each hazard independently, a consideration of relation and interaction between hazards-and probability of the occurrence of events leading to potential cascade effects-constitutes important work. This is particularly important in the Nepal Himalaya, where the occurrence of landslides due to a combination of natural (monsoonal rainfall and earthquakes) as well as anthropogenic processes (construction of roads) is prevalent.
Although our model was able to identify risky areas in the study area, it demonstrated many challenges. Typically, variables for such a social and physical data intensive model are available in various temporal and spatial scales across the units of analysis and also exist in diverse data formats (Wilhelmi and Morss 2013). In this model, social data were aggregated from the household scale and represented in discrete administrative boundaries (i.e., the ward level), while physical data (e.g., elevation, rainfall) exist upon continuous surfaces at different time scales. Additionally, local level geological data are not available for recent time periods and are measured in a different datum and coordinates than other recent spatial data, such as high resolution WorldView-3 images. Our model was also limited in other ways. For example, it could not consider the cascading effects of a hazard as such analyses require long term and often historical records of the type, magnitude, and intensity of hazards. Available historical records on hazards from the city office of Dharan report only the total sum of money sanctioned for disaster mitigation and/or compensation provided to households after a disaster in a particular year. While these data offer insight, they were not able to be effectively incorporated in the model and their temporal extent was too brief to extrapolate in a robust and meaningful way.
The outcomes of this study and methods herein have utility in urban land use planning, allocation of scarce resources, and, most importantly, informed risk and policy planning at the local levels. GIS-based models, such as the one presented here, are promising in data scarce regions such as Nepal and neighboring countries that make up the Himalayan belt. In general, these countries are isolated, relatively poor, confront numerous hazards, and do not have comprehensive recordkeeping systems. By overlapping the current administrative map, it is evident that the newly developing settlement areas are under medium to high risk.

Conclusion
This study offers a composite risk and vulnerability assessment model for the data scarce city of Dharan, Nepal, based on geospatial data and techniques in combination with social data. This was accomplished by assessing and combining risk from three natural hazards (landslides, floods, and earthquakes) and subsequently developing a multi-hazard risk model that collectively considers the three types of hazards. The multi-hazard risk model was then merged with a social vulnerability analysis to produce a composite risk and vulnerability map. The model is primarily based on publicly available remote sensing imagery and socioeconomic data collected through fieldwork and from key informants in the local government. This method can be replicated in other data scarce regions that are at risk to multiple hazards. Additionally, results from this approach can assist decision makers to better understand comprehensive risk and, consequently, to design more effective and spatiallytargeted policies to increase capacity and resilience.
Risk is dynamic as changes occur due to geophysical processes, human activities, and their complex interactions. Hence, risk assessment should be measured routinely to capture changes in the natural and social environments and thus to delineate high risk regions, vulnerable populations, and monitor changes over time. One of the major challenges to risk quantification is inadequate long-term, historical information on hazard occurrence. This can partially be addressed with the help of satellite imagery, but the process can be tedious and sometimes insufficient to accurately model risk. Still, with the assistance of geospatial data and techniques, a risk assessment of multiple hazards can be conducted and utilized to inform evidence-based decision making.