Spatial landslide risk assessment in a highly populated Rohingya refugee settlement area of Cox’s Bazar, Bangladesh

ABSTRACT

In recent times, landslides have become a major concern in the southeast part of Bangladesh. This study aims to develop comprehensive landslide risk mapping by applying the analytical hierarchy process (AHP) and geospatial techniques in Ukhiya and Teknaf Upazilas, a highly populated Rohingya Refugee Settlement area. To assess the landslide risk, we selected 12 influencing criteria of hazard, vulnerability and exposure under the relevant components of risk, such as precipitation intensity, landslide inventory, distance to fault line, stream density, distance to stream network, elevation, aspect, slope, geology, normalized difference vegetation index (NDVI), land use and land cover (LULC), and population density. These spatial criteria were weighted using the AHP model. The weighted overlay techniques were then used to produce the landslide risk map. The findings demonstrate that 2.19% of the total area is classified as a very high risk zone, and 12.74% is categorized as a high risk zone. Moderate risk areas cover 23.08% of the total area. The risk map is validated by the landslides inventory. The outcomes can be used by the concerned authorities to take the necessary steps to reduce the impact of landslides.

KEYWORDS:

• Landslide risk
• analytical hierarchy process
• GIS
• remote sensing
• Bangladesh

Introduction

Landslide is a motion of landmass (earth or debris) or downward movement of rocks on a slope. This movement can be triggered by various external factors such as rainfall, earthquakes, variations in water level, storm waves, or rapid erosion of streams that can cause shear stress to either increase or decrease quickly in slope-forming materials (Cruden 1991). In hilly terrain areas around the world, landslides have resulted in numerous deaths and significant economic losses. It killed over 18,000 people and affected about 4.8 million people worldwide from 1998 to 2017 (UNISDR 2018). In comparison to other continents, Asia is one of the most prone and vulnerable areas to landslides, with the largest number of adverse accidents and deaths. In Asia, 3,456 disasters were recorded between 1970 and 2019, resulting in the deaths of 975,778 people and economic losses of US$ 1,204 billion (WMO 2020). India, Nepal, China, and Japan are considered the most vulnerable regions to landslide disasters in this continent (Kirschbaum et al. 2010; Petley, Dunning, and Rosser 2005; Wu 2017), and Bangladesh is the most recent addition to this list (CDMP-II 2012). Bangladesh is extremely vulnerable to natural disasters like cyclones, floods and storm surges due to its geographic position, terrain features, abundance of rivers, coastline morphology, and monsoon climate (Moniruzzaman 2012).

Landslide disasters have been common in the past in the heavily settled area of Bangladesh. However, in recent times, many landslides have resulted from intensified human activities, including hill cutting for slums and suburban housing construction (Ahmed 2017). In the southeast part of Bangladesh, landslides have become a significant concern since 2000 (Alam 2020). For example, a series of devastating landslide incidents resulted in over 725 deaths in the Chittagong Hill Districts from 2000 to 2018 (Sultana 2020). The impacts of the intense weather incidents (Khan et al. 2020), demographic influence in flatlands, acceleration in urban growth, changes in hill slopes, deforestation, and a lack of cultural awareness have been classified as emerging challenges to mitigate the loss due to landslides at the national level (Alam 2020).

Recently, Cox’s Bazar District has been categorized as a highly susceptible area to landslides. More specifically, Ukhiya and Teknaf Upazilas¹ are highly susceptible due to the changes in hill slopes for the development activities to accommodate a large number of Rohingya refugees. Since 25 August 2017, 655,000 Rohingya refugees had arrived in Cox’s Bazar District, and the number reached nearly 860,000 at the end of the year (Imtiaz 2018). According to the Inter Sector Coordination Group (2019) weekly report, 12 different landslides and erosions affected 98 refugees from 14 to 20 June, and four different landslides and erosions affected 40 refugees from 20 to 27 June 2019. A continuous 14 in. of rainfall for 72 h in July 2019 resulted in 26 landslides with one causality and 4,500 homeless in the refugee camps of Cox’s Bazar, Bangladesh. Moreover, in September 2019, 16 separate floods affected nearly 20,000 refugees, including 2,000 refugees who had been displaced by new landslides (Jonathan 2019).

It is necessary to identify the landslide risk zones to mitigate the probable loss in the area that contains a large number of local inhabitants as well as Rohingya refugees. As a result, landslide risk mapping can be a solution to reduce the damages (Regmi et al. 2014). The interaction of certain dangers that make humans, entities, and ecosystems vulnerable and exposed are described as risks (Hoque et al. 2019). The risk of a damaging landslide in a spatial (“where”) and temporal (“when”) sense, as well as magnitude (“how large”) of the occurrence, can be described as a landslide hazard (Bui et al. 2011). Vulnerability, on the other hand, refers to the degree of damage that can be incurred by hazards to people's possessions, livelihoods, and the environment (Hoque et al. 2017). In this study, the produced risk map provides crucial resources to help in sustainable construction and land use planning decisions. Furthermore, the successful use of these maps will significantly minimize landslide damage risk and other cost effects.

Qualitative maps are generated by adding weights to all variable landslide maps that rely on the expertise of the researcher, while quantitative methods include strategies to examine the relation between landslide criteria based on various models (Park et al. 2013). Based on factors connected to the consequences of landslides and locations where historical landslides have occurred, both qualitative and quantitative approaches are tested (Ahmed 2015; Park et al. 2013). In landslide risk mapping, different methods are used with the key goal of designing landslide hazard zoning maps (Bathrellos et al. 2012; Highland and Bobrowsky 2008). Geographic Information System (GIS) and Remote Sensing (RS) techniques have recently gained popularity in landslide mapping, and are widely used for landslide susceptibility (Ahmed 2015; El Jazouli, Barakat, and Khellouk 2019; Tehrani and Husken 2019), hazard evaluation (Bachri and Shresta 2010), risk analysis (Pathak 2016), as well as landslide management (Kavzoglu, Sahin, and Colkesen 2014; Talaei 2014). The Analytical Hierarchy Process (AHP) developed by Thomas L. Saaty in 1970 has been used by many researchers to assess landslide risk worldwide (Bachri and Shresta 2010; Boroumandi, Khamehchiyan, and Nikudel 2015; El Jazouli, Barakat, and Khellouk 2019). It incorporates multiple criteria for the risk components to develop a detailed spatial risk map by adding the weight and ranking those criteria on a hierarchical structure (Roy and Blaschke 2015).

Landslide-related studies in the southeast part of Bangladesh are mainly concentrated on the areas of Bandarban and Rangamati. Because of their high elevated hilly areas and frequent landslides during the monsoon, researchers most of the time concentrate on those areas while studying landslides. However, they do not take Ukhiya and Teknaf Upazilas into account for an extensive landslide risk assessment. Most of the studies that have been done in Ukhiya and Teknaf Upazilas are focused on the camps of Rohingya refugees (Ahmed et al. 2020a; Ahmed, Firoze, and Rahman 2020b; Tehrani and Husken 2019). One more study was also conducted in Cox’s Bazar Municipality area, following the AHP method (Ahmed 2015). The study aims to develop a risk scenario for the Cox’s Bazar Municipality, integrating landslide hazard, vulnerability and exposure.

A massive influx of Rohingya refugees occurred in the south-eastern region, mostly in Ukhiya and Teknaf Upazilas in 2017. As a result, a large number of forests have been degraded, which accelerated many landslide conditioning factors. Because of that, there might be a devastating landslide which can result in loss of life and property on a larger scale. In this study, a variety of criteria have been used to produce the landslide risk map more accurately. No other previous studies concentrated both on the host communities and refugees at the same time. Whereas most of them focused only on the refugee camps (Kamal, Ahmed et al. 2022a; Kamal, Hossain et al. 2022b), this study assessed far beyond that. This paper is focused on developing a comprehensive landslide risk map emphasizing both the highly populated Rohingya refugee camps and the local communities. In addition to a wide variety of factors to get a more accurate output, we have considered a historical landslide inventory dataset, which will make the research more scientific. It can be helpful for the researchers who will conduct future research in this area to identify and validate the risk zones in the host communities and the refugee camps. This study will give a clear view to the concerned authorities to take the necessary steps to mitigate the loss due to landslides.

Materials and methods

In this study, geospatial techniques, a combination of RS and GIS, and multi-criteria weighting technique AHP were used to assess the landslide risk in Ukhiya and Teknaf Upazilas, focusing on the refugee camp area. To assess the landslide risk, 12 conditioning factors have been considered (Figure 1), such as precipitation intensity, landslide inventory, distance to fault line, stream density, distance to stream network, elevation, aspect, slope, geology, NDVI, LULC, and population density. All data layers were projected to the Universal Transverse Mercator (UTM) zone 46 North system, using the Datum of World Geodetic System (WGS)-1984.

Figure 1. Flowchart of landslide risk assessment method used in this research.

Risk is typically a combination of hazard, exposure, and vulnerability. Therefore, risk assessment studies must consider how risk elements such as hazard, exposure, and vulnerability interrelate with each other (Hoque et al. 2019). In this study, the following equation (Eq. 1) was used for risk assessment in the study area. (1) $\mathrm{Risk}=\mathrm{Hazard}\times \mathrm{Vulnerability}(V)\times \mathrm{Exposure}(E)$(1)

Study area

This study focuses on the landslide risk assessment in the highly populated Rohingya refugee settlement area of Ukhiya and Teknaf Upazilas (excluding Saint-Martin Island), which is in the southeast part of Cox’s Bazar District, Bangladesh. Ukhiya and Teknaf Upazilas cover an area of 261.80 and 388.66 km², respectively, and their population are 207,379 and 264,389, respectively (BBS 2011). The settlement area is located between 20.867°N to 21.283°N latitude and 92.30°E to 92.10°E longitude, bounded by the Naf River and the Arakan State of Myanmar on the east, the Bay of Bengal on the southwest, and Ramu Upazila of Bangladesh on the north (Figure 2).

Figure 2. Study area, showing the administrative boundary of Ukhiya and Teknaf Upazilas and the Rohingya refugee camps located in the area.

This study is conducted in five Unions (local administrative area) of Ukhiya Upazila named Haldiapalong, Jaliapalong, Palongkhali, Rajapalong, and Ratnapalong and five Unions of Teknaf Upazila named Baharchhara, Nhilla, Sabrang, Teknaf, and Whykong. Both Upazilas are bordered by tropical evergreen and semi-evergreen forests and other vegetation (Rahman, Islam, and Chowdhury 2019). This area is considered a highly landslide risk-prone area for several reasons, including its high-slope mountains, geological and geomorphological formation, and the presence of the Indian tectonic plate underneath the Burmese plate. Slope is a key contributing factor that can destabilize the land, and the average slope of the area is about 15 degrees to more than 20 degrees, which is higher than the anticlinal area (Hossain et al. 2023). Forest cover in the mountainous regions acts as a natural barrier against the landslide, and our study area is primarily characterized by its dense forests, with more than 25 percent of that region consisting of protected forest areas, natural reserves, and national parks (Hassan et al. 2018; Hossen, Hossain, and Uddin 2019). However, the major concern arose after the influx of refugees when almost 1284.84 hectares of vegetation cover were cleared in the Teknaf Upazila and 102.87 hectares were cleared in the Teknaf Wildlife Sanctuary alone (Imtiaz 2018). The existence of these natural features and anthropogenic changes after 2017 made this area a hotspot for landslide risk assessment to mitigate future loss.

Data set and sources

In this study, 12 different types of landslide conditioning factors were considered under hazard, vulnerability, and exposure to assess landslide risk. Data from a wide range of national and international organizations and multiple field investigations were used to map all the chosen criteria (Table 1).

Table 1. Data types and sources used to assess landslide risk.

Data type	Source	Period	Mapping output
Sentinel-2 (10 m resolution)	United States Geological Survey (USGS) Earth explorer	December 2018	LULC, NDVI
Digital Elevation Model (30 m resolution)	NASA Earthdata	November 2013	Elevation, slope, aspect, stream network, stream density
Geology	Geological Survey of Bangladesh	Not applicable	Geological formation of surface
Precipitation	Bangladesh Meteorological Department (BMD)	2018–2020	Precipitation intensity
Historical landslides	Field survey, different articles and newspapers	2009–2020	Landslide inventory
Fault line	Geological Survey of Bangladesh	Not applicable	Distance to fault line
Population	Humanitarian Data Exchange (Projected)	2020	Population density

Production of thematic data layers

A number of thematic data layers for precipitation intensity (Boroumandi, Khamehchiyan, and Nikudel 2015; Jaedicke et al. 2014; Torkashvand, Irani, and Sorur 2014), landslide inventory (El Jazouli, Barakat, and Khellouk 2019; Jaedicke et al. 2014), distance to fault line (Boroumandi, Khamehchiyan, and Nikudel 2015; El Jazouli, Barakat, and Khellouk 2019; Pourghasemi, Pradhan, and Gokceoglu 2012; Torkashvand, Irani, and Sorur 2014), stream density (Boroumandi, Khamehchiyan, and Nikudel 2015), distance to stream network (Pourghasemi, Pradhan, and Gokceoglu 2012), elevation (El Jazouli, Barakat, and Khellouk 2019; Pourghasemi, Pradhan, and Gokceoglu 2012; Torkashvand, Irani, and Sorur 2014), aspect (Boroumandi, Khamehchiyan, and Nikudel 2015; El Jazouli, Barakat, and Khellouk 2019; Pourghasemi, Pradhan, and Gokceoglu 2012), slope (Boroumandi, Khamehchiyan, and Nikudel 2015; El Jazouli, Barakat, and Khellouk 2019; Pourghasemi, Pradhan, and Gokceoglu 2012; Torkashvand, Irani, and Sorur 2014), geology (Jaedicke et al. 2014), NDVI (Ahmed 2015; Tehrani and Husken 2019), LULC (Boroumandi, Khamehchiyan, and Nikudel 2015; El Jazouli, Barakat, and Khellouk 2019; Jaedicke et al. 2014; Pourghasemi, Pradhan, and Gokceoglu 2012), and population density (Jaedicke et al. 2014) were produced. All these factors are responsible for the occurrence of landslides in the area. In recent times, many anthropogenic activities have also made the area susceptible to landslides. These criteria were selected for the risk assessment upon analyzing all the historical data on landslides in the area.

Precipitation intensity

Depending on the ground condition, including terrain, soil type, and land cover characteristics, precipitation can work as a triggering factor for the landslide (Naranjo 2007). Analyzing previous rainfall data, it is evident that Cox’s Bazar District experienced the highest amount of rainfall in the southern part of Bangladesh. In this study, a precipitation intensity map is produced from the historical averaged monthly rainfall dataset for the station of Cox’s Bazar and Teknaf in the year 2018–2020, which is collected from the Bangladesh Meteorological Department (BMD) (Figure 3a). The created raster layer demonstrates the area receiving low to high rainfall during this period. A high rate of precipitation in an area can destabilize the slope which can result in the landslide hazard.

Figure 3. Hazard mapping criteria layers (a) precipitation intensity, (b) landslide inventory, (c) distance from fault line, (d) stream density, and (e) distance from stream

Landslide inventory

Landslide inventory plays an important role in providing the basic details for any landslide zoning, which includes vulnerability, risk, and hazard zoning (Fell et al. 2008). Usually, a landslide inventory database comprises detailed information about the landslide and is collected through extensive field work and other sources (Gopinath et al. 2023). A total number of 58 landslide points were collected from historical landslide data sourced from newspapers, articles, and various organizations and later verified through extensive field surveys that were conducted between 2019 and 2020. Landslides at these 58 points were mainly caused by rainfall and hill cutting. Casualties had also been recorded for different landslides. The landslide inventory map shows the historical landslide points of the area (Figure 3b).

Distance from fault line

Landslides are common in geological fault areas because tectonic breaks deteriorate the resilience of the surrounding rock (Chen, Panahi, and Pourghasemi 2017). Fault lines that exist in the study area are mainly east-dipping thrust faults on the western limb, while some of them also have west-dipping thrust faults on the eastern limb. Fault line data were collected from the Geological Survey of Bangladesh. Euclidean distance method is applied by using the ArcMap software to measure the distance from the fault line (Figure 3c). The generated map shows the distance of the areas from the fault line. Areas located near the tectonic structures are more susceptible to landslides. These phenomena are aided by selective degradation and water flow along fault planes (Pourghasemi, Pradhan, and Gokceoglu 2012).

Stream density

Slope stability largely depends on a slope’s closeness to the drainage structures where streams may significantly make an impact on the stability by the erosion of the slopes or the saturation of the bottom part of the material, resulting in a rise in water levels (Gökceoglu and Aksoy 1996; Saha, Gupta, and Arora 2002). A drainage density map was produced using digital elevation model (DEM) data (Figure 3d). To create a stream density map, we extracted all the streams from the DEM using ArcMap software. First of all, the flow direction function from ArcGIS’s hydrology tool was used to ascertain the direction of water based on the slope derived from the DEM. Then, a flow accumulation operation was performed to determine the potential accumulation of water at each point to get the stream layer. After extracting all the streams, stream density was calculated by using the kernel density tool.

Distance from stream

The saturation degree of the materials on the slope is a vital variable that determines the slope’s stability (Yalcin and Bulut 2007). In this study, the steam network was developed from the DEM, and the Euclidean distance method was applied to measure the distance of places from the streams (Figure 3e). Areas with a lower distance from the stream tend to contribute more to the landslide events.

Elevation

The impact of elevation on landslides is largely determined by its effect on various biophysical factors and human activities (Adnan et al. 2020). Because of the topographic nature of the research region, which consists of various highly elevated areas that influence the likelihood of landslides, elevation was chosen as a significant determinant. Such topography plays an important role in landslide activity, where higher elevated areas are considered the most vulnerable zone than the lower elevated areas because of hill cutting and clearing of the forest (Ahmed 2021; Mind’je et al. 2020). In this research, an elevation map was produced from a 30 m × 30 m DEM (Figure 4a). The maximum elevation is found in the Teknaf peninsula with a height of more than 150 m above sea level.

Figure 4. Vulnerability and exposure mapping criteria layers: (a) elevation, (b) aspect, (c) slope, (d) geology, (e) NDVI, (f) LULC, and (g) population density.

Aspect

Aspect is known as a primary landslide-influencing component that is considered in many research related to landslides (Fernández et al. 1999; Saha et al. 2005). It plays a crucial role in determining how meteorological events, such as the intensity of rain and the amount of daylight, impact the susceptibility of landslides by influencing the morphological structure of the region. For example, slopes with certain aspects may receive more intense rainfall or sunlight, leading to faster saturation and increased landslide risk. Additionally, the slope’s filtering capability, which is influenced by various variables, also contributes to this susceptibility (Pourghasemi, Pradhan, and Gokceoglu 2012). In this research, using the aspect toolbox from ArcMap on the DEM, aspects are categorized into nine classes, i.e. flat (−1), north (0–22.5, 337.5–360), northeast (22.5–67.5), east (67.5–112.5), southeast (112.5–157.5), south (157.5–202.5), southwest (202.5–247.5), west (247.5–292.5), and northwest (292.5–337.5), to present the relationship between aspect and landslides (Figure 4b).

Slope

Slope degree is the key parameter for analyzing its stability (Lee and Min 2001). As there is a direct relation between slope degree and landslide, it is often used to prepare landslide susceptibility maps (Clerici et al. 2002; Lee 2005). Because of the development of settlements and the cutting of the hills, a lot of steepness is created in the area, which makes it highly vulnerable to landslide risk. As a result, a slope map is prepared from the DEM by using ArcGIS 10.5 software’s slope tool. Based on the slope map, the study area is divided into five categories to demonstrate how steep the places are in the study area (Figure 4c).

Geology

Since rock permeability varies based on the rock type, geology is a significant landslide conditioning element (Vojteková and Vojtek 2020). In this study, the geological data type of the study area was collected from the Geological Survey of Bangladesh. The surface geology map of the study area includes a total of eight classes, where csd, H₂O, QTdt, QTg, Tb, Tbb, Tt represent beach and dune sand, water body, Dupi Tila Formation, Girujan Clay, Bhuban Formation, Boka Bil Formation, and Tipam Sandstone, respectively (Figure 4d). The predominant geological composition of the study region consists of Tipam sandstone (Tt), Girujan clay (QTg), and Dupi Tila (QTdt) formations, which are predominantly composed of fluvial sandstone, clay, and beach soil. The data revealed that the majority of slope collapses took place in the residual soil and weathered silty sandstone layers (Kamal et al. 2022a; Kamal et al. 2022b). These geological layers have been rasterized and reclassified into classes based on their contributions to the slope instability.

Normalized difference vegetation index (NDVI)

Vegetation cover is one of the major factors for assessing the landslide risk. If the amount of vegetation cover decreases in the hilly area, then the area becomes more susceptible to landslides. To assess the amount of vegetation cover in the study area, we performed NDVI analysis by using 10 m × 10 m sentinel-2 satellite imagery (Figure 4e). This method helps to assess the amount of vegetation from multispectral satellite images (Gandhi et al. 2015). The magnitude of the NDVI index ranges from − 1 to 1. Values ranging from 0.2 to 0.4 represent sparse vegetation, 0.4 to 0.6 indicate moderately dense vegetation, and dense vegetation lies between 0.6 and 1 (Alex, Ramesh, and Hari 2017; Quader et al. 2021). An NDVI map is created by using the following formula (Eq. 2(2) $\mathrm{NDVI}={\displaystyle \frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}}$(2) ) to analyze the sentinel-2 imagery: (2) $\mathrm{NDVI}={\displaystyle \frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}}$(2) For sentinel-2 imagery, band 8 and band 4 represent the near-infrared (NIR) and red bands, respectively.

Land use and land cover (LULC)

Slope stability is mostly dependent on land use and vegetation cover (Ocakoglu, Gokceoglu, and Ercanoglu 2002). In this study, a hybrid classification method was used to classify seven different types of land cover (Table 2) in the study area using 10 m × 10 m sentinel-2 satellite imagery with the help of ERDAS Imagine 2015 (Figure 4f). Hybrid classification is the combination of both unsupervised and pixel-based supervised classification techniques (Tang and Pannell 2009). For this image, first of all, ISODATA (iterative self-organizing data analysis) classification was applied to find the clusters of spectrally related pixels. A number of 35 classes were chosen to segregate the image. As there were some errors in the unsupervised classification to identify the features as per location, the recode method was applied based on the pixel-based supervised classification. An area of interest (AOI) was delineated in the wrongly classified areas to recode them into the desired one. The final output demonstrates that the area is mostly covered by the forest.

Table 2. Number of classes for hybrid classification of the study area.

Class name	Description
Forest	Deep vegetation
Homestead vegetation	Residential, commercial, industrial, roadside, and mixed
Crop land	Agricultural field
Settlement	Residences and industries
Refugee camps	Kutupalong RC, Kutupalong Expantion, Hakimpara RC, Jamtoli RC, Unchiprang RC, Leda RC, and Nayapara RC
Barren land	Abandoned land
Water body	River, tributaries, ponds, lakes, reservoirs, and wetlands

Any image classification must include an accuracy assessment to determine the degree of precision of a classified image (Islam, Uddin, and Hossain 2021). Confusion matrix was used to calculate the overall accuracy of the image (Brandt et al. 2013). The overall accuracy is determined by the total number of correctly classified points divided by the total number of reference points (Islam, Uddin, and Hossain 2021). 100 points were used for the accuracy assessment. Those training samples were generated using stratified random sampling techniques. The overall accuracy of the classified image is 86%, and the kappa coefficient is 0.83. The kappa value above 0.80 represents strong agreement between the reference and the classified data (Negassa, Mallie, and Gemeda 2020).

Population density

The number of people living per sq. km is defined as the population density of that area which needs to be taken into consideration while determining the vulnerability. A higher amount of vulnerability lies in the area that contains a large population, as there is a high chance of greater damage to life, property, and services (Arrogante-Funes et al. 2021). In this study, a projected population dataset for the lowest administrative area for 2020 was used to produce the population density layer (Figure 4g). Then, population per sq. km was calculated by dividing the total population by the total area of that particular administrative unit. The layer was subsequently rasterized and reclassified into five groups by using ArcMap.

Standardization of criteria layers and alternative ranking

Each criterion was rated into five points, where 1 represents very low risk, and 5 indicates very high risk. Low, moderate, and high risk were identified by classes 2, 3, and 4, respectively. A raster layer was created for each criterion in a 30 m pixel to apply the weighted overlay technique, where each criterion was converted into a raster from the vector layer.

The following equation (Eq. 3) was used to create a standardized layer for each criterion, where linear scale transformation brought the entire alternative values into a common range from 0 to 1 (Hoque et al. 2019). (3) $p={\displaystyle \frac{x-min}{max-min}}$(3) where p represents standardized score, min and max refer to the minimum and maximum values of each dataset, respectively, and x = a value of the single cell in the dataset.

Analytical hierarchy process (AHP)

Several approaches for landslide mapping have been used in different studies, such as on-ground monitoring, geospatial approach, geomorphologic approach, factors overlay, statistic models, as well as geotechnical process model (Bachri and Shresta 2010; Zinck et al. 2001). This study used the analytical hierarchy process to weight the criteria for developing the risk map. Saaty’s AHP approach is based on multi-criteria assessment (MCA) techniques, which were developed in the mid-1970s (Ahmed 2015). The AHP method was used to rank a set of alternatives. First of all, the AHP hierarchy was created. Then, a pairwise comparison was conducted for each landslide-triggering factor with every other factor that we considered in this study. These comparisons were set out by specialized groups of people including faculties researching landslides and stakeholders working on landslide management. In this process, comparison matrices were formed for the influential elements, and the elements were evaluated pairwise following Saaty’s scale. On that scale, a score of 1 suggests that two criteria are equally significant, while a score of 9 indicates a significant disparity in importance. Finally, the AHP model considers the pairwise comparison as an input and generates the relative weights as output (Bachri and Shresta 2010; Sangchini et al. 2016). These scores, reflective of expert judgments, underpin the weight matrix developed for this study, as shown in Table 3. The application of these scores is consistent with the methodology reported by Swetha and Gopinath (2020), ensuring that our weighting approach aligns with established AHP practices.

Table 3. Weights and consistency ratios of criteria calculated from the pairwise comparison matrices.

Component	Criteria	Weight
Hazard	Precipitation intensity	0.04	Consistency ratio: 0.03
	Landslide inventory	0.08
	Distance to fault line	0.39
	Stream density	0.24
	Distance to stream network	0.25
Vulnerability and exposure	Elevation	0.09	Consistency ratio: 0.04
	Aspect	0.23
	Slope	0.34
	Geology	0.15
	NDVI	0.04
	LULC	0.07
	Population density	0.10

The consistency ratio (CR) is a consistency index used in AHP to calculate the degree of consistency (Reis et al. 2012). It is calculated using the following formula (Eq. 4): (4) $\mathrm{CR}=\mathrm{CI}/\mathrm{RI}$(4) where the random index (RI) indicates the randomly produced average consistency index, and the consistency index (CI) can be calculated as follows (Eq. 5): (5) $\mathrm{CI}=(\lambda max\text{–}n)/(n\text{–}1)$(5) where λmax means the greatest eigenvalue of the matrix and n indicates the matrix order (Malczewski 2010; Reis et al. 2012). The weights and consistency ratios of the criteria are presented in Table 3.

In determining the weights for the 12 criteria in the AHP model, each criterion was thoroughly assessed to determine its relative significance in contributing to landslide risk. Direct influence on slope stability was assessed by considering factors such as precipitation intensity, proximity to fault lines, and hydrological parameters. Additionally, the area’s characteristics like population density were assessed for their relation to the possible impact of landslides on human settlement and the surrounding environment.

Results and discussion

Hazard mapping

Landslide hazard refers to the probability of encountering a landslide that has the potential to cause harm within a certain timeframe, at a particular location, and of a specified scale (Kamal, Ahmed et al. 2022a; Kamal, Hossain et al. 2022b). The produced hazard map (Figure 5a) shows the spatial variation and degree of the landslide hazard in the study area. This spatial scenario was developed by overlaying the layers of distance to the stream network, stream density, precipitation, landslide frequency, and distance to the fault line. When there are no seismic events, heavy rainfall is often the main cause of landslides. However, certain factors, such as the saturation of hillslopes caused by previous rainfall, may also play a role (Emberson, Kirschbaum, and Stanley 2021). In the hazard map, 149.30 sq. km. (26.07%) of the area is considered a very high landslide hazard zone. On the other hand, 30.53 sq. km. (5.33%) of the area is considered a very low landslide hazard zone, where 89.64 sq. km. (15.65%), 138.27 sq. km. (24.14%), and 165.02 sq. km. (28.81%) of the area are identified as low, moderate, and high hazard zone respectively. Most of the very high hazard zone is in Teknaf Union, where a large number of hilly areas exist. Some parts of this zone are found near Inani and Shamlapur areas; and in the Teknaf wildlife sanctuary. High and moderate hazard zones cover most of the area along with Rajapalong, Court Bazar, and Bhalukia of Ukhiya in the north; Balukhali, Whaikhyang, Noapara, and Bardell refugee camps at the center; and Teknaf Union in the south. The low hazard zone covers the low land on the border of the study area. The main reasons for the high hazard of landslides are high stream density and access to the fault line.

Figure 5. Risk components maps along with the risk map of the study area; (a) hazard, (b) vulnerability and exposure, and (c) risk.

Vulnerability and exposure mapping

While vulnerability encompasses the human dimensions of disasters, such as social, economic, and cultural factors, exposure refers to the elements at risk of damage by potential hazards (Ahmed 2021). The spatial pattern of the vulnerability and exposure is derived from the GIS layers of slope, elevation, aspect, NDVI, geology, LULC, and population density. All of these factors have a substantial impact on the landslide vulnerability in the area. In this study, the produced vulnerability and exposure map (Figure 5b) illustrates that 41.06 sq. km. (7.22%) of the study area is very highly vulnerable to landslides, and 98.12 sq. km. (17.25%) of the area is of very low vulnerability. Teknaf Union, and part of Whaikhyang and Palongkhali Unions are identified as belonging to the very highly vulnerable zone to landslides. This map gives a comprehensive overview of the vulnerable and exposed communities living under the landslide risk zones. Some of the studies found that slopes with more than 4.6° and 15.7° have a frequency ratio value over 1 and 1.5, respectively (Rabby and Li 2020). Accordingly, these areas are highly vulnerable because of the high elevation, greater degree of slope, geology of the area, and reduction of vegetation cover. Besides, 163.14 sq. km. (28.68%) of the area is of low vulnerability, 160.92 sq. km. (28.29%) of the area is moderately vulnerable, and 105.59 sq. km. (18.56%) of the area is highly vulnerable to landslides. Also, the increasing population living along the hills makes them exposed to landslides. Uncontrolled and extensive hill cutting along with unplanned settlements are considered as some of the main reasons for landslides (Rahman, Ahmed, and Di 2017). The vulnerability and exposure map also shows the pattern in Teknaf area with large population density and settlement in the high elevated and steep slope of mountainous areas.

Risk mapping

Landslide risk refers to the occurring probability of a landslide and the resulting impacts and damages (Kamal, Ahmed et al. 2022a; Kamal, Hossain et al. 2022b). We simulated the landslide risk in the Ukhiya and Teknaf Upazilas using a weighted overlay analysis of all the parameters reflecting hazard, vulnerability and exposure. The produced risk map (Figure 5c) shows that 12.47 sq. km. (2.19%) of the study area falls under a very high risk zone and 72.46 sq. km. (12.74%) falls under the high risk zone for the landslide. A maximum of 131.26 sq. km. (23.08%) area has been identified as a moderate risk zone. 139.53 sq. km. (24.53%) area has been identified as a very low risk zone, and 213.01 sq. km. (37.45%) has fallen under a low risk zone for the landslide. Risk mapping in this area addresses the complexities and uncertainties associated with the quantitative assessment of landslide risk (Huang et al. 2022). Sites that are identified as highly risky are due to several criteria. As a result of a large number of Rohingya influx in Ukhiya and Teknaf Upazilas (shown in Figure 2), huge amounts of forests are degraded, and some parts of the study area have become risky for the landslide. High risk sites are found in central Teknaf, Teknaf wildlife sanctuary, and part of western Palongkhali. In the context of refugee camps, our study found Leda and Nayapara camps, where a large number of refugees are living, are in the highest risk zone.

The representation of the risk scenario in this analysis was developed based on the spatial overlap of the causative factors, i.e. the higher the spatial intersection of the positively associated contributing factors, the higher the landslide risk. The reasons why this area is mostly identified as a risk zone are forest clearing at a higher rate, a high level of precipitation, and the addition of a large number of refugees located in the area that reduce the resistance of the hilly area. During the monsoon season, this region receives a high amount of precipitation and becomes prone to landslides. IOM (2022) reported that more than 940 and 910 slope failure events were recorded in Kutupalong Registered Camp (KRC) alone between June 2019 and July 2021, respectively. At that time, an early warning system developed by Ahmed et al. (2020a) could have been applied to the high risk sites identified in this study to reduce the casualties. The current approach to landslide risk assessment can address the shortcomings of the existing techniques of evaluating landslide hazards based only on the physical factors without considering the incorporation of social vulnerability elements (Zhou et al. 2021). Similarly, we conducted our assessment considering the physical aspects along with the population density.

Correlation of the assessed risk with experimental impact

A total of 58 landslide points were collected to identify the deadliest landslide zones in the field. These incidents are considered for comparison with the assessed risk output of the study. According to the landslide inventory dataset, almost all the deadliest landslides took place from the central part to the southern part of the Teknaf Union. The assessed risk zone of this study strongly correlates with the inventory dataset and found the assessed area mostly falls either in the high risk zone or very high risk zone. Some of the landslide points from the inventory dataset do not align with the assessed risk area, as those points fall under the moderate risk zone and a few fall under the low risk zone. For example, this study designates the central-north part of Teknaf Union as a moderate risk zone, whereas there is a pretty high chance of a landslide based on the landslide inventory dataset. Almost 67.27% of the historical landslide episodes are located in the high or very high risk zones of the study, whereas 22.41% of the landslide episodes located in the moderate risk zones. Only 10.34% landslide events fall under the low and very low risk zones of the assessed result.

Conclusion

In this study, we prepared a comprehensive landslide risk map of Ukhiya and Teknaf Upazilas, a highly populated Rohingya Refugee Settlement area, located in the southeast part of Bangladesh. We used geospatial techniques to prepare and process data of various landslide triggering criteria, while multi-criteria weighting technique AHP was used to weight the criteria. The landslide map was validated using the landslide inventory data. The risk map illustrates high-risk areas located in the center of Teknaf, Teknaf wildlife sanctuary, and part of western Palongkhali, where some parts are close to the Rohingya refugee camps. Some of these sites had already experienced low to moderate level of landslides. The concerned authorities can be proactive in mitigating future losses in these areas by using the results of this study. This study will help them assess and identify the severe and moderate landslide-prone areas properly so they can take the necessary steps before any landslide disaster occurs.

Although we tried to utilize the available data sources and expert opinions properly, several limitations including data quality may affect the findings of the study. The collection of quality and up-to-date data required for the consideration of 12 criteria was challenging. Compared to using a 10-m resolution DEM, high-resolution topographic data such as LiDAR might add more value to the assessment. For LULC classification, high-resolution satellite images could have been better than sentinel-2. Historical landslide inventory data is also not that much available and authentic for the region. Furthermore, an assessment of mitigation capacity could make the result more fruitful, but again this was prevented by a lack of relevant data.

Despite all these limitations mentioned above, the results of the study are still very useful for developing strategies to mitigate the loss due to landslides in the area. Many forest covers that were degraded need to be planted again in the high risk zone so that they can make soils resist landslides. As many Rohingya are settled in these areas, they degrade the forest area at a faster rate to meet the needs for makeshift shelter and firewood. So, these people should be managed properly and resettled from the forest area. Then, the hill cutting should be strictly monitored so that the slopes do not become vulnerable to landslides. If all these anthropogenic causes can be maintained properly, it is expected that the loss of life and property due to landslides may be reduced. Overall, the study suggests that some of the crucial anthropogenic factors that trigger the landslide in the study area are deforestation and slope alteration that make the land unstable. Also, natural factors that are contributing through high precipitation, large number of stream flows, and high elevation. As a result, proper land use planning is required in the study area through some fast-growing species plantation, relocation of refugees and local inhabitants living from the high risk sites, and ensuring slop stabilization through grass planting or concrete slope retention walls.

Acknowledgments

We are thankful to the Bangladesh Meteorological Department, the Geological Survey of Bangladesh, Cox’s Bazar local administration, NASA Earthdata, and the United States Geological Survey for providing the necessary data that we required. We are also grateful to the anonymous reviewers for their valuable comments to improve the quality of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Raw data were generated at Jagannath University, Bangladesh. Derived data supporting the findings of this study are available from the corresponding author (Muhammad Al-Amin Hoque) on request.

Additional information

Funding

This research is funded by the Ministry of Science and Technology, Bangladesh under the National Science and Technology Fellowship, Fiscal Year: 2020–2021.

Notes

1 Upazila is a sub-unit of a district.