Landslide Susceptibility Mapping Using Bivariate Statistical Models and GIS in Chattagram District, Bangladesh

Landslide is one of the most devastating hazards in Chattagram Dsitrict and has become a recurrent phenomenon in this region. This study attempts to produce Landslide Susceptibility Map (LSM) for Chattagram District of Bangladesh by using five GIS based bivariate statistical models, namely the Frequency Ratio (FR), Shanon’s Entropy (SE), Weight of Evidence (WofE), Information Value (IV) and Certainty Factor (CF). Landslide Inventory (2001–2017) of Chittagong Hilly Areas database was used to measure the relationship between the previous landslides with the landslide conditioning factors. SRTM DEM and Landsat satellite images were collected from USGS and the geological data were collected from GSB to produce the thematic layer of conditioning factors. Sixteen landslide conditioning factors of Slope Aspect, Slope Angle, Geology, Elevation, Plan Curvature, Profile Curvature, General Curvature, Topographic Wetness Index, Stream Power Index, Sediment Transport Index, Topographic Roughness Index, Distance to Stream, Distance to Anticline, Distance to Fault, Distance to Road and NDVI were used. The Area Under Curve (AUC) was used for validation of the LSMs. The predictive rate of AUC for FR, SE, WofE, IV and CF were 76.11%, 70.11%, 78.93%, 76.57% and 80.43% respectively. CF model indicates 15.04% of areas are highly susceptible to landslide. All the models showed that the high elevated areas are more susceptible to landslide where the low-lying river basin areas have a low probability of landslide occurrence. The findings of this research will contribute to land use planning, management and hazard mitigation of the CHT region.

Remote sensing and Geographical Information Systems (GIS) along with statistical models plays a significant role in landslide hazard zonation which is very important to reduce loss and damage in an area (Dahal et al. 2008;Ozioko and Igwe 2020). There are many popular statistical methods frequently been using for landslide susceptibility mapping for case studies at local scale using GIS environment. These can produce more accurate and reliable results in the complex topographic environments (Youssef et. al. 2015;Youssef et al. 2016;Yang et al. 2019). Present study employed Frequency Ratio (FR), Weight of Evidence (WofE), Shanon's Entropy (SE), Information Value (IV) and Certainty Factor (CF) models in GIS environment for landslide susceptibility mapping. Still now these models are used extensively and many researcher had get reliable and acceptable result using these models (Park et al. 2012;Youssef et al. 2015;Ba et al. 2017;Wen et al. 2017;Passang and Kubicek 2018;Elmoulat and Brahim 2018;Nohani et al. 2019;Silalahi et al. 2019;Wang et al. 2019).
In this research thematic layers of landslide causative factors such as Slope Aspect, Slope Angle, Geology, Elevation, Plan Curvature, Profile Curvature, General Curvature, Topographic Wetness Index (TWI), Stream Power Index (SPI), Sediment Transport Index (STI), Topographic Roughness Index (TRI), Distance to Stream, Distance to Anticline, Distance to Fault, Distance to Road and Normalized Difference Vegetation Index (NDVI) were used and correlate with landslide events. Though landslides of this region got attention of many researchers but no attempt had been taken to implement Bivariate statistical models and analysing so many causative factors for landslide susceptibility mapping in this region. As the result of Bivariate Statistical models varies from region and model to model so that five models were employed to produce landslide susceptibility map (LSM) and the accuracy of these maps were calculated using Area Under Curve. So the objectives of the research goes to calculate the degree of relationship between landslide conditioning factors and previous landslides by using five GIS based bivariate statistical models, identify the areas vulnerable to landslide according to the models, and calculate the prediction rate of these models in the Chattagram District.

Study Area
Chattagram District is located in between 21°54' and 22°59' North latitudes and in between 91°17' and 92°13' East longitudes (Figure 1). The Chattagram Hill Tracts were uplifted during Miocene period as a result of collision between the Indian and Eurasian Plates. Soils of the hills are composed of Tertiary and Quaternery sediments and there are many anticlines and synclines exist in this hilly region. The surface strata of this region have easily-weatherable feldspar and unconsolidated sedimentary rocks (yellowish-brown to reddish-brown loams). Geologically these Hill ranges are underlain by rocks of Tipam and Surma formation. The hill ranges have low height below 300 meter but there are steep slopes which are disconnected by rills and gullys (Brammer 1996;Sarkar and Rashid 2013;Brammer 2016). Surface soils of the region are dominated by Cambisols (Moderately developed soils) and Gleysols (Groundwater affected soils) (FAO 1988). Climatologically CHT region belongs to tropical monsoon climate with mean yearly rainfall 2540 mm to 3810 mm which is the main triggering factor of landslide (Ahmed 2014).

Methodological Summary
Flow chart of the broad steps followed in this research for deriving results of landslide susceptibility assessment and the model used for validation (AUC) are given below (Figure 2).

Data Processing
In this section, authors describe data processing and preparation method of the landslide conditioning factors including the background theory where necessary. Data were processed by using ArcMap 10.5 software to generate the spatial database of the conditioning factors.
Slope aspect mainly affects the amounts of sunlight and rainfall. Again wind can negatively impact slope equilibrium and indirectly influence landslides (Demir et al. 2013;Wen et al. 2017). Slope aspects of the study area were derived from DEM images by using ArcMap 105 software. They are classified into nine groups (  The slope angle controls the mass movement and gravitational energy and is directly connected to the landslide occurrence (Wen et al. 2017). DEM images were employed to calculate the slope angle of the study area and they are categorized as 0-1.29°, 1.29°-3.44°, 3.44°-6.02°, 6.02°-8. 82°,.
Three types of curvature were calculated, namely general curvature profile curvature and plan curvature. The curvatures significantly control the distribution of stresses plan curvature and the shrinking and expansion of the saturated landslide materials according to the landslide motion direction (Carson and Kirkby 1972;Guo et al. 2015). As a result, increase in convexity and concavity of any slope curvature enhances the likelihood of landslide occurrence (Lee and Talib 2005;Mersha and Meten 2020). The slope curvatures were categorized into three classes namely: convex, concave and flat figure (Figure 3d, e, f).
High elevated mountains contain weathered rock and intermediate or close to the ground are covered by thin colluvium materials which are more susceptible to landslide (Youssef et al. 2015). The elevation raster of the study area was directly derived from DEM processing by ArcMap 105 software (Figure 3c). In the analysis of landslide susceptibility distance to lineaments (faults folds joints and fractures) is considered as an effective factor (Youssef et al. 2015). In this research distance to faults ( Figure 3h) and distance to anticline ( Figure 3i) were used.
Road construction and inappropriate cut slope and drainage from road adversely impact slope stability as a result landslide occurs along roads (Dahal et al. 2008). Raster data of distance to the road (Figure 3j) was computed by using Euclidean distance. Drainage is another major cause of landslides (Elmoulat and Brahim 2018). This casual factor was derived from DEM by using hydrology tool of ArcMap 10.5 software. DEM raster analysis was executed to generate  Figure 3p) were calculated by equations 1 and 2 respectively. The high TWI value negatively impacts the slope durability causes landslide (Beven and Kirkby 1979). SPI impacts the debris flow development and this factor is extensively applied to landslide susceptibility assessment (Tehrany et al. 2019).
where α is the local upslope area draining through a certain point per unit of contour length, and β is the local slope gradient in degrees.
where A S indicates the catchment area that is studied and β indicates the local slope gradient measured in degrees.
The Sediment Transport Index (STI) (Figure 3n) is another index of measuring erosive power of the overland flow (Moor et al. 1993). To calculate it by using Universal Soil Loss Equation instead of 1D slope length the 2D catchment area is used (Kumar and Singh 2016).
where, A S is the area (m2/m) of specific catchment and β is the measurement of slope gradient (degree).
The Topographic Roughness Index (TRI) (Figure 3o) is measured as an output of surface ratio provides an index of topographic roughness. It is a very effective tool used for assessing slope stability (Jenness 2013). TRI was calculated by the Equation 4, Normalized Difference Vegetation Index (NDVI) was analysed from Landsat 8 OLI image by using AcrMap 10.5 software and classified into five classes ( Figure 3l). General equation of NDVI analysis is NDVI= (NIR-RED)/ (NIR+RED) where NIR denotes Near Infrared band and RED indicates the red band of Landsat image.

Frequency Ratio (FR)
Frequency ratio is the possibility of occurrence of an event for a given characteristics (Bonham-Carter 2014). By calculating the bivariate relationship between the landslide and their conditioning factors it reveals the degree of relationship among them (Lee and Talib 2005). Equation (5) shows the FR algorithm, where FR is the frequency ratio of a class i of a conditioning factor j, Lc is the number of landslide points of a class of a conditioning factor, Lt is total landslide point of the study area, Af is the number of pixel of a class of a conditioning factor and At is the total number of pixel in that conditioning factor.

Shanon's Entropy (SE)
It calculates the bi-variate relationships among the component of a system with the degree of disorder. This method can separate the most significant impact factor from less significant one. In landslide study entropy of a landslide indicates the degree of influence of factors to occur landslide. By using the following Equations value of the parameter is computed (Arabameri et al. 2019), where, Eij is the ratio between number of landslide occurrence in a specific class and total pixel within that class.
where H j and H jmax are the entropy of a class of a conditioning factor, I j is the coefficient of information value, M j is the number of classes in each conditioning factor, and W j is the achieved weight value for the given factor. The value of Wj ranges 0-1. The entropy of W j close to 1 has higher probability to be inconsistent and imbalanced (Khosravi et al. 2016;Nohani et al. 2019).

Weight of Evidence (WofE)
After Bonham-Carter et al. (1988) used this method for gold exploration, still now it has been used for mineral exploration integrating with GIS and also had become popular in landslide susceptibility assessment (Elmoulat and Brahim 2018;Wang et al. 2019). Statistically, WofE can be expressed by Equation 13.
WofE model calculates the positive weight (W + ) and negative weight (W -) of the indicator where W + indicates the presence of the event and Windicates the absence of the event as follows (Bonham-Carter et al. 1994). The original equation was rewritten in pixel format to calculate weights of landslide as follows: where Npix1 indicates previously occurred landslide point is specific factor class, Npix2 indicates the existence and non-existence of landslide point in the factor map, Npix3 s the existence of predictive factor but non-existence of landslide in the factor class, Npix4 is the non-existence of both causative factor class and the landslide point. (10) Relative Frequency of Goods Relative Frequency of Bads * 100 Weight of contrast (C) is measured by substituting W + and W -. It is the degree of spatial association of causative factor class and the corresponding landslides. C is calculated by Equation 16.
The weighted value of C was assigned to create factor map (Equation 17) and LSM was prepared by combining all the factor map (Equation 18):

Information Value (IV)
The information value method can calculate the probability of landslide occurrence in an unbiased way as it is considered as an indirect statistical approach (Zezere 2002). In this method, the ratio of landslide point density in a class of a factor and the landslide density of the entire map is employed and the natural logarithm of this ratio is considered as the Information of Value of that parameter class (Bui et al. 2011). IV is calculated by using equation 19.
where, I i is the weight of a given class of a factor; Densclass is the density of landslide of a given class; Densmap is the landslide density of the factor; N pix S i is the number of landslides in a certain factor class; N pix N i is the number of pixels in a certain factor class; N is all of the landslide points in the database; and n is the total number of classes. Finally all the factor maps were combined to produce the LSM using Equation 20.
where, m = number of thematic layers.

Certainty Factor (CF)
Certainty factor is the function of probability (Sujatha et al. 2012) that is modified by Heckeman in 1986 from the original version. The modified equation is as follows (Equation 21): where PPa is the ratio between the number of landslide point and the number of pixel within class a, to 1 the better the accuracy of the model (Youssef et al. 2015;Nohani et al. 2019;Vakhshoori et al. 2019).

Frequency Ratio
Results of weights of all the conditioning factors and their subclasses obtained for Frequency Ratio model tabulated in Table 1. Probability ratio of each map was used with every subclass for LSM. Equation 6 used in GIS environment to prepare the susceptibility map as follows: Results of frequency ratio showed that SPI has the highest impact on landslide occurrence in the Chattagram District. Results of SPI indicated that all the landslide points are in one class (0 -26,376,673.13) and it got the highest probability ratio of 5.171. Halda and Karnafuli river basin crossed over the Chattagram District, where the Halda Basin covers a large area causing high SPI in the study area. As a result, a high SPI value generates and most of the landslides belonging to the higher class produced a high probability of landslide for SPI. The Second highest influencing factor was TWI, the probability ratio is 4.2. Low TWI value produced in the steep hill slopes of the study area; as a result, most of the landslides concentrated in low TWI increased the probability of TWI to landslide occurrence. Slope, distance to anticline and distance to fault also had significant influence in triggering landslides in Chattagram District, followed after SPI and TWI. The probability ratios of these three causative factors were above 3. The presence of anticline and fault increased the weathering process causes deformation in soil stability. In the study area, slope played a significant role in causing landslides. Among the slope aspects, the Southwest and PPs is the ratio between the total landslide point and all the pixel in the entire factor layer. The value of CF varies between 1 to -1. Positive value indicates the increasing certainty and negative value indicates decreasing certainty of landslide occurrence while the 0 value means initial probability is very similar to the conditional probability (Wang et al. 2019). Finally the LSM produced by combining all the factor map (Equation 22): where, n= number of factor map layers.

Model Verification
The AUC (Area Under Curve) is a popular method used to identify the model's accuracy which shows the success rate of the model. The AUC indicates the ability of a model for predicting future probability of landslide occurrence and non-occurrence. If AUC value ranges between 0.5 and 1 it is considered that the model can predict future landslide perfectly. AUC value of 1 indicates perfect prediction and the more closely the value

Shanon's Entropy
The results of Shanon's Entropy for the relationship of the effective factors with the occurrence of landslides are presented in Table 1. Shanon's weight Wj was calculated for all the factors. Wj value is the relative percentage of the importance of each map. Finally Equation 11 and GIS environment used to prepare susceptibility map according to following formula: Legend: Nsclass = landslide point in class, nClass = pixel in the class, FR = frequency ratio, RF = relative frequency, PR = probability ratio, Pij = ratio of FR, Wj = weights of shanon's entropy, C = weight of contrast, IV = weights of Information Value model, CF = weights of certainty factor This model showed slope aspect, general curvature, STI, distance to stream and NDVI had the highest impact on landslide occurrence and the Wj value of these variables are ≥ 0.07. Slope, elevation, geology, plan curvature, profile curvature, TRI and distance to road had Wj value > 0.06. TWI, distance to anticline and distance to fault had Wj value >0.05 and SPI got the lowest value of 0.038 in this model. Wj value is the relative percentage of the importance of each map. Shanon's Entropy considers the frequency of each subclass of a factor. As a result, it greatly reduces the unevenness among the causative factors. Because of this homogeneity importance of a parameter differs from that of the Frequency Ratio.

Weight of Evidence
Results of weights of evidence model were tabulated in Table 1 In terms of individual factors, Slope, Elevation and Geology had a strong positive effect on landslides and in terms of individual class slope degrees between 25º and 54 º degrees had a strong positive impact on the landslide. But in this subclass, both the number of class pixels and the number of landslides were smaller than other subclasses. As a result, a small number of landslides increased the weight both in the WofE and FR model. Geological subgroups of Tb (C = 2.0), Tbb (C = 2.1) and Tt (C = 1.2) were more unstable than other subgroups. After the previous three causative factors, STI, distance to stream and slope aspect had a relatively more positive relation with landslide occurrence. Distance to anticline the subclass of <4000 meters had a significant positive impact on the landslide. Results showed that in the case of distance to anticline, the subclass of 0-2000 meter had a weighted value of 1 and 2000-4000 meter had a weighted value of 2. Comparing and considering the three curvature, it was evident that flat curvature had no impact on landslide occurrence where concave curvature (C = 1.26 for plan, 0.31 for profile and 0.92 for general curvature) were relatively more susceptible to landslide than convex curvature (C = 0.34 for plan, 0.7 for profile and -0.131 for general curvature). In the case of NDVI, it was seen that a relatively high value showed a more strong positive impact on landslide occurrence which means in terms of the whole Chattagram District, landslide occurred not only in the populated area but in vegetated areas of the hills. But the most devastating landslides have occurred in the Chattagram Municipality area (Ahmed 2015;Sultana 2020). Subclass 3500-7400 of distance to fault line has a high weighted value of 2.2 and 95% of the landslide occurred in the distance between 0-7400 meter. WofE model calculates both positive and negative weight of a class based on the information of landslide where the influence of the class is also determined by the percentage of area of the whole factor. It is seen evidently in the negative values of the classes. Compared to FR method, WofE methods give negative weight to a class more precisely. Such as, TWI class (11.392 -14.518) and (14.519 -24.659) have a probability ratio of 0 in the FR model, but they get -1.57 and -1.22 in the WofE model (Table 1). These types of odd relationships are also seen in Slope Aspect, Geology, TWI, SPI, TRI, Distance to Anticline, Distance to Fault, Distance to Road and NDVI, where some classes have no landslides. A similar pattern is also seen in the IV and CF models. The weight of the classes produced in the WofE model largely resembles with CF model.

Information Value
The Information Value method was utilized to calculate the relationship between landslide location and the landslide conditioning factors. The results of weights of each class of landslide conditioning factors of this model tabulated in Table 1. The final LSM is produced by using GIS and Equation 22.
This model showed that slope degree was the most important factor in landslide occurrence. Like the weight of evidence model slope class of 25-54 showed a strong positive impact on the landslide and slopes above 3.5º were marked as more susceptible to the landslide. With increasing slope degree it has become unstable. In the case of geological factors 87% landslide occurred in four subgroups namely Tb, Tbb, Tt and Qdt and they had a weighted value of 0.80, 0.77, 0.41 and 0.43 respectively. Elevation (25) LSM IV =Aspect IV + Slope IV + Geology IV + Elevation IV + PlanCurvature IV + ProfileCurvature IV + GeneralCurvature IV + TWI IV + SPI IV + STI IV + TRI IV + Dis_Stream IV + Dis_Anticline IV + Dis_Fault IV + Dis_Road IV + NDVI IV the landslides occurred within the subclass of 0-500 meter but this subclass showed a negative impact to landslide occurrence. Comparing the landslide conditioning factors slope degree, elevation and geology had a more positive impact on landslide and STI and distance to stream and slope aspect had the secondhighest impact. The weight of the IV model largely resembles with WofE model. Both of the models use logarithmic transformation to calculate the weights of each class. In the case of Geology, like all other models, the priority weights of the classes vary in this model. It is because they cover a very small area in the factor. Tb, QTg, QTdt and Tbb covers 2.2%, 3.17%, 4.41% and 5.57% areas respectively (Table 1). This small difference in areas and the ratio landslide density produces different weights in different models.

Certainty Factor
Results of Certainty Factor were tabulated in Table 1. Certainty factor, Information value and weight of evidence method showed similarities in weighting values of landslide conditioning factor. Equation 27 and GIS finally used to prepare the map according to the following procedure: Considering the values of the slope aspect, the Southwest facing class had the most significant impact (0.597) on the landslide. Other important classes having a positive impact on landslides are sequentially South (0.405), Southeast (0.345) and North (0.156) facing slope. Among the factors slope was the most significant impact on the landslide. With the increase of slope angle, the impact of it to landslide was also increased. Such as classes of 3.44-6.022, 6.022-8.819, LSM CF =Aspect CF + Slope CF + Geology CF + Elevation CF + PlanCurvature CF + ProfileCurvature CF + GeneralCurvature CF + TWI CF + SPI CF + STI CF + TRI CF + Dis_Stream CF + Dis_Anticline CF + Dis_Fault CF + Dis_Road CF + NDVI CF indicated that landslides started to occur above 15 meters height and 60% of the landslides occurred within 15-50 meter range and 75% of landslides occurred within a range of 15-75 meter and no landslide occurred above 160 meters. It indicates hills of the intermediate position were more unstable and more prone to landslide as they were composed of loose material. In terms of distance to stream, 50% of 8. 819 -11.830, 11.830 -15.272, 15.272 -19.574, 19.574 -25.811 and 25.811 -54.849 had a weighted value of 0.310, 0.649, 0.771, 0.763, 0.785, 0.873 and 0.973 respectively. In the case of Geology Tb, Tbb and Tt had shown a positive impact on landslide and the rest of the class showed a negative impact on the landslide. After slope, another important factor was Elevation. Five subclasses out of seven had shown a positive impact on landslide and the class 1008.01-159 had the greatest impact (0.868) on landslide. Among curvature, the concave and convex slope had a positive impact on landslide and the convex slope had a relatively high impact on the landslide. In the case of SPI only one class showed no relation with landslide and the rest of the class showed a negative impact on the landslide. Among the other factors, STI and Distance to Stream showed a more positive impact on landslide and Distance to Road showed a very little impact as the only class (0-322.103) had a value of 0.244. In the case of NDVI the higher value showed a more positive impact and the medium to very low value showed a negative impact on the landslide. The weight of individual classes in the CF model largely resembles with FR method. But CF model gives negative weight to classes where the last limit of FR is 0. In measuring negative values, the CF method varies from the FR model. In the case of SPI, it is seen that the 0 weights of FR get -1 weight in CF (Table 1). In the case of Geology, weights of CF vary from FR (Table 1). FR measures the relationship in correspondence to the class area; as a result, small areas with a higher number of landslide gets high weight. CF method calculates the weight comparing the conditional probability of the class and the probability of the entire factor. So that priority class differs from FR and other models.
Finally five LSM were prepared for FR, SE, WofE, IV and CF models (Figure 4). The LSM of the study area were categorized into five classes according to landslide vulnerability level, namely Very Low, Low, Medium, High and Very High. Percentage of areas containing these classes illustrated Figure 6.

Model Verification
Model verification was completed by using Area Under Curve (AUC) techniques. Prediction rate for AUC of the models FR, SE, WofE, IV and CF are 76.11%, 70.11%, 78.93%, 76.57% and 80.43% respectively ( Figure 5). Considering the results of evaluation of the five models CF showed most accurate prediction rate for identification of landslide susceptible areas in the Chattagram District. The least accurate prediction in identification of landslide found for SE model. The second best prediction rate showed WofE model. The FR and IV models had AUC of 76.11% and 76.57% respectively showed similar prediction rate in the study area.

Discussion
In the present study CF (80.83%) model produced excellent prediction of the landslide occurrence (Figure 5). Other models also produced acceptable result and SE produced the lowest accuracy among the models. The CF model showed 15% areas are at very high risk for future landslide occurrence in the study area ( Figure 6 and Figure 4). Very high landslide susceptible area found    Figure 6). Two previous studies were conducted over the entire Chattagram District to produce a landslide susceptibility map (Mourin et al. 2018;Chakraborty 2019). Findings of this research are slightly dissimilar to the previous study conducted by Mourin et al. (2018) and Chakraborty (2019). Low and very low susceptible area calculated from the present study was 40-48% where in previous study marked 60%. Chattagram District belongs to a tropical monsoon climate. It receives a huge amount of rainfall during the monsoon and receives considerable torrential rain in pre and post-monsoon seasons. Generally, northfacing slopes are comparatively wetter because of the less solar radiation (Pham et al. 2018). Because of the geographical location of the Chattagram District, the Southeast, south, and southwest facing hill slopes receive most of the rainfall. The condition makes the soil more saturated than the soil of the north-facing slope. These slopes also receive a high amount of sunlight, and heat increases weathering process. Moist air also hits these slopes directly almost all the year. So these slope aspects face more weathering compared to the other slope aspect. As a result, southeast, south, and southwest slopes are more prone to landslides. Low elevated areas are more vulnerable to landslides in the Chattagram District. 61.8% of landslides occurred below 50-meter height, whereas 40% occurred between 15 and 50 meters. This range of elevation is more favorable for weathering process because of the same reason of aspect. In the CHT region, soil moisture decreases with increasing slope and elevation because of the faster runoff (Emran et al. 2018). This resulted in a significant amount of landslide occurrence below 15-degree slope and 50-meter elevation.
Most landslide susceptible geological units are Boka Bil formation (Tbb), Tipam Sandstone (Tt) and Bhuban Formation (Tb), having landslides of 32.5%, 28.9% and 13.9%, respectively. The general soil composition of this region contains sandstone, shale and siltstone. These types of soil are less stable as they have a greater capacity for water retention and have high permeability. The upper layer of Boka Bil formation consists of bluish-grey shales. These shales have irregular fractures and also contain a high amount of calcite. The middle part (75-350 feet) contains sandstones, shales, siltstones and sometimes fossiliferous conglomerates. The lower unit has bluish shales with shaly sandstone. This type of structure is volatile because of the concentration of sand in the middle part. The Tipam sandstone is composed of coarse and granular sandstone. Along with an interlayer of clay, Tipam sandstone also contains fossil wood and lignite coal. The soil is iron-rich. The upper layer of Bhuban Formation contains bedded sandstone covered by a layer of siltstone at the top and bottom. The middle part dominated by shale and silt and the lower unit also contains some fossil shales (Quazi 1986;Brammer 1996;Sarkar and Rashid 2013). These conditions make the region more vulnerable to landslides.
The convex and concave slope has the highest control on landslides like its regular pattern, as they can retain more water to saturate the soil faster. Generally, soils having a high TWI value contain greater moisture and, therefore, a higher probability of landslide initiation (Guo et al. 2015). Landslide initiation and gully erosion are triggered by the same topographical and geotechnical conditions (Aigbadon et al. 2021). Arabameri et al. (2019) showed that gully formation is most evident between TWI values of 4.8-10.8 where he found 84.8% of the gully in this range in his study area (Semnan of Iran). Park et al. (2012) reported 93.65% landslide occurrence between TWI values of 4.5-6.6 in Inje area of Korea. Guo et al. (2015) found 81.23% landslide below the TWI value of 9. In Chattagram District, 83.5% of landslides occurred in the TWI class of 3-6.8 and 96.4% landslides occurred below the TWI value of 8.9. According to landslide occurrence, areas with low TWI values are more erosion-prone and susceptible to landslides in the study area. This is because the high TWI values generated in the Halda basin located in the middle of the study area and low TWI values generated in the high elevated regions (Figure 3). All the landslides in the study area occurred in very low SPI values and Park et al. (2013) also reported the same result. Landslide density in the study area is concentrated in medium to high TRI indices. Nakileza and Nedala (2020) also report the same relationship. Low TWI, low SPI and high TRI values were generated in the hilly areas. The Halda basin and the other low elevated areas showed a reverse pattern for their flat topography.
Landslide gradually increases with increasing distance to stream in the study area where 92.8% of the landslides occurred below 920-meter distance to stream and 71% of them occurred below 700 meters. The landslides are distributed almost equally in the subclasses into this range. Fan et al. (2017) found that 94.41% of landslides occurred above 600 meters and Ozioko and Igwe (2020) found 95.4% landslide occurrence within a 500-meter distance to stream in their study area. Distance to stream does not always have significant control over landslides. The result of this research rejects the findings of Pham et al. (2018) about control of stream on landslide occurrence but resembles with results of Nohani et al. (2019). The distance to stream does not significantly control the landslide in the study area.
Low distance to anticline (<4000 meters) has highest landslide concentration 98.5%. In the study area, 98.5% of the landslide occurred in the first two classes of distance to anticline, which is within 4000 meters from the anticline. Low distance to anticline supports weathering of soil that helps in triggering a landslide. Distance to fault produces natural discontinuity makes favorable conditions to landslide occurrence (Nohani et al. 2019). Nohani et al. (2019) found 46% landslide occurrence below 400-meter distance from fault and the rest of the landslide occurred above 400 meters. In the stud area, 96.4% of landslides occurred below 7375 meters from the fault. It seems apparent that distance from anticline and fault has a significant influence on landslide from the result, but the effect is very low compared to the other study (Guo et al. 2015;Nohani et al. 2019).
Distance from road showed a regular pattern in triggering a landslide. In the study area, 73% of landslides occurred below 322 meters from the road. Road distance has a significant impact on landslide occurrence in the study area. The study area is the port city of the country and there is a growing urban area developed by cutting hills. Most of the roads are also constructed by cutting hill slopes and there is higher road density evident in the study area (Figure 3j). Guo et al. (2015) reported 69.6% landslides below 500-meter distance to stream and Nohani et al. found 67.6% landslides below 400 meters. Ozioko and Igwe (2020) reported a reverse pattern of landslide occurrence with the road distance where only 21.2% landslide found below 500-meter distance from the road. Hill cutting for road construction produces discontinuity regular slope patterns and makes the soil more porous along with road results in a higher concentration of landslide. Sometimes this discontinuity itself reduces the slope's shear strength, removing the slope's basement. This condition makes the slope unstable above the road to fall without or with very little triggering force. As a result, minimal amount of rainfall can produce large landslides in the study area. Sometimes road diverts the natural water flow from a large hillslope during rainfall and concentrates it into a common route. This excessive flow in a single course produces landslides.
Studies reported that low NDVI is highly prone to landslide (Youssef et al. 2015;Guo et al. 2015;Nohani et al. 2019). In the study area, NDVI showed a reverse pattern. Most of the landslides occurred in the densely vegetated area. 78.9% of landslides occurred between Medium to High NDVI class. It is because the landslide flow material is produced in the vegetated areas due to external disturbances. It also needs to consider that the whole study area is densely vegetated and only 12.27 percent area is covered by water bodies and settlements (Table 1 and Figure 3.l).
The FR models showed that among the causative factors, SPI, TWI, Distance to Fault, Distance to Anticline and Slope had the significant causing landslides in the study area. At the same time, the SE model showed higher weight on Distance to Stream, Aspect, STI and General Curvature. The CF models sequentially showed the highest priority to the subclasses of OTdt of Geology, Slope classes of 25.812 -54.849 and 19.574 -25.812, Elevation of 108.01 -159, Tb and Tbb of Geology, STI class of 378.08 -697.34, Slope of 15.272 -19.574. In the study area, slope cutting to build infrastructure like roads and houses in less steep slopes causes most landslides. There is also an impact on the characteristics of Geological subgroups.

Conclusion
In this research bivariate statistical models of FR, SE, WofE, IV and CF was employed integrating with GIS software (ArcMap 10.5) for landslide susceptibility mapping. The causative factors of slope degree, slope aspect, elevation, geology etc. have great influence on landslide occurrence in Chattagram District. The excellent prediction rate found for CF model of 80.43%. The prediction rates of other models are also acceptable. CF model indicates 18.17%, 24.67%, 24.85%, 17.26% and 15.04% areas are very low, low, high and very high vulnerable respectively. All the models showed that only the high elevated areas are more vulnerable to landslide but the low-lying river basin areas have low probability of landslide occurrence. LSM produced in this research will be helpful for policy makers and development planners to avoid more damaging areas suggested in the susceptible maps ( Figure 4). Forest resource management authority and environmental conservation unit will also be benefited from this research. If more precise and high resolution remote sensing data can be used accuracy of the model will be increased.