A replication of landslide hazard mapping at catchment scale

Abstract

Landslide hazard assessment is a primary tool to understand the basic characteristics of slopes that are prone to landslides, especially during extreme rainfall. In this study, weights-of-evidence modelling a bivariate statistical method, and a logistic regression model, a multivariate statistical method, were used for landslide hazard mapping in two catchments of the Siwaliks in the Nepal Himalaya. Two typical catchments, Charnath and Jalad of the Siwaliks in easternNepal, were selected for the landslide hazard mapping. Both modelling approaches were applied by considering 10 intrinsic factors and one extrinsic factor. Mainly digital elevation model-based causative factors and field data were used to prepare data layers of landslide causative factors. In many approaches formodelling of landslide hazard in GIS, the model validation process is always dependent, and landslide data, which are used to calculate a landslide hazard index (LHI), are applied for verification. However, in this study, the LHI was calculated in one catchment (Jalad) and the same index for a different class of causative factors was applied for another catchment (Charnath), and the LHI wasverified. The verification results were very promising, with an independent prediction rate of about 75%. This validates weights-of-evidence and logistic regression models for landslide hazard assessment in the Siwaliks Range of Nepal.

1. Introduction

Landslides are amongst the most disastrous natural hazards in tropical, subtropical and temperate climatic zones of the world. Potential sites that are particularly prone to landslides should therefore be identified in advance so as to reduce damage frompotential landslide disasters. Landslide hazard or susceptibility assessment is aprimary tool used to understand the basic characteristics of the slopes that are prone to failure, especially during extreme rainfall.

In general, landslide hazard can be depicted as the physical potential of the process to produce damage. The damage depends on the particular impact characteristics of the landslide, as well as the magnitude and frequency with which it occurs or is encountered. Fundamentally, landslide hazard can be characterized by statements of‘what’, ‘where’, ‘when’, ‘how strong’ and ‘how often’ (Crozier and Glade 2005, Glade and Crozier 2005). Unfortunately, the ability to predict landslide hazard with precision is always limited, and consequently landslide hazard predictions are generally described in terms of likelihoods and probabilities. In many landslide hazard mapping techniques, probabilities of failure are categorized and landslide hazard zonation maps are prepared. Landslide hazard zonation is defined as the mapping of areas with an equal probability of occurrence of landslides within a specified period of time (Varnes 1984). Landslide hazard zonation demonstrates the spatial and temporal occurrence of future landslides in terms of probability.

The majority of the published literature on landslide hazard mapping addresses landslide susceptibility mapping. There are numerous studies involving landslide hazard evaluation, and Guzzetti (2005) in particular has explored many cases of landslide hazard studies. Landslide hazard levels can be addressed qualitatively orquantitatively, as well as either directly or indirectly (Guzzetti et al. 1999). Qualitative methods always describe landslide hazard in descriptive terms, whereas quantitative methods use probability values to predict landslide hazard in the study area. Direct methods consist of the geomorphological mapping of landslide hazard, whereas indirect methods usually deal with identification and mapping of groups ofphysical factors which are directly or indirectly associated with the landslides. Such physical factors are called as intrinsic and extrinsic factors. In many cases, intrinsic (bedrock geology, geomorphology, soil depth, soil type, slope gradient, slope aspect, slope curvature, elevation, land use pattern, drainage pattern and so on) and extrinsic (rainfall, earthquakes and volcanoes) factors are used to determine landslide hazard in an area (Siddle et al. 1991, Dai et al. 2001, Çevik and Topal 2003, Dahal et al. 2008a). The extrinsic factors are site specific and possess temporal distribution. Moreover, they are difficult to estimate because of the lack of information about their spatial distribution. Hence, in landslide hazard assessment practice, the term ‘landslide susceptibility mapping’ is addressed without considering the extrinsic factors in determining the probability of occurrence of a landslide event (Dai et al. 2001, Dahal et al. 2008a). A region is considered to be susceptible to landslides when the terrain conditions at that site are similar to those in the region where a slide has occurred (van Westen 2000).

The integrated analysis of all intrinsic variables in relation to the spatial distribution of landslides has gained enormous success by the introduction of geographic information systems (GIS), the ideal tool for the analysis of parameters with a high degree of spatial variability. For landslide hazard assessment, the assumption is made that conditions which led in the past to landslides will also result in potential unstable conditions in the present and future. Thus, a landslide inventory map, differentiating the type, activity and dimensions of mappable landslides, provides primary data for landslide hazard or susceptibility zonation. The inventory map also needs to cover information about time-span-based landslide distribution as far as possible. When mapping intrinsic factors, emphasis should be given to the most relevant terrain parameters related to the occurrence of landslides. Generally, it is true that the selection of intrinsic parameters takes the nature of the study area and the data availability into account. However, in a GIS-based technique it is also necessary to be sure that any selected factor is functional (has a certain degree of affinity with previous occurrences of landslides), complete (is reasonably represented over all of the study area), non-uniform (remarkable spatial variation), measurable (can be expressed by nominal, ordinal, interval, ratio scales) and non-redundant, i.e. the outcome of selected factors should not account for double effects in the final result (van Westen 2000, Yalcin 2008).

Landslide hazard and susceptibility mapping are assessed through heuristic, statistical and deterministic approaches (Okimura and Kawatani 1987, Yin and Yan 1988, Gökceoglu and Aksoy 1996, Soeters and van Westen 1996, van Westen and Terlien 1996, Pachauri et al. 1998, van Westen 2000, Dai et al. 2001, van Westen et al. 2003, Süzen and Doyuran 2004, Zêzere et al. 2004, Saha et al. 2005, Dahal et al. 2008b, Prahdan et al. 2009, Chauhan et al. 2010, Poudyal et al. 2010, Pradhan and Lee 2010, Regmi et al. 2010a, Srivastava et al. 2010, Odzemir 2011).

The heuristic approach is a direct or semi-direct mapping methodology in which a direct relationship is established between the occurrence of landslides and the causative terrain parameters during the landslide inventory. Therefore, in this approach, the opinions of the experts are very important to estimate the potential of landslides based on the intrinsic factors of landslides included in the analysis. Similarly, assigning weight values and ratings to the variables is very subjective and the results are not reproducible. Deterministic approaches, however, are based on slope stability analyses, and are only applicable when the ground conditions are relatively homogeneous across the study area and the landslide types are known. The infinite slope stability model has been widely used to assess landslide susceptibility in deterministic approaches (Wu and Sidle 1995, Gökceoglu and Aksoy 1996, Terlien 1996, Dahal et al. 2008c), and such a stability model requires a high degree of simplification of the intrinsic variables. The statistical approach, on the other hand, is an indirect susceptibility mapping methodology which involves statistical determination of the combinations of variables that have led to landslide occurrence in the past. All possible intrinsic variables are entered into a GIS and crossed for their analysis with a landslide inventory map. Both bivariate and multivariate statistical approaches have been used for mapping landslide hazards (Siddle et al. 1991, Atkinson and Massari 1998, Guzzetti et al. 1999).

In most of cases of the statistical approach to landslide hazard mapping, either bivariate or multivariate, the replication of hazard indexes of causative parameters is always a problem. It is very hard to find the same intrinsic parameters in two different regions. To avoid these limitations, in this study, two catchments are evaluated for replication of landslide hazard mapping. In this study, the landslide hazard is evaluated in the GIS platform using weight-of-evidence modelling and logistic regression modelling. Two catchments (Jalad and Charnath) of the Sub-Himalayan Zone (Siwalik Group) of the Nepal Himalaya are selected as the ideal sites for the study, because both catchments have serious rainfall-induced landslides and slope degradation problems. They have geological and geomorphological settings typical of the region.

2. The study area

The study area is located in the south-eastern hills of the Nepal Himalaya and belongs to the Sub-Himalayan Zone or Siwaliks (figure 1) of the Himalaya.

Figure 1. Location of the study area and general geological map of Nepal (modified after Dahal and Hasegawa, 2008). Available in colour online.

Siwaliks, also called Churia in Nepal, is the youngest low-altitude mountain in the Himalayan chain. Siwaliks consists of relatively unconsolidated soft rock, namely mudstone, sandstone and loose conglomerate, and it comprises many vertical slopes and rugged terrain. The average slope and relief of terrain are high in the western regions compared with other regions of Nepal. Many parts of the Siwaliks are inhabited by a relatively large human population, with agriculture and livestock as the mainstays of livelihood. River valleys, comprising narrow tracts of flatland, are extensively used for cultivation and settlement. The population and livelihood activities are highly vulnerable to floods, debris torrents, debris slides and debris flows from hill slopes. Obviously, the pressure of the population and fast-growing development and economic activities in forest and cultivation is high in the entire Siwaliks. Land use and land cover distribution patterns show the predominance of forest in the hill slopes and cultivated land in the plain area. Considerable patches of cultivated areas also lie on areas of steep topography, flood and slope failure-prone areas. Over the last four decades there has been a tremendous change in land use and land cover, but this change has been less in the last decade (Ghimire et al. 2008) because of the conversion of forest areas into cultivated lands. The change was more pronounced in slopes of Siwaliks and inner river valleys. A significant expansion ofcultivated land in steep slopes and in flood-prone areas clearly indicates environmental degradation and livelihood vulnerability in the Siwaliks (Ghimire et al. 2008, Ghimire 2011).

Two typical catchments, Charnath and Jalad (figure 2), were selected for the study. Both catchments are located in the northern part of Dhanusa district of Nepal (figure 2). The Jalad and Charnath are two major river systems originating from theSiwaliks and they are highly disaster-prone rivers. These catchments have very fragile lithology and steep topography, predisposing the landscape to significant problems of geology- and geomorphology-related slope failures. The Charnath River is about 8 km east from the outlets of the Jalad River in Terai. The Jalad catchment in Churia is orientated in a north–south direction, whereas Charnath is in northwest–southeast orientation.

Figure 2. Location map of Jalad and Charnath catchments in Siwaliks of Dhanusa district, Nepal (map modified after UNDP Nepal 2010). Available in colour online.

2.1. Geological and geomorphological settings of selected catchments

The geology of the Nepal Himalaya has been divided into five major east–west trending tectonic units (figure 1). From south to north, these units are Terai Zone, Sub-Himalayan Zone (Siwaliks), Lesser Himalayan Zone, Higher Himalayan Zone and Tibetan-Tethys Zone. As mentioned earlier, the study area belongs to the Sub-Himalayan Zone or Siwaliks. The Siwaliks 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. According to three fold classifications of the Siwaliks, the area can be divided into three geological units (Burbank et al. 1996, Ulak and Nakayama 1999): Lower Siwaliks, Middle Siwaliks and Upper Siwaliks, as per the occurrence of rock types. Lower Siwaliks usually has thin to thickly bedded green to grey-coloured mudstone and siltstone, whereas Middle Siwaliks has thick to very thickly bedded sandstone. The Upper Siwaliks contains thick beds of loose and fragile conglomerates.

Factors such as excessive rainfall and human intervention are the main triggering agents of landslides in Siwaliks. Factors such as groundwater condition, river undercuttings and deforestation on slopes also facilitate slope failures. Similarly, Lower Siwaliks and Middle Siwaliks have problems stemming from alternating bedsof mudstones and sandstone. In such alternating bands, mudstone can flow when saturated with water, which results overhanging sandstone beds. Such overhang-jointed sandstone beds are easily disintegrated into blocks. Similar to the rest of Nepal, the rainfall within the Siwaliks is normally in the range of 2000–2500 mm per year (Ghimire et al. 2008). As a result, geological conditions and the climate render the Siwaliks highly susceptible to landslides.

In the both Jalad and Charnath catchments, the Siwaliks can be divided into twounits: Middle Siwaliks and Upper Siwaliks, on the basis of lithological composition (figure 3). Middle Siwaliks consists of alternating thickly bedded light yellow to light grey-coloured mudstone and light grey-coloured sandstone with subordinate pebbles. The proportion of sandstone and mudstone is equal on the southern part, while the proportion of sandstone increases towards upstream of both catchments.

Figure 3. Thickly bedded alternating bands of coarse and fine-grained sandstones in the Middle Siwaliks (left) and loose conglomerate in the Upper Siwaliks as seen in the upstream of the Jalad River (right). Available in colour online.

Land degradation is a common problem in the fragile and steep hill slopes of both catchments, where pressure from human intervention and grazing animal are extensively high. The transport routes are all used extensively for smuggling wood from the forest. As a result, the areas in the vicinity of the routes are exceedingly deforested. Debris slides and debris flows are the major land degradation processes in both catchments. Charnath catchment is relatively more densely populated than Jalad catchment. In Nepal, the landless population is high in the northern Lesser Himalayan Zone, and migration to the southern Siwalik Range isa common problem. Both catchments have the same problem of migration, andimproper agricultural practices make both catchments highly vulnerable to landslides.

The stream channels are narrowed and meandered, showing cliffs formed by bank cuttings on their concave banks. The river width is relatively higher upstream than downstream. The sediment thickness of the channel lag on the riverbed is not more than 10 m, while in some places the river flows over the bedrock. The main geological difference between Jalad and Charnath is on the coverage of Upper Siwaliks (loose conglomerate). Nearly half of the Jalad catchment has Upper Siwaliks, whereas Charnath has Upper Siwaliks only in the upper reaches. As a result, there is less gravel in the Charnath river bed in comparison with the Jalad river bed.

3. Goals and objectives of the study

The main goal of this study is to prepare landslide hazard maps of the selected catchments by replicating the landslide hazard index (LHI) of one catchment into another catchment. The specific objectives for this study are as follows:

•	To conduct field work to gather the relevant data and information regarding shallow landslides;
•	To identify similar causative intrinsic and extrinsic factors of landslides in both catchments;
•	To analyse the geological and geomorphological data and information for landslide hazard mapping;
•	To prepare weights-of-evidence and logistic regression models and apply them in the Jalad catchment for assessment of LHI;
•	To apply the same LHI of each causative factor obtained from the Jalad catchment to similar causative factors of the Charnath catchment and validate the models;
•	To evaluate the replication procedure of landslide hazard mapping at catchment scale;
•	To prepare a landslide hazard zonation map (scale 1:10,000) of both catchments based on the validation and comparative results of the weights-of-evidence and logistic regression models.

4. Methodology and materials

Although there are various methods for statistical approaches in the study of landslide hazard, in this study weights-of-evidence modelling and logistic regression modelling were used for the hazard analysis in the two selected catchments. First, Jalad catchment was analysed with weights-of-evidence and logistic regression models, and then the results were applied in the Charnath catchment. Details of methods and data used in the analysis are described in the following paragraphs.

4.1. Weights-of evidence model

In this study, both bivariate and multivariate approaches were selected to assess LHIs. In the bivariate statistical approach, the weights-of-evidence model was used to assess the LHI. The method uses the Bayesian probability model, and was originally developed for mineral potential assessment (Bonham-Carter et al. 1988, 1989, Agterberg 1992, Agterberg et al. 1993, Bonham-Carter 2002). Several authors have applied the method to mineral potential mapping using GIS (Emmanuel et al. 2000, Harris et al. 2000, Tangestani and Moore 2001, Carranza and Hale 2002). Cheng (2004) used the method to predict the location of flowing wells, and Daneshfar and Benn (2002) used weights-of-evidence to analyse spatial associations between faults and seismicity, whereas Zahiri et al. (2006) applied weights-of-evidence model for mapping cliff instabilities associated with mine subsidence. The method has also been applied to landslide susceptibility and hazard mapping (Lee et al. 2002, Van Westen et al. 2003, Lee and Choi 2004, Lee and Sambath 2006, Neuhäuser and Terhorst 2007, Dahal et al. 2008a, 2008b, Sharma and Kumar 2008, Regmi et al. 2010a).

A detailed description of the mathematical formulation of the method is available in Bonham-Carter (2002), and mathematical relationships for landslide hazard mapping are available in Dahal et al. (2008a,b). For landslide hazard modelling, themethod calculates the weight for each landslide causative factor based on the presence or absence of the landslides within the area. Therefore, historical landslide data are essential for weighting factors. The modelling procedure also relies on thefundamental assumption that future landslides will occur under conditions similar to those contributing to past landslides. It also assumes that causative factors for the mapped landslides remain constant over time. The method calculates the weight for each landslide predictive factor (F) based on the presence or absence of the landslides (L) within the area, as indicated in Bonham-Carter (2002) and modified by Dahal et al. (2008a) for landslide hazard assessment as follows:

A positive weight indicates that the causative factor is present at the landslide location and there is a positive correlation between presence of the causative factor and landslides. A negative weight indicates an absence of thecausative factor and shows the level of negative correlation. The difference between the two weights is known as the weight contrast, , and the magnitude of contrast reflects the overall spatial association between the causative factor and landslides. In weights-of-evidence model, the combination of causative factors assumes that the factors are conditionally independent of one another with respect to the landslides (Bonham-Carter 2002, Lee and Choi 2004). In this research also, the assumption is made that all landslides in the selected catchment area occur under the same combination of parameters, and that all sets of parameters are conditionally independent.

4.2. Logistic regression model

In this study, a logistic regression model was used to assess LHIs. Logistic regression can be used to determine the relation of landslide occurrence and the related factors (Guzzetti et al. 1999, Dai and Lee 2002, Ohlmacher and Davis 2003, Ayalew and Yamagishi 2005, Lee 2005, Zhu and Huang 2006, Akgün and Bulut 2007, Chen and Wang 2007, Chauhan et al. 2010). It is useful when the outcome variable or dependent variable is binary or dichotomous. The dependent variable for this analysis is the absence or presence of a landslide.

Considering n independent variables x ₁, x ₂, x ₃, … , x_n , affecting landslide occurrences, the vector X = (x ₁, x ₂, x ₃, … , x_n ) has been defined. In logistic regression analysis, the logit y is assumed as a linear combination of independent variables and is given as follows:

where b ₀ is the constant of the equation, and b ₁, b ₂, …, b_n are the coefficients of independent variables, x ₁, x ₂, x ₃, … , x_n . For landslide hazard assessment, the dependent variable is a binary variable, with values of 1 or 0, representing the presence or absence of landslides. Quantitatively, the relationship between the occurrence and its dependency on several variables can be expressed as follows (Hosmer and Lemeshow 2000):

where P is the estimated conditional probability of landslide occurrence and e is the constant 2.718. From equation (3) and equation (4), a relationship can be obtained in which the natural logarithm of the odds, log(P/(1 − P)), is linearly related with the independent variables as follows:

The goodness of fit of the model was tested with the Wald statics and Hosmer–Lemeshow test (Hosmer and Lemeshow 2000). By examining the sign of a dependent variable's coefficient estimate, the effect of that variable on the probability of landslide occurrence can be determined. These values were utilized to decide the final LHI of the Jalad catchment of the study area.

4.3. Data preparation

The main steps for landslide hazard mapping are data collection and the construction of a spatial database, from which relevant causative factors are extracted. This is followed by assessment of the landslide hazard using the relationship between landslides and causative factors, and the subsequent validation of results. A key feature of this method is that the probability of landslide occurrence will be comparable with observed landslides.

For the modelling in weights-of-evidence and logistic regression, a number of thematic data of causative factors have been identified, including 1) slope; 2) slope aspect; 3) distance to drainage; 4) wetness index; 5) sediment transport index; 6)relief; 7) distance to transportation route; 8) curvature; 9) geology; 10) land cover; and 11) mean annual rainfall. Topographic maps and aerial photographs provided by the Department of Survey, Government of Nepal, were considered as basic data sources for generating these layers. Field surveys were carried out for data collection and to prepare data layers of various factors. Landslide distribution maps of both catchments were also prepared in the field. These data sources were used to generate various thematic layers using the GIS software ILWIS 3.7. Brief descriptions of the preparation procedure of each data layer are provided here.

4.3.1. Landslide characteristics and inventory maps

A landslide inventory map is the simplest output of direct landslide mapping. It shows the location of discernible failures. Two landslide inventory maps for Jalad and Charnath catchments were prepared. However, a detailed landslide map was prepared for the Jalad catchment and only recent locations of landslides were mapped in the Charnath catchment. Only scars (main failure portion) were used to delineate the landslides. With the use of field data and satellite imaging, 426 landslide scars were delineated in the Jalad catchment. Satellite imaging PRISM (Panchromatic Remote-sensing Instrument for Stereo Mapping) has 2.5 m resolution and is very useful to identify landslide images even from visual interpretation. The landslide scars mapped in the field were again revised with PRISM data and final landslide inventory maps were prepared for both catchments. All landslides were used to assess LHI from both logistic regression and weights-of-evidence modelling. Finally, to check the predictive power in all possible aspects of both models, the LHI of each class of the causative factors estimated from the Jalad catchment was used for the same class of causative factors in the Charnath catchment. Moreover, the prediction rates of both models were evaluated from thelandslide inventory of the Charnath catchment. The inventory maps of both catchments are given in figure 4 and figure 5.

Figure 4. PRISM image of the Jalad catchment in the Siwaliks with 426 scars of failures. Available in colour online.

Figure 5. DEM of Charnath catchment in the Siwaliks with 101 scars of recent landslides. Available in colour online.

4.3.2. Geological map

Geology plays an important role in landslide susceptibility and hazard studies because different geological units have different susceptibilities to active geomorphological processes of the Himalaya (Pradhan et al. 2006). Sandstone, mudstone, siltstone and gravel conglomerate are the main rocks of the selected catchments, and they belong to the Lower and Upper Siwalik groups. For the geological map preparation, previous studies (DMG 2004, Ghimire et al. 2008) were consulted and geological boundaries of formations were checked during the field visit and modified as per the field observations.

4.3.3. Geomorphology related causative factors

A Digital Elevation Model (DEM) representing the terrain is key to generating various geomorphological parameters which influence the landslide activity in an area. Hence, a DEM was prepared with digital contour data (1:25,000 scale) obtained from the Department of Survey, Government of Nepal. The DEM of study area was prepared in 10 m × 10 m pixel size. From this DEM, geomorphological thematic data layers such as slope, aspect, relative relief, distance to drainage, sediment transport index and wetness index were prepared.

4.3.3.1. Slope and aspect

Slope data layer, an important parameter in slope stability considerations, comprises eight classes. This classification was based on the slope angles of more than 70 landslides measured in the field. Field measurement signified that most of the landslides occur at slope angle between 15° and 30°. Thus, in total seven classes, <5°, 5–10°, 10–15°, 15–20°, 20–25°, 25–30° and >30° were used to prepare the slope map. Aspect is referred to as the direction of maximum slope of the terrain surface. For the study area, it was divided into nine classes, namely, N, NE, E, SE, S, SW, W, NW and Flat. Both slope and aspect maps were prepared in ILWIS 3.7 from DEM.

4.3.3.2. Relief

Relative relief is another DEM-based derivative and is defined as the maximum height dispersion of a terrain normalized by its length or area. In this study, relative relief is computed as the difference between the maximum and minimum elevations within the given class of elevation. The relief was sliced into six classes at 50 m elevation difference.

4.3.3.3. Distance to drainage

In the study area, the appearance of landslide was frequent along the stream. Thus, the location of the landslide relative to the stream was considered as another geomorphology-related causative factor. Subsequently, adistance-to-drainage map was generated as per the field identification that slope failure may be more frequent along stream, due to groundwater movement towards the stream and to undercutting. In order to produce the map showing distance to drainage, the drainage segment map was rasterized and the distance to the drainage was calculated in metres. The resultant map was then sliced into six classes. The six classes were 0–20 m, 20–40 m, 40–60 m, 60–80 m, 80–100 m and >100 m.

4.3.3.4. Curvature, sediment transport index and wetness index

Following rainfall events, water flows from areas of convex curvature and accumulates in areas of concave curvature. This process is known as flow accumulation and is a measure of the land area that contributes surface water to an area where water can accumulate. This parameter was considered as relevant to this study because it defines the locations of water concentration after rainfall, and those locations are likely to have a high slope failure incidence.

Three slope morphology-related parameters, namely sediment transport index, wetness index and curvature, were considered for slope morphology-related parameter study. The sediment transport process is enhancing the slope failure mechanism extensively in the Siwaliks. The sediment transport index accounts for the effect of topography on erosion. The planimetric catchment area is used to calculate this index. The wetness index sets the catchment area in relation to the slope gradient. An idea of the spatial distribution of zones of saturation can be obtained from the wetness index (Moore et al. 1991). As per the value of the wetness index and sediment transport index, as well as histogram information and equal area distribution of cumulative curves, the values were classified into different categories and analysed with the landslide inventory map.

4.3.4. Distance to transport route

One of the controlling factors for the stability of slopes is road construction activity. However, in the study area, road access has been a prime factor in forest degradation. As a result, many landslides are prominent in the degraded forest area. Every day wood smuggling is increasing (figure 6), and along the transport route the forest is also declining. Thus, a distance-to-road map was generated as per the hypothesis that slope failure may be more frequent along roads and trails, owing to excessive forest declination along the roads and trails. In order to produce the map showing distance to transport route, the road and foot-trail segment map was rasterized and the distance to thetransport routes calculated in metres. The resultant map was then categorized into six classes as <200 m, 200–400 m, 400–600 m, 600–800 m, 800–1000 m and >1000 m.

Figure 6. Photographic evidence of smuggling of tree logs and firewood from the selected catchments. Available in colour online.

4.3.5. Land cover map

More than 90% of the study area is covered by forest. To interpret land cover of the area, Normalized Difference Vegetation Index (NDVI) was considered and the NDVI map was obtained from a Landsat TM satellite image acquired on 1 November 2005. The NDVI value was calculated using the NDVI = (IR − R)/(IR + R) formula, where IR is near-infrared band image and R is red band image of Landsat TM. The NDVI value denotes areas of vegetation in an image. The presence of dense green vegetation implies high NDVI values, due to the high concentration of chlorophyll resulting in a low reflectance in the red band, as well as due to the high stacking of leaves. Sparse vegetation, on the other hand, implies low NDVI values. Likewise, water and clouds have larger visible reflectance than near-infrared reflectance. Thus, these features yield negative index values. Rock and bare soil areas have similar reflectance in the two bands and result in vegetation indices near zero.

To prepare categorical thematic layers of NDVI for Jalad and Charnath catchments, the cumulative frequency curve of NDVI values was segmented into five classes representing equal distribution to yield five NDVI classes: 0–20%, 20–40%, 40–60%, 60–80% and 80–100%.

4.3.6. Rainfall map

Rainfall is an extrinsic variable in hazard analysis, and the spatial distribution of mean annual rainfall is commonly used in this statistical hazard analysis (Dahal 2008a, Regmi et al. 2010b). Rainfall stations around the study area were used to prepare a mean annual rainfall map. The spatial distribution of rainfall was calculated through the application of the inverse distance squared method in ILWIS 3.7. The resulting map was sliced to give a raster map with four classes of <1500 mm, 1500–1550 mm, 1500–1600 mm and >1600 mm for the Jalad catchment, and <1450 mm, 1450–1500 mm, 1500–1550 mm and >1550 mm for the Charnath catchment.

All the data layers such as slope, slope aspect, distance to drainage, wetness index, sediment transport index, relief, distance to transportation route, curvature, geology, NDVI and mean annual rainfall of both catchments are given in figure 7 and figure8.

Figure 7. Various data layers used in landslide hazard analysis in the Jalad catchment. Available in colour online.

Figure 8. Various data layers used in landslide hazard analysis in the Charnath catchment. Available in colour online.

5. Analysis and results

5.1. Weights-of-evidence

To evaluate the contribution of each factor towards landslide hazard, the landslide distribution data layer of the Jalad catchment have been compared with various thematic data layers. For this purpose, equations (1) and (2) were converted to equations (6) and (7) for use in the GIS platform,

where:

Npix ₁ is the number of pixels representing the presence of both potential landslide causative factor and landslides;

Npix ₂ is the number of pixels representing the presence of landslides and absence of potential landslide causative factor;

Npix ₃ is the number of pixels representing presence of potential landslide causative factor and absence of landslides;

Npix ₄ is number of pixels representing absence of both potential landslide causative factor and landslides.

All thematic maps of causative factors were stored in raster format with a pixel size of 10 m × 10 m, and they were combined with the landslide inventory map for the calculation of the positive and negative weights. The calculation procedure was written in the form of a script file in ILWIS 3.7, consisting of a series of GIS commands to support equations (6) and (7). Consequently, the weight of each class of causative factors was calculated (table 1).

Table 1. Weights-of-index values of various classes of causative parameters estimated from the Jalad watershed and replicated in the Charnath watershed.

Classes	LHI	Classes	LHI	Classes	LHI
Aspect		NDVI		Sediment transport index
Flat	−0.4421	<0.476 (20%)	0.6330	<1.5	−0.5639
North	0.0225	0.476 to 0.532 (40%)	0.1424	1.5 to 5	−0.1859
North East	−1.3963	0.532 to 0.565 (60%)	−0.1548	5 to 10	−0.0574
East	0.1149	0.565 to 0.595 (80%)	−0.2549	10 to 20	0.0826
South East	0.2600	>0.595 (100%)	−0.7976	20 to 40	0.3243
South	0.5443	Relief		>40	0.2609
South West	0.2491	<200 m	−2.2734	Distance to tranportation route
West	−0.4646	200 m to 250 m	0.5115	<200 m	−0.6888
North West	−0.8037	250 m to 300 m	−0.8372	200 m to 400 m	0.8066
Curvature		300 m to 350 m	−0.3458	400 m to 600 m	0.1560
Convex	0.4662	350 m to 400 m	0.5374	600 m to 800 m	−0.1275
Planar	−0.9954	>400 m	1.1512	>800 m	−0.4717
Concave	−0.4579	Slope		Wetness index
Distance to drainage		<5 degrees	−0.7958	<6.5	0.8230
0 m to 20 m	−1.2829	5 to 10 degrees	−0.1492	6.5 to 7.5	0.4102
20 m to 40 m	0.3846	10 to 15 degrees	0.2616	7.5 to 8.5	−0.3179
40 m to 60 m	0.3301	15 to 20 degrees	0.3715	8.5 to 9.5	−0.7865
60 m to 80 m	0.4372	20 to 25 degrees	0.6108	>9.5	−1.4191
80 m to 100 m	0.5190	25 to 30 degrees	0.3618	Rainfall
>100 m	0.2604	>30 degrees	1.1978	<1500 mm	−0.3175
Geology				1500 mm to 1550 mm	0.2118
Lower Middle Siwaliks	−0.1703			1550 mm to 1600 mm	−0.8334
Upper Middle Siwaliks	−0.8965			>1600 mm	0.5680
Upper Siwaliks	0.5541

Only the Jalad catchment was used in the calculation. The weight factors are either positive value or negative value. If the weight is positive, the factor is favourable for the occurrence of landslides, and if it is negative, the factor is not favourable for the occurrence of landslides. The weights were assigned to the classes of each thematic layer, and the resultant weighted thematic maps were overlaid and numerically added to produce a LHI map (figure 9). The weight integration approach used for the Jalad watershed is as follows:

Figure 9. Landslide hazard index map of the Jalad catchment after weights-of evidence modelling. Available in colour online.

In many cases of landslide hazard mapping, old landslide data are used to assess LHIs and new landslide data are used in validation of the model approach (van Westen 2003, Lee et al. 2007, Dahal et al. 2008a). However, in the case of the LHI map of the Jalad catchment, the accuracy of map was not known because all mappable landslides in the catchment were used to assess the LHI. To overcome this issue, in the present study, a different attempt has been made to validate the model. First of all, the weights of each class of causative factors obtained based on weights-of-evidence modelling of the Jalad catchment were replicated for the factors of the Charnath catchment. For this purpose, the class domain of thematic maps of all causative factors in the Charnath catchment were converted to value domain with respect to the weight value of each class obtained after weights-of-evidence modelling of the Jalad catchment (table 1). After the replication of weights of each class, weighted thematic maps of the Charnath catchment were prepared which were overlaid and numerically added to produce a LHI map of the Charnath catchment (figure 10). During replication, special consideration was made forNDVI and rainfall causative parameters. For NDVI, equal distribution classes (0–20%, 20–40%, 40–60%, 60–80% and 80–100%) of the Jalad catchment were considered as the same class of NDVI causative parameter in the Charnath catchment. For therainfall, hazard indexes of <1500 mm, 1500–1550 mm, 1500–1600 mm and >1600 mm estimated from the Jalad catchment were replicated forclasses <1450 mm, 1450–1500 mm, 1500–1550 mm and >1550 mm of the Charnath catchment. This assumption was made as per the fact that the consequences of higher range of rainfall for landsliding are same in both catchments.

Figure 10. Landslide hazard index map of the Charnath catchment after replication of weights of each class of causative factors from Jalad catchment. Available in colour online.

5.2. Logistic regression

When using the logistic regression model, the number of samples needed to create dependent variables is always an issue. The current literature has mainly established three kinds of practice (Zhu and Huang 2006). The first uses data from all over the study area, which leads to unequal proportions of landslide and non-landslide pixels (Guzzetti et al. 1999, Ohlmacher and Davis 2003). Usually, a large volume of data is taken in this method. The second approach uses all the landslide pixels and equal non-landslide pixels, which leads to a decrease in data number. The third practice uses all landslide pixels and an equal number of randomly selected pixels from areas free of landslides (Yesilnacar and Topal 2005). In this study, the first category, as defined by Guzzetti et al. (1999), was used for logistic regression modelling. As the dependent variable is dichotomous and the relationship between the dependent variable and independent variables is nonlinear, the logistic regression model was used on the causative factors. In this study the causative parameters were categorized into various classes and each parameter was a nominal variable. These variables were converted from nominal to numeric by coding.

The landslide density (Carrara 1983) used to transform a nominal variable to a numeric variable. This avoids the creation of an excessively high number of dummy variables. The following formula was used to calculate the densities of landslide in each class of causative factors:

where LD is landslide density, c_i is the area of ith class of a factor (such as slope class 10–15°), l_i is landslide area within ith class of a factor (such as slope 10–15°) and n is the total number of a certain factor. Then the parameter maps were overlaid with the landslide inventory map to calculate landslide densities based on equation (9). All the classes of causative factors were converted to numerical variable landslide density as shown in table 2. The domain of the landslide inventory map was changed from landslide present and landslide absent to numerical variables 1 (for landslide present) and 0 (for landslide absent). All spatial databases of causative parameters and landslide inventory were exported to statistical software and logistic regression analysis was performed for the Jalad catchment. In addition, logistic regression formulae were created as shown in equations (10).

Table 2. Landslide density in each class of causative factors of the Jalad watershed; the same landslide density was used for the Charnath watershed for logistic regression modelling.

Classes	LHI	Classes	LHI	Classes	LHI
Aspect		NDVI		Sediment transport index
Flat	0.0796	<0.476 (20%)	0.3239	<1.5	0.1017
North	0.1245	0.476 to 0.532 (40%)	0.2274	1.5 to 5	0.1427
North East	0.0324	0.532 to 0.565 (60%)	0.1797	5 to 10	0.1581
East	0.1340	0.565 to 0.595 (80%)	0.1655	10 to 20	0.1916
South East	0.1515	>0.595 (100%)	0.1036	20 to 40	0.1998
South	0.1890	Relief		>40	0.2061
South West	0.1497	<200 m	0.0151	Distance to tranportation route
West	0.0816	200 m to 250 m	0.2050	<200 m	0.1358
North West	0.0578	250 m to 300 m	0.0719	200 m to 400 m	0.3343
Curvature		300 m to 350 m	0.1050	400 m to 600 m	0.2281
Convex	0.5139	350 m to 400 m	0.2085	600 m to 800 m	0.1764
Planar	0.1665	>400 m	0.3945	>800 m	0.1255
Concave	0.3196	Slope		Wetness index
Distance to drainage		<5 degrees	0.0561	<6.5	0.3879
0 m to 20 m	0.0509	5 to 10 degrees	0.0878	6.5 to 7.5	0.2810
20 m to 40 m	0.1821	10 to 15 degrees	0.1223	7.5 to 8.5	0.1640
40 m to 60 m	0.1778	15 to 20 degrees	0.1359	8.5 to 9.5	0.1065
60 m to 80 m	0.1971	20 to 25 degrees	0.1711	>9.5	0.0606
80 m to 100 m	0.2189	25 to 30 degrees	0.1378	Rainfall
>100 m	0.1733	>30 degrees	0.2891	<1500 mm	0.1914
Geology				1500 mm to 1550 mm	0.2985
Lower Middle Siwaliks	0.3262			1550 mm to 1600 mm	0.1429
Upper Middle Siwaliks	0.1778			>1600 mm	0.3671
Upper Siwaliks	0.4960

All data are used in the logistic regression modelling and the overall model statistics of the regression conducted in this study are given in table 3. The Hosmer–Lemeshow test showed that the goodness-of-fit of the equation can be accepted because the significance of chi-square is 0.122. In the Hosmer–Lemeshow test, if significance of chi-square value is less than 0.05 the logistic regression model cannot be accepted (Zhu and Huang 2006, Chen and Wang 2007). The value of Cox and Snell R ² and Nagelkerke R ² showed that the independent variables can explain the dependent variable in an acceptaable way. Similarly, the Receiver Operating Characteristics (ROC) value of the model is 0.795. For a valid regression model, the ROC value should be greater than 0.5 (Hosmer and Lemeshow 2000). Significance ofthe Wald statistic for all causative factors is less than 0.05, and it confirms that independent variables have significantly predicted the outcome.

Table 3. Overall statistics of logistic regression model after using all causative parameters.

Hosmer–Lemeshow test			−2 log likelihood	Cox and Snell R ^2	Nagelkerke R ^2	ROC value
Chi-square	df	sig
12.726	8	0.122	10779.340	0.008	0.108	0.795
The Wald statistic and its significance for all factors
Factors	Wald	Sig.
Aspect	109.423	0.000
Curvature	4.076	0.044
Distance to drainage	144.621	0.000
NDVI	280.002	0.000
Rainfall	40.574	0.000
Relief	42.218	0.000
Slope	46.399	0.000
Geology	2.012	0.013
Sediment transport index	3.943	0.047
Distance to transportation routes	98.131	0.000
Wetness index	28.103	0.000

From the logistic regression modelling of the Jalad catchment, the following logistic regression equation was obtained:

The factors given in equation (10) represent the landslide density values.

For application of the model in the Charnath catchment, the value of landslide density of each class of all causative factors calculated from the Jalad catchment was replicated in the Charnath catchment, as shown in table 2, and equation (10) was applied to estimate the hazard probability of the Charnath catchment after logistic regression modelling. Here also, special consideration was made for NDVI and rainfall parameters. For NDVI, equal distribution classes (0–20%, 20–40%, 40–60%, 60–80% and 80–100%) of the Jalad catchment were considered as the same classes of NDVI parameter in the Charnath catchment. For the rainfall, landslide density of <1500 mm, 1500–1550 mm, 1500–1600 mm and >1600 mm estimated from the Jalad catchment were replicated for classes <1450 mm, 1450–1500 mm, 1500–1550 mm and >1550 mm of the Charnath catchment. Here also, this assumption was made as per the fact that consequences of higher range of rainfall for landsliding are same in both catchments. Moreover, the slight difference of rainfall value range (50 mm) in the rainfall classes certainly will not introduce a huge uncertainty in the result.

The distribution of the landslide hazard probability of the Charnath catchment is given in figure 11.

Figure 11. Landslide hazard index map of the Charnath catchment after replication of logistic regression model from the Jalad catchment. Available in colour online.

5.3. Validation

The LHI maps of the Charnath catchment presented in figure 10 (after weights-of-evidence modelling) and in figure 11 (after logistic regression modelling) were prepared without any consideration of previous landslides occurrences in the Charnath catchment. To understand its accuracy and to validate weights-of-evidence model and logistic regression modelling, these maps were used to estimate prediction rates. In the science of landslide hazard and susceptibility analysis, success rate and prediction rate are used to validate hazard or susceptibility indexes (Chung and Fabbri 1999, van Westen et al. 2003, Remondo et al. 2003, Lee 2004, Lee et al. 2007, Dahal 2008a, 2008b, Regmi et al. 2010a, von Ruette et al. 2011). The success rate indicates what percentage of all landslides occurs in the classes with the highest value of hazard. When old landslides are used for LHI calculation and new landslides are used for prediction, the calculated accuracy rate is called the prediction rate. In this study, instead of using new landslides in the same area (Jalad catchment), landslides from a different area (Charnath catchment) were selected for estimation of prediction rate.

A ROC curve approach is used to analyse the prediction accuracy of the proposed models. The ROC value is measure of success rate and prediction rate, usually obtained from ROC curve. These ROC curves give the area under the curve, and it is a measure of goodness of fit. The ROC curves of the Jalad catchment account the success rate of the models (figure 12), whereas ROC curves of the Charnath catchment (figure 13) depict prediction rate of the models. From the ROC curve, the success rate and prediction rate of weights-of-evidence model were 0.780 and 0.743, respectively. Similarly, the success rate and prediction rate of logistic regression model were 0.795 and 0.749, respectively.

Figure 12. ROC curve (area under the curve) of landslide hazard index estimated from the Jalad catchment after weights-of-evidence and logistic regression modelling. Available in colour online.

Figure 13. ROC curve (area under the curve) estimated from Charnath catchment after replication of weights-of-evidence and logistic regression models. Available in colour online.

5.4. Landslide hazard zonation maps

From the independent evaluation of accuracy of the weights-of-evidence and logistic regression models in the Siwaliks for landslide hazard assessment, it was found that the logistic regression model had slightly better accuracy than the weights-of-evidence model. Therefore, LHIs calculated from the logistic regression model were considered as a suitable model to prepare landslide hazard zonation maps of the selected catchments. Considering the prediction rate of logistic regression from ROC curve, five landslide hazard classes were defined: very low (<30% class of low to high LHI value), low (30–50% class of low to high LHI value), moderate (50–70% class of low to high LHI value), high (70–90% class of low to high LHI value) andvery high (>90% class of low to high LHI value, i.e. most higher LHI values). Next, landslide hazard zonation maps after logistic regression modelling of both catchments were prepared (figure 14).

Figure 14. Landslide hazard zonation maps of both catchments after logistic regression modelling. Available in colour online.

The resulting hazard zonation maps were again compared with the existing landslides of both catchments, and the results are given in figure 15. The results show that higher hazard classes correspond to 87% of the landslides in the Jalad catchment and more than 82% of landslides in the Charnath catchment.

Figure 15. Percentage coverage of existing landslides in different hazard zones of the Jalad and Charnath catchments. Available in colour online.

6. Conclusions

In this study, weights-of-evidence modelling and logistic regression modelling were used for landslide hazard mapping. Models were applied by considering 10 intrinsic factors. The thematic layers of all causative factors and existing slope failures were prepared in GIS (ILWIS 3.7). Mainly DEM-based data and field data were used to prepare data layers of causative factors. The conclusions of this study can be explicitly summarized as follows:

•	In many approaches of modelling of landslide hazard assessment in GIS, the model validation process is always dependent and the same landslide data used for analysis is considered in the verification. However, in this study, the model was verified independently, from landslides which were not used in the hazard analysis, and the results were promising with an independent prediction rate of about 0.75. This validates weights-of-evidence and logistic regression models for landslide hazard assessment in the slopes of Siwaliks of the Nepal Himalaya.
•	The evaluation of predictive factors used in this modelling suggests that the role of factors related to slope morphology, geology and human interventions (such as transportation routes) is found to be significant for landslide processes in the Siwaliks Range of Nepal. In particular, the Upper Siwaliks Group is very susceptible to landslides, and the transportation route to the forest area is enhancing forest degradation which ultimately contributes to landslide disasters in the study area.
•	In the science of landslide hazard mapping, replication of the hazard index or methodology is always a problem. In this study, the LHI estimated in one catchment was replicated to the other catchment, and the prediction rates of the models were estimated. This study confirms that replication of the landslide hazard analysis procedure is possible if the causative factors used in the models are similar in nature and origin.

Acknowledgments

We sincerely thank Mr. Bharat Prasad Pudasaini, Director General, the Department of Soil Conservation, Government of Nepal for his guidance and support for the research. We also acknowledge local community forest user groups of the study area for providing permission to enter the forest for investigation. Mr. Anjan Kumar Dahal is acknowledged for his technical support during the preparation of this paper. We thank Dr. Netra Regmi, Division of Earth and Ecosystem Sciences, Desert Research Institute, Reno, USA for his constructive comments. The authors are very grateful for the detailed review comments from three anonymous reviewers. The study has been partly funded by the Ministry of Forestry and Soil Conservation (Nepal), the Office of President of Nepal, CARE Nepal and Sasakawa Fund for Scientific Research, The Japan Science Society.