ML Rabonza (a), RP Felix (a,b), IJG Ortiz (a,b), IKA Alejandrino (a), DT Aquino (a), RC Eco (a,b), AMFA Lagmay (a,b)
(a) Nationwide Operational Assessment of Hazards, University of the Philippines, Diliman, Quezon City, Philippines
(b) National Institute of Geological Sciences, University of the Philippines, Diliman, Quezon City
Situated in the humid tropics, the Philippines will inevitably be a locus of climate-related disasters. Super Typhoon Haiyan, considered as one of the most powerful storms recorded in 2013, devastated the central Philippines region on 8 November 2013. In its wake, Haiyan left 6,190 fatalities, 28,626 injured and 1,785 missing, as well as damage amounting to more than USD 823 million. To mitigate damage from similar events in the future, it is imperative to characterize hazards associated with tropical cyclones such as those brought by Haiyan, with detailed studies of storm surges, landslides and floods. Although strong winds and powerful storm surges up 15-17 feet were the primary causes of damage, landslides studies are also vital in the rehabilitation of typhoon damaged areas. In order to delineate areas susceptible to rainfall induced shallow landslides and generate a worst-case scenario hazard map of the two provinces, Stability INdex MAPping (SINMAP) software was used over a 5-meter Interferometric Synthetic Aperture Radar (IFSAR)-derived digital elevation model (DEM) grid. SINMAP has as its theoretical basis in the infinite plane slope stability model. Topographic, soil-strength and hydrologic parameters (cohesion, angle of friction, bulk density and hydraulic conductivity) were used for each pixel of a given DEM grid to compute for the corresponding factor of safety. The landslide maps generated using SINMAP are found to be highly consistent with the landslide inventory derived from high-resolution satellite imagery dated 2003 to 2013. The landslide susceptibility classification found in the landslide hazard maps are useful to identify no-build, areas that can be built upon but with slope intervention and monitoring as well as places that are safe from shallow landslides. These maps complement the debris flow and structurally controlled landslide hazard maps that are also being prepared for rebuilding Haiyan’s devastated areas.
Typhoon Haiyan, a category-5 super typhoon, made landfall in central Philippines on 8 November 2013 leaving 6,190 fatalities, 28,626 injured and 1,785 missing, as well as damage amounting to more than USD 823 million . Haiyan recorded a maximum sustained winds of 235 km/h near its center and gustiness of up to 275 kph . Its powerful winds and heavy rains triggered widespread flooding and landslides leaving the central Philippines under a state of calamity. Although strong winds and storm surges up 15-17 feet were the primary causes of damage, landslides studies are also vital in the rehabilitation of typhoon damaged areas . Occurrences of shallow landslides depend on local terrain conditions such as slope-materials, topography, groundwater, and land cover in addition to triggering events, like torrential rainfalls and earthquakes which modify those characteristics and produce changes that cause slope instability . Slope instability is a geodynamic process that naturally shapes up the geomorphology of the earth; however, it becomes a major concern when those slopes affect the safety of people and property. Elevation of the ground water table due to pro- longed or heavy rainfall also naturally change the hydro-static and dynamic forces at the slope. This often initiate most of the slope failures in the Philippines. Existing slopes that have been stable for years can still experience significant movement when natural or man-induced conditions change their slope stability. Other natural conditions that can change slope stability are earthquakes, surface erosion, development of shrinkage or tension cracks followed by water intrusions, and bedrock weathering . With the growth of population, the demand for suitable and necessary infrastructure and other services increases. Therefore, in a country that is mostly hilly and mountainous, utilization of land on slopes is inevitable. It is therefore very important to map out unstable areas in order to ensure the safety of the people and delineate suitable areas for development.
The purpose of this study is to generate a shallow landslide susceptibility map for Samar and Leyte. The map will then be compared to a landslide inventory derived from satellite imagery dated 2003-2013. Unlike delineating hazard zones using historical data, maps generated using a deterministic slope stability GIS model, considers pixel-specific soil-strength parameters and hydraulic characteristics. Prior to field validation, the generated maps can provide a reasonable estimate of the hazard zones for shallow landslide susceptible areas. The shallow landslide susceptibility maps when complemented with debris flow and structurally controlled hazard maps are useful to identify no-build, areas that can be built upon but with slope intervention and monitoring as well as places that are safe from shallow landslides.
1.2. The SINMAP (Stability INdex MAPping) Approach to Stability Index Modelling
SINMAP is an ArcView GIS extension and an objective terrain stability mapping tool that complements other types of terrain stability mapping methods. This methodology presented in this paper is based upon the infinite slope stability model which balances the resisting components of friction and cohesion and the destabilizing components of gravity on a failure plane parallel to the ground surface. It is implemented to shallow translational landsliding phenomena controlled by shallow groundwater flow convergence. The slope stability theory does not apply to deep-seated instability including rotational slumps and deep earthflows . Slope failures occur frequently during or following period of heavy rainfall. The mechanism leading to rainfall-induced landslides can be summarized in general terms as follows: when rainwater infiltrates a soil profile that is initially in an unsaturated state, a decrease in negative pore pressure (or matric suction) occurs. This causes a decrease in the effective normal stress acting along the potential failure plane, which in turn diminishes the available shear strength to a point where equilibrium can no longer be sustained in the slope . Therefore, a slope stability modelling considering hydrologic wetness, soil-strength properties of soil and topography in a certain region can possibly predict zone of failure initiation in slopes. The input data required to implement the methodology include topographic slope, specific catchment area and parameters quantifying material properties (such as soil strength) and climate (hydrologic wetness parameter). The topographic variables are computed from digital elevation model (DEM) data. SINMAP does not require numerically precise input and accepts ranges of values to account uncertainty. Other input parameters are specified to SINMAP in terms of upper and lower bounds on the ranges they may take. The methods implemented in SINMAP software rely on grid-based data structure. Each of the input parameters is delineated on a numerical grid over the study area. The accuracy of the output is therefore heavily reliant on the accuracy of the input DEM .
The primary output of this modeling approach is a stability index (SI), the numerical representation of the stability of terrain and hence the possibility of landslide occurrence at each pixel in the study area. The indices are not to be interpreted as a numerically precise and are most appropriately interpreted as indications of ’relative hazard’.
The stability index is defined as the probability that a location is stable assuming uniform distributions of parameters over the uncertainty ranges. It does not predict that shallow translational slope movements will occur, but it forecasts that if they do, where they are more likely to initiate given the assumptions and input parameters used in the analysis. This value ranges between 0 (most unstable) and 1 (least unstable). Where the most conservative parameter ranges still results in stability, the stability index is defined as the factor of safety (ratio of stabilizing to destabilizing forces) at a location. This case yields a stability index value greater than 1. If the computed stability index value is less than 1, it is then defined as the probability that a location is stable given the range of the parameters used .
Table 1 presents the definition of terms of the SI based on broad stability classes. The selection of breakpoints (0.0, 0.5, 1, 1.25, 1.5) is subjective which requires interpretation. According to the model, regions that should not fail given the most conservative parameter used are classified as ’stable’, ’moderately stable’, and ’quasi-stable’. For this cases, the SI is defined as the factor of safety. Moreover, the terms ‘lower threshold and ‘upper threshold’ are used for regions where, as computed 29 by the model, the probability of instability is less than or greater than 50% respectively. To induce instability in these areas, external factors are not required. Failure may simply occur due to a combination of parameter values within the specified range. The term ’defended slope’ is used to classify regions where the model is not appropriate, or such slope are held in place by forces not represented in the model (e.g. bedrock outcrops, man-made slope protection) .
For mapping purposes in the Philippines, these 6 classes are reduced to 3 major hazard ratings with corresponding interpretations. Class 3 (1.25 > SI > 1.0) in Table 1 are classified as zones of ‘Low Susceptibility’. These regions must be built with slope intervention. Class 4 (1.0 > SI > 0.5) are classified as zones of ‘Moderate Susceptibility’. Such zones require slope intervention, protection and continuous monitoring. Class 5 (0.5 > SI > 0.0) are ‘No Build Zones’ and categorized as ‘High Susceptibility’ areas. The following color scheme is selected for each landslide hazard rating.
Background Multiple approaches are widely used to assess landslide hazards and slope stability :
- Field inspection using a check list to identify zones of landslide susceptibility
- Estimate future occurrences based on historical data
- Stability ranking based on criteria for slope, land form, lithology or geologic structure
- Probability analysis of failure based on slope stability models with stochastic hydrologic simulations
Each one is significant for certain applications. It is desirable to take full advantage of the fact that landslide source areas are, in general, strongly controlled by surface topography through shallow subsurface flow convergence, increased pore water pressure, reduction of shear strength and increased soil saturation. Recently, the availability of DEM data has prompted the development of methods that take advantage of the geographic information system (GIS) technology to calculate topographic attributes related to slope instability. GIS technology permits patterns of slope instability to be mapped at the scale of the DEM. The approach implemented in this paper combines steady state hydrologic concepts with the infinite slope stability model and incorporates grid-based DEM methodology. Parameter un- certainty is incorporated through the use of uniform probability distributions within upper and lower limits of the parameters. Moreover, in substitute to dynamic modeling over a range of rainfall events, the range of uncertainty of the hydrologic wetness parameter can be applied. In an approximate sense, this capability enables to predict the zones of landslide susceptibility without the additional input of weather data and analysis.
The variables theta and a are obtained from the DEM topography. Other parameters, C (cohesion), phi (soil angle of friction), R/T (recharge divided by transmissivity), and r (the ratio of water and soil density) are manually entered into the model. These are the more uncertain parameters and are set in terms of lower and upper bound values. The smallest C and tan phi together with the largest R/T define the most conservative (worst case) scenario within the assumed variability in the input parameters .
2. Geology of the Study Areas
Leyte Island has three major tectonostratigraphic terranes: Northwestern Leyte (part of Visayan Sea Basin), Leyte Central Highlands (has arc/ophiolite affinity), and northeastern Leyte (ophiolitic association in the Leyte Gulf). The northeastern segment of the Visayan Sea Basin is rep-resented by the sedimentary sequence found in northwestern Leyte. This has a north-north-west trending anticlinorium. Southwest, on the other hand, is represented by ophiolitic basement and Paleogene sedimentary rocks. Central Highland is dominantly underlain by igneous rocks. Its topmost layer is the Late Pliocene- Recent Leyte volcanic arc complex which is composed of andesitic volcanic cones and flows with minor basalt. This covers the old volcanic rocks of the Central Highland. The eastern Leyte has ophiolite complex overlain by three formations wherein the topmost layer is the Early Miocene Bagahuin formation which consists of conglomerate, sandstone, shale and fine tuffaceous sequences with intercalations of volcanic flows .
Samar Island has two stratigraphic groupings namely the Samar Block and the Leyte Gulf. Samar block constitutes the whole island except for southwestern Samar since it is part of the Leyte Gulf together with the islands of Bucas Grande, Dinagat, Homonhon and Siargao. Samar ophiolite constitutes the basement rocks overlain by different formations. Balo, Hagbay and Catbalogan formations are distributed in areas of Samar. Late Cretaceous Balo formation, found in Bagacay and San Jose de Buan, consists of limestones, conglomerate, sandstone, mudstone and shale. Middle Miocene Hagbay formation, distributed in Barrio Hagbay, San Jose de Buan, consists of reefal limestone and siltstone. Middle Miocene – Early Pliocene Catbalogan formation is found in Catbalogan, Loguilocan and Bassey. It is composed of marl, siltstone, sandstone, pebble and conglomerate .
3. Soil Type
3.1. Leyte Soil
Leyte province has an extensive cover of clay soil (Fig. 5). Clay loam, fine sand, hydrosol, sandy loam, silt loam and silt are also present but covers only relatively smaller portions of Leyte. Mountainous areas are located in places identified only in the soil map as rough mountainous land .
3.2. Samar Soil
The soil type in Samar varies from hydrosol, clay, clay loam, loam, loamy sand, sandy loam to beach sand. Majority of the area on the west is covered by clay loam (Fig. 6). Faraon clay type extends from north to south of Samar’s central area. However, on the eastern side, soil is only described as undifferentiated .
The following figure outlines the methodology for creating a landslide susceptibility map using the coupled hydrological and deterministic slope stability model (SINMAP).
4.1. Digital Elevation Model Grid Data
DEM data were acquired from the National Mapping and Resource Information Authority (NAMRIA) of the Philippines in September 2013. These are Interferometric Synthetic Aperture Radar (IFSAR) – derived digital terrain models (DTM) with 5-meter resolution and 0.5 meter vertical accuracy.
4.2. Geotechnical Data
Parameters not determined by the DEM, i.e. the geotechnical data, are considered more uncertain and are specified in terms of upper and lower boundary values . The following text rationalizes the specific input values and distributions used in this study.
The equation used to determine the dimensionless cohesion combines root and soil cohesion. Theoretically this is the ra-tio of the cohesive strength of the roots and soil relative to the weight of a saturated thickness of soil.
4.2.2. Internal Friction Angle
Internal angle of friction is the measure of the shear strength of soil due to friction. It can be determined in the laboratory using Direct Shear Strength or Triaxial Stress Test. Although no independent soil analysis was completed, the 31 to 37 degrees soil-friction angles used to realistic for the study area since the study requires that parameters be generalized over large areas with variation encompassing the properties of several different soil types. This is more realistic than providing a small range of input values.
4.2.3. T/R (Ratio of Transmissivity to Effective Recharge)
The ratio of transmissivity of the soil (square meters) to the effective recharge (m/hr), when multiplied by the sine of the slope, the T/R value can be interpreted as the length of the hillslope (in meters) necessary to develop saturation (Pack et. al, 1998b). The recharge rates used in this study have been derived from a lower and upper precipitation values or limits, i.e., 50 mm/day and 200 mm/day. 50 mm/d rate was chosen as minimum rate and 200 mm/day is used as a maximum, i.e., an extreme example of the precipitation that can produce shallow translational movement. For comparison, most municipalities affected by Haiyan experienced an accumulated rainfall of 60-100 millimeters over a three-hour period (Project NOAH, 2013).
The transmissivity rate (m2=hr) was calculated using basic equation
Where K is the hydraulic conductivity of the soil and b is the soil depth in meters. In the context of shallow landslide movement, soil thickness (b) of 1.5 meters is assumed to be uniform in the simulation. Data on Hydraulic conductivity is obtained from Hamazaki’s Database on Red-Yellow and Related Soils in the Philippines Part 2 Visayas and Mindanao Soils. Soil classifications (from the Soil Texture Map of Bureau of Soils and Water Management, 2013), descriptions, and test results, along with literature values for soil properties given in Hammond and others (1992) were used to constrain reasonable ranges of soil input parameters for the stability index modeling. Parameter values (primarily dimensionless cohesion, soil thickness, internal friction angle, and hydraulic conductivity) were then adjusted within reasonable ranges to maximize the number of slope movement locations per area. Based on the dominant soil classification in Samar and Leyte, the lower and upper limit values for the parameters are estimated to correspond to a worst case scenario simulation.
4.2.4. Landslide Inventory
A landslide inventory has been completed for the study areas using high- resolution satellite imagery. These landslide points were later used to calibrate and validate the model. The stability index map produced by SINMAP (6 classes) was further reclassified into 3 classes (Low, Medium, High) to obtain the landslide susceptibility map and was compared with the landslide inventory locations. The relationship between land use and landslide susceptibility classes was also obtained to locate the vulnerable areas.
5. Results and Discussion
Figure 10 is the SINMAP simulation of Leyte province. It shows that flat land areas are identified as stable and mountainous areas as zones with dominantly lower and upper thresholds. Stable zones account for 69.35% of Leyte’s total land area whereas unstable zones account for 22.75% of the area.
Marginal unstable zones covers 5.80% of the area. This type of zone becomes unstable when minor destabilizing forces are applied such as surface erosion, development of shrinkage or tension cracks followed by water instrusions, and bedrock weathering. The accuracy of the model is verified by overlaying the inventoried landslides over the simulated model. Table 3 shows the summary of the landslide inventories under dierent stability class. Under the unstable zones predicted by the model, 85.82% of the total inventoried landslides fall under it. Another 1.4% is predicted as part of marginal unstable zones. However, 12.68% of the inventoried landslides fall under stable zones and hence not predicted by the model. Figure 11 shows a close up view of an area with identified landslides. Three of them is located in areas identified as unstable and 1 in marginally unstable.
SINMAP model for Samar province is shown in figure 12. High hazard areas are mostly concentrated on the central mountainous part of the province. These high hazard areas together with areas classified as build only with slope intervention are considered as unstable areas which constitute 28.81% of Samar’s total area. Another 7.48% of the area is part of marginal instability zones. Stable areas which is dominantly found near the coastlines cover 66.72% of the total land mass area. Table 4 shows the accuracy of the model with respect with the landslide inventories. Landslide points found in predicted unstable zones comprise 75.49% of the total inventoried landslides. No landslide point falls under the marginal unstable zone but 24.51% of the points fall under the zone predicted as stable. A zoomed in view in an area in the northern part of Samar is shown in Figure 13 where most of the identified landslides are located in the predicted unstable zone.
In both simulated models, there are no landslide points identified in some areas regarded as mostly unstable. This is due to lack of high resolution images in those particular sites for landslide inventory to be done. On the other hand, those landslide points found in predicted stable areas can be explained by the fact that some of these are old landslides. There may be some cases when their location became more stable after the land-slide event so when recent DEMs are used, they are shown as already stable areas. Another reason for this is that the calibration parameters used has its own limitations.
The landslide maps generated using SINMAP are found to be 81.40% accurate with the landslide inventory derived from high-resolution satellite imagery dated 2003 to 2013. Accuracy can be further improved through field detailed geological and geotechnical assessment of the study areas with preference in areas having satellite images with lower resolution and areas with high count of landslide inventories. The landslide susceptibility classification found in the land- slide hazard maps are useful to identify no-build, areas that can be built upon but with slope intervention and monitoring as well as places that are safe from shallow landslides. These maps complement the debris flow and structurally-controlled landslide hazard maps that are also being prepared for rebuilding Haiyan’s devastated areas.
We thank the National Mapping and Resource Information Authority (NAMRIA) for the Digital Elevation Model data. Also to Bureau of Soil and Water Management (BWSM) for the soil maps.
 NDRRMC, Ndrrmc situation report on the effects of typhoon Yolanda, January 22, 2014 (6:00 a.m.) @ONLINE (Jan 2014). URL http://www.gov.ph/
 A. F.-P. KD Suarez, Typhoon Yolanda weakens as it exits PH @ONLINE (Nov 2013). URL http://www.rappler.com/
 R. Soeters, C. Van Westen, Slope stability: recognition, analysis and zonation, in: Landslides: investigation and mitigation, National Academy Press, Washington, D. C., 1996, pp. 129–177.
 J. F. L. S. K.M. Weerasinghe, H.V.M.P. Abeywickrema, Use of a deterministic slope stability predicting tool for landslide vulnerability assessment in Ratnapura area, Sri Lanka, Geo Informatics Center (GIC)of the Asian Institute of Technology (AIT), 2003.
 D. R. Montgomery, W. E. Dietrich, A physically based model for the topographic control on shallow landsliding, in: Water Resources Research, 1994, pp. 1153–1171.
 R. Orense, Slope failures triggered by heavy rainfall, Philippine Engineering Journal.
 C. G. A. P. R.T. Pack, D.G. Tarboton, SINMAP 2 – A Stability Index Approach to Terrain Stability Hazard Mapping, Utah State University (August 2005).
 R. Pack, D. Tartabon, C. Goodwin, Terrain stability mapping with SINMAP, Terratech Consulting Ltd., Salmon Arm, B.C., Canada, report number 4114-0 Edition, report and software available online: http://moose.cee.usu.edu/sinmap/sinmap.htm (1998).
 A. Witt, Using a geographic information system (GIS) to model slope instability and debris flow hazards in the french broad river watershed, North Carolina, Master’s thesis, North Carolina State University (2005).
 C. Hammond, D. Hall, S. Miller, P. Swetik, Level i stability analysis (LISA) documentation for version 2.0: General technical report int-285, Tech. rep., U. S. Department of Agriculture , Forest Service, Intermountain Research Station (1992).
 MGB, Geology of the Philippines, 2nd Edition, Mines and Geosciences Bureau, North Avenue, Quezon City, Philippines, 2010.