Himalayas of Nepal is prone to several mass movements and landslides every year . Most of the landslides occur in rainy season during heavy and continuous rainfall. Weak geological structures, steep and rugged topography, high altitudinal variation along with monsoon triggers give rise to high degree of fragility in the whole mountain system.
Landslide is defined as massive movement of materials (soils, rocks, organics, etc.) from up slope to down slope forming materials along surfaces of separation due to the effect of gravity . Landslides are characterized by movement of mass of rock, debris or earth down a slope  . The Himalayas are one of the most active mountains which have resulted from the collision of the Indian Plate and the Tibetan Plate some 50 million years ago in which the Nepal Himalayas is the longest and youngest mountain system of the world . In the context of Nepal, Landslides are the most occurrences natural hazards in the Nepalese mountain areas  and have become one of the main natural hazards among other natural hazards like earthquake and lighting . Many hill slopes in Nepal are situated on or adjacent to unstable slopes and old landslides, which are reactivated from time to time . In recent years, globally, landslide has been previously assumed hazard, both in terms of economic losses and injuries . Not only economic loss but also the large number of death is increasing every year in the world due to landslide   .
Besides the natural processes, anthropogenic activities such as improper land use, unplanned development activities such as construction of roads toe cutting and haphazard irrigation canals without proper protection measures in the vulnerable mountain belt further intensifies the risk in Nepal . Nepal can be divided into five geological zones from north to south, namely: Tibetan Tethys Zone, Higher Himalaya, Lesser Himalaya, Sub-Himalaya (Siwalik) and Terai Zone. Out of these zones, Terai is plane and therefore, landslide and associated hazards are not possible but areas near to the foot hills of Siwalik which are considered as Bhabar zone also suffer from debris flow as this is the first zone that gets the materials brought from the north for deposition because of sharp decrease in elevation . Siwalik (Churia) is made up of geologically very young sedimentary rocks such as mudstones, shale, sandstones, siltstones and conglomerates form the southernmost mountain range of the Himalaya and is bounded by MFT in the south and Main Boundary Thrust (MBT) in the north  . These rocks are soft, unconsolidated and easily disintegrable. So the geology in this range is very weak and fragile. Thus the Siwalik range is becoming more vulnerable and prone to slope failure and debris flow as well as flash floods in Nepal . Most of the landslide mass converts into the debris flow during intense and continuous rainfall and thus travels for long distance towards downstream  .
The study areas Chatara-Barahakshetra, which lies in Siwalik zone is one of the victim of landslides problems. Barahakshetra is one of the famous religious places of Hindu people. On average, 300 people visit daily in this area. The landslides block the road and transportation problem appears in every rainy season. Now, the Barahakshetra-Dhankuta highway is under construction along this area. No landslide related research has been done in this area. In order to minimize the impacts of landslides on human life and property, understanding the spatial distributions of landslides is necessary as it provides an opportunity to recognize the areas that may be susceptible to landslides in the future or not     because landslides tend to occur in areas where they have occurred in the past  .
Beside this, landslide inventory and susceptibility mapping have become the most useful tool in landslide hazard management which are based on Google earth images, GIS and remote sensing data of particular study site. Landslide inventories are the basis for assessing landslide susceptibility, hazard and risk  . They are very much essential for the susceptibility model that can predict landslide on the basis of past information whereas landslide susceptibility mapping shows degree of susceptibility of area to landslide occurrences. These maps can be generated based on the spatial prediction of landslides on the assumption that future landslides will occur under same conditions as in the past . In recent years, many landslide susceptibility maps have been generated in many regions all over the world using Geographic Information System (GIS). Presently, statistical approach is the most popular for landslide susceptibility assessment. Many methods have been applied using this approach in the past such as Frequency Ratio, Weights of Evidence, and Logistic Regression. Out of these methods, Frequency Ratio (FR) method is used widely for landslide susceptibility assessment with good performance . The main objective of the current study is to create landslide susceptibility map of Chatara-Barahakshetra area, a part of Siwalik zone of Nepal using the FR model for landslide hazard management. Most of the past study has been done in Metamorphic and igneous terrain globally. This study has been done in the youngest Mountain of the world called Siwalik hill and the terrain is quaternary sedimentary terrain. This study will be a database and reference for future research related to landslide susceptibility globally.
2. Materials and Methods
2.1. Study Area
Barahakshetra is a Hindu pilgrimage site in Barahakshetra Municipality of Sunsari District of province No.1 of Nepal. The study area falls in between the Sunsari and Udayapur district along some portion of Dhankuta district and lies in the eastern part of Nepal. The area is bounded between the latitudes N26˚58'0" and N26˚47'0" and longitudes E87˚18'0" and E87˚1'30" respectively occupying a total area of 145 sq. km (Figure 1).
The study site has the tropical climate where the rainy season starts from June and ends in September. The average temperature in the study area ranges from 19˚C to 31˚C. The nearby rainfall station of the study area is Chatara index no. 1316 according to which the area has annual rainfall of 1457.41 mm per year and the average rainfall is 2140.30 mm per annum . The Chatara-Barahakshetra falls within the watershed of the Koshi River basin. The Koshi River is the largest
Figure 1. Location map of Chatara-Barahakshetra.
trans-Himalayan river passing through Nepal. Sapta Koshi is the major river basin that drains across the study area. Agriculture including animal husbandry is an important economic activity in the study area.
2.2. Geological Setting
Geologically and tectonically, the Nepal Himalaya is divided into five major tectonic zones, namely, Terai (Indo-Gangetic plain), Sub Himalaya (Siwalik), Lesser Himalaya, Higher Himalaya, and Tibetan-Tethys Himalaya  . The study area mainly lies in Siwalik zone and some portion in lower part of Lesser Himalaya which is made up of geologically very young sedimentary rocks such as mudstones, shale, sandstones, siltstones and conglomerates. These rocks are soft, unconsolidated and easily disintegrable.
Tectonically the Siwalik of Chatara-Barahakshetra is bounded by the Himalayan Frontal Thrust (HFT) and the Main Boundary Thrust (MBT) in south and north respectively. The topography changes from gentle to steep, pressures ridge can also be observed. A major syncline can be observed near the Bramhananda Khola which is east-west trending and shows extension eastward. Numerous outcrop scale folds can be observed in the study area.
Geologically, study area comprises rocks of midland group and the recent deposits. 3-fold classifications show that Siwalik can be classified as Lower Siwalik, Middle Siwalik and Upper Siwalik which follow with Takure formation, Syangja formation and Seti formation. According to , the Lower Siwalik consists of irregularly laminated beds of fine-grained greenish sandstone and siltstone with mudstone. The Middle Siwalik is represented by fine to very coarse-grained sandstone, pebbly sandstone interbedded with mudstone and siltstone . The Upper Siwalik is comprised of conglomerate and boulder beds and subordinately sand and silt beds. The Mid-land Group consists of the Gondwana sub-group, Lakharpata subgroup and the Pokhara subgroup belongs to the Siwaliks rocks . Similarly, the recent deposit consists of alluvium, boulder, gravels, sands and clays.
2.3. Methods of Data Collection
To fulfill main objective i.e. preparation of landslide susceptibility map; spatial data collection was done at first step and different factors were extracted. After extraction of different factors, the relationship with the landslides was calculated and validated. Both primary and secondary data were collected to prepare the spatial inventory map and for landslide susceptibility map based on interpretation of the aerial photographs using Google earth, remote sensing data and GIS. For the production of the thematic data layer, various thematic data layers representing landslide conditioning factors were prepared.
2.4. Primary Data
Primary data were collected through field observation with the help of instruments like GPS and checklist. Three field surveys were conducted to collect data for the study purpose. The first preliminary field visit was conducted from 18th January, 2018 to 26th January, 2018 where the study area was briefly studied and idea regarding the study area was gathered to ease the second phase field visit. The second phase field visit was conducted from 16th march, 2018 to 22nd march, 2018 which was focused on primary data collection for landslide inventory.
2.5. Secondary Data
For secondary data, various data were collected from different sources. Topographic map of the study area was collected from Department of survey (DoS), Government of Nepal. The topo-sheeet number 268701D (CHATARA) with a map scale of 1:25,000, was used to study the topographical features and related information of study area. For the study of the geological condition of the research area, Geological Map was collected from Department of Mines and Geology (DMG), Government of Nepal. The geological map with a scale of 1:50,000 were used to study the geology of the study area. Google Earth Imagery was used first to mark landslide polygon in desk before verification of landslide. It was also used to draw the verified landslide polygon and for further analysis purpose. Historical time-based images within Google Earth were used for the analysis purpose. The different factor maps were prepared by using 12.5 m DEM downloaded from Alaska satellite facility.
Meteorological data collection included collection of 30 years monthly rainfall data from 1988 to 2018 from the nearby five rainfall station were selected The rainfall data was later analyzed through interpolation Inverse Distance Weighted (IDW) method within Arc GIS to determine rainfall within the study area.
2.6. Landslide Inventory Mapping
Landslide inventory is prepared based on the database prepared from our previous study . In addition to that, visual interpretation and digitization of landslides over satellite images available in Google Earth captured from 2004 to 2018 (October, 2018) was used to update the landslide inventory remotely. Landslide was marked in polygon shape and then imported into Arc GIS software for future analysis. For the detailed study, polygon shape landslide was converted into raster data type and projection was WGS 1984 UTM zone 45 N. A total of 382 landslides are identified and mapped by evaluating aerial photos in 1:25,000 scale with well supported by field survey and these landslides polygons were mapped for the inventory.
These landslide locations were validated from historical landslide reports, and field data. The landslide events were randomly divided into 267 (70%) for training and 115 (30%) for validating the landslide susceptibility map. The study has experienced many types of landslides, for example, debris flow, rock fall, debris slides etc.
2.7. Landslide Susceptibility Mapping
For Landslide susceptibility, Frequency Ratio Model (FR) was used as stated by . Landslide susceptibility and analysis involves a number of steps and methodologies. It considers those intrinsic factors which are responsible for making the area prone to landslide. For landslide susceptibility analysis, it requires base map and factor map. Landslide Inventory map was used as a Base map. For this research, the intrinsic factors such as slope, aspect, elevation, relative relief/hill shade, geology, distance from road and streams, curvature and land cover were considered along with extrinsic factor like rainfall. The factors were totally selected on the basis of their effectiveness and availability. The intrinsic factors are chosen after the long and detail field visit. Some factors like hill shade, distance from road is taken after reviewing the global literature where as rest of the factors were chosen after analyzing their roles to occurs landslide. Total of ten factors were used to generate the landslide susceptibility map. According to , the assessment of the local landslide areas should be convenient and relevant and the factors should be illustrative and effectively available. These all factors were used to prepare different factor maps. These factor maps and base maps were used to generate landslide susceptibility map for the study area. The rainfall, being one of the main triggering factor for mass movement, it is considered within the landslide susceptibility analysis of the study area. Every factors map like slope, slope aspect, Elevation, hill shade map and curvature was prepared by using the Digital Elevation Model (DEM) (Figures 2(a)-(e)). Similarly, land cover map was prepared from Landsat image (Figure 2(i)). Geological map was prepared by using the Geo referenced map of Department of Mines and Geology (Figure 2(f)). Drainage to distance map and distance to road map was prepared from the river map obtained from Department of Survey and Open Street Map (OSM)
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)
Figure 2. Landslide affecting factor maps. (a) slope angle map; (b) slope aspect map; (c) elevation map; (d) Hill shade; (e) Curvature map; (f) Geology map; (g) Distance to stream map; (h) Distance to road map; (i) Land use map; and (j) rainfall map.
road map (Figure 2(h), Figure 2(g)) respectively. For Rainfall map, IDW interpolation method was done (Figure 2(j)). Using these factor maps and base map, the landslide susceptibility map was prepared. In this study, the frequency ratio model was applied in order to evaluate quantitatively the landslide susceptibility of the area based on the observed spatial relationship between landslide locations and each predisposing factor. Before the analysis was applied, the original landslide map was divided randomly into two sub-set groups, one with 70% (267 cases) of the total landslides was used as an estimation group to produce the susceptibility map and the second with 30% (115 cases) was used for validation of the results. To evaluate the contribution of each factor towards landslide susceptibility, the landslide estimation group was overlaid with thematic data layers separately then the frequency ratio of each factor’s class was calculated in three steps. Firstly, area ratio for landslide occurrence and non-occurrence in each factor’s class was calculated, secondly, calculating the area ratio of each factor’s class to the total area of the factor. The calculation steps for an FR for a class of the landslide-influencing factors are below  (Equation (1))
where, A is the number of pixel with landslide for each factor, B is the number of total landslide in study area, C is the number of pixel in the class area of the factor, D is the number of total pixel in the study area and FR is the frequency ratio of a class for the factor.
The obtained ratio values using FR were assigned as weight values to the classes of each factor map to produce weighted factor thematic maps, which were overlaid and numerically added using the raster calculator to produce the Landslide Susceptibility Index (LSI) map  using Equation (2):
where, Fr = the frequency ratio of each factor’s type or range).
The calculated values of FR for each pixel in the LSI indicate the relative susceptibility to landslide occurrence. The higher pixels value of LSI are the higher susceptible to landslides and the lower pixels value are the lower susceptible.
Many methods can be employed for classification of landslide susceptibility indexes such as the equal interval, the natural break and the standard deviation . Out of these, the natural break method is the most widely used  thus it has been selected for classifying the landslide susceptibility indexes in this present study. Using this method, landslide susceptibility indexes were classified into 5 susceptible classes as: 1) Very low; 2) Low; 3) Moderate; 4) High; 5) Very high.
Map validation is the final step in landslide susceptibility analysis to validate the predicted results. There are a number of methods to validate LSI map such as Receiver Operating Characteristic (ROC)    and Area under Curve methodology using Success Rate Curve or Prediction Rate Curve. In the present study, the landslide susceptibility maps were verified using the success rate and prediction rate curves  -  to carry out a validation of the prediction results.
Success rate curve was obtained by plotting cumulative percentage area of training landslide data against cumulative percentage of susceptibility map area. Whereas, prediction rate curve was obtained by plotting cumulative percentage area of landslide not considered in the model against cumulative percentage of susceptibility map area. To obtain the success rate curve for each LSI map, the calculated index values of all pixels in maps were sorted in descending order. Then the ordered pixel values were categorized into 100 classes with 1% cumulative intervals. These classified maps were crossed with landslide inventory map. Then the success rate curve was made from cross table value. Similar procedure was also done for prediction rate curve. Then the accuracy assessment was carried out for the frequency ratio model using the area under curve (AUC) method   . The area under curve was calculated by Trapezoid rule.
3.1. Landslide Inventory Map
Landslide inventory was done by using the high resolution Google Earth Imageries. Total of 382 landslides were recognized. In the field survey, 115 landslides were verified. Out of 115 landslides, 8 landslides were found in Saptakoshi river sub drainage consisting Setikhola and Mahakalkhola, Bagundrekhola and Kothukhola of udayapur district, 27 were from Kosapakhola of Udayapur district, 20 were Patnalikhola, Baghkhola and Bhaguwakhola of Sunsari district, 22 were from Kokahakhola of Sunsaridistrict and remaining 38 from Bramanandkhola, Bahunekhola, Dhansarkhola, Dudhpanikhola, RaniKhola, Gunjamankhola, Chinnamastakhola, Gaurikhola and Karamkhola of Sunsaridistrict (Figure 3).
The study area has total area of 145.3 sq. km in which landslides covers 1.53 sq. km area (1.30%) of whole study area. In the area of 145.3 sq. km, the landslide density was found to be 14.73 per km sq. The landslide was distributed unevenly in the study area.
3.2. Landslide Susceptibility Mapping
The landslide susceptibility of the study area was identified using frequency ratio (FR) model. The calculation for each factor map was conducted based on FR model after reclassifying each factor map and overlaying them with the train landslides. The final landslide susceptibility map is obtained by the FR model (Figure 4 & Figure 5). Frequency ratios for the classes of each factor, which are shown in Table 1, were calculated by dividing the landslide occurrences ratio by the area ratio.
Table 1 shows that frequency ratio increased with increasing of slope gradient. At the slope gradient between 15˚ - 60˚ the frequency ratio were > 1 in the south west, south, southeast, west and North West aspect with landslide occurrence whereas negative correlation between other aspects. Correlations between elevation and landslide occurrence indicated that landslides are more likely to occur at the elevation range 240 - 840 m (Table 1) whereas as the elevation increases
Figure 3. Landslide inventory with drainage and hill shade map of study area.
Figure 4. Final weighted map of study area.
Figure 5. Landslide susceptibility map of Chatara-Barahakshetra area.
from 840 m, less FR values and landslide occurrence are observed. In the study area, relative relief was found to be varying between 0 to 254 m. FR value showed the higher correlation between landslide occurrence and elevation from 0 to 100 m along with elevation from 200 - 254 m (Table 1). In the case of curvature, FR is highest for concave slope (1.53), followed by convex slope (1.14), while this value is less than 1 for flat slopes (Table 1). For the lithology, it can be seen that Upper Siwalik contains the highest value of FR (2.83) indicating very high correlation with landsides, followed by Takure Formation (1.84) and Middle Siwalik (1.50) (Table 1) while other have low correlation. Analysis from FR value shows that landslide frequency decreases as the distance from roads increases. At a distance of < 400 m, the ratios were < 1 indicating a very weak relation with landslide occurrence, whereas at distances > 400 m, the ratios were > 1, indicating strong relationship and shows that it has a very less influence in landslide formation in the study area (Table 1). In the case of distance from stream, the distance in between 150 and 250 m shows the highest correlation with the landslides (Table 1). Most part of the study area is covered by forest and barren land as well as agriculture land shows positive correlation with the landslides (FR > 1). All other land use classes show negative correlation with the landslide formation (Table 1). In the case of relationship between precipitation and landslide occurrence, FR values showed positive correlation between rainfall and landslides occurrences from 1895 mm to 2059 mm.
The LSI value was calculated based on Equation (2). The landslide susceptibility map was prepared then after the combination of all these factors, and the
Table 1. Landslide conditioning factors and its Frequency Ratio values.
areas were categorized into five classes. The higher the LSI value, the greater is the susceptibility of the zone. The result showed that as very low, low, moderate, high and very high classes occupied 21.29%, 27.25%, 30.40%, 17.36%, and 3.68% respectively of the total study area (Figure 6). The performance of the FR model was evaluated using the area under curve (AUC) method In the case of success rate curve; the success rate curve constructed showed 72.55 percentage of the area lying under the curve with the training dataset (Figure 7). Similarly, the prediction rate curve prepared using the validation group of landslides showed 71.73 percentage of the area lying under the curve indicating that prediction ability of the Frequency Ratio model.
Figure 6. The percentages distribution of the susceptibility classes and landslide occurrence.
Figure 7. Area under curve shows the success rate curve of landslide susceptibility.
The landslide susceptibility index in the study area ranged between 3.56 and 40.99 which were classified into five major susceptibility classes. Very high susceptible (VHS) areas are found mostly in barren land, forest areas and the settlements areas. The study shows that very low, low, moderate, high and very high classes occupied 21.29%, 27.25%, 30.40%, 17.36%, and 3.68% respectively of the total study area. Present study on this context considered both intrinsic and extrinsic factors. Considering rainfall factor within the study can make the information more reliable to predict the landslide susceptible areas. It was observed that the slope, aspect, curvature, stream buffer, geology, land use along rainfall strongly influencing the landslide susceptibility of the area. After intense work in the field, nine intrinsic factors were chosen for the susceptibility analysis. As per the result, the landslides within the study area were contributed by the slope degree between >30˚ in southeast, south and south-west aspect.  showed increasing relationship of landslide density with slope angle up to 35˚ and is highest in 25˚ - 35˚ slope class and similar results were obtained by  and  in Mugling-Narayanghat road, Chitwan District of the Narayani Zone and Garuwa sub-basin, East Nepal respectively. Also, this has probably to do with the prevailing direction of the Monsoon storms that enter from the east and slowly move toward the west while producing a lot of precipitation on the south facing slope . Similarly, elevation between 240 to 640 m and surface structure (convex and concave)was also responsible for weathering process through the repeated dilation and thawing process due to exposed sun rays and rainfall in higher elevation which increases soil water pressures and reduces shear resistance .  found similar result in Chure Khola Catchment area of the Siwalik region, Central Nepal. Distance from drainage also has a clear influence on landslides as nearer distance from stream indicating a higher probability of landslides. The debris flow and sliding are high in Siwalik zone due to intense rainfall near to the stream  which was Similar to result obtained by  in Bhalubang-Shiwapur Area of Mid-Western Nepal and  in Tulsipur-Kapurkot road section, Nepal respectively. The weak relation of road distance contradicts the result obtained by . This may probably due to the ongoing road construction only in district highway; very less number of small landslides was found on road side. Beside this, the weak geology topography and rocks mainly quartize, phyllites, presences of sandstone, siltstone and mudstone creating high probability of landslide occurrences in the study area as the landslide number is higher in middle Siwalik, Upper Siwalik and Takure formation    found similar result in their study in Siwalik zone. The obtained landslide susceptibility maps also demonstrated that the large area is susceptible and prone to landslide.
For validation, both success rate and prediction rate curves were used by comparing the existing landslide locations with the landslide susceptibility maps. The success rate and prediction rate results were obtained by using the training dataset 70% (267 landslide locations) and validation dataset 30% (115 landslide location). Since the success rate method used the training landslide pixels that have already been used for building the landslide models, the success rate is not a suitable method for assessing the prediction capability of the models. However, the success rate method may help to determine how well the resulting landslide susceptibility maps have classified the areas of existing landslides  . The prediction rate explains how well the model and predictor variable predict the landslide. This method is already widely used as a measure of performance of a predictive rule   . In our study, the area under the curve showed the success rate of 72.55% and prediction rate of 71.73% for the frequency ratio model. In this study, both the success rate and prediction rate curve show almost similar result, and the models utilized in this study showed reasonably good accuracy in predicting the landslide susceptibility of the study area.
There are many studies where success rate curve and predictive rate curve were used to validate the result using frequency ratio model.
 applied frequency ratio, weights of evidence, and statistical index models in their comparison in landslide susceptibility mapping in Central Nepal Himalaya. The results depicted that the FR model had the best performance with a success rate of 76.8% and a predictive accuracy of 75.4%.  obtained 75% success accuracy and 70% prediction accuracy using FR model in a part of Uttarakhand Himalaya, India for landslide susceptibility assessment.
Another study carried out by  in Bhotang village development committee, Nepal using FR (Frequency Ration) and SI (Statistical Index) methods and AUC value for FR model was found to be 0.713 for validation, 0.725 for training and 0.722 for whole dataset. Similarly,  obtained 83.31% and 78.58% success accuracy and prediction accuracy for landslides susceptibility mapping along Bhalubang-Shiwapur Area of Mid-Western Nepal Using Frequency Ratio.
Our result is accurate and mostly similar with the previous study and FR model is found as suitable model for landslide susceptibility mapping in the Siwalik Hills.
In this study, remote sensing and GIS are effectively used to develop a landslide susceptibility map. Total of 382 landslide polygons were mapped from Google earth and by field verification. Most of the landslides were observed to be concentrated along the Main Boundary Thrust and top of the hill. Slope greater than 30˚, southeast, south and southwest aspect, elevation from 240 m to 640 m, Middle Siwalik, Upper Siwalik along with Takure Formation, barren land, intense rainfall and concave slopes are more susceptible for the landslide occurrence. The derived relationship between the landslide inventory and causative factors is used in the Frequency ratio technique to develop a landslide susceptibility map. The susceptibility map was classified into very low, low, moderate and high and very high susceptibility class. The result shows that 17.3% and 3.68% of the study area are identified as high and very high susceptible zone for landslide.
The validation results showed that the success rate curve with 72.55 percentage of the area lying under the curve and the prediction rate curve with 71.73 percentage of the area lying under the curve indicating that prediction ability of the Frequency Ratio model. These landslide susceptibility maps can be used as a planning tool by prioritizing areas for controlling the landslide effects. More than 71% success rates indicate that frequency ratio model is suitable model for the landslide susceptibility in the study area. During our study, we could not manage LIDAR and other latest technique. We just applied statistical method to landslide susceptibility study. We strongly recommend to the future researchers to prepare LIDAR based DEM of resolution less than 5 × 5.
We want to acknowledge University Grant Commission Nepal for providing research grant for this study having award number MRS-75/76-S&T-35. Similarly, our deepest gratitude goes to the local people and local government of Chatara-Barahakshetra area of Nepal.