Landslide hazard, susceptibility and risk assessment (HSRA) based on remote sensing and GIS data models: a case study of Muzaffarabad Pakistan

The notion of this research is based on the two devastating earthquake events that happened on October 8, 2005, and September 24, 2019, in the regions of Azad Kashmir and Muzaffarabad. This study aims to (i) identification of the susceptible zones where landslides can occur in the future; (ii) preparation of landslide inventory maps using vector data, satellite imagery, Shuttle Radar Topographic Mission (STRM) and Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER) DEM; (iii) implementation of Analytical Hierarchy Process (AHP) model using weighted overlay analysis (WOA). For this purpose, key factors such as land use, faults, slope, contours, soil, and seismology maps are used to develop a landslide hazard zonation map. The output landslide susceptibility map has four susceptibility levels such as low, medium, high, and very high vulnerable zones. The results indicated that a highly susceptible landslide zone is found in the northwestern part of Muzaffarabad, which is a metropolitan region. Moreover, there are 127 active landslides are identified and collectively about 9% of the study area is very highly susceptible to future landslides. Furthermore, research findings are helpful in tactful thinking for future infrastructure development, and ecological protection in high-susceptible landslide regions in Muzaffarabad. It also helps the Government to make strategies based on any specific zones on a priority basis to reduce the casualties and destruction in future landslide events.


Introduction
Landslides are significant geological hazards that are often difficult to predict, similar to other natural disasters such as floods, earthquakes, and typhoons because they are uncertain with potentially negative consequences. Landslide is a physical occurrence that involves an extensive range of ground motions, including rock falls, slope failure, and debris flow. It refers to the downhill sliding of large amounts of landmass over steep slopes, and it can be quick or gradual. Landslides occur due to slope instability and it includes all types of slope failures in earth materials caused by gravity (Ha et al. 2022;Ram and Gupta 2022). The slope uncertainty can be induced by various factors, for example, earthquakes, flash floods, volcanic eruptions, and excess precipitation. It can also occur due to human activity such as mining, grading, and terrain cutting. In general, a landslide can travel by 3.260 feet per second (Pudasaini and Krautblatter 2021). Landslides also affect the world's ecological settings (Sim et al. 2022), it surplus sediments that can damage streams and waterbodies. Moreover, Landslides can have an impact on key ecological factors such as temperature, soil, and water. Landslides cause forests deterioration, animal habitats vanish and remove productive soils from slopes in terms of their ecological properties. Earthquakes are the most common source of landslides (Yan et al. 2021). The Azad Kashmir region is extremely prone to landslides, which are mainly caused by earthquakes. Therefore, it is essential to research in this region.
Historically, numerous land sliding incidents happened in Azad Kashmir but two devastating events that occurred in the recent two decades caught the author's attention. The Kashmir earthquake occurred on October 8, 2005, and it is Extended author information available on the last page of the article recorded as the second most devastating earthquake in the history of Pakistan.It had an intensity of 7.6 M w magnitudes and triggered thousands of landslides, mostly rock falls and rock slides, in the central area near the cities of Muzaffarabad and Balakot, Pakistan (Huang et al. 2021). The land sliding induced by the Kashmir earthquake caused numerous casualties and destruction to infrastructure. The Mirpur, Azad Kashmir earthquake that happened on September 24, 2019, with a magnitude of 5.6 M w , also triggered several hundred landslides.
Presently, landslides can be monitored by enhanced mapping techniques and other systematic measures such as remote sensing satellite products and GIS expertise. Remote sensing and Geographic Information Systems (GIS) technologies with their excellent spatial data process efficiency have attracted great attention in natural disaster assessment applications (Osako 2021). In terms of spatial landslide analysis, GIS is an important tool to carry out remote sensing data input, handling, visualization, combination, query, and analysis. GIS is an essential tool for geohazards assessment. The significant capability of GIS is to interoperate the dynamic spatial data of ground elevation, soil properties, and slope for landslide hazard assessment (Mahmoody Vanolya and Jelokhani-Niaraki 2021). Remote sensing and GIS were together used for landslide susceptibility assessment and the construction of a landslide prediction model (Wang et al. 2021a, b;Zhang et al. 2020).
There are various assessment models discussed by previous researchers conferred to landslide hazard assessment such as the Weights-of-evidence model (Sifa et al. 2020) and the susceptibility model based on cumulative frequency. The Weights-of-evidence model is one of the bivariate methods to calculate the weight for each landslide predicting factor that is dependent on the existence or absence of the landslide within the defined area (Bopche and Rege 2022).
Furthermore, different methods such as AHP (Panchal and Shrivastava 2022), Heuristic, Probabilistic, simple statistical bivariate, regression, and deterministic models have been used by the researchers. In a variety of case studies, meaningful findings were obtained using AHP (Meghanadh et al. 2022;Moragues et al. 2021;Roccati et al. 2021) and Weights-of-evidence models with high accuracy and reliability .
Muzaffarabad is a landslide-susceptible region due to its rugged terrain, topography, climatic conditions, location on an earthquake fault zone, geology, and geomorphology. This research work is aimed to identify the level of vulnerability and the potential landslide hazards sites in Muzaffarabad that will be helpful for future land sliding events. An AHP-based analysis is being adopted to compute landslide susceptibilities and their impacts with the progress of prior research and ongoing initiatives.
For this purpose, the existing system and literature have been reviewed thoroughly. The scope of the study can be evaluated by (i) review of previous studies and related published articles (ii) interpretation of satellite images and geological maps (iii) preparation of landslide maps inventory (iv) evaluation of hazard and risk assessment of identified landslide zones (v) Modernized Landslide hazard susceptibility map for Muzaffarabad. This research work is significant as it contributes to the identification of the landslide hazard zone map of Muzaffarabad for future use. It also provides a methodology for implementation of the AHP model using Weighted overlay analysis in ArcMap spatial analyst and model builder tools.
2 Study area and data

Study area
Muzaffarabad is the capital city of Azad Kashmir, located at the north-eastern border of the Khyber Paktoon Khawa (KPK) province with a geographical location of 34. 3500°N and 73.4667°E. The district of Muzaffarabad is situated at the convergence of the Jhelum and Neelum rivers, originating from Indian-occupied Kashmir. Both 2005, Muzaffarabad Kashmir earthquake (Farooq and Akram 2021) and 2019, Kashmir earthquake (Khan et al. 2021;Papathanassiou et al. 2022) happened along the Jhelum fault line, which is well-known for its high seismic activity in the South-Asian region. This seismic barrier was formed by the movement of the Indian, Arabian, and Eurasian Plates. Northern Pakistan and India still have active fault lines that can cause displacement of crustal material and landslides in the future. The location of the study area on the map is shown in Fig. 1 below. In terms of precipitation and temperature settings, Muzaffarabad has a 48.93 average annual rainfall mm (inches), an average maximum temperature of 22°, and a minimum temperature of -3°. According to (Burgan and Aksoy 2022) it is essential to analyze the precipitation metrics that influence the underlying mechanism of landslides.

Dataset overview
The basic data needs are to prepare a variety of maps inventory for the landslide hazard zone map. But, some existing maps were available such as seismic, geology, topography and soil maps. These data sets and maps are compiled, digitized and re-classified according to the requirement of research work.
However, Landsat satellite imagery and multiple DEM are the key datasets used to produce land-use maps, slope maps, elevation maps, and contours maps. Furthermore, these maps are integrated with vector data (roads, fault lines) and precipitation data for an AHP-based synthesis, which is then utilized to create a final landslide hazard susceptible map. Following Table 1 is showing the respective platforms for all satellite datasets used in this research and is further explained below.

LANDSAT data
NASA/USGS provides the extensive and most continuous space-based Landsat satellite products worldwide (Ahmad et al. 2022). Especially, Landsat 5 data product is operated for almost 29 years (Hemati et al. 2021) and achieved a Guinness World Record for the longest continuously operated earth observation satellite. Landsat 8 is a cooperation of NASA and the US Geological Survey (USGS), and it includes instruments such as the Operational Land Imager (OLI) and the Thermal Infrared Sensor (TIRS). This research is being carried out using Landsat 8-OLI products for land-use/land cover change detection(Ghayour et al. 2021; Jamali 2021). The Landsat 8, imagery is being acquired from the USGS website on October 6, 2021, with ID: LC08_L1TP_150036_20211006_20211013_02_T1 and 5% cloud cover.

ASTER DEM
The Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER) is a 14-channel imaging instrument operated by NASA's Terra satellite since 1999. ASTER visible near-infrared (VNIR) instrument, has three bands and 15 m' instantaneous field of view (IFOV). ASTER GDEM is developed by stacking all individual cloud-masked scenes and non-cloud-masked scene DEM (Abrams et al. 2020). Then several algorithms are being applied to remove any abnormal data to finalize the end products.
Moreover, in August 2019, the ASTER Science Team released Version 3 of the global DEM (GDEM) with the help of stereo correlation of almost 1.8 million ASTER  image tiles. The DEM has 1 arc-second/30 m (Bakiev and Khasanov 2021) latitude and longitude postings. We used the ASTER DEM in TIF format with 34°Northing and 73°E asting for this research work. In addition, we obtained the Nasa LPDAAC collection, ASTER Level 1 T, using the Earth Explorer platform.

SRTM
The Shuttle Radar Topography Mission (SRTM), was launched on 11 February 2000 (Azam et al. 2021). STRM product covers digital elevation data for the earth's landmass between 60-degree north and 54-degree south latitude. It has a higher resolution in terms of elevation and is captured with a spatial resolution of 30 m (Carrera-Hernandez 2021). SRTM characterizes the world's highest quality, publicly accessible digital elevation models (DEMs). During an 11-day Space Shuttle flight, SRTM captures elevation data across 80 percent of the Earth's land area with a horizontal precision of 3 arcsec. In this study, we used USGS/NASA SRTM data of year 2021 which is acquired using row-51 and path-6 for Muzaffarabad area in Geo-tiff format.

Methodology
In this study, a landslide susceptibility model is developed by combining AHP and weighted-overlay analysis using model builder in ArcMap. This methodology presents a preliminary overview of potential landslide areas and is adopted to achieve the hazardous sites in the final landslide zonation map. These interpretations are based on slope angle, soil type, classified satellite images, seismic map, ASTER data and other GIS maps. Below Fig. 2a, b are showing the (i) landslide evaluation procedure, and (ii) overall workflow diagram and AHP implementation, for this research work.

Analytical hierarchy process (AHP) model
The AHP method was developed by Thomas L. Saaty (Sur et al. 2020) in the 1970s. It is a significant decision-making model, which has been applied in various research studies related to landslide susceptibility assessment (Du et al. 2020;Mao et al. 2022). The current study is achieved through the application of the AHP approach for fusing weights of the different factors for weighted overlay analysis. The AHP approach is commonly used in site selection and suitability analysis. It is founded on three principles: decomposition, comparative assessment, and priority synthesis. It provides essential knowledge to understand the procedures for occurring landslides and for preparing a susceptibility map.
Thus, the AHP method (Bahrami et al. 2021) is significant to generate maps in land-use planning and setting up safety measures. The analytical hierarchy process is a systematic technique for complex decisions and analyses based on criteria and alternatives. The Analytical Hierarchy Process model is implemented in the development of maps inventory by changing alternatives and criteria, after weighted overlay analysis for soil, slope, fault line, land use map, and other layers respectively.

Weighted overlay analysis
Weighted overlay analysis is one of the most common spatial GIS operations (Kumar et al. 2021), that integrates remote sensing data with other GIS (raster) data layers. The weighted overlay method (WOM) is available in ArcMap (spatial analyst) tool which helps to produce a susceptibility/prediction map (Saleem et al. 2020) based on user input. It is commonly used to solve multi-criteria problems such as optimal/hazardous site selection (Lau and Zawawi 2021;Senouci et al. 2021) or suitability modeling.
Many researchers have implemented weighted overlay analysis to produce a landslide susceptibility map (Devara et al. 2021). It is a method of creating an integrated analysis by applying a single scale of values to varied and distinct sources using visual or analytical operations (Tiwari and Ajmera 2021). In this research, the authors have designed, distributed the model into sub-models, and identify the input layers using weighted overlay analysis. The data layers are ranked with different weight scales for each criterion and assigned weights from 1 to 10, with 10 being the most and 1 being the least suitable.
Moreover, raster reclassification has been performed using the reclassify tool in ArcMap for layers such as slope, and aspect, then utilized as input for weighted overlay analysis. First, we assigned weights according to their importance and controlling factor then we determined all those layers using the AHP model in model builder. Table 2 is showing the eight targeted layers with their corresponding weight importance and resampling techniques. The cumulative weight of all input layers is usercontrolled. Following Eq. 1 (Ali et al. 2018) has been used on the backend to integrate all layers for weighted overlay analysis.
where Wi represents the weight of the ith factor, Sij represents the subclass weight of the jth factor and S is the final output map.

Potential landslide hazard zone map
All the above-reclassified maps are further used in Analytical Hierarchy Process by applying Arc Map Model builder. The weighted overlay analysis tool is used to assign weights to specific criteria i.e. slope of 45 Degree assigned the highest weight comparatively 15 and 30 degrees of Slope. Similar to slope the other parameters are also assigned weights concerning their hazardous factor and compositions. Figure 3 is illustrating the spatial data model (SDM) designed for the analysis. The Model Builder is a graphical environment that allows you to design a model to finish a complex geoprocessing task (Kahal et al. 2021;Polat 2021). Multiple tools such as buffer, join, plus, distance, raster calculate, reclassify and other mathematical operations can be added according to analysis requirements. Precipitation data is acquired from Pakistan metrological department from 2000 to 2020 and intergraded with spatial data layers. Annual rainfall data is interpolated to create a precipitation graph for each year and then the sum of overall rainfall is calculated. Following Fig. 4 is showing the sum of rainfall in Muzaffarabad city, with maximum precipitation in July and August.

Land use/land cover map
Land-use maps are usually made to reflect the possible uses of a ''unit'' of land, by human beings. The most common physical factors are vegetation, buildup area, water bodies, and slope. Generally, land-use/land cover can have a significant influence on landslide activity (Rabby et al. 2022). It can cause by human activities such as road construction and cultivation on steep slopes, forest logging and deforestation. In previous studies, experts used Maximum Likelihood (ML) supervised classifier for land-use classification using satellite imagery.  Thus, Landsat imagery is ortho-rectified and atmospherically corrected by using Geomatica software. The maximum likelihood classifier is frequently used to estimate learning parameters using training samples. As a result, maximum likelihood (Fang et al. 2020;Pham et al. 2021;Wang et al. 2021a, b) estimates may be generated for a wide range of estimating circumstances.
Then composite images are classified as shown in results Fig. 5e using Maximum Likelihood (ML) supervised classification in the ENVI (Sun et al. 2021) image processing software. The description of the classes is described in Table 3 which provides the details regarding the name of the classes and the number of signatures made for each category.

Digital elevation model and slope map
Digital elevation models (DEM) record information about factors that change over time, such as elevation and groundwater depth. DEMs are very useful to extract key parameters for susceptibility assessment (Brock et al. 2020;Chen et al. 2020;Long et al. 2022) such as elevation, slope, aspect, curvature, and watershed. A slope map depicts an area's topography, as well as an evaluation of topographic elements that have impacted and may continue to influence land development. SRTM and ASTER DEMs are used in this study to derive the slope maps as shown in Fig. 5d, f. SRTM DEM data is also used to create a contour map with an interval of 500 m. The majority portion of the study area lies between 1000 and 2500 m and the least area lies \ 500 m and a very small portion present at 4000 m of elevation as shown in Fig. 5h. Following Table 4 is presenting the composite area based slope-aspect and percentage.

Soil map
A soil map is a geographical description of soil types and soil characteristics (Hakim et al. 2022) such as soil pH, textures, organic matter, and depths of horizons, and is typically prepared through a soil survey. Soil maps are frequently used for the classification of land, spatial planning, agricultural development, environmental protection, and other soil-based projects. In past, soil maps indicate the overall distribution of soils types and are accompanied by a soil survey report. Now, digital soil mapping techniques are used to create innovative soil maps. The soil map of the study area is updated based on a soil map prepared by the Soil Survey of Pakistan. There are seven soil types are present in the study area, with a maximum count of 35,857 Rockland and 217 counts of Didal, explained in Table 5. The map was scanned, georeferenced and then the soil units were digitized as shown in Fig. 5b.

Seismic and fault line map
A seismic map is generally emphasizing the expected earthquake ground motions at the surface of the earth. These maps help in public safety consideration in terms of existing natural conditions and man-made structures. A seismic hazard map (Mori et al. 2020) helps to determine the relative velocity on a local, regional or national scale for any certain location. PMD, (Pakistan Meteorological Department, n.d.) is responsible for mapping the seismology and data is acquired from PMD for both seismic and fault line maps respectively. Fault Lines Map represents the lines that separate two tectonic plates of the Earth's crust. In this study, both faults and seismic maps are digitized using geological and seismic maps. Then multiple buffer polygons of 250 and 500 m are produced around the digitized fault lines shown in Fig. 5g. Then all the vectors layers are converted into raster to incorporate in the weighted overlay tool. Following Fig. 5 is presenting the overview of the satellite image, soil map, SRTM, ASTER, and their corresponding slopes (Figs. 6, 7).

Results
This study revealed that map inventories and susceptibility maps are important components of landslide identification and management. The factors such as rainfall, elevation, slope, erosion, geology, fault lines, and a variety of human causes directly or indirectly have an impact on landslides. The method used in this research for landslide hazard zonation is mostly based on the weights assigned to the parameters that cause a landslide. Therefore, in the weighted overlay analysis, indicators such as slope, aspect, elevation, soil map, seismic map, distance from faults, distance from highways, land-use, and geological map are used to generate the landslide susceptible map.
Based on the results, we found that the northwestern part of the study area is more susceptible to future landslides and about 127 active landslides are present in the study area. After creating a landslide inventory map to validate the results, the susceptibility map is compared with active landslides. Furthermore, the comparison map shown in Fig. 8 demonstrates the importance of the land-use indicators.
In addition to spatial extents of landslide vulnerability, it is observed that the very high and high hazard zones are mainly concentrated around the urban center of the study area. The majority of the activated landslides are found along roadways, which indicates that human interference, such as deforestation and construction, has an extensive influence on the stability of the slope. The risk index of future landslides in the study area is derived using Fig. 7   the final susceptibility map. The weights are assigned to each controlling factor using the AHP model. Table 6 shows the % area covered for susceptibility levels ranging from low to very high. The final output map Fig. 7 shows the area with low, medium, high, and very high susceptible levels. It is clear from the results that AHP and WOA are significant models for landslide hazard mapping and assessment. It also revealed that collectively 25% of the geographical distribution of the area lies in a very high susceptibility level that is highly prone to landslide occurrences. Therefore, the research on slope instability and susceptibility mapping are essential components of hazard management in decreasing the risk of future landslides. More advanced GIS and Remote sensing (artificial intelligence) based analysis can help more efficiently in understanding the mechanisms that cause landslides, as well as to prepare susceptibility maps.

Discussions
The AHP-based weighted overlay approach was employed in this research work to create a landslide susceptibility map for Muzaffarabad. A review of previously published articles aided us in determining the causative variables (Bounab et al. 2022) involved in the occurrence of landslides. Muzaffarabad was selected as a study area but analysis can be performed on all northern areas for the measurement of landslide hazard zonation and comparison to a large spatial extent. It can be discussed that map inventories and susceptibility maps are necessary factors for landslide identification and management. Landslides are difficult to forecast (Ahmed et al. 2021) since they are unpredictable and can have potential consequences similar to other natural disasters such as floods, earthquakes, and avalanches. For many years, governments and academic institutes throughout the world have tried to assess landslide hazard risk and map its spatial distribution. At present, landslides risk can be measured systematically and machine learning algorithms combined with conventional neural networks can help more efficiently (Aslam et al. 2021). It is also recommended that engineering geologists should determine the extent of possible unstable conditions, especially for very high susceptibility levels before future development. It also concludes that Muzaffarabad's deforestation and infrastructure developments are vulnerable to earthquake-caused landslides that can cause significant damage to both lives and properties.

Conclusions
This research work has applied the AHP approach using weighted overlay analysis to construct a landslide susceptibility map. The landslide hazards susceptible map reveals that 75% of the study area has low susceptibility, 11% has a medium, 5% has high and 9% area has very high susceptible zones respectively. It is concluded that collectively about 25% of the area is susceptible to future landslides in Muzaffarabad. The results also demonstrated that almost 127 active landslides are found in the study area. These ongoing landslides show a strong correspondence with the high and very high susceptibility class of the final map. Furthermore, the comparison map depicts that landslide zones with very high susceptibility are located in or near urban areas in the northwest direction of Muzaffarabad. The AHP methodology combined with weight overlay analysis has demonstrated itself to be a suitable tool for mapping landslide hazard zonation. The landslide susceptibility map created as part of the current research is an important step forward in the management of landslide   hazards in the Muzaffarabad region. Ecological factors such as water, temperature, and soil are also found to be strongly linked with landslides events. Landslides indirectly have an impact on soil erosion, water channels wedge, and temperature variations in the study area. Furthermore, the integration of some machine learning algorithms should be taken into account with the weighted overlay tool and AHP model for the advancement in future landslide mitigation. In the future, artificial intelligence and deep learning models must be adopted to map landslides more efficiently.
Funding It is declared that the authors have not received any funding for this work.

Declarations
Competing interests The authors have no conflict of interest.