Landslide susceptibility modelling in Nyahode and Buzi sub-catchments of Zimbabwe

Landslides are natural hazards that have severe and even disastrous consequences that include fatalities, diseases, construction and infrastructure failures and land damage (Lynett 2012; Uwihirwe et al. 2020). Landslides can lead to tree loss, forest fragmentation, land use/land cover changes to agricultural land, and produce loss and can exacerbate slope instability and collapse. Landslides can also affect the water quality of river channels and reservoirs by transporting large amounts of sediment, especially in landslide-prone areas. Landslides involve the downslope movements of debris, rocks, or earth material under the influence of the force of gravity. A landslide occurs when the pull of gravity exceeds the resistance force due to the destabilization of natural soil or rock slopes by natural and anthropogenic activities (Wubalem 2021). Slope destabilization is mostly the result of nonstop unsettling influence by human exercises such as the development of farming lands, deforestation, road cuttings, and poorly planned settlements among others (Shahabi & Hashim 2015; Elmoulat & Brahim 2018; Segoni et al. 2018). Landslides commonly occur on unstable slopes because of outside variables such as heavy precipitation, tremors, flooding, snow dissolving, stream disintegration, change in groundwater level, volcanic ejection, or any blend of these normal elements (Batar 2021). The recurrence, greatness, and volume of landslides are relied upon to increment by both natural and anthropogenic variables (Elmoulat & Brahim 2018). Besides, the recurrence, extent, and magnitude of landslides are further expanding because of climatic limits in delicate sloping or precipitous regions (Khan et al. 2019).

The vulnerability to landslides in developing countries is exacerbated by poor resilience, preparedness, and lack of financial support to deal with hydro-meteorological events (Lynn 2008). Actions to curb the effects of landslides in the area that is already failing or prone to rock and debris movement require identification of existing landslides, determination of the contribution of predominant causal factors, and development of landslide susceptibility maps (Mersha & Meten 2020). Landslide susceptibility is the probability that a landslide will occur in an area, depending on the nature of the topography (Tunas et al. 2020). Landslide susceptibility mapping is useful for the identification of areas that are at risk and for estimating the relative contribution of each factor to slope failure (Patricia & Zafra 2015).

Before the utilization of remote sensing data and the advancement of GIS tools, mapping susceptibility to landslides has been a difficult task because of lack data (Moresi et al. 2020). Although various approaches have been developed for landslide vulnerability mapping, they can generally be found in three categories qualitative, deterministic, and quantitative methodologies (Batar 2021). Studies utilizing these approaches of landslide susceptibility mapping have been done mostly at a regional scale (Devia et al. 2015; Elfriede et al. 2019). Few or no studies have been done at fine spatial scales to quantitatively map and model landslides in small catchment areas such as Nyahode and Buzi. Landslide susceptibility modelling at a fine spatial scale provides improved understanding of spatial variation of landslide causative factors at micro scale (Uwihirwe et al. 2020).

This study focused on Nyahode and Buzi sub-catchments located in eastern Zimbabwe in Save catchment and they drain to Mozambique, which was highly affected by Cyclone Idai 2019. As mentioned by Chatiza (2019), the frequency of hydro-meteorological hazards has increased in southern Africa, particularly in eastern Zimbabwe and western Mozambique. A series of cyclones have occurred in eastern Zimbabwe including cyclone Eline 2000, Cyclone Dineo 2017, Cyclone Idai 2019, and Cyclone Kenneth 2019. The 2019 Cyclone Idai highly affected Chimanimani and Chipinge districts (Chatiza 2019). It has been estimated that at least 634 died and more than 300 people remained missing 10 months after the event (Resilience 2020). This study focused on the Buzi sub-catchment where Chimanimani and Chipinge towns are located. The study area is characterized by settlements with high population density (Chipinge and Chimanimani town), intensive farming (banana cultivation), and steep mountainous terrain.

Despite the highly devastating landslides experienced in the Nyahode and Buzi sub-catchment, landslide susceptibility modelling has not been done for improved spatial resolution of 30 m and less, which is an imperative tool for landslide hazard management. Therefore, this research prepared a landslide inventory map using geospatial platforms in the data-scarce environment and evaluated landslide susceptibility using the geostatistical approach in a Geographical Information System environment. The frequency ratio model approach was implemented for developing a landslide susceptibility model. The same approach was implemented by several researchers, who include Batar (2021) in the Indian Himalayan region, Azeze (2020) in part of Abay Basin in, north-western Ethiopia, Wubalem (2021) in Uatzau Catchment area, north-western Ethiopia, Pirasteh & Li (2017) in the Iranian plateau and Elfriede et al. (2019) in Bogor, West Java, Indonesia. Landslide susceptibility zonation map is useful in policymaking and environmental impact assessment and landuse management (Chanza et al. 2020). It also provides comprehensive information that helps in reducing landslide vulnerability. Hence this study enhances the successful building of comprehensive resilience to hydro-meteorological events.

The study focused on the Nyahode and Buzi sub-catchments in Save Catchment district (see Figure 1). The area is located in the Eastern Highlands of Zimbabwe on the border of Zimbabwe and Mozambique in Manicaland province. The study area is bounded between xMin, yMin 32.5349, −20.4211 and xMax, yMax 33.0644, −19.6492. The study area covers an area of 2,375 km² with the lowest altitude of 554 m at the Rusitu river outlet and the highest altitude of 1,792 m in the Chimanimani Mountains. Nyahode, Haroni, Rusitu, and Buzi are the major rivers in the study area. Rusitu and Nyahode intersect at the Kopa area in Rusitu Valley. The river system forms dendritic and parallel drainage patterns in Nyahode and Buzi River catchments (see Figure 2). The rivers are characterized by vertical and horizontal erosion leading to the formation of deep vertical valleys susceptible to landslide phenomenon. The area is characterized by different topography which includes gullies, ridges, and valleys due to the high activity of vertical erosion.

The area receives high rainfall which ranges from 1,500 to 2,000 mm/year and falls under Zimbabwe region one, which is a unique highland climatic zone. High rainfall events experienced in the area can trigger landslides. The average temperature varies from 15 to 18 °C, with the risk of frost during the winter season. The rainfall season starts in October up to April when high rainfall is received between November and March. Favourable rainfall supports a dense network of rivers and streams, most of which drain mountains and perennial wetlands. Long years of denudation have resulted in deep and steep valley slopes that explain the high velocity of rivers that flow throughout the year. This means that most rivers become malignant during the rainy season with high erosion forces.

Agriculture is highly specialized and diversified, characterized by a range of horticultural products (mainly potatoes, bananas, mangoes, sugar cane, among others), corn production (mainly for personal use), and plantations (tea, coffee, and wood). The main source of income for the agricultural population in the study area is small-scale banana cultivation.

The study area consists of three main lithologies, which are dolerite, limestone, and quartzite. Limestone is more dominating lithology with high coverage in Buzi catchment. The area is highly characterized by red loam soils in Nyahode catchment. The soils have high silt content; hence they have high moisture-holding capacity.

The methodology for landslide susceptibility modelling followed in this study was:

• 1.

  Landslide inventory creation

• 2.

  Landslide conditioning factors selection and analysis

• 3.

  Landslide susceptibility model development

• 4.

  Landslide susceptibility zonation

• 5.

  Landslide susceptibility model validation

Figure 2 represents the methodological flow process of how a landslide susceptibility zonation map was deduced. The final result is the validated landslide susceptibility zonation map, which is a tool for hydrometeorological disaster risk assessment in the study area.

The inputs datasets used in the development landslide susceptible model are shown in Table 1. Data type source and format are summarized in Table 1.

Detailed fieldwork together with Google Earth satellite image interpretation was carried out from 5 to 25 July 2021 to create a landslide inventory map (Figure 3). Coordinates of the accessible landslides were taken using a Garmin GPS. To estimate the landslide size, a tape measure was used to measure landslide length, width, and depths. The slope angle was measured using a clinometer. For the inaccessible landslides, the coordinate closer as possible to the landslide was taken together with the direction of the landslide location from the taken point. Information on landslide location was obtained from residents and from google earth. Using Google Earth Pro coordinates from the fieldwork were loaded, which guides the digitization of inaccessible landslides. The same approach was used by Lupiano et al. (2019) and Hao et al. (2018).

The digitized landslide inventory was exported from Google Earth Pro as .kml file. The .kml files were converted into shapefile using QGIS 2.18 to perform geometry analysis in a GIS environment. Landslides, area, maximum length, minimum width, and the total area covered were calculated using geometry tools in QGIS 2.18. The landslide inventory shapefile was randomly divided into two groups: 75% for model development and 25% for model validation. The relationship between each selected causative factor and landslide inventory was established in the GIS environment.

Elmoulat & Brahim (2018) illustrate that for the delineation of areas that are susceptible to landslides it is important to consider geology. This is because each geological unit has a different impact on triggering geomorphological processes. The lithology map for this study was created by digitizing the scanned georeferenced existent 1:50,000 geological map of Zimbabwe. The study area consists of three lithologies, which are dolerite, limestone, and quartz. Manual digitization was performed in QGIS 2.18. A shapefile for lithology in the study area was converted to a raster file to perform a frequency ratio statistical test with a landslide inventory map. The dominating lithology in the study area is limestone followed by dolerites, and gabbro and quartz are the least dominating.

It has been recognized that landuse change is one of the most important factors influencing the occurrence of such landslides across the globe (Mersha & Meten 2020). Most of the anthropogenic activities such as intensive farming, cultivation on the steep slope, deforestation, settlement development, and road construction initiate slope instability. Vegetation cover has an impact on slope stability. The root system network determines shearing resistance and the rate of erosion hence determining slope stability. Sentinel 2A 10 m resolution satellite image was downloaded from Earth Explorer for July 17, 2021. The boundary shapefile of the Nyahode and Buzi sub-basins was used to clip bands four (red band), three (green band), and two (blue band). Supervised image classification was performed in ILWIS 3.3 Academic using the maximum likelihood algorithm for eight land use/land cover sample set classes. At least 30 pixels were assigned as training samples for each class. The image was classified into eight land use classes, which are: bare-field, bare-land, cultivated area, plantation, rock outcrop, settlement, shrubs, and grassland and water.

A confusion matrix was used to test the accuracy of image classification. Randomly selected ground control points of different classes from the field data and Google Earth were used for accurate assessment of the classified image. The points were rasterized and crossed with the classified image in an ILWIS GIS environment to obtain a confusion matrix table. As To improve land use/land cover classification accuracy assessment more points were assigned in the sample class. The final classified land use/land cover was crossed with landslide inventory for frequency ratio calculation in each land cover class.

The elevation map was calculated using the DEM SRTM 30 m resolution. Five elevation classes were created and crossed with landslide inventory map frequency ratio statistical calculation.

NIR = near infrared band

RED = red band

The NDVI values range from +1 to −2 where high values represent vegetation and low values represent barren rock and ice.

The river system is another conditioning factor used in landslide susceptibility for streams that are responsible for undercutting and eroding slopes. Elmoulat & Brahim (2018) noted that fluvial erosion results in triggering landslides across the globe. For drainage network extraction DEM hydro-processing was performed in ILWIS GIS software following the method used by Gumindoga et al. (2017) in the Upper Manyame catchment of Zimbabwe. Major rivers in the study area such as Buzi, Nyahode, Rusitu, Haroni, and their tributaries were extracted. These extracted rivers were validated using a 1:50,000 hydrological map of Zimbabwe. The Distance operation in ILWIS 3.3 was used for stream channel buffering. The extracted stream channel was used as the source for distance calculation. The distance map was sliced into six classes and crossed with inventory training data for frequency ratio statistical calculation.

The slope aspect is the direction of the maximum values changing from each cell to its neighbourhood in the output raster. Values in the output raster for each cell indicate the compass direction faced with the surface at that location. The direction is measured from 0 to 360 degrees in the clockwise direction.

The slope map in degrees ASPECTD was sliced using the group domain directions, as well as the corresponding degrees, which are described in Table 2.

The slope aspect map was crossed with a landslide inventory map for frequency ration calculation.

Ln = natural log

TANB = slope angle in radians

The undefined area obtained from TWI was replaced by 0.001. The TWI map was sliced into five classes and crossed with a landslide inventory map for frequency ratio calculation.

Water holding capacity is determined by physical soil characteristics. The landslide behaviour is mostly influenced by fluctuating water content and stresses in the unsaturated zone. Physical soil characteristics have a significant effect on the mechanical properties of the slip zone. Soil grids data for clay, sand, and silt content grams per kilogram was downloaded from the FAO website. The soil maps were sliced into five classes and crossed with landslide inventory raster map for frequency ratio calculation.

The histogram of developed landslide susceptibility map was used in the creation of four landslide zones.

The produced susceptibility map from the model was validated by the receiver operating characteristic curve (ROC curve) the approach used by Batar (2021). The area under the curve (AUC) of the ROC method represents the quality of the probabilistic model (its ability to predict the occurrence or non-occurrence of an event). An AUC value close to 1 indicates high accuracy, and an AUC value close to 0.5 indicates inaccuracy. Therefore, the greater AUC values show a higher accuracy performance of the produced susceptibility map. The success rate shows how many landslides areas, which are used in the model creation, are successfully captured by the susceptibility map to represent the model efficiency The AUC values obtained from the susceptibility maps and testing landslides show the model's success rate.

Figure 4 shows the developed landslide inventory map, which reveals a high occurrence of landslides observed in the northern part of the study area near Chimanimani town. The area was highly affected by Cyclone Idai in 2019, especially sites closer to Ngangu residential area. The total spatial extent of the mapped landslides inventory covers an area of about 12.1 km² constituted by more than 300 landslide polygons. The landslide width stretches from one metre up to around 40 metres and some occur in series across the slopes, for instance in the Chimanimani mountains. The occurrence of landslides in series may be a result of the vibration of earth, which triggers the occurrence of new landslides. The debris flow moves downslope up to a distance of 100 metres depending on the slope angle. The flow path of debris flow mostly follows stream channels were changes in stream morphometry and in some instances high distractive new flow paths not along the stream can be created. This was observed in Skyline mountains and areas along Nyahode river and Rusitu Valley.

In the development of the landslide susceptibility model, the FR was used to determine the level of correlation between the landslide inventory and the landslide conditioning factors. Landslide conditioning factors used in this research are slope degree, slope aspect, altitude, lithology, land use/land cover, distance from rivers, NDVI, TWI, and soil physical characteristics (silt, clay, and sand content). The results are shown in Figures 5 and 6 for each classified landslide causative factor, The relative frequency (RF) for each class within landslide causative factors is shown in Table 3 together with the overall prediction rate (PR) used as the weight of each causative factor in landslide susceptibility model.

As shown in Table 3 for the lithology conditioning factor, the quartzite lithology class has the highest correlations with landslide susceptibility with RF of 40. From the lithology map used, quartzite rock is located at a high elevation. For some filed observation in the study area along Rusitu Valley quartzite rocks are not balanced, hence highly susceptible to slope failure. Dolerites and gabbros lithology has a second relative frequency of 33. This lithology has the second predominance in the study area following limestone. Although dolerites and gabbros are located on steep slopes, they are well stabilized because of vertical and horizontal joints; hence they have a low relative frequency compared to quartzite lithology, which is not well balanced. There are several landslides located on limestone and shale lithology class; however, the RF is low. Although there is a high number of landslides located in the limestone lithology class it is the lowest landslide predictor in the study area. The limestone lithology class has the highest spatial extent compared to other lithologies in the study area (see Table 3) hence the number of landslides mapped is low concerning the surface area covered with other lithologies. Figure 6 (map A) shows that limestone lithology is highly dominant in low-lying and gentle slope areas with certain small coverage in high elevation (refer also to Table 3) for lithology conditioning map.

The most dominated land cover in the study area is shrubs and grassland and water is the least dominating class. As noticed in Table 3, the value of the RF for land use/land cover of bare land, rock outcrop, settlement, shrubs grassland, and plantation are greater than 10. The bare land class has the highest RF of 33, indicating high landslide probability in this land use/land cover class. Increased gully soil erosion in bare land resulting moisture dynamic and reduced shear strength of soil material and a high correlation with slope stability. Although rock outcrop and settlement land used land cover have moderate RF since these activities are being practised on gentle slopes, they influence landslide susceptibility prediction.

Landslide susceptibility decrease as the distance from the river increase due to the reduction in slope angle. At a distance of about 300 m–600 m RF starts to increase. RF is high again from 300 m to 100 m since it's the average distance to the river catchment divide. This is because of the effects of slope modification, gully erosion riverbank erosion, and river undercutting. The relation between TWI and landslide probabilities showed that the lower classes have the highest FR. The variations in slope angles observed in the Nyahode and Buzi catchments range from zero degrees, which correspond to the flatter slope, and 54 degrees and above for the steeper ones. As the slope angle increases, the RF increases; hence slope stability is depending on slope steepness. As the slope angle increases, the shear stress in the soil or other unconsolidated material generally increases. Gentle slopes are expected to have a low frequency of landslides because of the generally lower shear strength.

Table 3 shows the summary of the frequency ratio of a landslide of the class of each landslide causative factor and the overall PR of each causative factor.

Figure 7 shows the PR of each landslide causative map. The elevation map had the highest prediction factor, whereas distance from the river is the least predictor.

The final landslide susceptibility model for each pixel was obtained by summing overlapping pixels of the input conditioning factor multiplied by their PR (see Equation (8)). A landslide susceptibility map was created as shown on Figure 8.

The produced landslide susceptibility map was classified into four classes low, moderate, high, and very high-class (see Figure 10) The zonation of the landslide was done using the histogram of the continuous landslide susceptibility map. Nyahode Sub-catchment where Chimanimani town is located is dominated by a very high landslide zone and Buzi sub-catchment where Chipinge town is located is dominated by a low landslide susceptibility zone. This left the northern part of the study area Nyahode sub-catchment highly susceptible to landslides compared to the Buzi sub catchment. High landslide zones have been modelled in high elevation zones since it has high PR. Also, the model is highly sensitive to the causative factor with high PR for example Elevation, silt content, and slope degree. Figure 9 shows the spatial extent of the landslide zone increase with a decrease in landslide susceptibility.

The landslide susceptibility was divided into four classes as shown in Figure 10. The very high and high landslide susceptibility zones are highly dominated in the Nyahode river headwater. The area covered with these two has a high elevation steep slope and dip red loam soils. the drainage density in these is high which promotes erosion and slope instability. The low and moderate landslide susceptibility zone covers the highest area of the study area and is highly dominant in the Buzi catchment. The Buzi subzone has lower elevation with gentle slopes so the landslide conditioning factor classes have low prediction frequency in this part of the study area. Across the study area, the low landslide susceptibility zone is dominated by major rivers. However, some parts belong to the moderate landslide susceptibility zone. This is mostly the result of steep slopes in the valley, for example downstream on Nyahode river and along Rusitu valley.

The validation results in Figure 11 showed a frequency success of 0.80 and an FR PR of 0.82. These results reflect well on the landslide model for they are towards the perfect score value of the model.

This study produced a landslide susceptibility map based on statistical analysis, GIS methods, and remotely sensed together with field observed data. Landslides occur in a certain location within a given period and the intensity differs from place to place and time to time. The results of this study are important for the development of landslide early warning systems in light of the increased frequency of tropical cyclones in the southern African region. The main findings relating to each objective are:

• i.

  This study mapped landslide scars cover an area of about 12.1 km² where over 65% of the mapped landslides are located in Nyahode sub-catchment and others in Buzi sub-catchment A total of 385 polygon landslides were mapped in high elevated slopy areas and areas closer to the stream bank, mostly the northern part of the study area where Chimanimani town is located. The size of the landslide is highly related to the slope angle and landuse/landcover.

• ii.

  After assessing the performance of the landslide model with ROC curve the geostatistical approach shows reliability in landslide susceptibility modeling with AUC of 0.85 for frequency success rate and 0.82 for FR PR. Nyahode river sub-catchments are highly susceptible to landslides compared to Buzi river sub-catchments.

• iii.

  The final outputs of the landslide susceptibility map provide cheaper, faster, and more reliable results that are useful in land planning and management that reduce the impact of such hydrometeorological hazards. The results are useful decision-making supporting tools that can be used by government and non-governmental organizations in planning and enhancing resilience building to hydrometeorological hazards.

• 1.

  There is a need for geophysical surveys to assess slope stability in very high landslide susceptibility zones.

• 2.

  It is highly recommended that the master plans of Chimanimani and Chipinge towns should take note of landslide susceptibility areas provided in this study.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.