GIS-BASED LANDSLIDE VOLUME ESTIMATION USING DIGITAL TERRAIN MODELS DERIVED FROM LIDAR AND RADAR SYSTEMS: CASE OF PIDIGAN, ABRA

Southwest Monsoon (Habagat) and Typhoon Luis caused a deep-seated landslide that struck Sitio Kayang, Brgy. Immuli, Pidigan, Abra on August 15, 2018. Rainfall-induced deep-seated landslides displace partially at a time which necessitates the determination of remaining landslide volume along the slope. In this study, the potential landslide volume and mass transport were estimated using several remote sensing products, including SAR (Synthetic Aperture Radar) data and LiDAR-DTM (Light Detection and RangingDigital Terrain Model). The post-landslide DTM was generated using Sentinel-1 SAR data. The potential landslide volume and landslide failure surfaces were ascertained through the stability analysis using Scoops3D, while the mass transport volume was obtained from the preand post-landslide DTM. Results showed that the estimated total volume in all the landslide areas was 135,962 m. Meanwhile, the remaining landslide volume (i.e., difference between potential volume from pre-landslide event and volume of transported mass) yielded illogical values due to the derived large mass transport values. This blunder may be attributed to the generalization of the transported volume (due to Sentinel-1 DTM coarse resolution), and decorrelation due to vegetation cover. Overall, the LiDAR-DTM data delivered a high-resolution estimation of the potential landslide volume and proved to be useful for landslide application studies. Future studies may incorporate field data (e.g., geotechnical parameters, groundwater, landslide actual measurements) for more accurate performance of stability analysis and may best to utilize LiDAR-DTM in post-landslide volume computation for a more reliable estimation of mass transport and potentially remaining landslide volume.


Background
Landslides can cause destruction and damage to properties that may even lead to human fatalities. By definition, landslide is the movement of mass of rock, soil and debris down the slope under the influence of gravity (Cruden, 1991). This hazard may be induced by rainfall, earthquakes, volcanic eruption, by other natural phenomena or even man-made activities (Wieczorek, 1996).
On August 15, 2018, the combined effects of Southwest monsoon (Habagat) and Typhoon Luis caused a landslide that struck Sitio Kayang, Brgy. Immuli, Pidigan, Abra, forcing residents to evacuate from their homes (NDRRMC, 2018). According to DOST-PHIVOLCS (2018), a deep-seated landslide hazard was present in the area and was also assessed by DENR Mines and Geoscience Bureau (MGB) via HazardHunterPH (https://hazardhunter.georisk.gov.ph/) as having moderate landslide susceptibility with a possible landslide occurrence.
Remaining landslide volume along the slope of deep-seated landslides may pose hazards on the area when triggered by heavy rainfall and other natural or anthropogenic causes (Lin et. al, 2020). To assess the extent and mitigate the possible hazard, the volume of the possible remaining landslide must be determined in the area.

Objectives and Expected Output
This study aims to determine the potential volume of the 2018 landslide event in Sitio Kayang, Immuli, Pidigan, Abra and to obtain the remaining landslide volume (i.e., the difference between potential volume from pre-landslide event and volume of transported mass) in the study area using LiDAR and RADARderived DTMs. Possible damages to the nearby houses in the landslide hazard area was estimated as well.

Scope and Limitations
This study focused on landslide hazard areas situated in Sitio Kayang,Brgy. Immuli,Pidigan,Abra (17°31'24.47" N,120°36'16.71" E) as shown on the hazard map (see Fig. 2) which was requested from PHIVOLCS Dynaslope Project. This paper will solely focus on the use of LiDAR DTM data and publicly available remotely sensed data such as Sentinel-1 images. No further ground validations were performed due to mobility restrictions caused by the ongoing global pandemic.

Deep-Seated Landslide
Deep-seated landslides are gradual-moving landslides, rooted in bedrock that has wide coverage, and capable of destroying infrastructures and residential development (Washington State Department of Natural Resources, 2017). According to Cruden et al (1996) and Lin et al (2014), a deep-seated landslide may be described by the following topographic features: a crown with tension cracks, a main scarp or head scarp, a related slide blocks, a minor scarp, main body, transverse tension cracks, and a toe bulge or swelling area.
A landslide scarp is a steep and inclined exposed soil located at the head of the landslide where the ground surface failure occurred (Ohlmacher, 1999). An example of the deep-seated landslide captured by PHIVOLCS Dynaslope Project during their field surveys conducted from 2015 to 2018 in Pidigan, Abra is shown in Fig. 3.

LiDAR DTM for Landslide Applications
Light Detection and Ranging (LiDAR) is a remote sensing technique that measures distances in the Earth's surface or features by using light in the form of pulsed laser emitted by a laser scanner. The recorded pulses bouncing can produce dense, three-dimensional data on Earth's terrain and surface features, thus the generation of a high-resolution digital elevation or terrain model is made possible (NOAA, 2021).
LiDAR technique has been used in different landslide applications such as detection of mass movements, landslide hazard assessment and susceptibility mapping, landslide 3D modelling, and landslide monitoring (Jaboyedoff, 2010).
A study by Lin et al. (2020) made use of LiDAR DTMs to obtain possible landslide volume at the critical surface of a pre-landslide event in Basihlan River Basin, Taiwan, by determining topographic variations between pre-and post-landslide event and computing the remaining volume of deep-seated landslides that may be transported in the future.

Stability Analysis
Stability analysis is a static or dynamic process of examining and assessing the stability of a slope. Some of the physical models used in evaluating stability analysis are limit equilibrium models, empirical approaches for rock slopes, finite element models, and district element codes (Geoengineer, 2021).
Limit Equilibrium Models (LEM) initially established the failure surfaces then analysed the slope of the stability based on physical laws of equilibrium of forces on the predetermined failure surfaces through the calculation of the factor of safety (Ma, 2018). Factor of safety (FOS) is the ratio of the resisting strength of the material to the applied load or forces. It should always be greater than 1 to declare a stable slope condition (Bagtang et al, 2019), otherwise the slope will be considered as a critical slip surface. Factor of safety can be expressed by Eq.
(1) as given by Duncan (1996): (1) where c = cohesion ϕ = angle of internal friction of soil σ = normal stress on the slip surface τeq = shear stress required for equilibrium Some of the commonly used limit-equilibrium techniques are Ordinary Method of slices, Bishop's modified method, Janbu's simplified method, and Modified Swedish method (Duncan, 1996).
The Bishop's simplified method determines the normal force acting upon the potential failure surface by calculating the force equilibrium with respect to the vertical direction on the base of each slice (Reid, 2015). This technique will be used in this study for the evaluation of the stability of the ground in the study area. The soil parameters required for running a Bishop simplified technique are cohesion, angle of internal friction and unit weight. Cohesion describes how soil particles are attracted to each other or bonded at the molecular level (Molla, 2012). The angle of internal friction pertains to the friction shear resistance of soils together with the normal effective stress while unit weight is the weight per unit volume of soil (Koliji, 2018).

METHODOLOGY
The flowchart shown in Fig. 4 shows the general methodology of the study. Succeeding subsections will discuss the data and methods used in this study.

Data and Software
To identify the location of the landslide scarps in the study area, the landslide hazard zone maps, site photos, and landslide features were requested and provided by Dynaslope Project. The digital terrain model (DTM) for the pre-landslide incident was obtained from UP TCAGP, DREAM Commercial Services. Due to limited availability of resources, the soil type was identified using NAMRIA's Geoportal website. The soil parameter data was taken from Bagtang et al (2019). Sentinel-1 images dated March 3 and March 15, 2019 were downloaded in Copernicus Open Access Hub (https://scihub.copernicus.eu/) to generate the post-landslide DTM due to the unavailability of a recent LiDAR DTM.
The software used in this study are USGS Scoops3D for stability analysis, ESA SNAP and SNAPHU for the DTM generation using Sentinel-1 images, and Quantum GIS and ArcGIS for visualizing and mapping of the data.
Scoops3D is a computer program, developed by the U.S. Geological Survey (USGS), for slope stability analysis by using a digital elevation model or DEM. It uses a "3D method of columns limit-equilibrium analysis" to compute for the stability of the potential slope failures (landslides) with a spherical potential slip surface (Reid, 2015).

Delineation of Scarps
The initial identification of landslide scarps was done using archived high-resolution images in Google Earth since it was difficult to identify the landslide scarps in Sentinel-2 images alone. Additionally, the shaded relief from the LiDAR-derived DTM in different azimuth angles aided on the delineation of scarps.

Sentinel-1 to DTM
The workflow for the generation of DTM from Sentinel-1 images is shown in Fig. 6.  The two downloaded Sentinel-1 SAR products have a perpendicular baseline of 127 meters and a temporal baseline of 12 days. These images were imported in the TOPS Split tool to limit the bursts useful for the analysis and were employed with precise orbit information. Orbit-applied images were then coregistered through Back Geocoding and performed with Enhanced Spectral Diversity to improve the quality of coregistration. Consequently, the interferogram and coherence images were created from the stacked images which includes the phase variation due to the earth's curvature, topographic phase, and other components (Braun, 2020). To eliminate the seam lines between the bursts, TOPS Deburst was performed to the interferogram. A subset was applied to the previous output to trim the artefacts at the edge of the image. Multilooking and Goldstein Phase Filtering were then employed to optimize the output image and to reduce the random noise from the true backscatter signal (Braun, 2020). Using SNAPHU via SNAP, the altitude ambiguity was eliminated from the phase image to produce altitude variation in terms of phase between the two SAR images. To convert radian units of the phase variation to metric units, the Phase to Elevation tool was applied to the previous output. Finally, the DTM was applied with geometric correction through the Range Doppler Terrain Correction tool.

Stability Analysis using Scoops3D
Soil parameters were identified based on the soil type of the area due to unavailability of field data, the. According to the soil type layer from the NAMRIA's Geoportal, the soil in Sitio Kayang, Pidigan, Abra is made up of Bauang clay loam. A clay soil has a cohesion of 5-20 kPA, angle of internal friction of 5°-25° (Bagtang et al, 2019), and a unit weight of 18 kN/m 3 (Zhu, 2014.).
To perform the stability analysis, LiDAR DTM (in ASCII file) is imported in Scoops3D. The soil parameters used for the material properties in the area is shown on The generated products from the stability analysis are the raster data of the factor of safety (FOS), volume, slope, and other output files. Landslide susceptible areas were mapped out by isolating FOS values of less than 1.

Topographic Analysis
Sentinel-1 derived DTM was resampled from 13.5 meters to 1 meter to match the resolution of the LiDAR DTM. The terrain variation between the pre-and post-landslide was derived by obtaining the difference of the before and after landslide DTM rasters. The mass transport volume was then computed for each of the landslide areas as mapped from the stability analysis.

Potential Landslide Volume and Remaining Landslide Volume Computation
The potential landslide volume is estimated using outputs of the stability analysis. One of the output files is a new DTM in which all masses with a FOS less than one are removed from the LiDAR DTM. This raster was subtracted from the original DTM to determine the potential volume of the critical surfaces at the brink of landslide failure. Remaining landslide volume is then calculated by getting the difference of the potential landslide volume and the volume of transported material during the landslide in 2018.

Sentinel-1 Derived DTM
The SRTM-DTM and Sentinel-1-derived DTM are shown sideby-side in Fig. 9. Comparing the two DTMs, the terrain was properly represented in some portions of the Sentinel-1 DTM compared to SRTM-DTM, such as on the upper left and lower right region of the images. Moreover, the difference between the two images was mainly apparent on the central portion where widespread vegetation cover was present.

Potential Volume Estimation
The map of factor of safety as generated from the stability analysis is shown on Fig. 10. Areas marked in red are the potential landslide failure surfaces (FOS less than 1). Green and orange colored areas imply stable and marginally stable landslide susceptibility, based on the FOS categorization by Bagtang et al. (2019).  Some of the identified landslide susceptible areas (Fig. 10) agree well with the landslide scarp features in Fig. 8, as they matched well with the location of each critical surface. Landslide susceptible areas with established scarp presence were highlighted in orange in Table 3.

Landslide Susceptibility FOS
The potential landslide volume in landslide-susceptible areas is the landslide volume at the verge of surface failure (Fig. 11). From the extracted masses, the volume was computed using zonal statistics. Table 3 lists the volume, horizontal area, and centroid coordinates (UTM Zone 51N) of the landslide susceptible areas. The deep-seated landslide (ID=4) obtained a landslide volume of 22,893.89 m 3 . Overall, the total landslide volume is 135, 962 m 3 .   Table 3. Potential Volume.

Topographic Analysis and Remaining Landslide Computation
Figs. 12 and 13 show the DTM of the pre-landslide and postlandslide events. As expected, the LiDAR-DTM exhibited more detailed, smooth, and continuous variations in elevation compared to Sentinel-1-based DTM with a coarse pixel resolution (~ 13.5 m). The terrain changes between the two DTMs are displayed in Fig. 14. The deep-seated landslide area (ID=4) and two other landslide regions showed erosion marked in red color from the generated negative values. Fig. 15 shows the volume of the transported materials for each landslide area. The resulting difference image produced a discrete appearance due to the coarse resolution of the post-landslide DTM.

Cost Estimation of Damages to House Structures
From the digitized houses provided by Dynaslope Project, 24 houses were identified that were located within the landslide hazard zone. The estimated damage costs for a farmhouse structure are shown in Table 5 (Mag-aso, 2018).

Variables
Costs of Damage, PHP Costs of new materials for repair 65, 500.00 Labor costs 21, 350.00 Total 88, 850.00 Table 5. Total value of damage to farmhouse structures.
Since houses are mostly situated in agricultural lands, these houses were assumed to be constructed as farmhouses. Therefore, the estimated damage cost for house structures due to a landslide event is estimated to be Php 2,132,400.

CONCLUSION AND RECOMMENDATIONS
In this study, the landslide volume in Immuli, Pidigan was estimated using stability analysis of pre-landslide event DTM in Scoops3D. The LiDAR DTM data delivered a high-resolution estimation of the potential landslide volume in the study area.
Landslide susceptible areas (FOS is less than 1) were identified corresponding to landslide scarps as identified by PHIVOLCS Dynaslope. Meanwhile, the topographic change in the area was mapped and only few areas were observed to have experienced soil erosion. Factors such as the coarse resolution and decorrelation of the post-landslide event DTM affected the reliability of derived mass transport and remaining landslide volumes yielding unrealistic values.
For future studies, the use of high-resolution imagery or aerial images is suggested for exact delineation of landslide scarps. Also, elevation profiles and field data (e.g., geotechnical parameters, actual landslide measurements) can be incorporated for more precise results and validation of the results of stability analysis. The effect of groundwater in slope stability may also be explored since it influences soil strength and consequently, landslide occurrences. Lastly, it is advantageous to utilize LiDAR-derived DTM for post-landslide events for a more detailed representation of terrain variations and obtain reliable computation of mass transport volume and remaining landslide volume that may be transported in the future leading to another disaster event.