Evaluating the Effects of Digital Elevation Models in Landslide Susceptibility Mapping in Rangamati District, Bangladesh

Digital elevation models (DEMs) are the most obvious data sources in landslide susceptibility assessment. Many landslide casual factors are often generated from DEMs. Most studies on landslide susceptibility assessments rely on freely available DEMs. However, very little is known about the performance of different DEMs with varying spatial resolutions on the accurate assessment of landslide susceptibility. This study compared the performance of four different DEMs including 30 m Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM), 30–90 m Shuttle Radar Topographic Mission (SRTM), 12.5 m Advanced Land Observation Satellite (ALOS) Phased Array Type L band Synthetic Aperture Radar (PALSAR), and 25 m Survey of Bangladesh (SOB) DEM in landslide susceptibility assessment in the Rangamati district in Bangladesh. This study used three different landslide susceptibility assessment techniques: modified frequency ratio (bivariate model), logistic regression (multivariate model), and random forest (machine-learning model). This study explored two scenarios of landslide susceptibility assessment: using only DEM-derived causal factors and using both DEM-derived factors as well as other common factors. The success and prediction rate curves indicate that the SRTM DEM provides the highest accuracies for the bivariate model in both scenarios. Results also reveal that the ALOS PALSAR DEM shows the best performance in landslide susceptibility mapping using the logistics regression and the random forest models. A relatively finer resolution DEM, the SOB DEM, shows the lowest accuracies compared to other DEMs for all models and scenarios. It can also be noted that the performance of all DEMs except the SOB DEM is close (72%–84%) considering the success and prediction accuracies. Therefore, anyone of the three global DEMs: ASTER, SRTM, and ALOS PALSAR can be used for landslide susceptibility mapping in the study area.


Introduction
Local terrain conditions, including terrain relief, hydrology, geology, and land use, are crucial for the assessment of landslide susceptibility [1]. These local features are often termed as "causal factors" [2]. Identifying these causal factors is considered as the steppingstone of landslide susceptibility assessment. Landslide susceptibility represents the likelihood of landslide occurrence in an area. It assumes that future landslides may occur in previous landslide locations where the causal factors already created conducive environments for triggering landslides [2][3][4]. Several natural and anthropogenic

Landslide Inventory
Landslide inventory records different information, including the exact location, size, type, time of occurrence, causalities, trigger, and causes [2]. Rabby and Li [33] prepared and published [34] landslide inventory of the Chittagong hilly area, Bangladesh. They used participatory field mapping to prepare this inventory [41]. In our study, we selected 168 landslides from their inventory as it has been advised to use more than one method in landslide inventory preparation [42], we analyzed available Google Earth images on the Google Earth platform to map more landslides in the study area. We used the method proposed by Rabby and Li [33] for Google Earth mapping and mapped 93 landslides that occurred from January 2001 to January 2019. Therefore, in our study, we used a total of 261 (168 + 93) landslide locations. The mean size of the landslide was 274.2 m 2 . The smallest and the largest dimensions of the landslides were about 14.6 m 2 and 3422.4 m 2 , respectively.

Landslide Causal Factors
In our study, we used a total of 15 causal factors (Table 1) to produce landslide susceptibility maps (see Appendix A-C). Of these factors, seven factors were derived from DEM: elevation, slope,

Landslide Inventory
Landslide inventory records different information, including the exact location, size, type, time of occurrence, causalities, trigger, and causes [2]. Rabby and Li [33] prepared and published [34] landslide inventory of the Chittagong hilly area, Bangladesh. They used participatory field mapping to prepare this inventory [41]. In our study, we selected 168 landslides from their inventory as it has been advised to use more than one method in landslide inventory preparation [42], we analyzed available Google Earth images on the Google Earth platform to map more landslides in the study area. We used the method proposed by Rabby and Li [33] for Google Earth mapping and mapped 93 landslides that occurred from January 2001 to January 2019. Therefore, in our study, we used a total of 261 (168 + 93) landslide locations. The mean size of the landslide was 274.2 m 2 . The smallest and the largest dimensions of the landslides were about 14.6 m 2 and 3422.4 m 2 , respectively.

Landslide Causal Factors
In our study, we used a total of 15 causal factors (Table 1) to produce landslide susceptibility maps (see Appendices A-C). Of these factors, seven factors were derived from DEM: elevation, Remote Sens. 2020, 12, 2718 5 of 35 slope, plan curvature, profile curvature, topographic wetness index (TWI), stream power index (SPI), and aspect. As we are comparing the performance of different DEMs, we derived each of these factors from four different DEMs: ASTER, SRTM, ALOS PALSAR, and SOB. The rest of the eight factors namely land use/land cover, land use/land cover change, geology, distance from the road networks, distance from the fault lines, distance from the drainage networks, rainfall, and normalized difference vegetation index (NDVI) were collected from different datasets. Maps of different causal factors had different resolutions, but for the convenience of comparison, we kept the 30 m resolution as the standard for landslide susceptibility maps. In the following sub-section, we provide a brief overview of the causal factors that we used in this study. We classify these factors into several classes primarily using Jenks Natural Break method in ArcGIS 10.7, unless otherwise mentioned.

Elevation
A change in elevation can bring changes in geomorphology, vegetation, and rate of erosion in an area and thus alters the landslide susceptibility [8]. We derived elevation from ASTER, SRTM, ALOS PALSAR, and SOB DEMs ( Figure A1a-d of Appendix A) and divided them into five classes (see Table A1 of Appendix D).

Slope
The slope is one of the most critical factors of landslides. Generally, with the increase of slope, shear stress increases, and therefore landslide susceptibility increases [43,44]. Like elevation, we derived slopes from four different DEMs using the slope tool in ArcGIS 10.7 ( Figure A1e-h of Appendix A) and divided them into five classes using the Jenks natural break method (see Table A1 of Appendix D). We found different maximum slope values for different DEMs-ASTER (51.89º), SRTM (61.24º), ALOS PALSAR (65.36º), and SOB (46.4º). As these values were different for the same study area, the five classes of slopes (see Table A1 of Appendix D) were different.

Aspect
The direction of precipitation, sunlight, and wind depend on aspect, and therefore it has effects on the growth of vegetation, rate of erosion, and thickness of soil [45]. From the DEMs, four aspect maps ( Figure A1i-l of Appendix A) were prepared and were divided into ten classes (see Table A1 of Appendix D).

Plan Curvature and Profile Curvature
Curvature is the rate of change of slope over time in an area. We used the four different DEMs to produce four plan and profile curvature maps (Plan: Figure A2a-d; Profile: A2e-h of Appendix B). Profile and plan curvatures were divided into three classes: concave, convex, and flat. Among these classes, concave slopes are more prone to landslides because water cannot disperse equally on these slopes [5].

Topographic Wetness Index and Stream Power Index
Topographic wetness index (TWI) increases with the decrease of the slope; therefore, it is inversely related to landslide susceptibility [4,43]. Stream power index (SPI) represents the erosion power of streams. SPI is directly related to slopes; in a steeper slope, SPI will be higher, representing more erosion power while in a flat alluvial plain, SPI is low [45]. TWI and SPI maps were derived from four DEMs using Equation (1) and Equation (2) (TWI: Figure A2i where, A = Area of a specific catchment and α = Slope gradient of the specific area. We divided TWI and SPI into five classes (Table A1 of Appendix D)

Rainfall
The intensity and duration of rainfall controls the initiation of landslides [44]. We used the mean annual rainfall of five weather stations of Bangladesh Meteorological Department (BMD) to prepare the rainfall map using the Kriging interpolation method ( Figure A3a of Appendix C). We later divided it into five classes (Table A1 of Appendix D).

Distance-Based Causal Factors
Distance from the road networks, drainage networks, and fault lines were the three distance-based causal factors in this study. We used the Euclidean distance tool in ArcGIS 10.7 to derive the distance of landslides from the targeted features: road, drainage, and fault lines ( Figure A3b-d of Appendix C) and divided the distances into five classes (Table A1 of Appendix D). Distance from the road networks is one of the most critical factors. The undercutting of slopes during road construction and the vibration created by vehicles damage the slope stability [43]. Drainage network indicates the zone of erosion in an area, and erosion is indirectly linked with the landslide susceptibility [4]. Fault lines indicate the geomorphological discontinuity in an area. Near the fault lines, the shear strength of rock is minimum. Therefore, areas near to the fault lines are prone to landslides [44].

Normalized Difference Vegetation Index (NDVI)
NDVI indicates the growth of vegetation and biomass of an area [46]. Generally, the probability of the occurrence of the landslide on the naturally vegetated surface is lower than the bare lands [8]. We used Landsat 8 level 2 imagery of 11/10/2017 to prepare the NDVI map ( Figure A3e) and divide it into five classes (Table A1 of Appendix D).

Geology
The strength of rock and soil permeability depends on the geology of an area. Therefore, geology has an impact on landslide susceptibility [43]. In this study, we used the geological map ( Figure A3f) (1:100,000) of the Bangladesh Geological Survey (BGS). There are eight types of geologic formation found here: Dihing and Dupi tila formation; Bokabil formation; Bhuban formation; Tipam Sandstone; Valley alluvium and colluvium; Dihing formation, Girujan clay, and waterbodies.

Land Use land Cover
Land use/land cover and land use/land cover change are the two most crucial landslide causal factors in the study area [35,47]. Ahmed [31] and Rabby and Li [33] found that the rate of change of land use/land cover in our study area is high compared to other adjacent areas. In this study, we used two Landsat imageries to analyze the land use/land cover change (Landsat 5: date of Acquisition: 24/12/1998; Landsat 8: Acquisition Date: 29/11/2018). We used supervised maximum likelihood classifier to classify the 1998 and 2018 images into four land use classes: bare land, vegetation, built-up, and water bodies ( Figure A3g of Appendix C). We later employed post-classification change detection techniques to analyze the land use/land cover changes between 1998 and 2018 ( Figure A3h of Appendix C).

Methodology
To compare the effects of four different DEMs: ASTER, SRTM, ALOS PALSAR, and SOB on landslide susceptibility maps, we used a bivariate method: modified frequency ratio (MFR), a multivariate method: logistic regression (LR) and a machine learning method: random forest (RF). We assessed two scenarios: (a) considering only DEM-based seven causal factors in the models, and (b) considering all 15 causal factors, including the DEM-based factors, in the models. As we used three methods: MFR, LR, and RF on four different DEMs under two scenarios, therefore, the outcome would be twenty-four landslide susceptibility maps. For legibility, we used different acronyms in later sections (see Table 2).

Training and Validation Dataset
We divided the 261 landslide locations randomly into training (75%) and validation (25%) datasets. Bivariate MFR is a one-class classification method where non-landslide locations or absence of landslides are not required [48]. On the other hand, for multivariate LR and machine learning-based RF, the selection of non-landslide locations (pseudo-absence points) is essential [4]. Any place that does not have landslide can be considered as non-landslide. We randomly selected 261 non-landslide locations ( Figure 1) (pseudo absence points) from the study area [49]. We split these non-landslide locations into training (196) and validation (65) data sets. In total, we had 392 (196: landslide locations; 196: non-landslide locations) data points for training and 130 (65: landslide locations; 65: non-landslide locations) data points for testing the LR and RF models.

Modified Frequency Ratio (MFR)
MFR is an improved version of the widely used frequency ratio (FR) method [35,50]. Lee and Talib [51] proposed FR, which assesses the spatial relationship between the landslide locations and the landslide causal factors [52]. In the FR method, each of the causal factors must be divided into subclasses or categories; for example, in this study, we split slope into five categories using the Jenks natural break method (Table A1 of Appendix D). We calculated the FR values using Equation (3), and later these FR values of each of the subclasses of causal factors were used in the MFR model.
where FR ij = Frequency Ratio of jth Subclass of Factor i N ij = Total area of the landslide pixels within the jth subclass of factor i N = Total area of landslide pixels in the study area M i j = Total area of the pixels in the jth subclass of factor i M = Total area of the study area FR > 1 means association of landslides with that subclass. In other words, there is a probability of occurrences of landslides in that subclass. FR < 1 means no association [53].
For calculating the MFR, then we normalized the FRs using Equation (4).
where Rf ij = Relative frequency of jth subclass of factor i FR ij = Frequency ratio of the jth subclass of factor i FR i = Sum of the frequency ratios of factor i Later, we calculated the prediction rate (PR) using Equation (5). In the FR model, the overall contribution of causal factors to the occurrence of landslides is not measured. Only the subclass wise contribution is measured [6,50]. In MFR, we can measure the overall contribution because PR indicates overall association of a causal factor. The lowest value of PR is 1 and the higher the PR value the stronger is the association of causal factor with the landslides [6]. To calculate the landslide susceptibility ondex (LSI) and to produce the landslide susceptibility maps Equation (6) was used Rf ij × PR i (6) where LSI = Landslide susceptibility index Rf ij = Relative frequency of jth subclass of factor i PR i = Prediction rate of factor i For landslide susceptibility mapping, we used the reclassify tool in ArcGIS 10.7 to reclassify the categories of causal factors according to the Rf values. Later, we multiplied each of the reclassified raster layers with the prediction rates and summed up to produce the final landslide susceptibility maps.

Logistic Regression (LR)
Logistic regression (LR) is one of the most widely-used multivariate statistical methods in landslide susceptibility mapping [4,6,[54][55][56]]. An LR model predicts the presence of landslides using the binary landslide data (presence and absence of landslides or landslides and non-landslides) and their relationship with the landslide causal factors [56,57]. Here, landslide and non-landslide locations are dependent variables, and causal factors are independent variables. These independent variables can be numerical or categorical [4]. Equation (7) is used in LR model where Y = The presence of landslides X i = ith Causal factor ß i = Regression coefficient of the ith causal factor e = Error We used Equation (8) to determine the probability.
We used the R software environment for the forward stepwise LR method. We multiplied the raster layers of the statistically significant causal factors with the coefficients and summed up using Equation (7) in the R software environment. Finally, we used Equation (8) to produce the landslide susceptibility maps.

Random Forest Classification (RF)
The random forest method was developed by Breiman [58] and is an ensemble learning method [19]. Lately, the use of the RF method for landslide susceptibility mapping has increased due to its high performance in predicting landslide locations [37,59].
This method uses bootstrapping techniques to generate a bunch of classification trees based on subsets of observations [27]. There is high variance among the individual trees, and therefore classification based on a single tree is unstable and prone to overfitting [37]. Random forest is improved over commonly used tree-based methods, such as a decision tree or bagged tree because it decorrelates the trees. RF uses ensembles of trees and lets each tree define the class membership, and finally, the respective class is assigned based on the highest votes [27,37]. Since the bootstrapping method is used, a set of data is not used in the model training stage and this set of data is known as out-of-bag (OOB) [27]. These OOB data are used to calculate the mean decrease of accuracy and Gini coefficient [37]. The accuracy and Gini coefficient are used in variable selection and ranking [19]. They also provide the statistical weights or variable importance of each of the predictors used in the model [27].
There are several advantages of using RF methods, such as rescaling and transformation of data are not essential; missing data and outliers can be ignored [60]. Moreover, it can deal with both numerical and categorical data, and the use of a dummy variable is not required [37]. In this study, we developed the Random Forest model in the R software environment using the "randomForest" package [61].

Multicollinearity Diagnostics
In the LR model, multicollinearity can bring inaccuracies in variance and unsuitability in estimates. On the other hand, in the RF model, it can affect the variable importance [62]. Therefore, multicollinearity diagnostics: variance inflation factor (VIF) and tolerance were used before using the causal factors in LR and RF models. Since VIF values were <10 and tolerance were <0.3, causal factors were independent and were used in these two models [63]

Model Validation and Comparison
Success and prediction rate curves were used to test the performance of the susceptibility models. The training data set was used for the success rate curve, and the validation data set was used for the prediction rate curve [6]. "The area under the curve" (AUC) of success rate indicates how well the model fits the training data while the AUC of prediction rate suggests how well the model will predict future landslides [15,64]. AUC value ranges from 0-1 or 0%-100% and it can be grouped into the following categories: 0.50-0.60 (fail); 0.60-0.70 (poor); 0.70-0.80 (fair); 0.80-0.90 (good), and 0.90-1.00 (excellent) [53]. We also used two non-parametric tests: Friedman and Wilcoxon Signed Rank test to assess whether there are any significant differences in performances between the susceptibility [65][66][67][68]. Friedman's test is used for multiple comparisons. This test determines whether there is any significant difference in performance in multiple models [29], while the Wilcoxon Signed Rank test is used for pairwise comparison of susceptibility models and, therefore, can indicate which models are significantly different [68].
However, statistical performance assessments such as success and prediction rate curves cannot show the level of agreement among the models. Therefore, we used convergent validation through the coverage based cross-comparison [69,70]. The MFR gives landslide susceptibility index while LR and RF provide the probability of landslides. We reclassified them into five susceptibility zones: very low, low moderate, high, and very high using the Jenks natural break method. Later, we used the raster calculator in the ArcGIS platform to subtract the reclassified model from one another. The outcome can be any integer, but only "zero" will indicate the areas which were classified into the same susceptibility zones by two compared models. We calculated the percentage of area under this "zero" class and this percentage indicates the spatial convergence or agreement between two compared models [69].

Susceptibility Assessment Using Modified Frequency Ratio (MFR)
We found that the prediction rates (PRs) of the slope, aspect, and elevation derived from three DEMs (ASTER, SRTM, and ALOS PALSAR) are similar. The PRs for these three factors are around 2.25, 1.0, and 3.0, respectively (see Table A1 of Appendix D). Compared to these three DEMs, the SOB DEM showed different PRs. Our analysis revealed that for slope, aspect, and elevation, the SOB DEM-based PRs are 1.79, 2.37, and 2.41, respectively. We further found that with the increase of slope, the probability of landslide increases. The relatively safer zones are in areas below 8º slope where frequency ratio (FR) < 1. For ASTER, SRTM, and ALOS PALSAR DEMs, we found that the slope class 14-23º had the highest probability of landslides. However, for SOB DEM, the highest probability of landslide is in the 8-14º slope class. In the case of TWI, we found that ALOS PALSAR DEM has the highest PR (4.30) followed by ASTER DEM (PR = 3.19) and SOB DEM (PR = 1.61). For all the four DEMs, the probability of landslides was higher in areas where TWI is less than 6. The SPIs derived from four DEMs have lower PRs compared to other causal factors and the class-wise weight (FR values) showed the same sort of pattern. We further found that for plan curvature, ALOS PALSAR has the highest PR (4.30), while for profile curvature, SRTM has the highest PR (4.68).
In general, we found no specific pattern of FR and PR values for these seven factors. We observed that causal factors derived from either ALOS PALSAR or SRTM have higher PR values, and SOB has the lowest PR among the four DEMs. The causal factors derived from different DEMs did not have a significant impact on the FR and PR values of MFR. This finding is similar to the findings of Chang et al. [29]. Most of the topographic factors are the first derivative of the DEMs other than the TWI and SPI. Thus, for TWI, DEMs have more substantial impacts than other factors. It is because there can be a small difference among the DEMs, and when the second derivative is used, these become pronounced [29].
The class-wise FR values of the eight factors that are not derived from the DEMs are generally similar; however, the PR values are different (see Table A1 of Appendix D). The causal factors derived from the DEMs played a crucial role in determining the PR values for these eight factors. For ASTER (0.141), SRTM (0.134), and ALOS PALSAR (0.139), the lowest difference was between the maximum and minimum relative frequency of aspect. While for SOB (0.160), it was SPI. Since these values are slightly different, it affected the PR values.

Landslide Susceptibility Maps (MFR)
We produced Landslide Susceptibility Indices (LSIs) of the four MFR models using Equation (6) for DEM-based causal factors. The LSI of MFR_ASTER_DEM ranged from 994.6 to 8388.1. The LSIs for MFR_SRTM_DEM LSIs for MFR_SRTM_DEM; MFR_ALOS_DEM and MFR_SOB_DEM ranged from 591.4 to 9458.8 and 527.8 to 10056.5 and 638.8 to 6130.0, respectively. LSI does not have a unit. It is the product of relative frequency and prediction rate Equation (6), and both of these do not have units. The greater the LSI, the greater is the landslide susceptibility, and the smaller the LSI value, the lower is the susceptibility [35,44]. The ranges of LSIs indicate that MFR_SRTM_DEM and MFR_ALOS_DEM models had a comparatively broader range than the rest of the two models. This happened because of the variable FR and PR values. For ASTER DEM, the highest PR value was for plan curvature, while for other DEMs it was for profile curvature (see Table A1 of Appendix D). For all DEM-based factors, SOB had the lowest PRs among the four models and therefore it affected the LSIs. Later, we used the same Equation (6) to produce four MFR models based on 15 causal factors. The LSI of MFR_ASTER DEM ranged from 1613.00 to 20370.10. The LSIs for MFR_SRTM, MFR_ALOS, and SOB DEMs ranged from 1314.40 to 22300.34 and 1234.95 to 22180.24 and 1995.7 to 17316.9, respectively. The LSIs of MFR_SRTM and MFR_ALOS DEMs have a comparatively broader range than the rest of the two models. The highest PR value for ASTER was 6.25 and for SOB, it was 5.64 (Table A1 of Appendix D) for distance to the road network, while for SRTM and ALOS PALSAR, it was 6.61 and 6.38, respectively. For other causal factors (Table A1 of Appendix D), SOB had the lowest PRs among the four, and therefore LSI was the lowest. As mentioned earlier, the FR values varied for seven topographic factors derived from four different DEMs, and these factors had impacts on the PR of eight common factors. Ultimately these variations defined the LSI of the susceptibility maps.
Since different models had different LSIs, Rescale by Function tool in ArcGIS was used to normalize the LSIs into a 0.0-1.0 scale. Later, we used the Jenks natural break method to classify the normalized LSIs into five susceptibility zones: very low, low, moderate, high, and very high. Generally, the spatial appearances of the landslide susceptibility maps have similarities with the map of causal factors that have a higher contribution to landslides. In this study, this contribution is shown by PR and FR values (see Appendices A-C). We found that the spatial appearance of seven causal factors derived from SOB was different from the spatial appearance of seven causal factors derived from ASTER, SRTM, and ALOS PALSAR. ASTER, SRTM and ALOS PALSAR based susceptibility maps (Figure 2a-c) show a comparatively lesser percentage of area as very low or low susceptibility zones than the SOB-based landslide susceptibility map (Figure 2d). However, the SRTM based map shows comparatively more areas as high and very high susceptibility zones that the other maps. We found when all factors were considered, areas near to the road network are highly susceptible to landslides (Figure 3a-d), because the PR values of the distance from the road networks were the highest.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 35 significant factors. Therefore, the susceptibility map took the shape of the map of these two factors (see Appendix A). When all (15) causal factors were used, the spatial appearance of LR_SOB (Figure 3h) was different from the other three maps (Figure 3e-g). LR_SOB map was influenced by the distance from the fault lines, distance from the road networks, and land use/land cover. Although aspect was a significant factor, the coefficient value of aspect was similar to other significant factors, and distance from the fault lines (ß = 1.07) (Table A2 of Appendix E) had a higher coefficient value than aspect. Therefore, most of the study area was classified as low or very low susceptibility zones. In the LR_ASTER model (Figure 3e), the slope had the highest coefficient (ß = 0.31) (Table A2 of Appendix E), and therefore, areas with steeper slopes were classified as high or very high susceptibility zones. But in LR_SRTM and LR_ALOS models, the slope had comparatively lower coefficient values than ASTER. As a result, some areas were classified as moderate susceptibility zones in these two maps (Figure 3f-g). In LR_SRTM and LR_ALOS susceptibility maps, common factors such as distance from the road networks and fault lines, land use/land cover, and land use/land cover change did not have higher coefficient values than the DEM-based causal factors. That is why, unlike LR_SOB, the spatial appearance of the susceptibility maps did not follow the appearance of the maps of common factors.

Susceptibility Assessment using Random Forest (RF)
RF_ASTER_DEM, RF_SRTM_DEM, and RF_ALOS_DEM models detected (Figure 4a) slope and RF_SOB_DEM model identified the aspect as the most critical factor. When all 15 factors were used in the models, RF_ASTER and RF_SRTM (Figure 4b) detected slope and the rest of the two models detected distance from the road network as the most important causal factor. For RF_ASTER, RF_SRTM, and RF_ALOS, DEM-based factors such as elevation, TWI, and aspect had higher

Susceptibility Assessment using Logistic Regression (LR)
For the DEM-based factors, the LR model detected two to five statistically significant factors (see Table A2 of Appendix E). Elevation and slope were the two common statistically significant causal factors for ASTER, SRTM, and ALOS PALSAR based models. Since the ALOS PALSAR based model, the highest number of causal factors was chosen, DEM had the highest impact on the landslide susceptibility map. Odds ratio (Table A2 of Appendix E) shows that slope was the most important factor of landslides for ASTER and SRTM based models, while aspect came out as the most crucial factor for ALOS PALSAR and SOB based models.
When all (15) causal factors were used, we found a total of four to eight statistically significant causal factors (see Table A2 of Appendix E). Slope, elevation, SPI, and aspect were the significant DEM-based causal factors for LR_SRTM and LR_ALOS based models. For LR_SOB, aspect was the statistically significant DEM-based causal factor. When 15 causal factors were used, the model assessed the interaction of DEM-based factors with the common eight factors. Therefore, when only the DEM-based causal factors were used, some factors came out statistically significant (e.g., SPI for LR_ASTER_DEM). When all the factors were used in the LR_ASTER based model, SPI came out insignificant. For LR_SRTM and LR_ALOS, three DEM-based causal factors were detected as statistically significant; therefore, for these two landslide susceptibility maps, DEM would have more impact than the LR_ASTER and LR_SOB susceptibility maps.

Landslide Susceptibility Maps (LR)
We used the Jenks natural break method to classify the probability of landslides into five zones: very low, low, moderate, high, and very high. The spatial appearances of the LR_SOB_DEM model ( Figure 2h) have a different appearance than the other three maps. LR_SRTM_DEM (Figure 2f) and LR_ALOS_DEM (Figure 2g) have an almost identical spatial appearance. While in LR_ASTER_DEM models, the slope had a higher coefficient (ß = 0.26) (Table A2 of Appendix E) than the SRTM (ß = 0.14) and ALOS PALSAR ((ß = 0.05). Therefore, in the LR_ASTER_DEM map, the areas with steeper slopes, mainly in the mid-north and mid-south of the study area, were classified as high to very highly susceptible. In LR_SRTM_DEM and LR_ALOS_DEM maps, these areas were classified either as moderate or high susceptibility zones. In LR_SOB_DEM, only elevation and aspect were two significant factors. Therefore, the susceptibility map took the shape of the map of these two factors (see Appendix A).
When all (15) causal factors were used, the spatial appearance of LR_SOB (Figure 3h) was different from the other three maps (Figure 3e-g). LR_SOB map was influenced by the distance from the fault lines, distance from the road networks, and land use/land cover. Although aspect was a significant factor, the coefficient value of aspect was similar to other significant factors, and distance from the fault lines (ß = 1.07) (Table A2 of Appendix E) had a higher coefficient value than aspect. Therefore, most of the study area was classified as low or very low susceptibility zones. In the LR_ASTER model (Figure 3e), the slope had the highest coefficient (ß = 0.31) (Table A2 of Appendix E), and therefore, areas with steeper slopes were classified as high or very high susceptibility zones. But in LR_SRTM and LR_ALOS models, the slope had comparatively lower coefficient values than ASTER. As a result, some areas were classified as moderate susceptibility zones in these two maps (Figure 3f-g). In LR_SRTM and LR_ALOS susceptibility maps, common factors such as distance from the road networks and fault lines, land use/land cover, and land use/land cover change did not have higher coefficient values than the DEM-based causal factors. That is why, unlike LR_SOB, the spatial appearance of the susceptibility maps did not follow the appearance of the maps of common factors.

Susceptibility Assessment using Random Forest (RF)
RF_ASTER_DEM, RF_SRTM_DEM, and RF_ALOS_DEM models detected (Figure 4a) slope and RF_SOB_DEM model identified the aspect as the most critical factor. When all 15 factors were used in the models, RF_ASTER and RF_SRTM (Figure 4b) detected slope and the rest of the two models detected distance from the road network as the most important causal factor. For RF_ASTER, RF_SRTM, and RF_ALOS, DEM-based factors such as elevation, TWI, and aspect had higher importance in the models than the common factors. But in the RF_SOB model DEM-based factors had less importance than the common factors. There is no similarity among the models in detecting the importance of DEM-based causal factors. For example, In RF_SOB, the slope was ranked as one of the least important factors, but for other models, it was ranked as the most important factor (Figure 4b).
RF_ASTER_DEM, the difference of variable importance between slope and other factors was comparatively higher than the other models, the effect of slope on the susceptibility map was visible.
RF_ASTER; RF_SRTM, RF_ALOS; and RF_SOB models (Figure 3i-l), spatial appearances were different from each other. Since in the RF_SOB model, distance from the road networks was the most crucial factor, areas near to roads were classified as high or very high susceptibility zones. Distance from the road network was not ranked as the most critical factor in the other three models.   (15) factors used in the models. Pl=plan curvature; Pr= profile curvature; LC= land use/land cover; LLC= land use/land cover change; DF= distance from the fault lines; DD= distance from the drainage networks; DR= distance from the road networks.

Landslide Susceptibility Maps (RF)
Like MFR and LR, we used the same method to classify the probability of landslides into five susceptibility zones. The spatial appearance of the RF_SOB_DEM susceptibility map (Figure 2l) was different from the susceptibility maps of the other three models (Figure 2i-k). For, RF_ASTER_DEM (Figure 2i) areas in the mid-north to mid-south were classified as high or very high susceptibility zones. While in RF_SOB_DEM and RF_ALOS_DEM the same areas were classified as moderate susceptibility zones. In RF_ASTER_DEM, slope was the most critical factor. Similarly, in RF_SRTM_DEM and RF_ALOS_DEM slope was the most crucial factor, but in these two models, the contribution of the slope (Figure 4a) to the model is lesser than the RF_ASTER_DEM model. In RF_ASTER_DEM, the difference of variable importance between slope and other factors was comparatively higher than the other models, the effect of slope on the susceptibility map was visible.
RF_ASTER; RF_SRTM, RF_ALOS; and RF_SOB models (Figure 3i-l), spatial appearances were different from each other. Since in the RF_SOB model, distance from the road networks was the most crucial factor, areas near to roads were classified as high or very high susceptibility zones. Distance from the road network was not ranked as the most critical factor in the other three models.

DEM-Based Causal Factors
When only DEM-based seven causal factors were used for MFR models, among all DEMs the MFR_SRTM_DEM model gave the superior performance for both success (AUC = 80.73%) and prediction (AUC = 77.37%) rate curves (Figure 5a,b). The MFR_SOB_DEM model appears to perform the weakest in assessing success and prediction. The AUCs of success rate curves (Figure 5a  For LR, LR_ALOS_DEM outperformed the other three models (Figure 5c,d). LR_SOB_DEM presented the weakest performance among the four models and thus fell under the fail category. The AUCs of success and prediction rates (Figure 5c,d) show that the other three models are under the fair category. For RF models, we got similar results as the LR model. RF_ALOS_DEM outperformed other models, and RF_SOB_DEM was the least accurate model. RF_SRTM_ALOS_DEM and RF_SRTM_DEM gave an almost similar performance.

All Causal Factors
We found that when all 15 causal factors are used, different models showed variable performances. For the MFR model, the use of 15 causal factors decreases the predictive performance on an average by 5% of landslide susceptibility maps based on ASTER, SRTM, and ALOS PALSAR (see Figure 6a,b).
However, for the SOB based model, it increased the performance by 3%. It indicates that for bivariate models, DEM-based causal factors can give better prediction performance and use of non-DEM-based factors can reduce the accuracy. Inclusion of more DEM-based causal factors and more landslide locations may increase the accuracy of the models in the study area.   For the LR model, the use of 15 factors increased the accuracy of the model. For three global DEMs success rates increased by around 5.0% but for SOB it increased by 21.7%. On the other hand, prediction rates showed the same trend as for three global DEMs the increase of performance was around 3.5% but for SOB it was 17.0%. It proves that the use of common factors increased the accuracy substantially for SOB DEM. Like the LR model, for RF models, the use of 15 causal factors increased the accuracy of the model. For ASTER and SRTM the increase of the success rate was around 7%. For ALOS PALSAR the success rate increased by 12.4%. Here again, SOB had the highest increase (22.6%) in success rate. The increase in prediction rate was not as high as the success rate. For three global DEMs prediction rates increased by around 2%-3% and for SOB the prediction rate increased by 10.1%. Machine learning algorithms such as the random forest learn the behavior or the training data. Therefore, the increase of success rate due to inclusion of new variables was high. Since machine learning algorithms learn the behavior of the training data it fails to predict the validation or unknown data [60]. Therefore, in our study, the increase of prediction rate is around 50% lower than the increase of prediction rate for RF models.

Spatial Comparison of Landslide Susceptibility Maps
Spatial convergence indicates how much area is classified into same susceptibility zones. When seven DEM-based factors were used in the MFR model, MFR_SOB_DEM had 30% of spatial convergence while the remaining DEMs had around 40% (Table 3). As we discussed before, the landslide susceptibility maps of MFR_SOB_DEM has a different spatial appearance (Figure 2d) than the rest of the three susceptibility maps and these results (Table 2) support the previous discussion of our study. For the LR models, Table 3 shows a similar trend. Spatial convergences of the LR_ASTER_DEM, LR_SRTM_DEM, and LR_ALOS_DEM were around 44%, while the LR_SOB_DEM showed approximately 19% of spatial convergence. For the RF models, Table 3 shows the similarities with the findings of MFR and LR. When all factors were considered for modeling, spatial convergence between the DEMs (Table 3) increased around 40% for MFR models. While for the LR and RF models, the spatial convergence was approximately 25% and 12%, respectively. In the MFR model, all causal factors were used, while in the LR model, significant causal factors were used and in the RF model, 2-3 causal factors, for example, profile and plan curvatures had comparatively low or no variable importance in the model.
The results of Friedman tests (Table 4) show that in both the scenarios (a. seven DEM-based causal factors used, and b. 15 causal factors used) for MFR and RF models P < 0.05. It means at least one of the landslide susceptibilities models had significantly different performance than the rest of the models. While for LR models, when seven DEM-based causal factors were used, at least one of the models was statistically different in performance than the rest of the models. When all factors were used in LR models, there was no statistically significant difference in performance between the models. Wilcoxon signed-rank test conducted the pairwise comparison. The results (Table 4) show that in scenario two, the performances of landslide susceptibility maps produced using SRTM and ALOS PALSAR did not have a statistically significant difference. Other than that, all the performances of the MFR based landslide susceptibility maps were statistically (α = 0.008 after Bonferroni correction) different from each other. When seven causal factors were used for LR models, the performance of SOB based models was significantly different (Table 5) from all other models. But when all factors are used, these differences become insignificant. It indicates that eight common factors overshadow the effect of SOB based causal factors. In the MFR model, it did not happen since it did not consider the interaction of causal factors. In the case of RF models, RF_ALOS_DEM (Table 5) was statistically different from the other models. But when all causal factors were used the difference of performance became insignificant for RF_SRTM model. RF used a more complex algorithm than the LR and MFR models. Therefore, the Wilcoxon Signed-Rank test gave different results for RF models than MFR and LR models.

Discussion
This paper evaluates the suitability of three available global DEMs: ASTER, SRTM, and ALOS PALSAR and a local DEM: SOB for landslide susceptibility mapping in Rangamati district, Bangladesh. Causal factors derived from ASTER and ALOS PALSAR DEM have been used in landslide susceptibility mapping in different parts of the Chittagong hilly areas, Bangladesh [7,9,[31][32][33][34][35]. Since the study areas of these studies were different, we could not compare them to find out which DEM gives the best accuracy in the prediction of landslide susceptibility [36]. Our study showed that three global DEMs outperformed the local DEM in all three landslide modeling scenarios.
In the first scenario, only DEM-based causal factors were used in modeling and for MFR models, MFR_SRTM_DEM outperformed the other three models. The difference of AUCs of both the success and prediction rate curves between _SRTM_DEM and ALOS_DEM was around 3%, indicating a similarity in predictions. For the processing of ALOS PALSAR DEM, SRTM GL1 data is used for radiometric correction [71]. Therefore, the quality of ALOS PALSAR DEM depends on the quality of SRTM DEM. SOB DEM-based MFR did not show a good performance for the study area and the prediction performance can be improved when a more representative landslide inventory with more landslide locations is used. We found that the ASTER DEM-based MFR model showed a weaker performance than the other two open-source global DEMs. Our result is consistent with other studies that utilized ASTER DEM [29,72,73]. It may happen because ASTER DEM contains many artifacts such as the presence of peaks in the flat terrain, and it ultimately affects the landslide susceptibility map [29]. The poor performance of the SOB DEM can be attributed to the interpolation methods that were used to extrapolate elevations from the available spot heights in the hilly parts of Bangladesh [74]. It affected the accuracy and quality of DEM in the hilly parts of Bangladesh. On the other hand, for LR and RF, ALOS PALSAR based models outperformed the rest of the models. Here, again the difference between LR_ALOS_DEM and LR_SRTM_DEM was low for success and prediction. In all cases, SOB based models gave the worst performance and causal factors derived from SOB DEM cannot explain the landslide susceptibility of the study area. For example: In LR_SOB_DEM, the LR model used two significant causal factors: elevation and aspect (Table A2 of Appendix E). The low coefficient (ß) values of these two factors indicate that these two causal factors cannot adequately explain the landslide susceptibility of the study area.
In the second, scenario, for MFR models, the use of 15 causal factors increased the prediction accuracy for MFR_SOB. But for the other three models, it reduced accuracy. It indicates that causal factors derived from three global DEMs were capable enough to explain the landslide susceptibility of the study area. MFR is a bivariate model, and it does not consider the interaction of the causal factors. Moreover, unlike LR and RF models, it does not require non-landslide (pseudo absence point) in modeling [44]. When all causal factors were used, PRs of some of the common causal factors were comparatively higher than the PRs of the DEM-based causal factors. For example, PRs of distance from the road networks were around 6.00 (Table A1 of Appendix D), while PRs of the DEM-based causal factors ranged from 1.00-4.43. Therefore, the higher PRs of the common causal factors overshadowed the DEM-based factors [6,35,44]. On the other hand, the quality of the SOB DEM was poor and was not capable of explaining the landslide susceptibility of the study area. That is why, in the MFR_SOB model, the prediction accuracy increased, but the increase was very low. For LR and RF models, in second scenarios, the prediction accuracy increased by 2%-3% for global DEM-based models, while for SOB based models, it exceeded 10%. Both LR and RF models consider the interaction of causal factors, and therefore, in these two models, the effects of common factors on prediction performance were revealed better than the MFR models.
The findings of this study have similarities with other research where the suitability of different global and local DEMs was evaluated for landslide susceptibility mapping [26][27][28]. In most of the studies, ASTER DEM-based landslide susceptibility maps were outperformed by the other global DEM-based landslide susceptibility maps. On the other hand, local DEM-based landslide susceptibility maps have better prediction accuracy than the global DEMs [1,29]. In these studies, local DEMs were mainly light detection and ranging (LiDAR)-based DEMs, which generally have very high spatial resolution compared to global DEMs [29]. SOB DEM was prepared under the project titles "Improvement of Digital Mapping System of Survey of Bangladesh", where the main aim was to prepare a 1:5000 scale DEM of major cities in Bangladesh. This SOB project prepared DEM for the hilly areas of Bangladesh, including our study area using the local spot heights. Different interpolation methods were used to prepare 25m DEM from these local spot heights [73]. Therefore, in hilly areas of Bangladesh, the quality of SOB DEM is not good enough, and the application of this DEM in geomorphological studies such as landslide susceptibility mapping can give questionable results similar to what we got in our study. Our study also pointed out the importance of the preparation of very high-resolution DEMs using LiDAR for the hilly areas of Bangladesh. It will help in various geomorphological studies, including landslide susceptibility mapping. As we did not find a substantial difference among global DEMs, any global DEM can be utilized for landslide susceptibility mapping in Bangladesh in the absence of very high-resolution DEMs.
For an in-depth study, we utilized non-parametric tests. Chang et al. [30] used non-parametric tests to detect a significant difference in performance in landslide susceptibility maps prepared using different DEM-derived causal factors. Their study did not find any significant difference in performance for two machine learning methods: RF and support vector machines. In our research, we did not see any significant difference in performance for RF and LR based models. But for the MFR model, we found a significant difference in performance. It indicates that the effect of DEM-based causal factors on the performance of bivariate models is more than on the multivariate and machine learning models. Thus, we suggest using multivariate and machine learning methods and any one of the global DEMs for landslide susceptibility mapping in Bangladesh.

Conclusions
This paper assesses the effects of the DEM-derived causal factors on the landslide susceptibility maps produced using the bivariate (e.g., MFR), multivariate (e.g., LR), and machine learning (e.g., RF) models. In this study, we tested two scenarios: a. susceptibility assessment with only seven DEM-based causal factors; b. inclusion of other 8 causal factors along with DEM-derived factors. The success and prediction rate curves showed SRTM DEM outperformed under an MFR model in both scenarios. Our analysis revealed that ALOS PALSAR DEM provided the best prediction accuracy while using both LR and RF models in landslide susceptibility mapping. For all models and scenarios, the SOB DEM does not perform well compared to other DEMs.
The prediction accuracies of landslide susceptibility mapping using only DEM-derived factors is low compared to the utilization of all casual factors. Although SOB DEM has a poor performance in susceptibility assessment with only a DEM-derived factor, the accuracy is significantly improved when other non-DEM-based factors are added. Therefore, we argue that causal factors derived from the SOB DEM have limited influence in landslide susceptibility assessment for the study area. Besides, for the LR and RF models, the use of all the causal factors increased the prediction performance. Spatial convergence analysis showed that three global DEM-based models have similar accuracies and performed far better than SOB DEM-based models. Therefore, we recommend that SOB DEM should not be used for landslide susceptibility assessment in Bangladesh. Although ASTER DEM-based models showed the weakest performance among three global DEMs, the difference of performance was negligible. Therefore, we recommend any one of these three global DEMs can be utilized for landslide susceptibility assessment in Bangladesh.
Our study also highlights that, for the MFR model, DEM had the highest impact on the accuracy of landslide susceptibility assessment as Wilcoxon rank tests showed that the performance of susceptibility maps was significantly different. For the LR and RF models, the effect of DEM was less significant. We contend that while using the bivariate model, we must be careful about the quality of DEM. Two scenarios in this study helped to understand the impact of DEMs in landslide susceptibility assessment. Moreover, we used multiple metrics to evaluate the accuracies of susceptibility assessment including AUC, spatial convergence analysis, and non-parametric test. Although the use of one performance measure is a common practice, the utilization of various measures helped this research to understand the impacts of DEM and makes the conclusion robust.
The main limitation of our study is that we cannot ascertain if adding more DEM-derived factors such as terrain roughness could improve the result. The impacts of DEM on landslide susceptibility maps are not universal and may vary from place to place. Therefore, we cannot conclude that a specific DEM can be better than others. Advanced machine learning and deep learning methods can be used to check whether these algorithms can reduce the differences in prediction performance of landslide susceptibility maps prepared using different DEM-derived causal factors.