Pr oj e c tNOAH Ope nFi l eRe por t s

Vol. 1 (2013), pp. 10-21, ISSN 2362 7409

Landslide hazard mapping of Leyte Province, Philippines using remote sensing, numerical modeling, and GIS RN Ecoa,b , DT Aquinoa , AMF Lagmaya,b,∗, I Alejandrinoa , MK Alemaniaa , AA Bonusa,b , CM Escapea,b , R Felixa,b , PK Ferrera,b , RC Gacusana,b , J Galanga , F Llanesa,b , PK Luzona,1 , M Magcamita , KR Montalboa,b , J Obriquea , IJ Ortiza,1 , C Quinaa , M Rabonzaa , V Realinoa,b , JM Sabadoa,b , H Ariolaa , JJ Sulapasa,b , S Salvosaa , NL Timbasa,c a Nationwide

Operational Assessment of Hazards, Department of Science and Technology, Philippines Institute of Geological Sciences, University of the Philippines- Diliman c Agricultural Systems Cluster, University of the Philippines, Los Banos

b National

Abstract When Super Typhoon Haiyan struck the Philippines in November 2013, the province of Leyte was among the most heavily damaged. Detailed hazard maps of Leyte’s cities and municipalities–including storm surge, flood, and landslide hazard maps– could have mitigated the damage by providing crucial information on the hazard extent. However, high-resolution hazard maps were not available to communities prior to Haiyan. This study seeks to generate detailed 1:10,000 scale hazard maps for Leyte Province to aid in future disaster prevention and mitigation efforts. Techniques in remote sensing, numerical modeling, and GIS were combined using a 5-meter IfSAR-derived Digital Terrain Model (DTM) to rapidly generate high-resolution community-level landslide hazard maps for the entire province. All of the 45 landslides identified in Leyte plot within the shallow and deep-seated landslide hazard areas mapped. The detailed landslide hazard maps are useful for identifying safe and unsafe areas for rehabilitation and reconstruction efforts for Haiyan devastated areas, as well as for community disaster planning and preparedness efforts. Keywords: Haiyan, Yolanda, Leyte, landslide, GIS, remote sensing, COLTOP, Conefall, debris flow, Flo-2D, hazard mapping, Matterocking, SINMAP

1. Introduction

dead, 28,689 injured and 1,601 missing. It also caused about PhP89.6 billion worth of damages to property, agriculture, and infrastructure (NDRRMC, 2013). One of the hardest hit provinces was Leyte. Unfortunately, the only available hazard maps prior to the disaster were drawn to 1:50,000 scale. It is insufficient to capture the details at the community level and therefore not useful at all for detailed planning and emergency response situations. Large-scale detailed hazard mapping can be tedious and time-consuming to create and even more so in regions with extensive rugged terrain as found in Leyte. Given the urgency of the reconstruction and rehabilitation efforts of the Haiyan-affected areas, downscaled hazard maps that can be produced with the least time possible are needed to help communities identify safer areas for rehabilitation. More importantly, the next event may come soon, before maps from traditional field survey are finished. This may yet result in another preventable disaster, if only such hazard maps were produced and disseminated more quickly. In this study, we first provide a general overview of the geophysical setting of the province. We then use an integrated approach to generate a provincial landslide hazard map. Several techniques were used, namely analysis of satellite imagery to identify recent landslides, deterministic numerical modeling to generate slope stability indices, structurally controlled landslide susceptibility and debris flow hazard maps, and geomorphic analysis using GIS tools to identify and delineate alluvial

Typhoons are a common occurrence in the Philippines, with an average of about 20 typhoons entering the Philippine area of responsibility each year (David et al., 2013). However, the continued heating of the oceans, which leads to the formation of low-pressure areas, is likely to increase the incidence of extreme weather events such as intense rainfall, tropical cyclones, storm surges, landslides and other hydro-meteorological hazards (IPCC, 2007). The effects of future climate change will almost certainly increase the risks faced by communities located in these hazard-prone areas. Such heightened risks were made even more glaring when Supertyphoon Haiyan (Philippine name Yolanda) made arrived in the Philippines on 7 November 2013 at 0400H local time, packing maximum sustained winds of 235 kph and gusts of 275 kph. It made initial landfall in Guiuan, Eastern Samar at that time, then by 0700H of the same day it was traversing through Tacloban, Leyte. After passing through several provinces in the Central Philippines, Haiyan exited the Philippine Area of Responsibility on 9 November (Fig.1a). The National Disaster Risk Reduction Management Council reported 6,300 people I Published online on 24 October 2014 at http://blog.noah.dost.gov.ph/2014/10/24/leyte-landslide-assessment/ ∗ Corresponding author Email address: [email protected] (AMF Lagmay)

10

fans. This method takes full advantage of the 5-meter resolution IfSAR-derived digital terrain model provided by the National Mapping and Resource Information Authority (NAMRIA), the national mapping agency of the Philippines. We then compare the results of the numerical models with the landslide inventory to see if the results of the models are consistent with the remote sensing observations.

1.2.2. Drainage Patterns The drainage system is largely influenced by the geology and topography of the province having a volcanic central highland flanked on both sides by large plains. The Leyte Valley is drained into the sea by the Binahaan, Daguitan, Guinaroma, and Marabang Rivers. The headwater of Binahaan River is Lake Danao. The accordant drainage patterns are generally either subparallel or dendritic. Most of the streams on the central part are offset along the traces of the Philippine Fault. Subparallel streams on the western side drains toward the Visayan Sea. Other areas has generally dendritic patterns following the slope of the terrain and implying a homogeneous underlying rock material. The rivers are usually short and intermittent and drainage is excessive on cultivated hills. Radial drainage patterns are also evident especially on volcanic edifices. Streams on the eastern side merge with the Pagsangahan River emptying into Ormoc Bay. Steep slopes in the area contribute to its highly unstable state and high channel failure potential. The transition from a high to a moderate inclination of side slopes along its channel indicate a continuous and progressive downflow of water towards Ormoc City. 1.3. Geologic Setting 1.3.1. Regional Geology Northwestern Leyte is part of the Visayan Sea Basin’s northeastern portion (MGB, 2004). The Malitbog Ophiolite forms the basement rocks of western Leyte, is an almost complete ophiolite sequence that includes serpentinized harzburgite, minor dunite, cumulate and isotropic gabbro, diabase dike complex, pillow basalts, and pelagic sedimentary rocks . Overlying this sequence are conglomerates, sandstones, siltstones, tuff, and limestones. Occurrence of a turbidite sequence composed of sandstone, siltstone, and mudstone was also described (Florendo, 1984; Jurgan, 1980; Corby et al., 1951a; Pilac, 1965).

Figure 1: a) Track of Supertyphoon Haiyan (Track from IBTrACS https://www.ncdc.noaa.gov/ibtracs). b) Physiographic and geographic setting of Leyte Province

1.1. Geographic Setting Leyte is one of the largest provinces in the Eastern Visayas Region. It is bounded to the north by the province of Biliran separated by the Biliran Strait and bounded to the south by Southern Leyte. To the west, across the Camotes Sea lies the island of Cebu. To the east and northeast is the San Juanico Strait that separates Leyte with the island of Samar (Figure 1B). The latter is the easternmost island in the Visayas and forms the barrier between the Visayan Sea and the western portion of the Philippine Sea. Southeast of Leyte lies the Leyte Gulf. 1.2. Geomorphology 1.2.1. Relief and Topography The province of Leyte is generally mountainous and broken by steep slopes, especially at the central portion of the island. These are expressions of a chain of volcanoes trending NNWSSE, straddling the trace of the Philippine Fault. This also gives the island well-expressed linear depressions, including elongated lakes and displaced or offset rivers. In the Tongonan area, a volcanic cone is observed to be displaced left-laterally (Lagmay et al., 2003), and deep escarpments along sag ponds and pressure ridges are also seen to align with the orientation of the fault (Aurelio et al., 1994). Where there are no mountain ranges, the topography of the island is generally flat to gently rolling. To the east and west of these are narrow plains that transition to hilly and rolling terrains going to the coasts on both sides of the island (Figure 1B).

Figure 2: Active faults and geologic map of Leyte (Edited from: http://www.phivolcs.dost.gov.ph). (R) to (PG1) are sedimentary rocks while (BC) is for pre-Jurassic basement complex of undifferentiated amphibolites, mics schists, quartzofeldspars and phylliste-slates. (NI) to (UC) are intrusive rocks while (QVP) to (N1) are volcanic rocks. Map modified from the Phivolcs active faults map and MGB geological map.

Central Leyte is primarily composed of igneous rocks owing to the presence of the volcanic chain (Figure 2A). Diorite, andesite, basalt, dacite, and pyroclastic rocks have been 11

identified in this area. In addition, occurrences of sedimentary rocks such as conglomerate, sandstone, shale, and limestone have been noted. Intruding and covering older volcanics are late Pliocene to Recent andesitic volcanic cones and flows, later named the Leyte Volcanic Arc Complex (MGB, 2004; Aguilar, 1995; Cabantog and Escalada, 1989; Pilac, 1965). Eastern Leyte is host to the Cretaceous-aged Tacloban Ophiolite, characterized by serpentinized harzburgites, dunites, gabbro, sheeted dikes, basalts, pillow basalts, and pelagic sedimentary rocks (MMAJ-JICA, 1990; Balce and Cabantog, 1998). Early Miocene sedimentary rocks, composed of conglomerates, sandstones, shales and fine tuffaceous sequences with intercalations of volcanic flows, are found found overlying older volcanic rocks. Middle Miocene limestone interbeds have also been noted to occur within the area, along with sandstones and marly tuffaceous shales with basal conglomerates, and is dated Late Miocene to Pliocene (MGB, 2004; Pilac, 1965).

land, well-drained flat lowland, and well-drained rolling uplands (Barrera et al., 1954). The poorly drained soils in Leyte are found in the low-lying areas around Palo, Tanauan, and Alangalang (Pawing Series); around the municipality of Palo, and extending to the inner regions at the foot of the central Cordilleras near Jaro and Dagami, Alangalang and Barugo to the north, and in Abuyog to the south (Palo Series). Both the Pawing and Palo series developed from recent alluvial deposits and exhibit minimal soil profile development (Barrera et al., 1954). The moderately drained flat lowland soils belong to the soil series of San Manuel, Umingan, Dagami, Mandawe, and Bantog series. These soils are found in Babatngon, Palompon, Hilongos, in the alluvial soil formation on Western Leyte from Ormoc to Baybay, in the municipality of Dagami, on the riverbanks along the Hilongos and Bangerahan Rivers, and in the municipality of Villaba. These soils, which are mostly cultivated for agriculture, developed from alluvial deposits. They 1.3.2. Tectonic Setting exhibit moderate soil profile development and level to nearly Leyte Province is traversed by the Central Leyte Fault (PHIVOLCS, level topography (Barrera et al., 1954). 2008), a segment of the Philippine Fault (Figure 2A). The PhilipUnder well-drained flat lowland soil are the Obando, Uminpine Fault is a 1,200-km-long left-lateral strike slip fault which gan, and Taal series. These series developed from marine decuts across the Philippines. The fault was identified using aerial position. These are found on the eastern shores of Leyte, occuphotographs because the trace is evident in the morphology pying the narrow coast from Palo to Abuyog. In addition, this of landforms (Kimura et al., 1968; Aurelio, 1992). The fault soil has low organic matter content and a texture that ranges is typically manifested in fault scarps, elongated depressions, from coarse sand to fine sandy loam allowing water to freely sag ponds and compressive ridges. Mountain range offsets, percolate through the profile (Barrera et al., 1954). drainage system patterns also present evidence that the fault The well-drained rolling upland soils are divided into nonmovement is left-lateral although a vertical component has been calcareous and the calcareous soils. The noncalcareous soils noted as well (Pinet and Stephan, 1990; Barcelona, 1981; Nakata include the Guimaras, Tacloban, Guimbalon, Luisiana, Palomet al., 1977; Allen, 1962; Rutland, 1968). There are 2 to 3 wellpon, Malitbog, and Maasin soil series. These soils are found in defined branches with curvilinear trace convex to the northeast, the foothills of the eastern slope of the range running from Palo and striking N25-30W (MGB, 2004). West of the Philippine to Babatngon, in the low range of hills and mountains between Fault in Leyte, the trends range from 28 to 42 degrees, while Palo and Babatngon, in Ormoc and in between Capoocan and east of the Philippine Fault in Leyte, the trends have a 298Dagami, in the upper area of Ormoc and in Barugo, in the mudegree to 318-degree range (Lagmay et al., 2003). nicipalities of Villaba and Palompon, in the small area south of Folds are also notable in Leyte Island, particularly in the Baybay, and in the interior regions from Inopacan to Malitbog. Batang Formation, the Taog Formation, and the Babatngon Schist. Another soil series found in Leyte is the Lugo series, which is The Batang Formation, composed of a sedimentary sequence, a calcareous soil. This soil is found in the municipalities of San is observed to be composed of two thrust slabs that are eroIsidro, Calubian, and in some areas of Villaba (Barrera et al., sional remnants of an antiform having a northeast trend (Flo1954). rendo, 1987). Taog Formation in western Leyte is also a folded sequence of sedimentary rocks (Corby et al., 1951b), while Ba2. Methodology batngon Schist has highly folded, well-foliated, and crenulated schists with schistose planes verging steeply in an east-west di2.1. Landslide Inventory rection (Pilac, 1965). The landslide inventory map produced for the province of According to PHIVOLCS (2008), there are 18 volcanoes Leyte was created through the use of high-resolution satellite within the province, although only one is potentially active– imagery. Sets of satellite images covering the study area were Cancajanag Volcano–and the rest are inactive (Figure 2a). Caninspected for landslides, which were then delineated as polycajanag is a 0.9-kilometer high dome complex, composed mainly gons and marked with points using ArcGIS. Data such as the of andesite. There is an active volcano in the southeastern-most total area of landslides (total damaged area in the province) tip of the island, Cabalian Volcano, however it is outside the and the total area of landslides with respect to the area of the province. province were also computed (Fig. 3). Satellite imagery used include SPOT, Quickbird, and World1.4. Soil Series view. SPOT imagery has a spatial resolution of 5 to 10 meters. The soils of Leyte Province (Figure 2B) are categorized as poorly drained flat lowland, moderately drained flat low12

hw = height of water R = recharge r = water density (ρw ) to soil density (ρ s ) ratio T = soil transmissivity = soil hydraulic conductivity x h φ = soil internal angle of friction θ = slope A FS <1 indicates unstable conditions, whereas a FS >1 indicates stable conditions given the assumptions and parameters based on the user input. SINMAP then assigns a stability based on the computed factors of safety. The stability zones are assigned the following relative hazard rankings: High with 0
Figure 3: Landslide inventory process flow.

Quickbird has a spatial resolution of 0.6 to 2.4 meters. Worldview has a spatial resolution of 0.5 meters. The satellite images used are from the years 2002 to 2012. The dates of the satellite images, however, are not uniform throughout the study area. Satellite images were selected based on their availability in Google Earth, a freeware program. Landslides detected in the satellite image were marked with ‘dots’ and delineated as polygons. These were saved as vector files so that they could be easily accessed and analyzed using GIS programs. The types of landslides that were easily identified in satellite images–and that were mainly considered in this inventory–were flows and slides. Both the total area of landslides in the study area and the total area of landslides with respect to the area of the province were computed. 2.2. Shallow landslide susceptibility The shallow landslide hazard map was produced using SINMAP (Stability Index MAPping) software, an ArcGIS 9.3 extension developed by Pack et al. (1998). SINMAP computes a factor of safety using the infinite slope model (Pack et al., 1998; Hammond et al., 1992) and based on the input data for hydrologic and geotechnical data of the soil and topographic data (Figures 4) for each pixel on a digital elevation model. The factor of safety (FS) is a dimensionless number that represents the ratio of the stabilizing forces to destabilizing forces at a location. Applying the infinite slope model, it can be expressed by h   i a C + cosθ 1 − min TR sinθ , 1 r tanφ FS = (1) sinθ

Figure 4: Process flow of shallow landslide susceptibility mapping using SINMAP

The model input parameters include upper and lower bound values for recharge to shallow groundwater system, soil transmissivity, and other soil properties (i.e., unit weight, depth of landslide, effective internal friction angle, and effective cohesion). SINMAP randomly samples the bounded input parameter values using a uniform probability distribution to account for the variability and uncertainty inherent to the natural system. Soil classifications, descriptions, previous test results (constant head permeability test and soil classification from the Bureau of Soils and Water Management and direct shear test from Philippine Geo-Analytics) along with literature values for soil properties (Hammond et al., 1992) were used to obtain reasonable ranges of soil input parameters for modeling the stability index.

where: FS = factor of safety a = topographic catchment area C = dimensionless cohesion = (Cr + C s )/hρ s g Cr = root cohesion; C s = soil cohesion h = soil thickness; ρ = soil density; g = gravity constant 13

To simulate the landslide susceptibility of the area, information on geology, soil type and elevation were acquired. The soil texture is generally clayey (Table 1). In this process, though landslide inventory may be used as an input, it is only used as an assessment for the accuracy of the model. Parameters 2*T/R Lower Upper 2*Cohesion Lower Upper 2*Phi Lower Upper soil

2.3.2. Rockslide zone and runout propagation zone identification Matterocking and Conefall were the software used in mapping structurally controlled landslide prone areas. Matterocking locates unstable slopes that fall under specified failure conditions (Wyllie and Mah, 2004). For planar sliding to occur; (1) slope >dip of discontinuity, (2) orientation of slope is within ± 20◦ of the discontinuity’s dip direction , and (3) slope >critical 45◦ angle of failure (Baillifard et al., 2003). Wedge failure could occur if (1) slope >dip of the line of intersection of two discontinuities, (2) normal of the slope is within ± 10◦ the direction of two discontinuities’ angle bisector, and (3) slope >critical 45◦ angle of failure. From the results of the Lineament Pattern Analysis and Coltop 3D, the most common discontinuity sets are keyed into the software, simulating both planar and wedge failures. Conefall is a supplementary software to Matterocking used to roughly estimate the extent of rock falls given a source point. The source points used were the unstable slopes from Matterocking. The structural landslide propagation extent in Conefall is designed for application of a cone slope as boundary condition. The energy of the rockmass to travel is proportional to the difference in height of the cone slope and DTM expressed by the equation.

Clay 40 400 0.606 1 18 26 1746

Table 1: Calibration parameters used as input to SINMAP.

2.3. Structurally-controlled landslide susceptibility 2.3.1. Lineament Pattern Analysis Lineament patterns are identified through morphologic interpretation of topographic data. The Coltop3D software is used as a part of the validation process. Coltop3D is a pseudo 3Dimensional topographic analysis software. It can simultaneously represent the slope aspect and slope angle of a digital terrain model. It provides a unique color representation for a given combination of dip and dip direction using classical Hue Saturation Value with a Schmidt-Lambert projected stereonet. Using Coltop3D, trends could be isolated and plotted as Rose Diagrams to explain the frequency of the given orientation. The resulting rose diagrams are compared to existing literature on the geology of the area, or if possible through measurements obtained in the field.

4h =

v2 (x) 2g

(2)

The 20◦ cone slope value is based on the data gathered on the actual landslide event last 2006 in Guinsaugon, Southern Leyte.

Figure 5: a) Lineament analysis of Leyte Province c) Coltop and Rose Diagrams of Leyte Province.

14

where i, j are the flow directions, pif d the susceptibility proportion in direction i, tan βi the slope gradient between the central cell and the cell in direction i, and x the variable exponent. The spreading is similar to the multiple flow direction when x=1, and as x → ∞, the divergence results into a single flow direction (Horton, 2013). According to Claessens et al. (2005), on the basis of field and laboratory measurements, the suggested value of the exponent for debris flows is 4. A factor dh, which is a change of height of the central cell, is also introduced in the model that will alter the gradient values to allow smoothing of DEM roughness and production of a more consistent spreading (Horton et al., 2013). For high resolution data, a dh value of 2 m is suggested. Inertial algorithm used for implementing the flow proportion in direction i according to the persistence is proportional to weights of the angle between previous direction and the direction from the central cell to the cell i 7. This distribution avoids backward propagation to save computing time. The implemented weights used are the values suggested by (Horton et al., 2013) which have of the spreading values of 1, 0.8, 0.4, 0, and 0 for the angles 0, 45, 90, 135, and 180, respectively.

2.4. Alluvial Fan Mapping Alluvial fans were initially delineated in ArcGIS 10, a mapping software, using a 5-meter IfSAR-derived digital elevation model (DEM). Using ArcGIS, watersheds and tributaries were extracted, and contour maps were generated from the DEM. The contour maps were used to manually delineate the fanning formation located directly downstream of the watersheds and then converted to shapefiles. These shapefiles were then overlaid on top of aerial and satellite imagery to verify if the location of the fan apex coincides with the location of the watershed and tributaries in these imageries (Figure 6). Both aerial and satellite imageries were also used to check for communities and built-up areas lying within the alluvial fans to create a hazard map. This hazard map serves as a guide for on-site fan validation to measure slopes and to characterize alluvial deposits through the analysis of outcrops.

Figure 6: Process flow of alluvial fan delineation using IfSAR

2.5. Debris Flow Simulation Figure 7: Spreading of susceptibility value to the neighboring cells (Horton et. al., 2013).

Debris flows are one of the most destructive mass wasting phenomena with sediments carried by a finer matrix of speeds ranging from 2 to 40 kilometers per hour. According to Rickenmann and Zimmermann (1993) and Takahashi (1981), critical combination of terrain slope, water input, and sediment availability is significant for the initiation of a debris flow. In this study, the generation of debris flow was done using Flow-R, a routing software used extensively in mud and debris flow simulations. IfSAR-derived digital terrain model and its corresponding flow accumulation data, obtained using a GIS software, are used for the simulation. The modified Holmgren (1994) algorithm was used for the spreading assessment of the flow direction. It adds a parameter to the multiple flow direction algorithm as an exponent x allowing control of the divergence: (tan βi ) x pif d = P8 { tan β > 0x ∈ [1; +∞[ x j=1 (tan β j )

Calculation of the runout distance follows the simplified friction-limited model (SFLM). This model is based on the maximum possible runout distance, which is characterized by a minimum travel angle that connects the source area to the farthest point reached by the debris flow (Horton et. al., 2013). A value of 3c irc was chosen for travel angle as it best fits the event that happened during the 2012 New Bataan, Compostela Valley debris flow event during Typhoon Bopha (Pablo). With an estimated maximum velocity of debris flow in New Bataan being 17 m s−1 , it is safe to choose a limit of 20 m s−1 for the simulations in this study.

(3)

15

3. Results 3.1. Landslide Inventory Landslide inventory for the province of Leyte yielded a total of 139 points (Table 2), each of which was identified and 74 of which had been delineated to calculate the damaged area. The total damage area calculated for Leyte Province was 211,749 square meters. The landslides within the province were noticeably concentrated in two regions: a western group and an eastern group, separated by the plains within the municipalities of Kananga, Ormoc City and Leyte (Figure 5). Most of the hazards occurred at the elevated region of the province which, in turn, coincides with the path of the Central Leyte Fault. The city of Ormoc in Leyte exhibited the highest damaged land cover at 67,444 square meters (Table 2). Municipality

No. of Landslides

Abuyog Alangalang Albuera Dagami Hindang Kananga Leyte Mahaplag Matalom Merida Ormoc City Tacloban City

5 2 2 1 4 2 12 9 1 6 26 4

Damaged areas (DA) in sq.m. 41,569.5 10,546.9 3,081.4 1,861.4 3,400.9 41,040.5 1,170.2 16,096.5 3,473.5 10,306.3 67,444.6 11,757.5

Figure 8: Example of a shallow landslide in the municipality of Mahaplag (A) and a deep-seated landslide within Ormoc City (B) that were identified in the inventory.

DA/Total DA [%]

they mobilize the bedrock. Also, they are usually more than 10 m in depth and have very prominent head and clearly defined main body (Lin et al., 2013). Therefore, the huge landslides wherein the bedrock had been exposed were classified as deep-seated landslides while small landslides that only involved movement of soil and weathered materials were classified as shallow landslides (Fig. 8).

19.6 5.0 1.5 0.9 1.6 19.4 0.6 7.6 1.6 4.9 31.9 5.6

Table 2: Inventory of landslides per municipality in Leyte

Most of the landslides that were identified were located in elevated regions especially beside rivers that cut through mountainsides. In the municipalities of San Isidro and Tabango, most of the landslides were located on rolling hills that had been modified for farming and/or had experienced vegetation loss. On the other hand, at the municipalities of Merida and Tacloban, landslides were observed to occur near roadcuts. This is a very common pattern found in landslide inventories all throughout the Philippines. Due to the roadcuts susceptibility for slope instability, road blockage due to landslides during typhoons is a recurring hazard. The recognition of deep-seated versus shallow landslides using satellite images had been the subject of numerous studies because the principles behind them are modern, powerful tools for landslide inventory mapping (Guzzetti et al., 2012). The difference between shallow and deep-seated landslides, is that fresh shallow landslides could be easily delineated using color differences between the vegetation and the landslide itself and they are less likely to be reactivated, therefore they pose less risk. On the other hand, deep-seated landslides are difficult to recognize especially if they occur at densely forested areas. They could also be reactivated during intense events wherein 16

Figure 10: Shallow Landslide Susceptibility Map of Leyte

Figure 9: Landslide Inventory for Leyte Province Figure 11: Landslide points within hazard maps produced from different models

3.2. Shallow Landslides A summary of the results of the SINMAP model shows that the total land of Leyte, with an area of 5,562.1km2 has 2.8% high landslide susceptibility: 30.5% under moderate and 8.5% under low landslide susceptibility. High susceptibility areas are found mostly on the central mountain range area. Unstable areas are concentrated on the western part of Villaba, San Isidro and Merida which are dominantly hilly. Low to moderate susceptibility are observed on moderately inclined areas of the whole province. Flat areas did not register as susceptible to landslides (Fig. 10). Validation of the SINAMP model results were done by: 1. Adding and comparing the locations of landslide inventories to the predicted instability with areas of actual instability and; 2. Evaluating the accuracy of the SINMAP model results. The accuracy of the SINMAP simulation was evaluated by overlaying the landslide points over the model and counting the number of landslides that fell within each susceptibility class (High, Moderate, Low). All 139 landslide points in Leyte based on the landslide inventory, are found in unstable areas as determined by SINMAP model. The majority of the points fell within the moderate (41.7%) and high susceptibility (58.3%) zones (Figure 11). There were no landslides found in low susceptibility and stable areas.

3.3. Structurally-controlled Landslides 3.3.1. Lineament Pattern Analysis The Philippine fault traverses the central mountain ranges of Leyte province. Along this line are conjugate shears and corresponding splays. Manually delineated lineaments of Leyte DTM represent these morphological features (Fig. 12). Location of the prominent trends are the slopes considered as discontinuities controlling structural failure. Stereoplots of the identified sets show observable similarity to a study of 44 established sinistral strike slip fault’s rose diagram (Fig. 13) (Abbassi and Shabanian, 1999).

17

Figure 14: Area of potential rockslide zones in Southern Leyte using Matterocking Figure 12: Lineament Analysis for Southern Leyte

3.3.2. Potential Rockslide Zones and Propagation Failure zones generated by Matterocking are classified as highly susceptible areas for structurally controlled landslide (Fig. 14). High activity of the Philippine fault segment in the province of 0.55cm/year (Cole et al., 1989) up to 3.5cm/year (Duquesnoy, 1997) tracks the line of high susceptibility in the mountain ranges of the province. Generally steep slopes and presence of discontinuities are the main factors that would cause instability as specified in the model. It is the case of wedge failure that could trigger most rockslides. This is due to the persistence of intersecting discontinuity sets. A total of 5 percent of the province’s area is classified as potentially unstable slope. Large portion of the municipalities of Leyte, Capoocan, Carigara, Jaro, Ormoc, Dagami, Albuena, Burauen, La Paz, Baybay City, Javier, Inopacan, Mahaplag, Abuyog, Hindang, Babatngon and Tacloban are susceptible areas. These slopes are used as sourcepoints to propagate in the digital terrain model. Rock mass propagation extent generated by Conefall (Fig. 15) have increased the susceptible areas due to structural failure to 21 percent of total area of the province.

Figure 13: Different rose diagrams for a segment of the Philippine fault and some parts of Leyte

18

Figure 16: a) Alluvial fans of Leyte Province with debris flow simulation and fan polygons overlaid on b) Google satellite and c) Bing aerial imagery. Only the fans which coincide with the debris flow simulation were included in the hazard map.

sent since the simulations represent scenarios and not events that have already transpired. Based on simulation, ten out of the 32 alluvial fans generated debris flows. The cities and municipalities that may be affected by debris flow hazards, indicated in red, are Albuera, Baybay, Burauen, Capoocan, and La Paz. Alluvial fans without simulated debris flow might be due to geomorphological setting. Gentler slopes, compared to steeper slopes, does not have the capacity to accommodate debris flows. Furthermore, the parameters used during simulations might be inappropriate for such setting.

Figure 15: Combined map of potential rock slide zone sources and propagation extent for Southern Leyte

3.4. Alluvial Fan and Debris Flow Hazards Thirty-two alluvial fans were identified in Leyte Province. The widest fan identified was located in Brgy. Taghuyan, Burauwen and covered a total area of 24.66 sq km. This was followed by the fan identified in Brgy Culasian, Capoocan with an area of 14.96 sq km, and by the Ormoc Bajada identified in Ormoc City and Albuera covering 14.54 sq km. Other fans were located in the municipalities of Babatngon, Tacloban City, Capoocan, Leyte, Ormoc City, Merida, Isabel, Albuera, Baybay, Abuyog, Lapaz, Burauwen, and Julita. The alluvial fan shapefiles together with the watershed extracted from a DEM, and the 100-year return period rainfall data of Leyte, were used to simulate debris flows from the selected catchments using FLO-2D (Figure 16). The resulting debris flow hazard map indicated similar risks in the same areas. While the debris overlaps the alluvial fan, it is safer to assume a debris flow will most likely not travel beyond the boundaries of the alluvial fan. The simulation used a digital elevation model with the lowest resolution as higher resolution images were still undergoing simulation. The map indicates that in the event of a debris flow, the areas covered in red will most probably be buried in 1 meter or above (at least 3 feet) of debris. Areas in orange pose a lower threat, with a debris depth of 0.2 meter to 1 meter (from 7 inches to below 3 feet). It is also important to note that the results of the simulation are continuously being validated in the field. Field assessment is also necessary to ensure the accuracy of the debris flow boundaries. However, the presence of debris flows can be ab-

4. Discussion The inherent complexity in the mobilization mechanisms of landslides limits the scope of what different numerical models can capture. For instance, SINMAP gives a good approximation of the overall extent of possible landslide hazard areas in a given area. However the model only takes into account susceptibility to shallow translational slides. Consequently, many of the landslides in the inventory fall in the moderate susceptibility level. Matterocking, on the other hand, targets only susceptibility induced by the presence of geologic structures, in this case planar and wedge failures. Susceptibility to shallow translational slides are thus not taken into account in the results. This is reflected as well when we cross-reference the landslide points from the inventory with the results from Matterocking. Numerous points fall outside the hazard areas identified by Matterocking. These are most likely shallow translational slides caused by soil instability. Debris flows have significantly different mechanism than slides. In addition, these can occur in relatively lower slopes than slides, often, but not always, in alluvial fans. It must be 19

number of landslides that fell within each susceptibility class (High, Moderate, Low). Among the 45 landslide points in Leyte, none were within the stable zones. The majority of the points fell within the medium (55%) and high susceptibility (45%) zones. Using a model that identifies landslide susceptible areas with high accuracy and at the same time maximizing safe and habitable areas is beneficial to the development plans of the local government. This is critical especially for the rehabilitation work for Haiyan-affected communities because it is essential that reconstruction be carefully planned to prevent future disasters in the area. The landslide maps that were provided are fit to the requirements of the rehabilitation process. It should be further improved for better accuracy by conducting more detailed geological and geotechnical assessments, and using more accurate LiDAR topographic data.

Figure 17: Unstable Slopes Map of Leyte on Project NOAH website

5. Conclusion The results derived from the high-resolution shallow landslide, debris flow, alluvial fan and deep-seated landslide models enabled the identification of areas that are safe from landslide hazards. Landslide hazard maps, when combined with other hazard maps such as flood and storm surge maps, can be used to identify locations that are favorable for development. Critical facilities such as hospitals and evacuation centers should be situated in areas that are accessible and least likely to be affected by hazards. The use of these maps can facilitate the visualization and deeper understanding of possible disaster scenarios and serve as a guide in hazard preparedness and mitigation.

Figure 18: Landslide Hazard Map of Leyte on Project NOAH website

noted, however, that there are numerous hazards that can occur within alluvial fans, including debris flows, hyper-concentrated flows, and sheet and flash flooding. To address these limitations, we combine the results from SINMAP, Matterocking, and Conefall to better depict the level of hazard present in the area. This takes advantage the strengths of each model, at the same time complementing what the others cannot capture. For alluvial fans, we merge these with the debris flow hazard maps to visualize where in the alluvial fans are most exposed to debris flow hazards. The results of SINMAP were combined with Matterocking and Conefall to create generalized landslide susceptibility maps. The SINMAP-Matterocking combination was labeled “Unstable Slopes Map” whereas the SINMAP-Conefall combination was labeled “Landslide Hazard Map.” The Unstable Slopes Map depicted areas that may potentially fail given extreme circumstances (Figure 17). By calculating the potential runout of the slope failures through Conefall, we included another aspect of landslide hazard mapping that cannot be captured by mapping of unstable areas. The areas covered by the runout were automatically classified as “high hazard” on the map (Figure 18). These data were then uploaded to the publicly accessible Project NOAH portal (http://beta.noah.dost.gov.ph) and where they could be overlain on different basemaps and zoomed down to the community level. To validate the accuracy of the maps, we compared them with the landslides identified in the landslide inventory. The accuracy of the SINMAP simulation was verified by overlaying the landslide points over the model and counting the 20

References

Jurgan, H., 1980. Geology and hydrocarbon prospects of the Visayan Basin. Appendix III: Younger Tertiary carbonate rocks. Kimura, T., Tokuyama, A., Gonzales, B., and Andal, R. A., 1968. Geological Structures in Tayabas Isthmu District, Philippines. Continental Geology and Paleontology Southeast Asia, , 49, 156–178. Lagmay, A. M. F. A., Tengonciang, A. M. P., Marcos, H. V., and Pascua, C. S., 2003. A structural model guide for geothermal exploration in Ancestral Mount Bao, Leyte, Philippines. 122, 133–141. Lin, C., Tseng, C., Tseng, Y., Fei, L., Hsieh, Y., and Tarolli, P., 2013. Recognition of large scale deep-seated landslides in forest areas of Taiwan using high resolution topography. Journal of Asian Earth Sciences, , 62, 389–400. MGB, 2004. Geology of the Philippines, volume 1. Mines and Geosciences Bureau. MMAJ-JICA, 1990. Mineral deposits and tectonics of two contrasting geologic environments. Technical report, MMAJ-JICA. Terminal Report. Nakata, T., Sangawa, A., and Hirano, S. C., 1977. A report on tectonic Lanlandf along the Philippine Fault Zone in Northern Luzon, Philippines. In Nakata, ed., Geography, volume 27 of 7, 69–93. Sendai, Japan. Science Reports, Tohoku University. NDRRMC, 2013. NDRRMC Update: SitRep No. 31 Effects of typhoon ”Yolanda” (Haiyan). Technical report, National Disaster Risk Reduction and Management Council. Pack, R., Tarboton, D., and Goodwin, C. N., 1998. The SINMAP Approach to Terrain Stability Mapping. 8th Congress of the International Association of Engineering Geology. PHIVOLCS, 2008. Active and Inactive Volcanoes. Website. URL http://www.phivolcs.dost.gov.ph. Pilac, J., 1965. Geology of Northern Leyte. Bureau of Mines. Pinet, N. and Stephan, J. F., 1990. The Philippine wrench fault system in the Ilocos Foothills, Northwestern Luzon, Philippines. Tectonophysics, 183, 14, 207–224. Rickenmann, D. and Zimmermann, M., 1993. The 1987 debris flows in Switzerland: documentation and analysis. Geomorphology, , 8, 175–189. Rutland, R. W. R., 1968. A tectonic study of the Philippine Fault Zone. Geological Society of London, , 123, 295–325. Takahashi, T., 1981. Estimation of potential debris flows and their hazardous zones: Soft countermeasures for a disaster. Journal of Natural Disaster Science, , 3, 57–89. Wyllie, D. and Mah, C., 2004. Rock Slope Engineering: Civil and Mining. Spon Press, 4th edition.

Abbassi, M. and Shabanian, E., 1999. Evolution of the stress field in Tehran region during the Quarternary. 3rd International Conference on Seismology and Earthquake Engineering (SEE3). Aguilar, Y. O. M., 1995. Neogene Stratigraphy and paleontology of the Calubin Peninsula, NW Leyte, Philippines. Master’s thesis, University of Tsukuba, Japan. Allen, C. R., 1962. Circum-Pacific faulting in the Philippines-Taiwan region. Journal of Geophysical Research. Aurelio, M., 1992. Tectonique du segment central de la faille philippine: etude structurale, cinematique et evolution geodynamique.. Ph.D. thesis, Universite Pierre et Marie Curie, Academie de Paris, Universite de Paris. Aurelio, M. A., Huchon, P., Barrier, E., and Gaulon, E., 1994. Displacement rates along the Philippine Fault estimated from slip-vectors and regional kinematics. Geological Society of the Philippines, 2, 49, 65–77. Baillifard, F., Jaboyedoff, M., and Sartori, M., 2003. Rockfall hazard mapping along a mountainous road in Switzerland using a GIS-based parameter rating approach. Natural Hazards and Earth System Science, , 3 (5), 435–442. Balce, C. and Cabantog, A., 1998. Tacloban Ophiolite Massif - Basement of northeast Leyte physiographic province. Barcelona, B. M., 1981. The nature of the faults in the Philippine fault zone. In Barcelona, ed., Geology of Southeast Asia, 267–275. Proc. 4th Regional Conference on Geology. Barrera, A., Aristorenas, I., and Tiongson, J. A., 1954. Soil Survey of Leyte Province. Department of Agriculture and Natural Resources Bureau of Soils. Cabantog, A. V. and Escalada, P. P., 1989. Geology of Tacloban Ophiolite. Technical report, Bureau of Mines and Geosciences Bureau. Claessens, L., Heuvelink, G. B. M., Schoorl, J. M., and Veldkamp, A., 2005. DEM resolution effects on shallow landslide hazard and soil redistribution modelling. Earth Surf. Proc. Land., , 30, 461–477. Cole, J., McCabe, R., Moriarty, T., Malicse, J., Delfin, F., Tebar, H., and Ferrer, H., 1989. A Preliminary Neogene Paleomagnetic data set from Leyte and its relation to motion on the Philippine fault. Technophys, , 168, 205–221. Corby, G., Kleinpell, R., Popenoe, W., Merchant, R., William, H., Teves, J., Grey, R., Daleon, B., Mamaclay, F., Villongco, A., Herrera, M., Guillen, J., Hollister, J., Johnson, H., Billings, M., Fryxell, E., Taylor, E., Nelson, C., Birch, D., Reed, R., and Marquez, R., 1951a. Geology and Oil Possibilities of the Philippines, volume 1. Bureau of Mines. Corby, G., Kleinpell, R., Popenoe, W., Merchant, R., William, H., Teves, J., Grey, R., Daleon, B., Mamaclay, F., Villongco, A., Herrera, M., Guillen, J., Hollister, J., Johnson, H., Billings, M., Fryxell, E., Taylor, E., Nelson, C., Birch, D., Reed, R., and Marquez, R., 1951b. Geology and Oil Possibilities of the Philippines, volume 1. Bureau of Mines. David, C., Racoma, B., Gonzales, J., and Clutario, M., 2013. A Manifestation of Climate Change? A Look at Typhoon Yolanda in Relation to the Historical Tropical Cyclone Archive. Science Diliman, 25, 2. Duquesnoy, B., 1997. Contributions de la geodesie a‘ letude de grands decrochements actifs associes a‘ des zones de subduction a‘ convergence oblique. Ph.D. thesis, University of Paris XI, Orsay. Florendo, F., 1984. Geologic map of Maasin Peninsula, Leyte Island. Regional office, 7. Florendo, F., 1987. The tectonic framework and the Cretaceous to Cenozoic evolution of the east-central Philippines. Master’s thesis, University of Tulsa Oklahoma. Guzzetti, F., Mondini, A., Cardinali, M., Fiorucci, F., Santangelo, M., and Chang, K., 2012. Landslide inventory maps: New tools for an old problem. Earth-Science Reviews, , 112, 42–66. Hammond, C., Hall, D., Miller, S., and Swetik, P., 1992. Level I Stability Analysis (LISA) Documentation for Version 2.0. Technical report, USDA Forest Service Intermountain Research Station. General Technical Report INT-285. Holmgren, P., 1994. Multiple flow direction algorithms for runoff modellingin grid based elevation models: An empirical evaluation. Hydrol. Process., , 8, 327–334. Horton, P., Jaboyedoff, M., Rudaz, B., and Zimmermann, M., 2013. Flow-R, a model for susceptibility mapping of debris flows and other gravitational hazards at a regional scale. Natural Hazards and Earth System Sciences, , 13, 869–885. IPCC, 2007. In Climate Change 2007: The Physical Science Basis: Contribution of Working Group 1 to the IPCC Fourth Assessment Report. Technical report, Intergovernmental Panel on Climate Change.

21