Modeling landslide susceptibility and landslide volume per geomorphologic landform unit

The analysis of landslide susceptibility and landslide volumes in landforms can provide information for planning disaster management in an area. Landslide susceptibility per landform unit and the potential contribution of material delivered from each unit was calculated for a 105 km 2 watershed on the south ank of Pico de Orizaba volcano, Mexico. The landslide susceptibility is calculated from the area and frequency of landslides. The volume is obtained from detailed geometric values of shallow landslides in order to establish an empirical relationship that takes the form of a power law from which the potential volume of all shallow landslides is calculated for the watershed. The study shows that most of the landslides are on volcanic landform units; however, the landslides in sedimentary units contribute more sediments per square kilometer. It also shows that landform units can be used to explain the predisposition and variability of landslide sediment production for a large and complex geological watershed. how landform units this paper is to provide a standardized method in a GIS system to evaluate landslide susceptibility per geomorphologic landform unit and the sediment production by shallow landslides. The study divides the watershed into 15 geomorphologic landform units. Each landform unit is derived by aerial photo interpretation, geological eldwork, and morphometric parameters (hypsometry, drainage density, vertical erosion, and slope). For each landform unit, the area and the spatial-temporal distribution of landslides are used to determine the landslide susceptibility. The potential total landslide volume in the watershed is obtained by establishing a power law model for shallow landslides. undergoing change by landslides, local land use, and environmental conditions. The main aim of this study was to model and assess landslides susceptibility and landslide volume in terms of landform units. The approach uses and adapts the LHZP of WSDNR to map landslides, and to dene landform units. It also uses aerial photographs and DEMs acquired with UAV surveys and geomorphological evaluation to establish current and reconstructed pre-slide topography of shallow landslides, from which a power law is established to calculate potential total landslide volume.


Introduction
In mountainous volcanic terrains, shallow and deep-seated landslides represent a major natural hazard for human settlements and their economic activities (Macías, 2005;Capra and Lugo-Hubp, 2006;Miyabuchi et al., 2015). Such is the case for Pico de Orizaba volcano (5675 m a.s.l.), the highest dormant stratovolcano in Mexico. This study refers in more detail to shallow landslides, including debris ows and debris slides, that affect the Río Chiquito-Barranca del Muerto watershed on the south ank of Pico de Orizaba. In the area, steep slopes, heavy rainfall, and changes in land use from forests to subsistence agriculture contribute to the potential mobilization of landslides and debris ows of weathered volcanic and epiclastic deposits that cover and/or coexist with sedimentary deposits. This creates a dangerous situation for inhabitants of cities on the southern anks of the volcano. In 2003 and 2011, landslides delivered material to the stream in the watershed; this increased the destructive power of the debris ows, which damaged human settlements and caused fatalities and disrupted economic activity in Balastrera, a town at the mouth of the watershed (Rodríguez et al., 2011;Legorreta-Paulín et al., 2013). Despite the importance of assessing such processes, there are few landslide inventory maps and few evaluations of the sediment budget and sediment availability connected to landslides in the area.
Worldwide studies of landslide susceptibility in an area consider the distribution, frequency, and size of landslides in relation to homogeneous areas, called landform units (Washington State Department of Natural Resources (WSDNR), Forest Practices Division, 2006;Hervás and Bobrowsky, 2009;Jasiewicz and Stepinski, 2013). The use of these landform units allows a better understanding of landslide susceptibility and the volume of material displaced by a landslide in a watershed with complex geology and geomorphology (Hungr et al., 2008;Broeckx et al., 2016). Techniques for calculating the landslide volume may be based on the geometric shape (conical, elliptical, circular, planar, irregular, etc.) of the current physical form of a landslide (rockfalls, deep-seated landslide, debris ow, etc.), and on analytical methods for resolution (geometric calculation, differential calculus, etc.) (Šíma and Seidlová, 2014;Tang et al., 2018). The physical form of the landslide and landform units are obtained by using detailed Digital Elevation Models (DEMs) through the use of differential GPS (Travelletti, and Malet, 2012;Young, 2015) and imaging equipment (digital cameras, multispectral sensors, Light Detecting and Ranging (LiDAR) instruments and Synthetic Aperture Radars (SAR)) carried by airplanes or unmanned aerial vehicles (UAVs) (Rango and Laliberte, 2010;UNEP, 2013;Chen et al., 2014). To determine the volume of landslide materials in large affected areas, empirical relationships are established by use of a power-law function that links geometric measurements of landslide area to landslide volume (Kalderon-Asael et al., 2008;Guzzetti et al., 2009;Wenkey et al., 2011).
In Mexico, landslide susceptibility has been addressed by use of Geographic Information Systems (GIS) with multiple scopes and scales (Capra and Lugo-Hubp, 2006;Oliva-González et al., 2014;Morales, 2016;CENAPRED, 2017). However, few works refer to landslide susceptibility per landform unit (Legorreta-Paulín 2013; Alanis-Anaya, 2017) and landslide volume (Capra and Lugo-Hubp, 2006;Muñoz-Salinas et al., 2009). For Pico de Orizaba, geological studies and computer simulations with GIS and remote sensing have produced lahar hazard maps for the stream systems (Sheridan et al., 2001;Hubbard et al., 2007). An overview of hazard assessment of landslide and low-concentration ows triggered by climatic and geological conditions in Veracruz State included the present study area (Rodríguez et al., 2011). A map of landslide susceptibility within the watershed was created by use of multiple logistic regression (MLR) and SINMAP (Legorreta-Paulín et al., 2013). Also, landform units have been de ned by following the Landslide Hazard Zonation Protocol (LHZP) of the WSDNR, Forest Practices Division 2006(Legorreta-Paulín et al., 2014. No work has hitherto modeled the landslide sediment production per landform unit. As a result, there is a lack of knowledge on how landform units differ in the intensity with which they trigger landslide sediment production. To address this de ciency, the aim of this paper is to provide a standardized method in a GIS system to evaluate landslide susceptibility per geomorphologic landform unit and the sediment production by shallow landslides. The study divides the watershed into 15 geomorphologic landform units. Each landform unit is derived by aerial photo interpretation, geological eldwork, and morphometric parameters (hypsometry, drainage density, vertical erosion, and slope). For each landform unit, the area and the spatial-temporal distribution of landslides are used to determine the landslide susceptibility. The potential total landslide volume in the watershed is obtained by establishing a power law model for shallow landslides.

Study Area
The study area lies on the south ank of the volcano Pico de Orizaba, which rises at the junction of Puebla and Veracruz states (Fig. 1). Pico de Orizaba is a dormant stratovolcano in the eastern part of the Neogene-Quaternary Transmexican Volcanic Belt physiographic province (Carrasco and Gómez-Tuena, 1997;Macias, 2005). The volcano shows complex phases of construction and destruction during its evolution; these include two well de ned volcanic collapses that have been attributed to tectonic faults and hydrothermal alteration (Carrasco and Gómez-Tuena, 1997;Hubbard et al., 2007).
Its lavas, pyroclastic ows, debris avalanches, fall deposits, and lahars overlie a weathered and highly dissected, folded and faulted, Cretaceous basement of limestone and shale. The sedimentary basement (with NW-SE trend) and preexisting NE-SW and N-S trending faults (presumed active) and fractures could potentially trigger shallow landslides (Concha-Dimas et al., 2005). The present study refers in more detail to the Rio Chiquito-Barranca del Muerto watershed on the south ank of Pico de Orizaba and to shallow landslides triggered by rain events. The area covers approximately 105 km 2 of rugged terrain, ranging from 1280 to 5675 m a.s.l, and with the slope angle generally increasing with elevation (from ~ 5° to ~ 65°). Rio Chiquito-Barranca del Muerto is a fourth-order watershed (Strahler's method) with a radial sub-parallel and dendritic pattern which ows into the Gulf of Mexico. Andesitic and dacitic Tertiary and Quaternary volcanic rocks form the main rock types in the study area followed by Cretaceous limestone and shales. Sedimentary landform units represent almost 30% of the total watershed area. Rainfall in the study area is concentrated mainly in the boreal summer (from May to November). Rainfall is mostly due to the formation of low-pressure systems (hurricanes) in the Paci c region and the in uence of westerly winds that causes clouds and precipitation over Pico de Orizaba. The maximum mean annual rainfall in a period of ~ 20 years is ~ 1100 mm, mostly concentrated in the mountain area at > 4000 m a.s.l. (CLICOM, 2014). The minimum precipitation is recorded in the lower parts of the watershed where annual rainfall is ~ 927 mm/yr (CLICOM, 2014). The contribution of rainfall from hurricanes is important because some of these events affect the hillslopes of Pico de Orizaba (García, 2004). The study area is naturally prone to landslides due to the erosive effects of high precipitation on steep and weathered volcanic and sedimentary terrains, exacerbated by deforestation and agricultural activities. The soils in the area are inherently weak and further lead to instability. The observed soils in volcanic terrain are Ocric Andosol, Eutric Cambisol, Lithosol, Eutric Regosol, Chromic Vertisol and Rendzina type for the sedimentary terrain (Alanis-Anaya, 2017). Within a period of nine years (2003 to 2012) changes in land use have been signi cant. Forested areas (2.64 km 2 ) have been lost in favor of urban, pasture, and agricultural areas. Deforestation and modi cation of the slopes in favor of human activities have left the main river and its tributaries with a buffer of only 50 100 m of forest along the hillsides (Alanis-Anaya, 2017). Such massive destruction of the forest has led not only to restriction of the natural vegetation to the gullies, but also to accelerated erosive and gravitational processes in sedimentary and volcanic deposits. For example, in 1974 when the government authorized felling of forest to establish agricultural lands along the western hills of Pico de Orizaba, a large gully (3 m wide x 350 m long x 4 m deep) developed as a result of heavy rains. In 1982, the gully size increased to 1.7 km long x 40 m deep (Werner, 1996).

Method
We obtained the landslide susceptibility and landslide volume per geomorphologic landform unit in the study area by the following four steps.
Step 1: Incidence of gravitational processes was mapped; a landslide inventory map was derived from preexisting landslide maps, reports, eld work, and interpretation of remote sensing data.
Step 2: A landform map was produced; a morphometric and morphographic analysis was conducted to de ne homogeneous landform units.
Step 3: A landform units susceptibility map was generated according to the LHZP of WADNR Forest Practices Division, 2006.
Step 4: An empirical relationship of area-volume was established in the form of a power law to estimate the potential landslide volume per landform unit. The implementation of each step was as follows.
The inventory (step 1) was obtained through compiling data from an existing landslide geo-database and with a landslide map that identi ed 442 gravitational processes (Legorreta-Paulín et al., 2013). This landslide inventory map uses the LHZP (WADNR Forest Practices Division, 2006) adapted from Cruden and Varnes (1996) and Wieczorek (1984) to classify landslides into shallow landslides (slope failures de ned by a failure surface within the forest rooting zone), debris ows (a shallow landslide that forms from the disaggregation of materials on a steep slope, involving the rapid movement of the soil and regolith over bedrock), debris slide ows (a shallow landslide that ows within a channel formed either by the valley walls of a low-order tributary or by levees of its own making), deep-seated landslides (a landslide in which most of the area of the slide plane or failure zone lies below the maximum rooting depth of forest trees), earth ows (a type of slow-moving, deep-seated landslide), and rockfalls (individual blocks or rock that become detached from a steep slope and descend through the air by falling, bouncing, or rolling before coming to rest on gentler slopes) (WADNR Forest Practices Division, 2006). The same approach was followed to develop a new landslide inventory (Fig. 2). It incorporates 165 newly identi ed landslides mapped on the basis of analog aerial photographs (scale 1:20,000 from 1994 and 1:10,000 from 2008) and panchromatic SPOT satellite imagery (from 2003 and 2012), and by local eld surveys. Landslide head scarp, evacuation zone, and deposit (these last two were not always visible) were mapped on transparent plastic and transferred into the GIS. The interpretation with aerial photographs for landslide mapping was aided by eld work. Three sessions of eld work per year were conducted from 2012 to 2017 along the main and secondary rivers, for data validation, landslide and landform mapping, veri cation, and landslide volume measurements. In the case of the SPOT satellite imagery, interpretation was made directly on the screen and then stored in the GIS. For the mapping, landslides that were substantially reforested and estimated to have been covered by vegetation for more than 22 years (the span of this study, from the 1994 aerial photographs to 2016) were not included in the inventory or in the evaluation of volume. Nor were landslides used whose area was less than 0.2 hectares (minimum mapping size for the expected nal publication at 1:50,000 scale).
The landform units map (step 2) was the result of combination and interpretation of thematic maps according to geomorphological criteria and eld evaluation (Demeck, 1972). The method for mapping landform units was adapted after Legorreta-Paulín et al. (2014). To de ne landform units, we used more speci c criteria that included photointerpretation, eld geology, land use, hypsometry, slope, drainage density, and the vertical erosion (Fig. 3). Drainage density is de ned as the sum of all talweg lengths per unit area. It was calculated by dividing the study area into 1 km 2 grids, and the total talweg length is obtained for each square kilometer (Simonov, 1985, as cited in Zamorano-Orozco, 1990; Lugo-Hubp, 1991). Vertical erosion was calculated by dividing the study area into 1 km 2 grids, and the maximum vertical distance was obtained for each square kilometer. The obtained values were interpolated. The map represents the average vertical height between the talweg and the watershed (Simonov, 1985, as cited in Zamorano-Orozco, 1990Lugo-Hubp, 1991).
Landform units were de ned by geomorphologic rules adopted by the Washington State Department of Natural Resources to address landslide hazards (WSDNR, Forest Practices Division, 2006). This process used GIS overlay of the topographic map, digital aerial photographs, geology, land use and morphometric and morphographic maps to delineate summits, runoffs, and edges between landforms. Field recognition with ninety-six control points per landform unit was used to verify the units. Field veri cation was based on geomorphology, geology, slope instability, and land use changes that de ne the units. The veri cation points were assigned by using a strati ed and random sampling according to the proportion of surface that covers each class, to avoid errors of subsampling or oversampling.
To produce a landslide susceptibility map (step 3) we followed the LHZP of WSDNR (WSDNR, Forest Practices Division, 2006). The guidelines of the LHZP are based on the experience of the Washington Landslide Hazard Zonation team in analyses of susceptibility of landforms along watersheds (Lingley, 2004;Wegmann, 2004). A semiquantitative approach was used to determine the landslide susceptibility in each landform unit. First, the number of landslides and their accumulative area in ha per landform unit was obtained by using the GIS landslide inventory map. Both values, the number of landslides and their accumulative area per landform unit were divided by the period of time spanned in the study (in our case 22 years), and the landform unit area (in acres). These values were then heuristically multiplied by one million and rounded for easier interpretation (Powell, 2007) to obtain the landslide frequency rate (LFR) and the landslide area rate (LAR). The LFR represents the number of landslide events registered in the reference landform unit in a time period, and the LAR represents the spatial area (in acres) impacted by landslides in the reference landform unit in a time period. After the LFR and LAR values had been determined, each was classi ed into Low, Moderate, High or Very High rating (Table 1). The qualitative LFR and LAR values were then entered into a matrix (Table 2) in order to determine the overall susceptibility rating for each landform (WSDNR, Forest Practices Division, 2006). For instance, the landform known as the pyroclastic ramp unit has a LFR value of 815 (15 landslides / 815 acres / 22 years * 1000000) and a LAR value of 258 (4.6 acres / 815 acres / 22 years * 1000000). Each of these values is classi ed as High (Table 1). According to the juxtaposition of these two qualitative values into Table 2, this landform has a high landslide susceptibility.
For the landslide susceptibility analysis, the uvial valley landform unit is not taken into account because this landform unit is contained in the remaining 14 landform units. Instead, the area of uvial valleys that affects each landform unit is given in the Results section.  Finally, we estimated the landslide volume per unit (step 4). Sixty-one landslides (~ 10%) from the new landslide inventory (encompassing 607 landslides) were veri ed in the eld for calculation of their area and volume. The area and volume of the landslide head scarp were measured with tape, stadia rod, laser range nder, differential GPS, and a UAV. We excluded the volume and area of the deposits and the slip evacuation area. Tape, stadia rod, and laser range nder were used to measure head scarp depths along transects by the procedures suggested by Gardiner and Dackombe (1983) and Compton (1985). Geometrical measurement of the 61 landslides used a Phantom 4 Pro quadcopter with a camera with 20 megapixels resolution and a CMOS sensor of 1". Flight planning used the Pix4D capture program. We used a photogrammetric capture on a grid ight pattern. With this ight pattern, aerial images were acquired with an average ight height of 100 m and at a speed of 1.5 m / s. The front and side overlaps between aerial photographs was 80%; this was because the actual ight was not absolutely straight owing to the wind and the potential GPS error. All images were acquired with the optical axis vertical, since the objective of the project was only to generate the DEM. Along with the ight planning, some ground control points (GCPs) were de ned in the area, and these served to georeference the aerial photographs. We used a dual-frequency differential GPS (with horizontal accuracy of 0.5mm and vertical accuracy of 0.10mm) to provide better georeferencing for the GCPs. Owing to the cost and effort required in remote areas with di cult access, only between 1 (90% of the sampled landslides) and 4 (one at each extreme corner of the landslides) GPS ground control points were taken to georeference the photographs during the post-processing.
The post-processing of the images used Agisoft PhotoScan photogrammetric software. We followed a standard postprocessing work ow adapted after Ouédraogo et al. (2014). The work ow in PhotoScan was as follows. 1) Selection and addition of all captured photos for each landslide; inspection of the general coverage and image sharpness quality to exclude low-quality images from the software geometric analysis. 2) Alignment of the loaded aerial photographs; PhotoScan employs the scale-invariant feature transform algorithm (SIFT) to nd matching points between overlapping images, to estimate camera position for each photograph, and to build sparse point cloud values. 3) Addition of markers that represent GCP were added (Root Mean Square Error (RMSE) range from ~ 1 m to ~ 0.5 m). 4) Construction of a cloud of dense point values; the high-accuracy option is used to densify the sparce points in the PhotoScan's work ow menu. 5) reclassi cation of the dense point cloud values within the PhotoScan environment by automatic division of all points into ground and the other elements (vegetation, buildings etc.). 6) creation of a DEM, and 7) creation of an orthophotograph (PhotoScan, 2013;Ouédraogo et al., 2014). DEMs at 50 cm per pixel were obtained.
In ArcMap, the current sliding topography and the reconstructed pre-sliding topography were created. This pre-sliding topography was set heuristically by eld observation and reconstruction of the original topography from the current sliding topography generated by the UAV. The nodes from polygons that represent the current landslide headscarp boundary were used to create a point le which in turn was used to extract the elevation data from the DEM for the topography of the landslide. The nodes were interpolated to create the interpolated pre-sliding topography. No internal topographic relics were observed in the eld that could be used to improve representation. It was considered that the reconstructed pre-sliding topography in general gives an adequate representation of topographic geometry, subsequently modi ed to varying degrees by gravitational and erosive processes. If exogenous forms such as valleys or sliding surfaces were removed from the topographic map, an approximation to the original relief could be obtained by interpolations (Lugo-Hubp, 1991).
The landslide volume was calculated from the algebraic difference between the reconstructed pre-sliding topography and the actual topography (Fig. 4), multiplied by the area of the pixel.
Selection for sampling landslides was based on eld reconnaissance of the landslide process and the accessibility to the sites. Data limitation included lack of sampling; also, statistical validity was precluded by the wide dispersion in a scatter plot of area against volume for deep-seated landslides (Fig. 5). Owing to such limitations, the empirical relationship between area and volume was calculated only for shallow landslides (51 of the 61 sampled landslides). From these 51 shallow landslides (including debris ows and debris slides), a validation subsample of 10 was used to test the power law model. The remaining 41 shallow landslides with narrow scatter plot dispersion (Fig. 5) were used to determine the empirical relationship. The area (A) and volume (V) of these 41 shallow landslides were used to determine a power law of V = 0.93239514 * A 1.0480 . With the landslide inventory map and related geo-database, the planimetric area of the remaining 465 shallow landslides was obtained. The planimetric area was the input for the empirical relationship areavolume in the power law. Finally, we added the calculated volumes of all shallow landslides to estimate the potential total volume that landslides can contribute per landform unit to the stream system.

Results
The 607 landslides in the inventory could be assigned to 6 types: 295 shallow landslides (48.6%), 126 debris ows (20.8%), 85 debris slides (14%), 77 rockfalls (12.7%), 19 deep-seated landslides (3.1%), and 5 earth ows (0.8%). In this watershed, the abundance and types of landslides are determined by geomorphological conditions in conjunction with land use change: 405 (66.7%) are in volcanic terrain from which 249 (61.48%) landslides are located in pyroclastic deposits where the land use has changed from forest activities to agriculture and 202 (33.3%) are in sedimentary terrain (mostly in folded and faulted Cretaceous limestones covered by broad-leaf and coniferous forest); of the 202 landslides in the sedimentary soils, 123 (60.9%) are in areas that have experience land use changes.
During the assessment of the watershed, 15 homogeneous landform units were identi ed as well as the landslide susceptibility associated with each one and their contribution in terms of volume ( Fig. 6 and Table 3). These landform units are: 1) The relict of volcanic collapse. This is an eroded, fractured, and weathered andesitic collapse structure of the ancestral Pico de Orizaba volcano (0.25 Ma) (Höskuldsson, 1992;Macías, 2005). No landslides have been mapped in this landform unit. It is not affected by the uvial valley landform unit since it has high permeability.
2) The historical lava ows (1545, 1566 and 1613). These andesitic and dacitic lava ows are associated with the end of the effusive period of the present Pico de Orizaba volcano (Höskuldsson, 1992;Macías, 2005). No landslides have been mapped in this landform unit. Seasonal uvial valleys affect 3.4% of its area.
3) The lava ows (Upper Pleistocene). These andesitic and dacitic lava ows (0.09-0.15 Ma) are associated with the second constructional stage of Pico de Orizaba volcano (Höskuldsson, 1992;Macías, 2005). This unit includes 2.1% of the mapped landslides, which contribute with 2.9% of the sediments that are delivered to the stream system. The uvial landform unit dissects this unit and covers 11.4% of its area. Five of its 13 landslides are along the stream system.
4) The lava ow hillsides from Sierra Negra (Pleistocene) are andesitic and dacitic lava ows from the Sierra Negra volcano. These ows were emplaced before those of the ancestral Pico de Orizaba volcano (Höskuldsson, 1992;Macías, 2005). This unit includes 0.3% of the mapped landslides, which contribute 0.4% of the sediments delivered to the stream system. The uvial landform unit dissects this unit and covers 11.4% of its area. Of the two landslides in this unit, one lies along the stream system.
5) The pyroclastic hillslopes unit is correlated with the eruptive sequence of ash and blocks from Loma Grande (8200 ± 120 year B.P.), that have been re-deposited along the highest regions of the watershed (Höskuldsson, 1992;Macías, 2005). No landslides have been mapped in this landform unit. Fluvial valleys affect 20.9% of its area.
6) The pyroclastic ramp landform unit is a volcanic depositional surface with slope less than 15° (Lugo, 2011). The unit is formed by andesitic-basaltic lava that is overlaid by andesitic blocks, lava, breccia, lahars, and tephra (Lower Pleistocene) (Höskuldsson, 1992;Macías, 2005). It includes 2.5% of the mapped landslides, which contribute 2.4% of the sediments delivered to the stream system. The uvial landform unit dissects this unit and covers 14.1% of its area. Five of its 15 landslides are along the stream system.
7) The undifferentiated pyroclastic ramp landform unit is formed by diverse deposits (pyroclastic, scoria, and pumice ows) originating from the successive eruptive phases of Pico de Orizaba volcano (Höskuldsson, 1992;Macías, 2005). This landform unit includes 58.6% of the mapped landslides, which contribute 53% of the sediments delivered to the stream system. The uvial landform unit dissects this unit and covers 18.3% of its area. Of its 356 landslides, 266 are along the stream system.
8) The volcanic plateau landform unit is formed by a sequence of pyroclastic ows with andesitic-dacitic blocks and ashes dated between the years 1860 and 1910 (Höskuldsson, 1992;Macías, 2005). It includes 2% of the mapped landslides, which contribute 1.3% of the sediments delivered to the stream system. The uvial landform unit dissects this unit and covers 10.9% of its area. Ten of its 12 landslides are along the stream system. 9) The basaltic hillslopes landform unit is formed by a sequence of Pliocene and Early Pleistocene basaltic lava ows associated with small cinder cones (Höskuldsson, 1992;Macías, 2005). It includes 1.2% of the mapped landslides, which contribute 0.4% of the sediments delivered to the stream system. The uvial landform unit dissects this unit and covers 10.8% of its area. Six of its 7 landslides are along the stream system.
10) The folded mountains (Upper Cretaceous limestone-shale) landform unit is formed by a sequence of limestones and shales that were formed from Upper Cretaceous reefs (Höskuldsson, 1992;Macías, 2005;Alanis-Anaya, 2017). It includes 18.3% of the mapped landslides, which contribute 13.9% of the sediments delivered to the stream system. The uvial landform unit dissects this unit and covers 13% of its area. Of its 111 landslides, 41 are along the stream system.
11) The folded hillslopes (Lower Cretaceous limestone-dolomite) unit comprises a sequence of limestones and dolomite that were formed from the Lower Cretaceous reefs (Höskuldsson, 1992;Macías, 2005; Alanis-Anaya, 2017). It includes 15% of the mapped landslides, which contribute 26.5% of the sediments delivered to the stream system. The uvial landform unit cuts across this unit and covers 11% of its area. Of its 91 landslides, 19 are along the stream system.
12) The piedmont (3-6°) comprises a sequence of reworked deposits along the base of the highest region of the watershed. No landslide has been mapped in this landform unit. The uvial landform unit cuts across this unit and covers 23% of its area (Alanis-Anaya, 2017).
13) The piedmont (< 3°) comprises a sequence of reworked sedimentary and volcanic deposits along the base of the highest region of the watershed. No landslide has been mapped in this landform unit. The uvial landform unit cuts across this unit and covers 6.1% of its area (Alanis-Anaya, 2017).
14) The uvioglacial ramp unit comprises a sequence of moraine and reworked uvial deposits. No landslide has been mapped in this landform unit. The uvial landform unit cuts across this unit and covers 13.3% of its area (Alanis-Anaya, 2017).
15) The uvial valleys unit is distributed throughout the watershed and is found in 13 of the other 14 landform units (Fig. 5). It is characterized by intermittent rivers with a dendritic pattern that form uvial deposits. This landform unit covers an area of 14.78 km 2 . If it had been taken into account for the landslide susceptibility analysis, the landform would have shown very high landslide susceptibility. The landform unit would cover 14% of the watershed area and would include 58.2% (353) of the mapped landslides, which would contribute 48.7% of the sediments. This landform unit is seriously affected by logging, leaving a small buffer of 50-100 m of vegetation along the hillsides (Alanis-Anaya, 2017).
In the watershed, ~ 74% of streams are rst-order with an average junction angle (the upstream angle between the tangent lines of two intersect streams) of ~ 69°, 21% second-order with a junction angle of ~ 56°, and 4% third-order with 53°. In volcanic terrain most streams follow a nearly N-S radial subparallel pattern whereas sedimentary terrains have a dendritic and subparallel pattern in which second and third-order streams follow a NE-SW direction, along structures and joints.
Landslide susceptibility is very high for one volcanic and two sedimentary landform units (Fig. 6, Table 3). Volcanic landform units cover 54.4 km 2 , of which 47.4% has very high landslide susceptibility (undifferentiated pyroclastic ramp), 14.6% has high susceptibility, 3.3% has moderate susceptibility, and 34.7% has low susceptibility. The remaining 51 km 2 is covered by sedimentary landform units; of this area, 60.2% has very high landslide susceptibility and the remaining 39.8% has low susceptibility.
Dispersion of data in the scatter plot shows that 97.8% (R 2 ) of the volume is explained by the area for shallow landslides, whereas only 48.1% (R 2 ) is explained for deep-seated landslides. That leads to a potential total shallow landslide volume of 1,037,954.2 m 3 .

Discussion
The present study shows that in the study area, some two-thirds of landslides occur in volcanic terrain, and that instability is highest in steep sedimentary followed by steep volcanic hills and valleys. These results agree with landslide distribution and landslide susceptibility tendencies reported with more quantitative techniques (eg. Legorreta-Paulín et al., 2013). The present study de nes 15 landform units to represent volcanic and sedimentary landscapes units. Our analysis shows that one volcanic landform unit and two sedimentary units contribute the most sediment per square kilometer and have the highest landslide susceptibility. Among the landforms, the uvial valley unit covers 14% of the watershed and encompasses more than a half of the landslides. However, since this landform unit is embedded in the other landform units, it is not taken into account for the landslide susceptibility analysis. Also, the result suggests that the combined effect of land use changes, the geological and geomorphological conditions, and the size of the study area covered by very highly susceptible landform units can explain the predisposition and variability of landslide sediment production. For example, among the sedimentary terrains, the folded mountains (Upper Cretaceous limestone-shale) and the folded hillslopes (Lower Cretaceous limestone-dolomite) units had the highest landslide susceptibility. Both have high-density drainage, vertical erosion, and clay rich soils, but there is a difference in the amount of sediment delivered, due in part to having been exposed to different lengths of weathering conditions. This is re ected in the fact that 26.5% of the sediments delivered to the stream system come from the folded hillslopes unit compared to 13.9% of the sediments delivered from the folded mountain landform unit. Another example is the undifferentiated pyroclastic ramp landform unit. In table 3, the landslide density (ratio between number of landslides and the landform area) shows that this landform unit has the highest number of landslides per kilometer (13.8). This landform unit, originally covered by forest, has fertile weathered and unconsolidated soils that have favored a change to large agricultural areas; however, unsuitable treatment of the remaining forest has triggered the high incidence of gravitational processes and erosion. The volume density or normalized volume (ratio between volume delivered by landslides and the landform area) shows that this landform unit contributes 29.6% of the total volume of sediments in the watershed, or 21,322.8 m 3 /km 2 ; of this contribution, almost 75% is derived from landslides along the uvial valley that cuts across this pyroclastic ramp landform. Even though the area of the undifferentiated pyroclastic ramp is ~ 4 and 2 times larger than that of the sedimentary Cretaceous mountains and hillslopes respectively, the sedimentary units rival the volcanic areas in number of landslides per square kilometer. The contribution of sediments from the folded hillslopes (Lower Cretaceous limestonedolomite) unit is the highest in the watershed (24,572.9 m 3 /km 2 or 34.1%), and the folded Cretaceous mountains (Lower Cretaceous limestone-shale) unit contributes 6921 m 3 /km 2 or 9.6%; hence, the sedimentary landform units contribute 43.7% of the sediments by volume whereas the volcanic landform units contribute 56.3%.
It is stressed that calculated sediments delivered per landform unit re ect only a potential sediment volume, because not all of them are incorporated into the streams in one event. During eldwork, it was observed that some landslide deposits were eroded by more than one subsequent rainfall event. Displacement distance and volume storage is determined by local conditions of precipitation, stream slope, landslide volume, and junction angle. Although local conditions are not considered in this analysis, it is observed that in volcanic landforms at high elevations (> 4000 m a.s.l,) with steep topography and resistant materials (volcanic plateau landform and basaltic hillslopes landform unit) the small rst-order streams have larger junction angles (> 65°); as the vertical erosion decreases, streams are nearly parallel downstream with small junction angles. Streams along volcanic landforms have a sub-parallel pattern that cuts deep into loose material across lithologic contacts with a N-S slope trend. In these volcanic areas, landslide deposits are remobilized more frequently. Along the folded Cretaceous sedimentary units, streams show a dendritic NE-SW pattern that may follow fractures parallel to normal faults that affect Cretaceous rocks (Concha-Dimas et al., 2005).
Our analyses also reveal that although no landslides were mapped in 5 volcanic and 3 sedimentary landform units, some units have the potential to develop shallow landslides, deep-seated landslides, and rockfalls due to their steep gradients and loose weathered material.
Although there is a more generic and rapid approach to the de nition and classi cation of landform units based on hillslope gradient and shape, rock type, and landslide density (WADNR Forest Practices Division, 2006;Legorreta-Paulín et al., 2014), the geomorphological method used here was preferable; its focus on geomorphological detail allowed speci c volcanic and sedimentary terrains of Pico de Orizaba volcano to be characterized. However, it may not be possible to compare these landform units with other volcanic areas in Mexico, as can be achieved with the generic approach. We acknowledge that the production of morphometric maps is time consuming and requires expert knowledge in geomorphology. Also, that the subjective evaluation of eld conditions by the researcher will in uence the delineation of landform borders.
Nevertheless, the resulting landslide susceptibility classes allow comparison between landforms because their semiquantitative value is normalized by the landform size. The landform units and the susceptibility classi cation can also be related with the volume delivered and may be useful in the quanti cation, assessment, and modeling of debris ows in the near feature.
In this study, landslide geometry and volume of 61 landslides were determined by eldwork and uses of a UAV. These measurements made it possible to establish an empirical area-volume relationship, expressed as a power law to calculate the volume of landslides for the entire landslide inventory. The power law obtained is in general agreement with existing relationships published in the literature (Martin et al., 2012;Guzzetti et al., 2009); however, several caveats and practical challenges need to be considered for further studies. First, generation of the DEMs used terrain data collected from the UAV that were not error free to generate the DEMs. Few GCPs were obtained for georeferencing photographs, and this is a source of potential errors in the volume calculation. Our results show that error for georeferencing photographs led to a RMSE of < 1 m. This value varies under different landslide sizes and conditions. Topographic displacement errors could lead to miscalculation of volumes, especially for long-term landslide studies where current and pre-existing topography are surveyed in the eld.
In our case the pre-existing topography was derived (reconstructed) exclusively from the current topography and not from two speci c periods of time. This decreases the potential topographic displacement. Owing to the abruptness of the terrain, only one landslide was measured with differential GPS and a GPS rover. A partial result showed that the volume calculated with the UAV overpredicted ~ 4% more than by using differential GPS. Second, cleaning the point cloud from non-topographic elements (trees, bushes, etc.) is also a source of potential errors in the volume calculation. This was partially resolved by using specialized programs of photogrammetry. However, there is still no robust methodology for processing images derived from UAVs. The choice of software is a crucial element in the methodology, and this has a signi cant in uence on the resulting accuracy.
The power law relationship was used to calculate the volume of the 465 landslides not selected for eld measurements.
The total landslide volume was calculated for each landform unit (excluding the uvial valley landform unit) and for the entire watershed (Table 3 and Fig. 6). Overall, the landslides represented a potential contribution of 1,037,954.2 m 3 of material to the main valley. The validation subsample shows that the calculated power law sub-predicts the total volume by 2.3%. On the other hand, use of the 61 sampled landslides (including deep-seated and rockfall) obtains a power law of V = 0.70242506 * A 1.1263 , and this leads to a total volume of 1,547,656.6 m 3 ; however, owing to the high dispersion of data in the scatter plot, only 71% (R 2 ) of the volume is explained by the area of landslides. Our study is in agreement with Broeckx et al. (2016) in showing a signi cant correlation between landslide sediment production and landslide susceptibility. However, we consider that the landslide volume calculation underpredicts the volumes delivered because the power law was calculated and applied only for shallow landslides and it does not include deep-seated slides, rockfalls, and earth ows or potential landslide areas predicted by a model (e.g. MLR model); hence, the calculated volume of sediment is well below the total potential volume that landslides can deliver to the main stream in extraordinary precipitation events. This is in agreement with previous work that shows that the relation area-volume differs signi cantly for very large landslides (Martin et al., 2012;Guzzetti et al., 2009). A large number of samples of rockfalls, deep-seated slides and earth ows are needed to improve the power law and decrease the sub-prediction of landslide volumes.

Conclusions
For large, remote watersheds that are di cult to reach, it is di cult to model landslide susceptibility and the potential volume of material that gravitational processes can deliver, since the topography is continuously undergoing change by landslides, local land use, and environmental conditions. The main aim of this study was to model and assess landslides susceptibility and landslide volume in terms of landform units. The approach uses and adapts the LHZP of WSDNR to map landslides, and to de ne landform units. It also uses aerial photographs and DEMs acquired with UAV surveys and geomorphological evaluation to establish current and reconstructed pre-slide topography of shallow landslides, from which a power law is established to calculate potential total landslide volume.  General procedure to create the DEM and calculate the landslide volume. A) DEM generated with the UAV topographic survey on a landslide. From the inventory map, the headscarp polygon is extracted and overlaid to the DEM (red line). Then the nodes (black dots) from the polygon are extracted to create a point shape le (green dots) which in turn is used to obtain the elevation data from the DEM. B) DEM of the reconstructed pre-slip topography generated by interpolation of point shape le. C) Altimetric changes between pre-and post-sliding surface. GIS algebraic comparison between the two DEMs indicates the landslide altimetric changes. D) Landslide volume map. The landslide altimetric change is then multiplied by the pixel size to derive the volume.