Anatomy of the Naga City Landslide and Comparison With Historical Debris Avalanches and Analog Models

Lagmay, Alfredo Mahar Francisco; Escape, Carmille Marie; Ybañez, Audrei Anne; Suarez, John Kenneth; Cuaresma, Genaro

doi:10.3389/feart.2020.00312

ORIGINAL RESEARCH article

Front. Earth Sci., 11 August 2020

Sec. Structural Geology and Tectonics

Volume 8 - 2020 | https://doi.org/10.3389/feart.2020.00312

This article is part of the Research Topic Rock Avalanches View all 14 articles

Anatomy of the Naga City Landslide and Comparison With Historical Debris Avalanches and Analog Models

$\nAlfredo Mahar Francisco Lagmay,$ Alfredo Mahar Francisco Lagmay^1,2^*

Carmille Marie Escape^1,2

Audrei Anne Ybañez¹

John Kenneth Suarez¹

Genaro Cuaresma¹

¹UP Resilience Institute, University of the Philippines, Quezon City, Philippines
²National Institute of Geological Sciences, University of the Philippines, Quezon City, Philippines

Debris avalanches pose some of the most destructive geologic hazards that threaten both urban and rural populations around the world. On 20 September 2018, villages in Naga City, Cebu, Philippines, were devastated by a landslide that claimed 78 lives with 6 missing, joining other catastrophic landslides in the country like the 1628 Iriga and the 2006 Guinsaugon debris avalanches. Understanding the mechanism of these gargantuan landslides and their correct nomenclature are useful for hazard prevention and mitigation. In this study, we compare the deposit characteristics of the Naga City landslide with analog models and well-known historical debris avalanche events/deposits in the Philippines to understand factors that led to the landslide disaster in Naga City. Physical characteristics obtained from aerial and satellite imagery, ground surveys, recorded footage, borehole data, and lithologic maps provided a detailed dataset for analyzing the conditions that led to the mass movement and the observed characteristics of the Naga landslide deposits. Comparison with analog models of hummock formation and the description of historical debris avalanche deposits show striking similarities, which were used to demonstrate that the Naga landslide was a Rockslide-Debris Avalanche. The equations of Corominas (1996) and Dade and Huppert (1998) for long-runout rockfalls support this analysis. The Naga landslide event is an example of a well-documented debris avalanche, complete with all the characteristics of this type of rapid mass movement. It is consistent with the descriptions found in the literature with respect to its deposit features and mechanical behavior as defined by laboratory models and empirically-derived equations. This study helps us understand historical and future long-runout debris avalanches in order for scientists and authorities to find ways to save lives. Unfortunately, there was lack of appropriate hazards assessment on the site, which had warnings in the form of the development of fractures at the headscarp of the landslide, a month prior to the disaster.

1. Introduction

On 20 September 2018, a massive landslide devastated Naga City, Cebu. The Naga City Landslide claimed the lives of 78 villagers and injured 18 while six people remained missing and are presumed dead. The majority of the fatalities were recovered at the landslide toe about 1.2 km from the 200 meter-high headscarp. This was due to the unexpectedly large landslide volume and unusually long-runout, which surprised villagers at the landslide toe (Figure 1). Despite early warnings from developing fractures near the headscarp a month prior to the disaster, no action was taken by villagers at the landslide toe in Sitio Sindulan, Barangay Tinaan and those in an adjacent area at the base of a cliff in Barangay Naalad (Figure 1). Only residents in Sitio Tagaytay, situated near the headscarp of the landslide, were evacuated by authorities the night before the mountain collapsed.

FIGURE 1

Figure 1. Aerial image of the Naga landslide taken a month after the incident. Numerous houses in Sitio Sindulan, Barangay Tinaan and the village situated near the base of the cliff in Barangay Naalad were buried by the landslide debris.

Discernible from the closed-circuit television (CCTV) footage that captured the event (see Supplementary Video 1, BJMP-NAGA, 2018) was mass movement that started as a translational slide. It also showed a landslide mass that initially moved as one intact block with no significant backward regression. The landslide lasted for a minute and traveled with a maximum velocity of 72 km/h. Fractures observed a month prior to the catastrophe (Lagmay, 2018; MGB, 2018; Catane et al., 2019) suggest that the sliding block started to move some time before the main collapse event, which happened in the early morning at 5:45 a.m. on 20 September 2018. The sliding mass moving as a whole, meant a larger volume was available, favoring increased runout of the landslide. Once in motion, the front of the sliding mass accelerated, stretching the limestone body to create a long-runout (L) relative to the collapse height (H) with a calculated H/L ratio of 0.17, which is a characteristic value of volcanic and non-volcanic debris avalanches (Ui, 1983; Siebert, 1984; Ui et al., 1986; de Vries and Delcamp, 2015).

Debris avalanches are catastrophic, large scale, mass wasting events with a fast moving body that can travel a long way relative to the collapse height, regardless of whether they are volcanic or non-volcanic in origin (de Vries and Delcamp, 2015). Though more commonly associated with volcanic mass wasting events, the use of the term debris avalanche also applies to non-volcanic landslides if they meet the criteria of having an amphitheater, hummocks, torevas, megablocks, jigsaw puzzle features (Ui, 1983) and a longer-runout compared to landslides. Debris avalanche deposits, irrespective of whether volcanic or non-volcanic in origin, exhibit similar characteristics (Ui, 1989). Well-known examples for volcanic debris avalanches are Mt. St Helens (Glicken, 1996), Jocotitlan (Siebe et al., 1992), Bezymiannyi (Belousov and Belousova, 1998) and Shiveluch (Ponomareva et al., 1998) volcanoes, whereas examples for non-volcanic debris avalanches are Guinsaugon (Lagmay et al., 2008; Futalan et al., 2010), Mt. Meager (Guthrie et al., 2012; Roberti et al., 2017), Blackhawk (Johnson, 1978; Ui et al., 2000), Sherman (McSaveney, 1978; Ui et al., 2000) and Luanshibao (Wang et al., 2018). Ui (1983) also demonstrated the similarity of the H/L ratio for volcanic and non-volcanic events and proposed that both types have similar mode of transportation related to gravitational sliding due to slope instability. Although a debris flow is another type of mass wasting event that has a long-runout, it can be differentiated from debris avalanches in the sense that their mobility is primarily controlled by the presence of water. The Naga City landslide was relatively dry and did not form debris flows at the distal portion of the deposits. Due to the morphological characteristics of the deposit field, the long-runout, the absence of debris flow deposits, and the initial slide movement of rock units, the Naga City landslide was classified as a rockslide-debris avalanche. Hungr et al. (2001) suggested that the distinction between rock avalanches and debris avalanches is gradational and subjective especially because there is a certain difficulty in separating welded materials (rock) from uncemented granular deposits (debris). The Carcar and Barili formations which comprised the collapsed Naga landslide material, are largely uncemented, soft to semi-hard and can crumble. The collapse started as a rockslide, but because of the extension and acceleration that leads to fragmentation, forms the block facies, matrix facies and other features that are distinctive of a debris avalanche, hence the classification.

Due to the high number of fatalities associated with the long reach of the landslide, it is important to analyze the conditions that contributed to the unusual runout of this event. To do this, we compared the Naga City disaster with two other known local debris avalanche events: The 1628 Iriga volcanic debris avalanche (Aguila et al., 1986; Paguican, 2012; Minimo and Lagmay, 2016) and the 2006 Guinsaugon debris avalanche (Lagmay et al., 2006; Evans et al., 2007; Catane et al., 2008). In this paper, we revisit the analysis of the landslide event and present new detailed field investigation. We explore the relationship of the volume with the runout length of the Naga landslide and compare the deposits with the worldwide dataset of debris avalanches (Corominas, 1996; Dade and Huppert, 1998) and existing analog models (Paguican et al., 2014). This is used to determine the emplacement mechanism of the landslide, which led to the fatal disaster in Naga City.

2. Geologic Setting

The City of Naga is located at the southeastern coast of Cebu Islands, Philippines. Cebu is part of an island group, which along with Panay, Negros, Bohol, Leyte, and Samar comprise the Central Philippines region that has a common geologic history (Deng et al., 2015). These islands, including their adjacent sedimentary basins, namely Iloilo and Visayan basins (Aurelio and Peña, 2002), are underlain by a basement complex composed of Cretaceous to Eocene igneous and metamorphic rocks and Cenozoic volcanic and sedimentary units (Santos-Ynigo, 1951; Dimalanta et al., 2006; Deng et al., 2015). The region lies within the Philippine Mobile Belt (PMB), a deforming and seismically active zone, which is bound by subduction zones of opposite polarities (Gervasio, 1967; Lagmay et al., 2009). West of Panay is the east-dipping Early to Middle Miocene Negros Trench, whereas east of Samar is the west-dipping Pliocene Philippine Trench. The Philippine Fault Zone traverses this part of the Philippine Archipelago along Leyte (Allen, 1962) (Figure 2).

FIGURE 2

Figure 2. The PMB is found within the convergent zone between the Philippine Sea Plate and the Eurasian Plate (adapted from Lagmay et al., 2009). The Visayas Region contains the Iloilo Basin and the Visayas Basin, wherein Cebu is located (adapted from Dimalanta et al., 2006). The City of Naga is located in central Cebu which is underlain by a basement complex of igneous, volcanic, sedimentary, and metamorphic rocks. Nearby is the active NE-SW trending Central Cebu Fault (PHIVOLCS) as well as numerous unnamed faults (MGB, 1983).

The oldest rock formation in Cebu is the Jurassic Tunglob Schist. It is overlain by the Cretaceous to Paleocene age Mananga Group, which consists of limestone, clastic sedimentary rocks, andesitic to basaltic pyroclastics and lava, calcareous mudstone, conglomerate and sandstone (Aurelio and Peña, 2010). Unconformably overlying the Mananga Group are a series of unconformable sedimentary and mostly calcareous formations that range in age from Late Eocene to Plio-Pleistocene age. Intruding into the Mananga Group at places is the Lutopan Diorite. In other areas of Cebu, the late Miocene Bulacao Andesite occurs as intrusive breccia and extrusive deposits of porphyritic andesite. Serpentinized ultramafic and mafic rocks occur as diapiric intrusions along major faults that cut across Cebu (Balce, 1977; Aurelio and Peña, 2010). The two youngest formations in the area of Cebu where the landslide took place and the subject of interest in this study, are the Late Miocene to Early Pliocene Barili and Plio-Pleistocene Carcar Formations, which are both calcareous in composition (Corby, 1951; Aurelio and Peña, 2010).

The nearest identified potentially active fault system from the Naga landslide is the Cebu Fault System, which is a northeast-trending fault system composed of two major structures: The Central and the South Cebu faults (PHIVOLCS, 2016). The Central Cebu Fault passes through Naga City and is located 5.5 km west of the landslide. Other unnamed faults that were previously mapped (MGB, 1983) are consistent with the northeast trend of the Central Cebu Fault. Two earthquakes near the Naga Landslide area were recorded in 2018 by the Philippine Seismic Network with magnitudes 3.0 and 3.4, respectively. Both of these earthquakes were less than 33 km deep with the epicenters located within 3 km of the Naga landslide deposits (Figure 2).

3. Methodology

Satellite data from Planetscope and other aerial images were used to analyze the pre-event (8 September 2018) and post-event (21 September 2018) conditions of the Naga landslide. These orthorectified images contained 4 multispectral bands (blue, green, red, near-infrared), with a resolution of 3 m. The pre- and post-event satellite images were compared to identify the extent of the landslide deposit field. The change analysis was done a day after the landslide event and was used during the search and rescue phase of the disaster to identify buried houses at the distal portion of the debris field (GMA News and Public Affairs, 2018).

For better analysis, satellite images were augmented by crowd-sourced drone photo and video footages. The drone images were processed using digital photogrammetry to generate Digital Elevation Models (DEMs). In succeeding field surveys, a DJI Mavic 2 Pro drone equipped with a Hasselblad L1D-20c camera with a field of view (FOV) of about 77 deg, aperture of f/2.8-f/11 and shooting range of a minimum of 1 m was used to fill in gaps of the initial DEMs. The latter drone surveys were focused in the headscarp area, which was extremely difficult to access due to dangerous and harsh terrain created by the landslide event. Point clouds were created and transformed into a DEM and combined with the crowd-sourced data to create a post-event DEM of the entire landslide area. A 1 × 1 m Lidar pre-event DEM was used for change analysis of the topography. The pre- and post-event DEMs, with the same resolution, were used for the volume calculation.

Field data collection was conducted from November to December 2018, and in January and June 2019 to investigate the landslide area and vicinities. Lithologies, geological structures and morphology of the area were mapped. The consecutive field surveys were conducted to characterize the landslide deposit in detail.

Drill core data collected in the year 2010 was obtained as a supplement to field activities and were used for the generation of a 3-dimensional geologic model of the area. Contact relationships derived from the core data were used as the guiding base for 3D modeling in the Leapfrog Geo software. Core logs were assimilated across a directed core line which utilizes an implicit model to define boundaries between stratigraphic units and geologic structures. These surfaces were generated using a triangulated irregular network in between known data points and were projected along the lithologic boundaries.

The pre- and post-landslide high-resolution DEMs were used to compute the traveled horizontal distance (L) over a vertical height difference (H) to determine the angle of reach. Known as the Fahrböschung or Heim ratio (Heim, 1932), the reach angle was calculated from the high-resolution DEMs to demonstrate the efficiency of landslide motion. According to the analysis of Corominas (1996) on a global landslide dataset, movements showing the lowest angles of reach attain the farthest horizontal distance in relation to fall-height of the landslide. The mobility plot of Corominas (1996) showing H/L vs. volume was used in the analysis to determine the regression limits and confidence interval in the classification of the Naga landslide. The angle of reach was computed based on Equation (1):

\begin{array}{l} log (\frac{H}{L}) = - 0.105 log v o l - 0.012 r^{2} = 0.763 & (1) \end{array}

where H is the vertical height difference and L is the horizontal projection of the distance. To be classified as a debris avalanche with a 95% mean confidence interval, the regression equation of Corominas (1996) requires a range of limits of −0.8607 and −0.6419, whereas for a translational slide, the range is from −0.7302 to −0.5454.

We also used the formula of Dade and Huppert (1998) to determine the area overrun by the landslide which has been demonstrated to be proportional to the potential energy of the debris mass. The long-runout scaling was computed using Equation (2):

\begin{array}{l} A = {(λ)}^{\frac{1}{3}} {(\frac{g M H}{τ})}^{\frac{2}{3}} & (2) \end{array}

where A is the area overrun by the landslide, ${(g M H)}^{\frac{2}{3}}$ is the potential energy of the debris mass before failure, τ is the resisting shear stress and λ is the geometry parameter of the landslide. According to Dade and Huppert (1998), the magnitude of shear stress of resistance (τ) and ratio of landslide width and length (λ) limit the runout extent and is proportional to the area of the landslide footprint.

However, the relationship between the Fahrböschung angle and landslide volume is not as straightforward as it seems (Lucas et al., 2014). Therefore, the Naga landslide was further examined relative to two well-known debris avalanche deposits in the Philippines. These are the 1628 Iriga (Aguila et al., 1986; Paguican et al., 2012; Minimo and Lagmay, 2016) and the 2006 Guinsaugon (Lagmay et al., 2006; Evans et al., 2007; Catane et al., 2008) landslides, which also had long reaches relative to their collapse height. Additionally, a recent study of the 2018 Naga landslide classifying it as a low-angle translational block slide (Catane et al., 2019) was included in the comparative analysis. Landslide parameters used in the comparative analysis include debris volume, area covered, vertical height, and horizontal distance for each of these landslide events as described in the literature. Where other parameters are absent (i.e., Fahrböschung angle), these were measured in maps or computed using given associated values included in the respective publications. Through calculations using the Corominas (1996) and Dade and Huppert (1998) equations, the deposits of the three landslides were compared to determine their similarities and differences, if any. The comparative analysis was made to better understand the Naga landslide and its nomenclature. The presence of geomorphic and structural features characteristic of debris avalanche deposits which includes: (1) megablock structures; (2) jigsaw puzzle effects; (3) hummocks; and (4) an amphitheater at the source (Ui, 1983; Siebert, 1984; Andrade and de Vries, 2010; Davies et al., 2010; de Vries and Delcamp, 2015), was critical in the analysis.

Scaled analog models that investigated geomorphic features, in particular hummock formation, were also used to characterize the deposit features observed in the Naga landslide deposit and to interpret their formation. For example, hummocks, which were demonstrated through laboratory models to form by extension of large blocks during the collapse event were compared to the pinnacle hummocks pervasive in the Naga landslide debris field. Analog model structures associated with hummock formation, such as normal faults, horsts and grabens were compared with those found in the Naga landslide deposits.

The Naga City landslide's morphology, geology, structures, runout behavior (according to Corominas, 1996; Dade and Huppert, 1998), comparison with known debris avalanches of the Philippines and analog models for large-scale volcanic and non-volcanic collapses of Paguican et al. (2014) were then integrated in the analysis. This was done to gain insights on the kinematics and dynamics of debris avalanches and advance our understanding on their long-runout behavior to prevent or mitigate their impacts in future events.

4. Results

4.1. Morphology

Planetscope satellite images taken a day after the landslide event reveal a striking land cover change (Figures 3A,B). The areas that exhibit the biggest difference in surface conditions are those within the scarp area at the elevated regions of the quarry site and the two lobes of the landslide deposit in Barangay Naalad and Sitio Sindulan, Barangay Tinaan (Figure 3). The region near the headscarp was vegetated as clearly seen from the pre-disaster satellite image (Figure 3A). After the collapse and mass movement, the underlying limestone was largely exposed and is seen in the post-disaster satellite imagery with high albedo. At the distal end of the landslide deposit in Barangays Tinaan and Naalad, 37 houses were buried in debris (Figure 3C). These residential areas at the toe of the landslide were buried up to 10 m as measured in the field and the DEMs. The Naga River was dammed as well. Quarried areas that were covered by the September 2018 Naga landslide deposit event did not show any significant change in NDVI values in pre- and post-disaster Planetscope imageries.

FIGURE 3

Figure 3. Planetscope satellite images of the Naga landslide area (A) before and (B) after the Naga disaster (C) close up of the Lobe 1 with building footprints (OpenStreetMap contributors, 2017). There were 37 houses buried in Sublobes A and B. (D) Pre-event cross-section from Google Earth show the “topographic high” that disrupted/diverged the flow of the materials which caused the formation of Lobe 1 and Lobe 2. A natural cliff also caused the formation of Sublobe A and B. Another notable feature of the event is the (E) remnants of the road that collapsed and formed a seemingly linear feature.

There are two lobes at the distal end of the landslide debris field. The first is located in the northern side (Lobe 1) whereas the second is in the southern side (Lobe 2). In between these lobes, is a topographic high that diverted flow toward two directions (Figures 3B,D). The topographic high has a maximum pre-event elevation of 120 masl whereas the areas which eventually became Lobe 1 and Lobe 2 only have maximum pre-event elevations of 55 and 80 masl, respectively. This elevated portion acted as a barrier along the landslide path and prevented the axial part of the landslide to equally spread further downslope. The diversion of flow effectively shortened the landslide runout to 830 m in the axial portion and caused the formation of a two-lobed landslide deposit field. The pre-event, frontal plane cross-section of the landslide area (Figure 3D) shows the topographic high that caused the diversion of materials. Furthermore, a natural cliff with an elevation of 60 masl dropping to 30 masl eventually caused the formation of sublobes A and B (Figure 3C). The collapse of this cliff and the overflow of landslide material from the top section can clearly be seen from the CCTV (see Supplementary Video 1, BJMP-NAGA, 2018). It may look like a waterfall or some dewatering has happened, but the deposits are dry and debris flows did not form. A significant amount of dust clouds was seen generated during mass movement, indicating relatively dry material. Debris flows would have also formed if there was a significant amount of water present.

The post-event DEM used for detailed analysis of the Naga Landslide has a resolution of 33 × 33 cm. This post-event DEM was resampled to match the 1 × 1 m resolution of the pre-event LiDAR DEM to calculate the 11,000,000 m³ volume of the Naga landslide. Cross-sections derived from the DEM, show prominent hummocks and rotated toreva blocks (Figure 4). Toreva blocks are commonly composed of one or more blocks that slide and retain the original stratigraphic sequence, whereas hummocks are mound features composed of block material in the surface with a form of conical shape and have a height of tens to hundreds of meters (Ui, 1983; Stoopes and Sheridan, 1992). Hummock sizes in the Naga deposit field decrease in dimension away from the headscarp, consistent with the descriptions in other debris avalanche deposits (Reiche, 1937; Ui, 1983; Crandell et al., 1989; Thompson et al., 2010; de Vries and Delcamp, 2015). The high-resolution DEM also reveals that the highest elevation of the failure is at 255 masl, whereas the elevation at the toe in Lobes 1 and 2 is 50 and 70 masl, respectively. In terms of maximum distance traveled from the headscarp, the landslide is measured at 1.1 km at Lobe 1 and 1.2 km at Lobe 2.

FIGURE 4

Figure 4. Cross sections of the Naga landslide based on a DEM generated from drone images.

Noticeable in the debris field are large tilted blocks found at the medial to proximal section of the 1.2 km landslide (Figure 5). The tilted blocks or toreva blocks manifest as rows of broken material with their long axis perpendicular to the direction of the landslide. The largest blocks within the debris field have widths ranging from 120 to 350 m and traveled up to 220 m. Careful analysis of a linear feature (Figures 3B,E) that appears to divide the torevas in the proximal collapse zone of the landslide, reveal a road that subsided and whose parts were preserved as a large, intact landslide block (Figure 3E). The torevas at both sides of the road are found to be contiguous upon closer examination and indicate a single collapse event.

FIGURE 5

Figure 5. (A) View looking north. Rotated toreva blocks at the proximal portion of the debris field (UP-CENVI, 2018). (B) View looking south-southeast. Horst and Graben structures and pinnacle hummocks formed due to faulting as the slide was extending.

Pinnacle hummocks were also identified at the medial section of the Naga debris field. The sizes of the hummocks vary but are as large as 15 m in height and 10 m in width at its base. These hummocks appear as the remains of highly stretched blocks that have developed normal fault structures with horst and graben structures (Figure 5). High-standing hummocks are horst structures with adjacent grabens separated from them.

The distal portion of the deposit field is comprised of smaller blocks compared to the proximal and distal portions (Figure 6). The range of sizes of the blocks are 5 to 15 m and are surrounded by non-graded angular finer-sized fragments. These fragments vary from clay (<1/256 mm) to larger than gravel sizes (2–64 mm). In Lobe 1, the largest blocks are up to 6 m in size whereas in Lobe 2, blocks are up to 15 m.

FIGURE 6

Figure 6. (A) Sandstone/siltstone in the debris field of the landslide. (B) A drill hole mark is evident in a boulder. (C) Sandstone and limestone deposit. (D) The contact between the limestone and the siltsone is clearly seen in a boulder in the debris field.

Large mobilized intact blocks (MIB) are present in the proximal and medial portion of the landslide. Relatively smaller blocks compared to the MIBs, but still up to several meters in size, are scattered in the proximal and medial portions of the debris field with their sizes generally decreasing away from the headscarp. These blocks, no matter their size, commonly exhibit jigsaw cracks and are surrounded by a finer-grained matrix composed of fragmented limestone (Figure 7).

FIGURE 7

Figure 7. (A–F) Boulders and blocks with jigsaw cracks.

Other notable morphological features include; a linear headscarp that extends for about 734 m and largely intact MIBs that had been translated over an average distance of 225 m without being overturned. Houses and trees on top of these large blocks remain standing.

4.2. Geology

The landslide area is composed of two geological formations as seen from the exposed amphitheater walls of the collapsed scarp. The lower part is the Barili Formation of Late Miocene to Early Pliocene age (Mines and Geosciences Bureau, 1981), which consists of a lower limestone member and upper marl member. The limestone is hard, light brown, coralline, locally porous or sandy, richly fossiliferous, whereas the marl is generally brown, slightly sandy poorly bedded and, fossiliferous with thin limestone interbeds (Del Rosario et al., 2005). In general, the bedded outcrops seen at the fringes of the landslide were classified as part of this formation. The Barili Formation is overlain by the poorly bedded to massive Plio-Pleistocene Carcar Formation, which according to the literature is composed of more coralline limestone and partly dolomitic.

The exposed calcareous formations within the landslide zone are more than 50 m thick and are underlain by a 3 m-thick sandstone/siltstone bed. This sandstone/siltstone bed, which is most likely the Marl component of the Barili Formation, is distinguished from the buff-colored limestone and appears in the quarry site of the Apo Tenement at 60 m elevation. Boulders of siltstone are also found in some areas of the debris field. The contact between the limestone and the sandstone/siltstone is clearly seen in some of the boulders within the avalanche deposit (Figure 6D).

4.3. Core Data

From April to May 2010, 12 boreholes within the Mineral Production Sharing Agreement (MPSA) 286 and 137 tenements of the Apo Land and Quarry Corporation were drilled and logged. Out of these, nine boreholes within the landslide area were used in this study. The elevations of these boreholes range from −0.5 to 242.7 masl with a drill depth of 40 to 121.5 m. The boreholes show limestone as the thickest and topmost lithology along the length of each core topped by a thin layer of topsoil. The limestone core logs are described as ranging from hard to soft with predominance of semi-hard to soft limestone descriptions. Sandstone with thickness ranging from 1.2 to 42.6 m underlie the limestone. Other lithologies described in the corelogs include lime, black clay, black shale, blackstone, gray and black pozzolan, green sands and soft clay. The Leapfrog plots reveal dipping beds with limestone above a less coherent layer of sandstone (Figure 8). The average dip angle of these layers is 7.65 deg with an average dip direction of N87.16°E toward the sea and in the direction of the landslide. This dipping bed and deposition plane between the sandstone and limestone beds is identified as the sliding plane where the landslide moved.

FIGURE 8

Figure 8. Nine borehole data were modeled to show the lithology of the study area. A dipping layer of limestone lays on top of a sandstone layer.

4.4. Structural Geology of the Landslide

The headscarp of the landslide is planar. Numerous measurements in different areas from top to bottom and north to south of this approximately 70 m-high and 734 m-wide planar structure show a northeast strike direction and a dip angle in the range of 58–90 deg. Slickenlines and slickensides (Figure 9) are also present on this plane with striations that have a rake angle of 69.9 deg. Striations with similar slickenlines are observed everywhere in this planar landslide head wall.

FIGURE 9

Figure 9. A N66°E fault plane at the landslide head showing slickensides and slickenlines with a 69.9° rake angle. Oblique slip movement is suggested by the direction and orientation of the slickensides.

About 30 m at the back of the landslide scarp, faults were also observed. Generally, northeast-trending fractures are also present in the north and southern margins of the headscarp whereas numerous northeast-trending tension cracks were observed on top of the head wall that continue to widen (Figure 10). The entire area of Naga City has numerous faults that are mainly oriented northwest or northeast directions (Figure 11). These faults, however, maybe of various ages given the range of rock types with various ages they cut through. Of particular interest are the set of fractures that have a general northeast direction, which correspond to the northward projection of the strike of the Naga landslide head wall (Figure 12).

FIGURE 10

Figure 10. (A–C) Faults at the margins of the headscarp. (D–F) Faults in areas farther from the landslide (please see Figure 11).

FIGURE 11

Figure 11. Stereoplots and rose diagram plots of structural measurements in Naga City. Figures 10A–C found in rose diagram A. Figure 10D in rose diagram B. Figure 10E in rose diagram E. Figure 10F in rose diagram D.

FIGURE 12

Figure 12. Landslide head wall with fault measurements. The image is an orthomosaic map overlain on a DEM image.

4.5. Quantitative Classification of the Naga Landslide

The H/L ratio of the Naga landslide yields a Fahrböschung angle of 9 deg. Using this angle and the measured volume of 11,000,000 m³ from DEMs, the Naga landslide deposit plots in the field of debris avalanches together with 71 similarly classified events out of the 204 landslides in the global dataset used by Corominas (1996) (Figure 13). Based on the regression equations for each type of landslide, the 2018 Naga landslide event falls under the category of debris avalanches with a 95% confidence interval.

FIGURE 13

Figure 13. Modified plot of debris flow and debris avalanche volume vs. tangent of reach angle for 71 events from Corominas (1996) and additional data. Labels correspond to the obstacles and topographic constraints of the path: f, dense forest: l, open forest; s, scree deposit; b, bends; d, deflections; t, toe thickening in a fan; h, channeling; w, opposing wall; u, unobstructed with free lateral extension. Black line represents the regression equation derived from the dataset whereas the red dashed lines show its corresponding upper and lower limits. Additional data are included to show how the Naga Landslide fits in the global context of debris and rock avalanches. References are as follows: (1) Catane et al., 2019, (2) Evans et al., 2007, (3) Lagmay et al., 2008, (4) Paguican et al., 2012, (5) Johnson, 1978, (6) Guthrie et al., 2012, (7) Naranjo and Francis, 1987, (8) Voight et al., 1983, (9) McSaveney, 1978, (10) Wang et al., 2018 whereas DA and RA stand for Debris Avalanche and Rock Avalanche, respectively.

The range of shear strength of resistance for terrestrial avalanches is from 10 to 100 kPa (Dade and Huppert, 1998). According to the empirically-derived equations of Dade and Huppert (1998) and given the range of shear strength for avalanches, the Naga landslide with a measured collapse height of 200 m and volume of 11,000,000 m³, classifies as a long-runout rockfall if its depositional area falls within 551,571–2,456,460 m². Our calculations of the area covered by the Naga landslide yields a value of 770,723 m², which is within this range. This power law relationship has a 2/3 exponent and is a best-fit regression line with a 95% confidence interval for 76 long-runout rockfalls or rock avalanches that were studied by Dade and Huppert (1998).

4.6. Comparative Analysis With Well-Known Philippine Debris Avalanches

The Naga landslide has a Fahrböschung angle of 9 deg (this study), which is less than the 16 deg Fahrböschung angle reported by Catane et al. (2019). The measured volume is 11,000,000 m³ and is much lower than the 27,000,000 m³ volume reported by Catane et al. (2019) (Table 1). Our estimate of the volume is based on the subtraction and cut and fill calculations of the pre-event and post-event high resolution DEMs. On the other hand, the volume of the Naga landslide reported by Catane et al. (2019) is from field estimates and terrain data derived from a 2013 Interferometric Synthetic Aperture Radar (InSAR), satellite images, and drone surveys.

TABLE 1

Table 1. Computational analysis of Naga Landslide and other known Debris avalanches.

For the Guinsaugon debris avalanche, Evans et al. (2007) reported a Fahrböschung angle of 12 deg and a volume of 15,000,000 m³, whereas Lagmay et al. (2008), published a volume of 15,000,000 m³ and a 10 deg Fahrböschung angle based on their reported H/L ratio. Lastly, the Iriga debris avalanche (DAD2) has a reported H/L ratio equivalent to a Fahrböschung angle of 6 deg (Paguican et al., 2012) and a volume of 1,500,000,000 m³ (Aguila et al., 1986). Based on the reported volumes, H/L ratios and equivalent Fahrböschung angles, the Iriga, Guinsaugon, and Naga landslides, all fall within the debris avalanches category with a 95% confidence interval according to the Corominas (1996) classification with the exception of the parameters used in the report of Catane et al. (2019).

In terms of the equation used by Dade and Huppert (1998), the area covered by the Guinsaugon landslide as reported by Evans et al. (2007) and Lagmay et al. (2008) is within the range of the shear strength of resistance for the prediction of extent of runout for debris avalanches. Similarly, the area of the 1628 Iriga debris avalanche falls within the range of possible shear strength values of debris avalanches and consistent with the power-law relationship between the area and potential energy for long-runout rockfalls or rock avalanches. On the other hand, the area reported by Catane et al. (2019), when calculated using the equation of Dade and Huppert (1998), yields a value of 102 kPa, which is out of the range of shear strength values for debris avalanches.

4.7. Comparison With Analog Models

Scaled analog models were used by Paguican et al. (2014) to study hummock formation and explore their importance in understanding landslide kinematics and dynamics. These models have been used to characterize hummocks in terms of their evolution, spatial distribution, and internal structure from slide initiation to final stop. The models were designed to replicate large-scale volcanic collapses but are also relevant to non-volcanic settings.

The analog model structures of Paguican et al. (2014), in particular the rotated torevas and hummocks, are consistent with field observations of the Naga landslide and suggest a general brittle slide emplacement (Figure 14). The sliding block is composed of brittle limestone, which when extended during transport, forms hummocky structures. Inter-hummocks or those in between hummocks are more broken, finer-grained matrix facies derived from the excessive extension during transport (de Vries and Delcamp, 2015). Hummocks with faults that formed horst and graben structures were preserved within the proximal to medial portions of the Naga landslide deposit field. Thus, the existence of these features imply the presence of sufficient cohesion of the landslide material, such that they remained intact despite the downward movement and spreading.

FIGURE 14

Figure 14. Comparison of analog models of Paguican et al. (2014) with the Naga landslide hummocks and toreva blocks. (A) The different stages of avalanche emplacement showing an interpretation of hummock formation and (B) degree faulting.

5. Discussion

5.1. Landslide Impact

The Naga landslide was an unfortunate event that caused fatalities because of its unexpected long-runout and vastly underestimated impact. Destabilizing conditions and possible triggers that culminated into a massive landslide have been discussed by Catane et al. (2019). They conclude that there was no apparent trigger for the landslide, citing minimal rainfall and no earthquakes immediately prior or during the slide, even though a M_s 3.0 tectonic earthquake occurred on 26 February 2018 about 3 km away. However, the failure was attributed to a marginally stable slope along a low-angle surface. According to Catane et al. (2019), failure could have been due to progressive weakening of the slope mass or further modification and disturbance of the slope. We agree that there is no apparent trigger but note that there was a post-landslide tectonic earthquake with M_s 3.4 on 21 October 2018, which happened 1.5 km away from the landslide area (Figure 2). These may be related to the faults identified in the immediate area of the landslide as both are shallow earthquakes with depths less than 33 km (Figure 11).

From the viewpoint of disaster prevention and mitigation, the manifestation of structures ranging from hairline fractures to several centimeter-wide cracks, which developed a month prior to the event, is very important as it represented a clear warning sign. Due to these telling events, residents close to the headscarp were evacuated by authorities a day before the landslide occurred but the more populated communities one kilometer downslope were not (MGB, 2018). The massive translation of part of the quarried mountain stretched the sliding body, which accelerated to create a long-runout landslide that buried houses far from its source relative to its collapse height. Had the residents in the distal areas been evacuated along with the residents in the areas near the headscarp, then unnecessary deaths could have been avoided.

A review of the 2018-2022 Local Disaster Risk Reduction and Management (DRRM) Plan of the City Government of Naga, certified by the Office of Civil Defense Region 7 based on their formal review along with the Technical Working Group composed of DRRM-mandated agencies, showed that the areas near the headscarp were highly susceptible to landslides while the distal areas buried by the Naga landslide were classified to have low susceptibility (CDRRMO of Naga City, Cebu, 2018). This may have been the basis for the evacuation of the highland areas proximal to the headscarp but not the hard hit lowland areas about 1.2 km away from source of the Naga landslide, which were mapped to have low susceptibility to landslide hazards. This demonstrates clearly that landslide hazard maps are very important as basis for disaster prevention efforts. Such landslide susceptibility maps need to reflect the appropriate understanding of the kinematics of landslides, in particular debris avalanches.

5.2. Landslide Classification

Various geometric parameters of the Naga event measured using satellite imageries and DEMs derived from drone aerial photos made it possible to characterize the morphology of the Naga Landslide and calculate the H/L ratio, volume, Fahrböschung angle and the involved resisting shear stress during transport and emplacement. Based on these mentioned parameters and along with the description of the structural and geomorphic features of the landslide deposit, we were able to categorize the Naga landslide as a debris avalanche.

In general, there is another type of landslide that generates low H/L ratios. These are debris flows which also have a long-runout. Although lahar (mud flow and debris flow) deposits have textures and internal structures similar to the matrix facies of a debris avalanche deposit, debris flow or lahar deposits do not contain debris avalanche blocks which exhibit three-dimensional jigsaw puzzle features and preserved intact primary stratigraphy in hummocks (Sigurdsson et al., 2015). Large boulders within a debris flow deposit are generally surrounded by finer-grained material and concentrate toward the upper surface of the deposit, forming reverse grading. There are also no steep cliffs that form at the distal and lateral edges of a lahar deposit (Ui, 1989).

Several elements that characterize volcanic or non-volcanic debris avalanches are found within the Naga landslide debris field. These features include the presence of an amphitheatre wall (linear headscarp), hummocks, jigsaw cracks, and a long-runout (Ui, 1983; Siebert, 1984; Ui et al., 1986; Andrade and de Vries, 2010; Davies et al., 2010; de Vries and Delcamp, 2015). The Naga landslide is notable for its linear head wall, which is distinct from a horseshoe-shaped amphitheatre commonly found in volcanic debris avalanches (Mt. Galunggung, Mt. St.Helens Siebert, 1984 and Mt. Iriga Paguican et al., 2010).

As for the non-conical shape of the Naga landslide source, we attribute the linear headscarp to a northeast-trending fault, which is similar in orientation to the Central Cebu Fault and one of the principal orientations of fractures measured within Naga City. The large nearly vertical planar feature comprising the head wall of the Naga landslide has slickenlines and slickensides with a consistent rake angle of about 70 deg northeast. Together with other normal faults and thrust faults found at the uncollapsed back- and side-margins of the headscarp (Figure 12) we interpret this planar head wall as an oblique strike-slip fault, which served as a discontinuity and one of the planes of failure of the landslide. The other discontinuity, which acted as the sliding plane, is the interface between the limestone and the underlying sandy to silty sedimentary strata (Figure 8).

All of these can be used to classify the event as a debris avalanche. However, to acknowledge the initial movement which is a translational rockslide, the Naga City landslide is more specifically classified as a Rockslide-Debris Avalanche.

5.3. Excessive Runout

The long-runout characteristic of the Naga landslide relative to the collapse height which was checked using the empirically-derived equations of Corominas (1996) and Dade and Huppert (1998) reveal a classification fit for debris avalanches. The 9 deg Fahrböschung angle and volume plot used by Corominas (1996) indicates relative mobility of the Naga landslide and falls within the range found in debris avalanches. We also calculated a shear stress value (τ) of 34 kPa, which is in the range prescribed for long-runout landslides and consistent with the description for a debris avalanche in terms of excess distance traveled. Such long-runout events, according to Dade and Huppert (1998), happen in both terrestrial and extraterrestrial environments and should have overall resisting shear stress values (τ) ranging from 10 to 100 kPa.

5.4. Comparison With Known Philippine Debris Avalanches

The Naga debris landslide, in terms of its morphology, field deposit description and runout, was compared with the Iriga and Guinsaugon debris avalanches. The collapse of Iriga volcano resulted in two main debris avalanche deposits in the southwest and southeast. This has been reportedly caused by a non-volcanic trigger and was not accompanied by an eruption (Aguila et al., 1986). The deposits are characterized by the presence of an amphitheatre crater, torevas, hummocks of intact conglomerate, sand, and clay units, jigsaw cracked blocks and long runout and cover wide areas in low, waterlogged plains. The presence of intact conglomerates derived from the base of the volcano indicates a very deep failure plane. The younger debris avalanche deposit features discrete hummocks made of ignimbrite (Paguican, 2012).

The Guinsaugon rock slide-debris avalanche slid along intersecting fault planes and joints of the Philippine Fault (Catane et al., 2008; Lagmay et al., 2008) and resulted in a 4.1 km long and 1.52 km wide deposit characterized by a rock slide which transformed into a debris avalanche and consequent debris flows. The debris avalanche deposit comprised of pointed conical hummocks (pinnacle hummocks) and jigsaw-cracked blocks surrounded by a matrix of granular material which was described as a mix of sand and soil from the collapsed mass, whereas the low-lying area of the deposit was reported to have numerous pressure ridges in the southern part, which was interpreted as debris flows (Catane et al., 2008).

The proportion of cohesive material has been observed as an important factor in determining runout length, lateral spread, shape and orientation of individual hummocks in debris avalanches (Vallance and Scott, 1997; Zernack et al., 2009; Paguican, 2012). More cohesive material is associated with pinnacle hummocks and shorter runout debris avalanches compared to the generally circular-based and flat-topped hummocks formed in longer runout debris avalanches (Paguican et al., 2014). Pinnacle or conical hummocks, as called by many authors, are present in both the Guinsaugon and Naga Landslides but are more pronounced in the latter case. The numerous pinnacle hummocks found in the debris field of the Naga debris avalanche indicate cohesive material. Limestone material from the Naga landslide is more cohesive than those of the Guinsaugon and Iriga debris avalanches and is most likely the primary reason for the dominance of pinnacle hummocks in the Naga landslide debris field. Conical hummocks are also found in other debris avalanche deposits aside from these local events. They have been described in Guinsaugon (Catane et al., 2008), Mt. St Helens (Glicken, 1996), Jocotitlan (Siebe et al., 1992), Bezymiannyi (Belousov and Belousova, 1998) and Shiveluch (Ponomareva et al., 1998) volcanoes.

In addition, recalculation of the Fahrböschung angle of the Naga landslide puts it in the same class of landslide as the Iriga and Guinsaugon debris avalanches. The long-runout of the Naga event is likewise confirmed by the equation of Dade and Huppert (1998) on long-runout rockfalls, which is characterized by shear stress values indicative of excess mobility and consistent with calculations made for the Iriga and Guinsaugon debris avalanches.

5.5. Analog Laboratory Models

The Naga landslide is a mass wasting phenomenon that constitutes a catastrophic geologic hazard. However, it is incompletely understood in terms of its kinematics and dynamics. In particular, the physical basis for the extent of runout remains poorly understood. There are many hypotheses to explain the excessive travel distance of such long-runout landslides, which can have H/L ratios of 0.6 for small events, but can be as low as 0.1 for large events with volumes of several cubic kilometers (Heim, 1932; Erismann, 1979; de Vries and Delcamp, 2015). These include: (1) elastic release of fracture energy (Davies and McSaveney, 2009); (2) granular fluidization (Okura et al., 2000; Manzella, 2008; Pastor et al., 2009); (3) trapped air (Shreve, 1968); (4) water pore pressure (Iverson et al., 2000; Manzella, 2008; Pastor et al., 2009); (5) vibrations (Wang et al., 2010); or (6) sudden loss of strength as the material breaks (Quinn et al., 2011; Hungr et al., 2014).

The Naga landslide offers a unique opportunity to test relatively recent literature on the mechanism for the development of long-runout landslides or debris avalanches. Geological fieldwork and DEMs from drone surveys of the Naga debris field allows for a comparison with analog models, where a block of sand material slides down a plane and stretches to create an avalanche. These scaled models allow a sliding body to stretch and lengthen its lobes for increased runout. Downslope movement stops when the resistance of the material is greater than its depth which provides the force for the motion. These laboratory replications reveal that hummocks which form during avalanche events, are morphological expressions of brittle layer deformation due to the spreading. These features are remains of tilted and rotated blocks, whose morphology and distribution depends on the material properties, such as cohesion and viscosity of the sliding layer. Furthermore, its shape, size and density of occurrence can change depending on subsequent spreading, breakup or merging due to a change or restriction in topography (Paguican et al., 2014).

Based on the hummock categorization proposed in these models, the toreva blocks found at the proximal and medial zones of the Naga landslide are transverse in orientation to the direction of the landslide. Torevas, which are actually elongated hummocks are described in the experiments as first-order landslide material, formed during the initial stages of spreading. In the medial portion are pinnacle hummocks. Due to continuous stretching, hummocks will proceed to disaggregate especially at the front portion of the avalanche. In the distal portion of the Naga landslides where there were no topographic barriers (i.e., Lobes 1 and 2 areas), the landslide continued to accelerate and stretch, further disaggregating the larger blocks into fragments.

The transport mechanism for the Naga debris avalanche requires low basal shear resistance to have formed the observed features in the deposit field. Extension during transport produced faulted blocks, including horst and grabens (Figure 5). Greater extension leaves behind the blocks, which are pinnacle in shape in 2D view and pointed conical shapes in 3D view. But the calcareous Carcar and Barili formations, which are sedimentary rocks can also be comminuted. The finer granular matrix surrounding the jigsaw-cracked blocks and those found in inter-hummock areas, are the broken (clastic) Carcar and Barili formation sedimentary rocks, which underwent extensive stretching. The hummocky topography of the Naga landslide, therefore, reflects the dynamics of emplacement, particularly the lower basal friction which becomes smaller for larger debris avalanche deposits because of the volume effect. Because of the presence of hummocks, much of the shear of the moving mass is interpreted to have been concentrated at the basal portion (Davies et al., 2010; de Vries and Delcamp, 2015), which are difficult to find in the deposit field.

Following the analog model sequence, the Naga landslide started as a simple translational slide, where failure along a south-southwest dipping bed of siltstone allowed a large part of the mountain or block of limestone to move down. Features resulting from the extending sliding mass were classified as torevas, megablocks and hummocks with the larger fragments more prevalent in the proximal and medial portion of the debris field. In terms of deposit facies classification, there are the toreva-block, matrix, mixed and basal facies (de Vries and Delcamp, 2015). The debris-avalanche blocks observed in the field were large rocks, sometimes mega blocks that preserve the structure of the source. These megablock features are characteristically found at the proximal to medial area of the deposit (Godoy et al., 2017). A matrix of fine granular rocks surrounds the blocks and is composed of non-graded materials ranging from clay to larger than gravel-sized sediments. The matrix are also present in between hummocks where spreading was much more extensive leading to the formation of more granular collapse material.

Often when water escapes, landslides convert into debris flows. However, this does not appear to have happened in the Naga event as the limestone did not hold much water content and any drained rapidly. Furthermore, although it rained lightly in the morning of the disaster and the weeks prior to the landslide event, there was insignificant rainfall to have caused debris flows to form. While this study relied heavily on comparison with analog models, future work may explore and focus on numerical models to gain new perspectives on the initiation and runout of the landslide.

6. Conclusions

Analysis of satellite images, ground surveys, review of video footage, processing of borehole data, and identification of structural discontinuities reveal details on the conditions that culminated into rapid mass movement and the formation of the Naga Landslide. Comparison of the observations of the Naga Landslide deposits with known examples of debris avalanches show striking similarities. In particular, there is consistency with the descriptions found in the literature in terms of the presence of the following features: (1) Amphitheater crater; (2) Hummocks; (3) Jigsaw cracks in blocks; (4) Megablocks and (5) Long-runout (Ui, 1983; Siebert, 1984; Andrade and de Vries, 2010; Davies et al., 2010; de Vries and Delcamp, 2015). In addition, there is also consistency of the characteristics of the Naga landslide with empirically-derived equations describing debris avalanches and long-runout rockfalls (Corominas, 1996; Dade and Huppert, 1998).

The anatomy of the Naga landslide as described in this work and its comparison with analog models of Paguican et al. (2014) allows a description of the emplacement mechanism of the landslide. Following the analog model sequence, the Naga landslide started as a translational slide when a large block of limestone comprising the mountain slipped along a southwest dipping bed of sandstone/siltstone. Once in motion, the front of the sliding block then accelerated, further stretching the limestone body creating the debris avalanche. Jigsaw-cracked blocks surrounded by a non-graded matrix indicate an en masse flow with a main body with a low shear stress and an underlying sliding boundary with higher shear stress.

This study demonstrates that a debris avalanche and not a simple translational slide devastated villages in Naga City, Cebu on 20 September 2018. This type of landslide with excessive runout is known as one of the most destructive geologic hazards. It claimed the lives of 78 Filipinos with 6 missing and now joins other catastrophic landslides in the Philippines, including the 1628 Iriga and 2006 Guinsaugon debris avalanches. The Naga landslide was used to understand the emplacement mechanism of debris avalanches to advance the knowledge on how to prepare against such hazards. By comparing the deposits of the catastrophic collapse with analog models and well-known debris avalanches, we were also able to provide proper nomenclature, essential in the understanding of factors that led to the landslide disaster in Naga City. The results of this work are an important step to understand Rockslide-Debris Avalanches, necessary for future hazard assessment and risk mitigation. Warnings a month before the catastrophe in the form of hairline fractures that progressed to centimenters-wide cracks saved people living near the headscarp of the landslide. Unfortunately, the long-runout potential of the landslide was neither anticipated nor understood. Such understanding of long-runout events, which this study advances, is crucial in hazards and risk assessment.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

AL is the head of the investigative group for this study while the rest of the authors are members. CE, AY, and JS contributed in data gathering, modeling, analysis, and writing different parts of the manuscript. GC contributed in the calculations and computational results used in the study.

Funding

This study was funded by the UP Resilience Institute under the Office of the President of the University of the Philippines.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors express their gratitude to the UP NOAH Researchers, especially Ms. Francesca Llanes, Ms. Feye Andal, Ms. Dianne Bencito, Ms. Krichi Cabacaba, Mr. Prince Goting, Mr. Marc Moises, Ms. Patricia Delmendo, Ms. Camille Cuadra, Ms. Girlie David, Ms. Lia Gonzalo, Mr. Joshua Abaygar, Ms. Trizia Muldong, Mr. Michael Ramos, Mr. Niño Joshua Tangaan and Mr. Archie Amba for their contributions to the entire report on the 2018 Naga Landslide event regarding information on topics not included or not the focus of this publication. We also would like to thank Dr. Benjamin van Wyk de Vries, Dr. Gioachino Roberti, Dr. Fabio Luca Bonali, and Dr. Malcolm Whitworth for virtual discussions and feedback on the analysis of the Naga Landslide. Last but not least, the authors also acknowledge the Center for Environmental Informatics (CENVI) for the additional drone image datasets used in this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.00312/full#supplementary-material

References

Aguila, L., Newhall, C., Miller, C., and Listanco, E. (1986). Reconnaissance geology of a large debris avalanche from Iriga volcano, Philippines. Philipp. J. Volcanol. 3, 54–72.

Google Scholar

Allen, C. R. (1962). Circum-Pacific faulting in the Philippines-Taiwan region. J. Geophys. Res. 67, 4795–4812. doi: 10.1029/JZ067i012p04795