Shallow Landslide Susceptibility Mapping for Selected Areas in the Philippines Severely Affected by Super Typhoon Haiyan

ML Rabonzaa,∗, RP Felix, a,b, IJG Ortiz, a,b, IKA Alejandrino a, DT Aquinoa, RC Eco, a,b, AMFA Lagmay a,b

a Nationwide Operational Assessment of Hazards, Department of Science and Technology, Philippines
b National Institute of Geological Sciences, University of the Philippines, Diliman, Quezon City

Abstract

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.

1. Introduction
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 [1]. Haiyan recorded a maximum sustained winds of 235 km/h near its center and gustiness of up to 275 kph [2]. 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 [2]. 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 charateristics and produce changes that cause slope instability [3]. Slope instability is a geo-dynamic process that naturally shapes up the geomorphology of the earth; however, it becomes a major concern when those slopes a ect 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 [4]. 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.

s1

1.1. Objectives
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.

s2

1.2. The SINMAP (Stability INdex MAPping) Approach to Stability Index Modelling SINMAP is an ArcView GIS extension and an objective ter- rain 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 [5]. 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 matricsuction) occurs. This causes a decrease in the e ective 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 [6]. 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 [7].

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 [8].

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) [7].

For mapping purposes in the Philippines, these 6 classes are reduced to 3 major hazard ratings with corresponding interpre-tations. Class 3 (1.25 ¿ SI 1.0) in Table 1are 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 interven-tion, protection and continuous monitoring. Class 5 (0.5 ¿ SI ¿ 0.0) are ‘No Build Zones’ and categorized as ‘High Suscep-tibility’ areas. The following color scheme is selected for each landslide hazard rating.

s3

Background Multiple approaches are widely used to assess landslide hazards and slope stability [5]:
1. Field inspection using a check list to identify zones of landslide susceptibility
2. Estimate future occurrences based on historical data
3. Stability ranking based on criteria for slope, land form, lithology or geologic structure
4. 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.

s4

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 [8].

2. Geology of the Study Areas

2.1. Leyte

Leyte Island has three major tectonostratigraphic terranes: Northwestern Leyte (part of Visayan Sea Basin), Leyte Central Highlands (has arc/ophiolite anity), 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 tu aceous sequences with intercalations of volcanic flows [11].

2.2. Samar

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 di erent 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 lime- stones, 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 [11].

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 [11].

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 typeextends from north to south of Samar’s central area. However, on the eastern side, soil is only described as undi erentiated [11].

s5

4. Methodology

The following figure outlines the methodology for creating a landslide susceptibility map using the coupled hydrological and deterministic slope stability model (SINMAP).

s6

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.

s7

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 [8]. The following text rationalizes the specific input values and distributions used in this study.

4.2.1. Cohesion
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.

s8

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 di erent soil types. This is more realistic than providing a small range of input values.

4.2.3. T/R (Ratio of Transmissivity to E ective Recharge)
The ratio of transmissivity of the soil (m2) to the e ective 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 200mm/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 a ected 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

T=Kb

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

5.1. Leyte
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.

s10

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 di erent 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.

s11

5.2. Samar
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.

s12

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.

Conclusion

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.

s13

Acknowledgements

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.

References
[1] NDRRMC, Ndrrmc situation report on the e ects of typhoon yolanda,
january 22, 2014 (6:00 a.m.) @ONLINE (Jan 2014).
URLhttp://www.gov.ph/

[2] A. F.-P. KD Suarez, Typhoon yolanda weakens as it exits ph @ONLINE
(Nov 2013).
URLhttp://www.rappler.com/

[3] R. Soeters, C. VanWesten, Slope stability: recognition, analysis and
zonation, in: Lanslides: investigation and mitigation, National Academy
Press, Washington, D. C., 1996, pp. 129–177.

[4] J. F. L. S. K.M. Weerasinghe, H.V.M.P. Abeywickrema, Use of a deter-
ministic slope stability predicting tool for landslide vulnerability assess-
ment in ratnapura area, sri lanka, Geo Informatics Center (GIC)of the
Asian Institute of Technology (AIT), 2003.

[5] D. R. Montgomery, W. E. Dietrich, A physically based model for the to-
pographic control on shallow landsliding, in: Water Resources Research,
1994, pp. 1153–1171.

[6] R. Orense, Slope failures triggered by heavy rainfall, Philippine Engineer-
ing Journal.
[7] C. G. A. P. R.T. Pack, D.G. Tarboton, SINMAP 2 – A Stability Index Ap-
proach to Terrain Stability Hazard Mapping, Utah State University (Au-
gust 2005).

[8] R. Pack, D. Tartabon, C. Goodwin, Terrain stability mapping with
SINMAP, Terratech Consulting Ltd., Salmon Arm, B.C., Canada, re-
port number 4114-0 Edition, report and software available online:
http://moose.cee.usu.edu/sinmap/sinmap.htm (1998).

[9] A. Witt, Using a geographic information system (gis) to model slope in-
stability and debris flow hazards in the french broad river watershed, north
carolina, Master’s thesis, North Carolina State University (2005).

[10] 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 Re-
search Station (1992).

[11] MGB, Geology of the Philippines, 2nd Edition, Mines and Geosciences
Bureau, North Avenue, Quezon City, Philippines, 2010.

Leave a Reply

Your email address will not be published. Required fields are marked *

two × 5 =