Landslide Hazard Zonation and Evaluation Around Debre Markos Town, NW Ethiopia—A GIS-Based Statistical Approach


 The main goal of this research was to perform a landslide hazard zonation and evaluation around Debre Markos town, North West Ethiopia, found about 300 km from the capital city Addis Ababa. To achieve the aim, a GIS-based probabilistic statistical technique was used to rate the governing factors, followed by geoprocessing in the GIS setting to produce the landslide hazard zonation map. In this research, eight internal causative and external triggering factors were selected: slope material (lithology and soil mass), elevation, aspect, slope, land use land cover, curvature, distance to fault, and distance to drainage. Data were collected from field mapping, secondary maps, and digital elevation models. Systematic and detailed fieldwork had been done for image interpretation and inventory mapping. Accordingly, the past landslides map of the research area was prepared. All influencing factors were statistically analyzed to determine their relationship to previous landslides. The results revealed that 17.15% (40.60 km2), 25.53% (60.45 km2), 28.04% (66.39 km2), 18.93% (44.83 km2), and 10.36% (24.54 km2) of the research area falls under no hazard, low hazard, moderate hazard, high hazard, and very high hazard respectively. The validation of the landslide hazard zonation map reveals that 1%, 2%, 3%, and 94% of past landslides fall in no hazard zone, low hazard, moderate hazard zone, and high hazard or very high hazard zones respectively. The validation of the landslide hazard zonation map thus, it has been adequately demonstrated that the adopted approach has produced acceptable results. The defined hazard zones can practically be utilized for land management and infrastructure construction in the study area.

Introduction geometry, and groun dwater are the main intrinsic parameters, while excessive rainfall, earthquakes, and man-made activities can be considered as the main causative factors for landslides.
Landslide hazard zoning (LHZ) was carried out according to various approaches. Such approaches can be categorized into three main classes; deterministic method, statistical approaches, and expert evaluation method (Karaman et al., 2006b;Fall et a., 2006;Dai and Lee, 2001;Du et al., 2017;Westen, 1993;Westen et al., 2003). Selection of approaches for landslide hazard zoning depends on the size of the analysis to be performed, the total coverage area, expertise, and skill set of evaluators, the geological or geomorphic parameters or methods used to produce parameter data (Fall et al., 2006;Carrara et al., 1992;Ermias et al., 2017).
The main aim of this research was to produce a landslide hazard zonation (LHZ) map for the area. Thus, a GIS-based landslide hazard evaluation and zonation method was used to produce the LHZ map.

Study area
The current research area is situated in the east-central region of the North-Western Ethiopian Plateau (NWP), in the East Gojjam district of the Amhara Regional State. The area covers small parts of Gozamin, Aneded, and Baso Liben districts as shown in Fig.1. It covers an area of nearly 236.8 km 2 . The study area can be accessed via Addis Ababa to Debre Werek main Asphalt road which is about 300 Km and then 27 km gravel road towards Yejube town. The area is generally not easily accessible by vehicles, and there are no standard routes to reach throughout all parts of the research area. The type of soil and its landscape is di cult to go through various regions of the research area.

Regional geology
The research area is situated in the Northwestern Plateau, which is distinguished by the geology of the Abay basin. The basin comprises of Precambrian basement, Paleozoic and Mesozoic sedimentary rocks, Tertiary to Quaternary volcanic rocks, and Quaternary deposits (Russo et al.,1994;Assefa, 1991;Lebenie and Bussert, 2009;Lebenie, 2010;Gani et al, 2008;Gani and Abdelsalam, 2006;Mogessie et al., 2002;Ahmed, 2008Ahmed, , 2009Ahmed, , 2010Ismail and Abdelsalam, 2012;Braathen et al. 2001). According to the regional map prepared by Tefera et al. (1996), the present area comprises different geological formations including Tarmaber Gussa, Ashange formation, Adigrat formation. The Tarmaber Gussa formation is distinguished by alkaline to transitional basalts that often form shield volcanoes and contain minor trachyte and phonolite ows, as well as Late Proterozoic Ultrama c rock composed of serpentinite, peridotite, dunite, and talc schists. The Ashangi formation is distinguished by strongly weathered alkaline and transitional basalt ows with unusual tuff intercalation, which is often tilted (contains Akobo basalts of SW Ethiopia). Adigrat formation is characterized by Triassic-Middle Jurassic sandstone.

Local geology
Page 4/26 Field surveys and secondary data sources were used to map the geology of the present research area. It consists of basalt, sandstone, alluvial deposits, colluvial deposits, and residual soils. The lithologies are visible throughout road cuts and on natural hillsides. The basalt unit is affected by a slightly to complete degree of weathering.

Methodology
This research was divided into three major phases: data collection, data processing and manipulation, analysis, and the development of an LHZ chart. The data collection phase of the study focuses on collecting existing topographical and geological maps, as well as conducting eldwork to visually examine the lithologic units and develop a landslide inventory map.
The probabilistic statistical approach was used for landslide hazard zonation. This approach can provide practical interactions between different triggers and landslide distribution in the past and present. It consists of a landslide inventory mapping accompanied by the production of a predictive hazard model based on eight triggering factors and their relationships with previous landslides (Hamza & Raghuvanshi, 2017). The causative factors are; slope material (soil and rock), slope, aspect, elevation, ground water surface traces, curvature, distance to streams, and distance to roads. Finally, the landslide zonation map was produced based on the relative in uence of various causative factors.
The statistical approach is used within each causative factor map and its parameter map classes to determine the densities of landslides (Wahono, 2010;Hamza & Raghuvanshi, 2016). Later, depending on the class distribution and landslide density, corresponding weights can be developed (Hamza & Raghuvanshi, 2017). The quantitative prediction was made using the density ratio of past landslides with the various factor classes. This helps to rate the various causatives that are responsible for the reoccurrence of landslides with similar conditions in the area (Shahabi et al., 2013;Hamza & Raghuvanshi, 2016). As a result, factor maps were integrated with the derived weights to obtain a landslide hazard zonation map.
The system followed in this research is based on the expectation that the causatives responsible were considered. The number of grid cells and a hazard index value can be used to determine these causative factors. This can lead directly to a landslide hazard assessment. The ratio of landslides that 'did' occur to landslides that 'did not occur' is the hazard index value of all factors (Chimidi et al., 2017). Later, based on the values of hazard indexes, equivalent weight was assigned for all factors. Then, using all six responsible factor layers, generic trial combinations were developed. Thus, landslide hazard evaluation was computed using the optimal combination of internal causative and external triggering factors. The actual data obtained from previous landslides was used to verify the developed landslide hazard zonation map of the research area.

Data collection and analysis
In this research, primary and secondary data were used. Secondary data sources include DEM, satellite images, topographical maps, and meteorological data. The primary goal of the eld investigation was to collect relevant information on the previous landslide and to validate the different factor maps developed.
The controlling factors including slope angle, slope aspect, and elevation were computed using a DEM obtained from the ASTER satellite. Lithological information was obtained from the geological map of Ethiopia (scale 1: 2.000,000), while Landsat + ETM data can be used to process land use land cover and soil cover.

Landslide inventory mapping
A landslide inventory map is important to landslide hazard mapping. It is frequently accepted that the conditions that caused previous landslides can cause landslides again if they occur elsewhere in the given area (Raghuvanshi et al.,2015;Lan et al., 2004;Dai et al., 2002).
In the present research, past landslide and no past landslide areas are delineated through extensive eld works. Accordingly, the Landslide inventory map was created by overlaying eld GPS failure point data on the areal photo and DEM_30m obtained from the USGS data source by using Arc GIS 10.5 (Fig.2). All maps in this study were produced at a scale of 1:50,000. During eld work, most past landslides were observed following fault structures (Fig.2). Besides, there are also a signi cant number of landslide activities following rives and valleys.

Evaluation of causative factors
The frequency of landslides is affected by the interaction of several factors, including topographic, hydrological, and geological factors; therefore, selecting the causative factors is regarded as a critical issue in LHZ mapping (Dou et al., 2015) In this research, based on eld observation and contribution to previous landslides, eight in uential causative factors were selected. These include lithology and soil mass, slope, aspect, elevation, curvature, land use land cover, and distance from the stream, and distance from lineaments/fault. These factors were identi ed based on eld observations of previous landslides and their contribution to slope instability. The lithological map of the research area was updated from the Geological Map produced by Tefera et al. (1996) and through extensive eld investigations (Fig. 3A). Past landslide distribution shows that 10.5% of the landslides occurred in Adigrat formation,1.8% of the landslides occurred in Ashangi formation, 1.6% of the landslides occurred in Tarmaber formation, 23.7% of the landslides occurred in alluvial soil, 51.5% of the landslides occurred in colluvial soil, 10.9% of the landslides occurred in residual soil.

Elevation
Higher elevation areas have a higher possibility of landslide occurrence than lower elevation areas (Kailas et al., 2018). Lower elevation terrain is usually gentle, and a high water table is needed to trigger slope failure (Kailas et al., 2018).

Aspect
The direction of the slope in respect to the sunlight has a huge effect on its structural and chemical composition of the soil slope ( Accordingly, 0.02% of landsides occurred in Flat surfaces, 6.5% of landsides occurred in slopes facing North, 6.1% of landsides occurred in slopes facing Northeast,3% of landsides occurred in slopes facing East, 3.6 % Southeast, 6.5% of landsides occurred in slopes facing South,13.7% of landsides occurred in slopes facing Southwest, 24.1% of landsides occurred in slopes facing West, 27.1% of landsides occurred in slopes facing Northwest and the remaining 9.4% of landsides occurred in slopes facing North direction.
The distribution of previous landslides in the area concerning aspect shows that the majority of landslides concentrated on slopes inclined toward the West and Northwest directions.

Land Use Land Cover
The strength of slope materials against sliding and the regulation of slope water content are affected by land use land cover (Bchari et al., 2019). Furthermore, plant roots strengthen the slope and are commonly regarded as reinforcements (Bchari et al., 2019). Land cover absorbs soil water and decreases the probability of a landslide (Bchari et al., 2019).
Secondary data and the Landsat image collected from the USGS data source were used to develop a land use land cover map. Then, the map was modi ed and updated based on visual eld observation, and the nal map was produced (Fig. 3E) in Arc GIS. Later, to assign SSEP ratings for the individual facet, facet wise percentage area coverage of land use -the land cover was generated by geoprocessing in the GIS environment.
The occurrence of landslides within different land use land cover categories reveals (Fig. 3E) that 1.4% of landslides occurred within the forest, 0.8% of landslides occurred within grazing land, 22.6% of landslides occurred within bush land, 24.2% of landslides occurred within barren land and the remaining 51% were recorded within cultivated land. DEM was used to develop a curvature map of the research area and this was divided into three main classes: concave (0), zero curvature ( at (0)), and positive curvature (convex (>0)) (Fig. 3F). The slope's concavity is shown by a negative curvature value, while the slope's convexity is showed by a positive curvature value. High convexity and concavity cause drainage concentration across space, resulting in slope saturation and instability (Mandal & Mandal, 2018). Accordingly, in the present study, 49.2% of landsides occurred in concave slope (<0), 1.1% of landsides occurred in at slope (0) and 49.6% of landsides occurred in convex slope (>0).

Distance to drainage
The saturation degree of the slope's material is a signi cant parameter that governs its stability (Anis et al., 2019;Yalcin, 2008). This other important factor in slope stability is the slope's proximity to drainage systems (Yalcin, 2008). River water is considered the main source of soil moisture (Nohani et al., 2019). Proximity to drainage lines of intensive gully erosion is a crucial factor controlling the occurrence of landslides (Foumelis et al., 2018) Similar to other causative factors, DEM was used to create a drainage network, and this was computed in ArcGIS using Euclidean distance. Then, the distance to drainage was classi ed into ve classes; 0-298m, 298-604m, 604-902m, 902-1249m, and 1249-2055m (Fig. 3G).

Distance to fault
Geological structures including folds, joints, bedding, and shear zones faults have a greater impact on slope stability (Michael & Samanta, 2016). The proximity of a slope to these features has a signi cant impact on its stability, increasing the probability of landslides occurring (Michael & Samanta, 2016). The number of joints in a rock reduces its strength, which increases the distance to faults/lineaments. The fault map was created using a Google Earth image, a geological map, and DEM 30m data from the USGS. The distance to fault was generated by Euclidean distance in ArcGIS, then reclassi ed into six classes; 0-500m, 500-2000m, 2000-2500m, 2500-5000m, 5000-6500m, and > 6500m (Fig. 3H). The distribution of landslides from past landslide inventory data revealed that 77.9%, 20%, 0.9%, 0.1%, 1%, and 0.1% of the landslides occurred in the fault distance class between 0-500m, 500-2000m, 2000-2500m, 2500--5000m, 5000-6500m and greater than 6500m respectively.

Landslide hazard evaluation
In this research, a bivariate statistical system was applied to develop an LHZ map. The densities of each factor were used to assign every factor. These were statistically integrated to produce the LHZ map of the research area (Mandal & Mandal, 2018;Chimidi et al., 2017;Suzen and Doyuran, 2004;Hamza & Raghuvanshi, 2017;Chimidi et al., 2017). The bivariate analyses are important to quantify the abundances of the landslide within each percentile class (Süzen & Doyuran, 2004). In each causative factor class, hazard index value' was used, obtained from the ratio of landslides that 'did' occur to landslides that 'didn't occur within each factor class. Thus, hazard index values were used to determine hazard classes and the weight assigned to each factor.
The quantitative relationships were determined by processing landslide inventory maps and all factor maps in a GIS environment. Thus, The total pixel count was calculated for each factor class where the landslide occurred and did not occur as presented in Tables 2 and 3. Using ArcGIS's raster calculator, the total pixel count for the whole research area was 263,112, while the 80, 276-pixel count was covered by a landslide.

Data preparation
Landslide locations were investigated using extensive eld surveying and satellite image interpretation. The whole existing landslide was demarcated using GPS along the fringe, and then polygons were generated using Google Earth images. Finally, using polygons, a landslide boundaries map was developed. Later, this map was used for landslide data analysis in a GIS environment. Moreover, vector to raster conversion was performed to generate a 30 x 30 m pixel raster data set of past landslides. In addition, an Arc map was used to create a spatial database on the internal causative and external triggering factors such as slop material (lithology and soil mass), elevation, slope, aspect, curvature, land use land cover, distance from the fault/lineament, and distance from drainage. Various data layers that were used to develop landslide hazard zonation maps and for statistical analysis were presented in Table  1.

Landslide Hazard Zonation(LHZ)
In this research, slope material (lithology and soil mass), slope, aspect, elevation, curvature, land use land cover, distance from drainage, and distance from the fault were selected as causative factors. It was believed that these causal factors were most likely to account for the landslides in the area.
The hazard map for each of the eight causal factors was prepared using a weight equal to '1.' (Table 3).
The possible explanation for assigning equal weight to all different internal causative and external triggering factors was that the relative contribution of each parameter to slope failure can not be quanti ed. Therefore, it was assumed that each of these factors would have a relative contribution to slope failure. Using the ''raster calculator" tool and ratings, a landslide hazard zonation map of the area was produced.

Results And Discussion
The produced landslide hazard map (Fig. 4), showed that 17.15% (40.60 km2), 25.53% (60.45 km2), 28.04% (66.39 km2), 18.93% (44.83 km2) 10.36% (24.54 km2) area fall in the no hazard, low, moderate, high and very high hazard classes respectively. Consequently, a glance at Fig. 4 reveals that very high hazard (VHH) areas generally fall in southwestern parts of the areas, while southern and southwestern parts are demarcated as high hazard (HH) zones, with a scarce distribution in the eastern and northern part. In addition, moderate hazard (MH) zones are widely distributed in the eastern and northern regions, while Low hazard (LH) are distributed in the northwestern region of the research area. More dispersed Low hazards areas are also found in the central and eastern regions. No-hazard (NH) areas are mostly clustered in the north, with a sparse distribution in the west and east.
The prevalence of colluvial and alluvial soils, relatively gentle slopes (14-25), and elevation class 1821-1954m are all factors that contribute to the wide distribution of the Very high hazard(VHH) zone and high hazard (HH) zone in the southwestern area. Most likely landslides occur on slopes with an angle of 14-25, according to previous landslide inventory data, with 82 % of landslides occurring only on this slope class. The majority of prone slopes are covered with alluvial and colluvium soil. furthermore, nearly 68% of past landslides were distributed following the elevation class 1821-1954m. This indicates that this elevation class is the most vulnerable to slope failure.

LHZ map validation
The landslide hazard map was validated with the help of past landslide locations. This was computed by overlaying landslide location data and landslide hazard map produced in this research. According to the overlay analysis (Fig. 5), 70% of past landslides were found in the very high hazard zone and 24% are in the high hazard zone. Therefore, 94 % of past landslide locations were satisfactorily matched with the LHZ map. The remaining 3% of landslides are found in the moderate hazard zones, 2% are found in the low hazard zones, and 1% are found in the no hazard zones. Furthermore, the produced landslide hazard map was physically validated in the eld at several regions, most notably in the very high hazard and high hazard zones. Landslide activities were observed in areas that de ed very high hazard and high hazard zones.
Thus, it is reasonable to assume this landslide hazard map (LHZ) map adequately demarcated the spatial and potential landslide hazard zones of the research area.

Conclusion
The southwestern region of Debre Markos town including Gozamin, Aneded, and Baso Liben districts affected by several landslides every year resulting in damages to infrastructures and propertiesEight governing factors were prepared for this study: lithology and soil mass, slope, elevation, aspect, curvature, land use land cover, distance from lineaments/fault, and distance from drainage. An inventory map of landslides was created using a digital elevation model (DEM), satellite images, and extent eld investigation. Using the inventory map, the statistical relationship between governing factors and landslides was investigated, as well as their probable contribution to landslides in the area. In this research, the governing factors were rated using a GIS-based statistical and probability method, and the landslide hazard map was generated using a customized raster calculation. In this study, landslide hazard evaluation showed that 17.15% (40.60 km2), 25.53% (60.45 km2), 28.04% (66.39 km2), 18.93% (44.83 km2) 10.36% (24.54 km2) area fall in the no hazard, low, moderate, high and very high hazard classes respectively.  status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. Inventory Mapping. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. Landslide hazard zonation map. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.