Intergration of GIS Using GEOSTAtistical INterpolation Techniques (Kriging) (GEOSTAINT-K) in Deterministic Models for Landslide Susceptibility Analysis (LSA) at Kota Kinabalu, Sabah, Malaysia

A practical application for landslide susceptibility analysis (LSA) based on GEOSTAtistical INterpolation Techniques (Kriging) (GEOSTAINT-K) for a deterministic model was used to calculate the factor of safety (FOS) and failure probabilities for the area of Kota Kinabalu, Sabah. In this paper, the LSA value can be expressed by a FOS value, which is the ratio of forces that make the slope fail and those that prevent the slope from failing. A geotechnical engineering properties data base has been developed on the basis of a series of parameter maps such as effective cohesion (C’), unit weight of soil (), depth of failure surface (Z), height of ground water table (Zw), Zw/Z dimensionless (m), unit weight of water (w), slope surface inclination (β) and effective angle of shearing resistance (). Taking into consideration the cause of the landslide, identified as groundwater change, two scenarios of landslide activity were studied. Scenario 1 considered the minimum groundwater level recorded corresponding to the actual situation of the most recent landslide while Scenario 2 considered the reverse. A simple method (infinite slope model) for error propagation was used to calculate the variance of the FOS and the probability that will be less than 1 for each pixel. The highest probability value of the various scenarios was selected for each pixel and final LSA 1 (scenario 1) and LSA 2 (scenario 2) maps were constructed. The validation between the examined LSA 1 and LSA 2 maps and the results of the landslide distribution map (LDM) were evaluated. This deterministic model had higher prediction accuracy. The prediction accuracy was 81 % and 85 %, respectively. In general for both factors, the LSA 2 map showed higher accuracy compared to the LSA 1 map. The resulting LSA map can be used by local administrators or developers to locate areas prone to landslides, determine the land use suitability and organize more detailed analysis of the “hot spot” areas identified.


Introduction
Landslide is an issue which is still under debate in the newspapers or any electronic media, especially in Kota Kinabalu, Sabah, Malaysia.The effect of the pressure of development activities that lead to rapid cut works or reclamation of slopes for road construction, infrastructure development and construction of dwellings or buildings is more widespread and has spread to hilly areas with a large population.The location of the study area is surrounded by the Crocker and Trusmadi Ranges and has a complex geological background, which also give a negative outlook for any exploration activities and further land development planning.
Landslides are amongst the most damaging natural hazards in the world.The term ''zonation'' in a general sense implies a division of the land into areas and their classification according to degrees of actual or potential landslide hazard or susceptibility (Varnes, 1984).The purpose of landslide susceptibility analysis (LSA) is to highlight the regional distribution of potentially unstable slopes based on a detailed study of the factors responsible for landslide (Aleotti and Chowdhury, 1999;Ayalew et al., 2005).LSA is defined as a quantitative or qualitative assessment of the classification, volume (or area) and spatial distribution of landslides which exist or potentially may occur in an area (International Society of Soil Mechanics and Geotechnical Engineering (ISSMGE), www.engmath.dal.ca/tc32).Susceptibility may also include a description of the velocity and intensity of the existing or potential landsliding.
In the literature, various approaches to deterministic models for LSA have been developed.Some popular deterministic models for LSA are the infinite slope, SHAllow Landsliding STABility (SHALSTAB), Stability INdex MAPping (SINMAP), Transient response, Transient Rainfall Infiltration and Grid-based Regional Slope-stability analysis (TRIGRS), etc.The infinite slope model is a static stability model in which local stability conditions are determined by the means of the local equilibrium along a potential slip surface.Other models couple the infinite slope stability model with more or less complex rainfall infiltration models (Okimura and Kawatani, 1986;van Westen, 1993;Montgomery and Dietrich, 1994;Dietrich et al., 1995;Terlien et al., 1995;van Westen and Terlien, 1996;Dymond et al., 1999;and Crosta and Dal Negro, 2003).SHAllow Landsliding STABility (SHALSTAB) model is a model that combines a hydrologic model (O'Loughlin, 1986) with an infinite slope stability equation, the Mohr-Coulomb failure law (Bolt et al., 1975), for the prediction of slope instabilities based upon the minimum amount of steady-state rainfall required to trigger landsliding (Montgomery and Dietrich, 1994).It required inputs are obtained from a DEM that is widely available within the U.S. and a few representative values of geotechnical parameters, such as soil bulk density, internal angle of friction and water table depth.This model calculates pore pressures for steady-state saturated water flow parallel to the slope plane.Pack et al. (1998) andZaitchik et al. (2003) combined a slope stability model (Stability INdex MAPping, SINMAP) with a steady-state hydrology model in selected watersheds of northern Vancouver Island, British Columbia and in the central highlands of Honduras, respectively.The main difference between these two models is that SHALSTAB assumes zero soil cohesion because of the spatial and temporal heterogeneity of soil cohesion (and therefore the difficulty in obtaining values) and because assuming a zero cohesion value results in the most conservative estimate of slope instability (Dietrich et al., 2001).However, new versions of the model do allow for the inclusion of soil cohesion.Great attention should be paid to the accuracy and variability associated with the input parameters.Transient response model developed by Iverson (2000) uses unsaturated flow to calculate pore pressures and vertical flow.The International Institute for Aerospace Survey and Earth Sciences (ITC) has developed a GIS called the Integrated Land and Water Information System (ILWIS) that has modules incorporated in the GIS for deterministic instability zonation (Van Westen, 1997a).The Level I Stability Analysis (LISA) prepared for the U.S. Forest Service by Hammond et al. (1992) uses average estimates for geotechnical parameters in their model.Similar examples of regional modelling and prediction of shallow landslides using a transient rainfall infiltration model in combination with slope stability calculation (Transient Rainfall Infiltration and Grid-based Regional Slope-stability analysis; TRIGRS) were applied for the Seattle area, Washington State, USA (Baum et al., 2005) and the Umbria Region, central Italy (Salciarini et al., 2006).The TRIGRS model predicts a larger area of instability than the area that actually failed, mainly due to uncertainty in soil thickness, local variation in soil properties, and Digital Elevation Model (DEM) errors (Chang and Kim, 2004;Dahal et al., 2008;Safaei et al., 2010;2011)

Materials and Methods
GEOSTAtistical INterpolation Techniques (Kriging) (GEOSTAINT-K) utilizes the statistical properties of the measured points.The purpose is to determine the probability of certain variables occurring over an area where identifying every possible location would be impossible.The approach used in GEOSTAINT-K method uses the concept of a combination of deterministic models and geostatistical interpolation.A deterministic-based approach, a landslide stability model, is advantageous over other approaches on allowing us to calculate the quantitative value, the factor of safety.The factor of safety represents the ratio between the sheer stress for failure and the sheer stress for stability.The variables in the equation were identified by many landslide models.
Many methods are associated with geostatistics, but they are all in the kriging family (Cressie, 1988;1990;Rivoirard, 1994& Stein, 1999).Ordinary, simple, universal, probability, indicator, and disjunctive krigings along with their counterparts in cokriging are available in Geostatistical Analyst (ArcGIS software).However, the simple kriging method was used in this research because it can use either semi-variograms or covariances (autocorrection), transformations and allow for measurement errors.Not only do these methods create prediction and error surfaces, but they can also produce probability and quantile output maps depending on the requirement.
GEOSTAINT-K method is divided into two distinct tasks: quantifying the spatial structure of the data and producing a prediction.Quantifying the spatial data structure, known as variography, is fitting a spatial-dependence model to the data.To make a prediction for an unknown value for a specific location, kriging uses the fitted model from variography, the spatial data configuration, and the values of the measured sample points around the prediction location.GEOSTAINT-K is a moderately quick interpolation method that can be exact or smoothed depending on the measurement error model.It is very flexible and allows the user to investigate graphs of spatial autocorrelation.It uses statistical models that allow a variety of map outputs including predictions, prediction standard errors, standard error of indicators, and probability (Cressie, 1988;1990;Rivoirard, 1994& Stein, 1999).The flexibility of these methods require a lot of decision making and assumes the data comes from a stationary stochastic process, which is a collection of random variables that are ordered in space and/or time.
In this paper, the LSA degree can be expressed by the factor of safety (FOS), which is the ratio between the forces that make the slope fail and those that prevent the slope from failing.F-values larger than 1 indicate stable conditions, and vice-versa (Table 1).At F=1 the slope is at the point of failure.Many different models exist for the calculation of FOS (as described above).In this study one of the simplest models, the so-called infinite slope model was used.This two dimensional model describes the stability of slopes with an infinitely large failure plane.It can be used in a GIS, as the calculation can be done on a pixel basis.In GIS environment, the grid form data of a layer consists of variable layers and calculates corresponding variable values using map algebra function.The pixels in the parameter maps can be considered as homogeneous units.For example, a grid value has 20 in the slope layer, and then the grid has uniformly 20°.The effect of the neighbouring pixels is not considered, and the model can be used to calculate the stability of each individual pixel, resulting in a hazard map of FOS.The FOS is calculated according the following formula (Brunsden and Prior, 1979;van Westen and Terlien, 1996): in which: c' = effective cohesion (kPa= kN/m 2 ), γ = unit weight of soil (kN/m 3 ), z = depth of failure surface below the surface (m), zw = height of groundwater table above failure surface (m), m = zw/z (dimensionless), γw = unit weight of water (kN/m 3 ), β = slope surface inclination (°) and ' = effective angle of shearing resistance (°).
Thus, the GEOSTAINT-K can be used either on profiles as well as on pixels, as shown in Figure 1.The entire analysis requires firstly the preparation of the data-base.Since this is not available for the Kota Kinabalu area, the process was simplified.

Phase 1: Landslide Hazard Identification Phase (LHIP)
Phase 1 as depicted in Figure 1 is the preliminary stages in landslide hazard assessment consisting of the landslide hazard identification phase (LHIP).LHIP involves three (3) main types of research namely desk, field and laboratory studies.The desk studies involved detailed aerial photograph interpretation and satellite images (using Erdas V.9.2 software) analyses, extensive literature review and secondary data collation.All of these sources were analysed and reclassified to get an idea or preliminary information about the landslide distribution and historical data aspects in the study area.The product from these desk studies were used to establish a landslide incidents background data-base and to generate the "Landslide Distribution Map" (LDM) for the study area by using Arc GIS V.9.3 software.The field studies in LHIP involved sampling of rocks and soils, engineering geology mapping, and observation of landslide hazard characterization parameters and extracting a digital elevation model (DEM).For the laboratory studies in LHIP, all samples of rocks and soils obtained from the field were analyzed and evaluated for their engineering properties in accordance to the standards recommended by the ISRM (1979;1979b& 1985) and BS1377-1990 (Method of Test for Soils for Civil Engineering Purposes) such as the direct shear test for rock testing and the classification of grain size, Atterberg limit, and triaxial test (Consolidated isotropically undrained, CIU) for soil testing.After all the laboratory studies are completed, several LHIP thematic maps were produced such as effective cohesion (c'), unit weight of soil (), depth of failure surface (Z), height of ground water table (Zw), Zw/Z dimensionless (m), unit weight of water (w), slope surface inclination (β) and effective angle of shearing resistance ().

Phase 2: Landslide Hazard Assessment Phase (LHAP)
In Phase 2, taking into consideration the cause of the landslide, identified as groundwater change, two scenarios of landslide activity were studied.Scenario 1 considered the minimum groundwater level recorded corresponding to the actual situation of the most recent landslide while Scenario 2 was vice-versa.A simple method for error propagation was used to calculate the variance of the FOS based on Equation (1) and the probability that will be less than 1 for each pixel.The slope stability calculation was carried out by a combination of the input parametric maps with the GIS operations using a grid base.Finally, the resultant LSA maps were compared and validated with the data from LDM.The highest probability value of the various scenarios was selected for each pixel and the final LSA map constructed.As the FOS in the final LSA map is assigned as a single value in every cell of a raster map, it is convenient to perform a reclassification of the calculated values (Table 1).

Geological Background
The Geological map portraying information about the main geological units is shown in Figure 2. The different rock compositions and textures affect slope instability, influencing strength, permeability, and susceptibility to chemical and physical weathering of the rock masses.This map also represents the structural setting of the study area.Features such as sequence and type of layering, lithologic changes, planes, joints, faults and folds are accountable for LSA.The exposed rocks in the study area and its surroundings vary in type and age, from Late Eocene-Early Miocene sandstone and shale of the Crocker Formation to Young Alluvial sediments which are still being deposited (Rodeano et al., 2006).The sandstone-siltstone-shale unit is defined by an interbedded sandstone and shale with occasional siltstone.The thickness of the individual beds ranges from 2 cm to 130 cm.The sandstone is normally fine to very fine-grained and highly fractured while the shale layers are sheared.The shale unit is generally composed of red and grey types of shale.The grey variety is occasionally calcareous.This alternating sequence is commonly interbedded with siltstone or very fine grained sandstone.The shale comprises about 12% of the total volume of the Crocker Formation.The sandstone composition is dominated by quartz with subordinate amounts of feldspars and chloritized, illitized or silicified lithic fragments.Calcareous fractions are rare.These are poorly sorted and well compacted with the pores filled by fine grained detritus or squeezed lithoclasts resulting in very low to nil primary porosity.The sandstone unit is characterized by very low to nil porosity but moderate to high secondary permeability.It is defined by its great thickness, medium to very coarse-grained and sometime pebbly.Thin shale or siltstone bed between 3 to 40 cm thicknesses occurs between the thick sandstone beds.The argillaceous beds are frequently the site of shearing while the sandstone beds are the site of fracturing or jointing.The alluvium is restricted to the low-land areas.It mainly represents unconsolidated alluvial sediment on river terraces and flood plains composed of unsorted to well-sorted, sand, silt and clay of varying proportions which were derived from upstream bed rocks.They occur in irregular lenses varying in form and thickness.Towards the coastal area, the alluvium becomes finer-grained and interbedded with argillaceous deltaic and marine strata.The alluvium may also consist of a very thin layer of organic matter.The alluvium sediment is soft, compressible and may be prone to settlement.

Landslide Distribution
Landslide distribution is very useful information to study the physical changes and the latest geological assessment of their vulnerability.It quantifies all the information for any complex phenomenon.In these studies, landslides were classified in terms of types, materials involved, estimated volumes and velocities, degrees of activity, and return periods; and distinctions are made between source and depositional areas.Approximately about 2,119 landslide locations have been identified through an extensive review of literature, aerial photograph interpretations and field studies (Figure 4).The landslides were classified into several types; 20% of fall, 30% of translational, 20% of rotational, 15% of flow, and 15% of complex (a combination of fall, slide and/or flow).In terms of landslide scale, the study area consists of small (<50 m 3 ) (20%), medium (50 m -500 m 3 ) (60%) and large (> 500 m 3 ) (20%).

Digital Elevation Model (DEM)
A digital elevation model (DEM) of the slope conditions provided by raster datasets on morphometric features (altitude, internal relief, slope angle, aspect, longitudinal and transverse slope curvature and slope roughness) and on hydrologic parameters (watershed area, drainage density, drainage network order, channel length, etc.) were automatically extracted from the DEM (Figure 5).In addition, the slope angle is also considered as an index of slope stability caused by the presence of a digital elevation model (DEM) which is evaluated numerically and is illustrated by the spatial analysis (Yalcin and Bulut, 2007).In terms of slope gradient, the results suggest that 48.37% of the area can be categorized as 0 o -5 o , 28.45% as a 6 o -15 o , 22.41% as 16 o -30 o , 0.75% as 31 o -60 o and 0.01% in excess of 60 o .Areas with slope angles in excess of 30 o represent a very steep slope segment in the study area where the steeper a slope, the higher the LSA value.

Geotechnical Engineering Properties
Geotechnical engineering properties are closely related to the mechanical behaviour of soil and rock, and their equilibrium.Since different geotechnical engineering properties have different LSL values they are very important in providing data for susceptibility studies.For this reason it is essential to group the geotechnical engineering properties properly (Carrara et al., 1991;Dai et al., 2001;Cevik and Topal, 2003;etc.).Soil and rock shear strength is an important engineering property.It is a fundamental property that governs the stability of natural and constructed slopes.It is not a unique value, but is strongly influenced by loading, unloading, and especially by water content.The shear strength is basically described as a function of normal stress on the slip surface (σ), cohesion (c), and internal angle of friction (φ).The relationship of these properties to other characteristics of natural soil has been given by Terzaghi and Peck (1967), and Fredlund and Rahardjo (1993).Geotechnical engineering properties of seventy two (72) soil samples indicated that the soil materials mainly consist of poorly graded to well graded materials of clayey, silty to sandy soils, which are characterized by low to high plasticity, effective cohesion (c') ranges 5.98 kPa to 38.00 kPa, unit weight of soil () ranges from 14.85 kN/m 3 to 23.14 kN/m 3 , depth of failure surface (Z) ranges from 3.96 m to 14.26 m, height of ground water table (Zw) ranges from 1.97 m to 5.19 m, Zw/Z dimensionless (m) ranges from 4.92 to 12.63, effective angle of shearing resistance () ranges from 3.88 o to 20.54 o , tan ranges from 0.08 o to 0.33 o , slope surface inclination (β) ranges from 42.65 o to 65.53 o , sinβ and cosβ ranges from -1.00 o to 1.00 o ; and cos 2 β ranges from 2.26e-013 -1.00 (Figure 6).

Landslide Susceptibility Analysis (LSA) Maps
After all the thematic maps in Phase 1 (LHIP) were produced, all of these parameters were compiled and computed according to the equation ( 1) and reanalysed together with LDM to generate two new thematic maps for different scenarios known as "Landslide Susceptibility Analysis (LSA 1 and LSA 2) maps" (Figures 7 and 8).
The resultant LSA 1 and LSA 2 maps provided a relative assessment of the landslide area using the FOS value.The resultant LSA 1 map suggests that 38% of the area can be categorised as Very Low Susceptibility (VLS), 14% as Low Susceptibility (LS), 28% as Moderate Susceptibility (MS), 6% as High Susceptibility (HS), 6% as Very High Susceptibility (VHS) and 8% as Extremely High Susceptibility (EHS).While the LSA 2 map suggests that areas of VLS to LS decreased to 33% and 12%; and areas of MS to EHS moderately increased to 29% of MS, 8% of HS, 8% of VHS and 10% of EHS.A comparison of LSA1 and LSA 2 indicate that β and Zw parameter factors have the highest influence on landslide instability.This is evidenced by the changes in the percentages of MS to EHS, which increased about 1 % to 2 % and decreased from 2 % to 5 % in the LS and VLS.In general, VLS to MS areas (FOS = > 1.00) refer to the stable conditions with flat to moderately steep slopes with pasture and this area is highly recommended for any future development planned (Table 2).In contrast, HS to EHS areas (FOS = < 1.00) represent unstable conditions with steep slope segments.HS to VHS areas are basically not recommended to be developed due to geological, hydrological and geotechnical constraints.However, if there is no alternative or the developer or the local authorities really want to develop this area, some procedures to be observed are as stated in Table 2. EHS areas are strictly not recommended to be developed and should have provisions, and suitable non-structural works planning control as shown in Table 2 are recommended.

Verification of the Landslide Susceptibility Map
For validation of landslide susceptibility calculation models, two basic assumptions are needed: (1) landslides are related to spatial information, such as topography and geology, and (2) future landslides will be precipitated by a specific impact factor such as rainfall or earthquake (Brabb, 1984;Chung & Fabbri, 1999;Varnes 1978).In this study, the two assumptions are satisfied because the landslides were related to the spatial information, and the landslides were precipitated by heavy rainfall in the study area.The LSA results were validated using known landslide locations.Verification was performed by comparing the known landslide location data (Figure 4) with the LSA maps (Figures 7 & 8).The comparison is shown in Figure 9 for the minimum (LSA 1) and maximum (LSA 2) groundwater level models.Figure 9 illustrates how well the estimators perform with respect to the landslides used in constructing those estimators.To obtain the relative ranks for each prediction pattern, the calculated index values of all cells in the study area were sorted in descending order.The success rate validation results were divided into 100 classes with accumulated 1% intervals, according to the landslide susceptibility index value.As a result, considering all the factors used in the study area, the 90-100% (10%) class, with the highest possibility of landslide, contains 38% and 62% of the landslide grid cells in success rate using the minimum (LSA 1) and maximum (LSA 2) groundwater level models and continuously until the calculation of 0% -100% (100%) of 100% of the total area in the rates obtained in this model, respectively.To compare the result quantitatively, the areas under the curve were re-calculated with the total area assumed to be 1, which means perfect prediction accuracy (Lee & Dan, 2005;Hyun et al., 2010).Hence, the area under a curve can be used to assess the prediction accuracy qualitatively.The area ratios were 0.810 and 0.847, which indicate prediction accuracy of 81% and 85 %, respectively.Overall the cases where both factors and maximum (LSA 2) groundwater level model was used showed a higher accuracy than cases where the minimum (LSA 1) groundwater level model used.

Conclusions
In light of the available information, the following conclusions may be drawn from the results of this study: GEOSTAINT-K is a revolutionary technology that provides a dynamic environment with a wide variety of tools and a friendly wizard interface to explore data, analyze anomalies, and optimally display an interpolated surface with associated uncertainties.It gives the user the power to fully understand the qualitative and quantitative aspects of their data.By providing the user with the freedom to statistically predict and model situations and incorporating a powerful visualization tool, GEOSTAINT-K effectively bridges the gap between deterministic -geostatistics models and GIS.
The benefit of GEOSTAINT-K is to provide insight and options for decision-making in practical problems.The benefits include: It encourages a rational, systematic approach in assessing the safety of slopes, and a framework to put uncertainties and engineering judgment into a system and allows comparison of hazard and risk for different slopes.The GEOSTAINT-K gave us a quantitative value for landslide susceptibility level, which was theoretically proven.
LSA allows collection, management, analysis and dissemination of a large amount of data for extensive areas.
All of these actions, based on continuous scientific and technologic research, with a strong multidisciplinary component and the involvement of local, regional, and interregional authorities, allow effective regional land-use planning.
This deterministic model has high prediction accuracy.The prediction accuracy is 81 % and 85 %, respectively.In general for both factors, LSA 2 map showed higher accuracy compared to cases of the LSA 1 map.In an iterative process the optimization of the LSA value are positive as long as the quality of the input data is good and sufficient knowledge exists on the relation of the occurrence of the triggering mechanisms in relation to the occurrence of landslides.
GIS geospatial technology capability of LSA provides a valuable tool for gaining susceptibility level estimates at the regional scale.This result highlights the importance of the potential effects of landslides in the study area.
The resulting LSA can be used by local administration or developers to locate areas prone to landslide, determine the land use suitability of the area, organize more detailed analysis in the identified "hot spot" areas and manage the impact of landslide disasters that may affect the regional economy (loss and damage to property) or welfare of the community (deaths and loss of homes) (risky areas).

Table 3 .
Calculation of the ratio of the average area under the curve