Landslide damming hazard susceptibility maps: a new GIS-based procedure for risk management

A complete landslide dam hazard management incorporates two assessment phases: the damming probability and the breach hazard. A prompt evaluation of the dam stability is crucial during the emergency to mitigate its consequences, but a reliable risk assessment can be realized only after the event has occurred, when the available time is very short. Therefore, it is necessary to develop tools able to help in mapping the spatial probability of damming over large areas for land-use planning, in order to better constrain consequence analysis and risk scenarios for setting up mitigation measures. In this work, a semi-automated GIS-based mapping methodology, based on a statistical correlation of morphometric parameters described by a morphological index, is proposed to spatially assess the likelihood of a river obstruction by landslide damming through two main mechanisms: the reactivation of existing landslides and the formation of new landslides. The two mapping methods (damming predisposition and damming probability) were used on a test area, the Arno River basin in Italy. The Eastern part of the basin resulted as the most susceptible to damming events in the whole basin. These are the highest mountain ridges in the basin (about 1600 m a.s.l.), characterized by calcareous, arenaceous, and marl lithology. The results are confirmed by the high concentration of the known historical landslide dams in the area according to existing inventories.


Introduction
Landslide dams are geomorphologic processes caused by the interference between landslides and river channels. River obstructions are common in mountain regions, where they can cause serious hazards such as upstream backwater formation, catastrophic downstream flooding, channel instability, changes in the river bed dynamics, and triggering of secondary landslides with a cascading effect (Swanson et al. 1986;Costa and Schuster 1988;Casagli and Ermini 1999). Rising impounded water and anomalous flood waves, consequent to dam breach, can cause catastrophic consequences at significant distance in up-and downstream areas, respectively (King et al. 1989;Dai et al. 2005;Chen and Chang 2016). As most of the human activities and population are located in valley floors, consequences of a downstream flooding have significant economic, social, and environmental impacts. Moreover, climate change and overpopulation worsen considerably the impact of natural hazards, such as landslide dams. Population can suffer casualties and restoration costs are often substantial and of two categories: direct (e.g., safety measures and infrastructure rebuilding) and indirect, more difficult to estimate (e.g., damage caused to industrial and agricultural productivity or loss in real estate value delayed in time). It is widely recognized (Costa and Schuster 1988;Ermini and Casagli 2003;Tacconi Stefanelli et al. 2015) that most of landslide dams last a short period of time.
About 40% of the dams collapse within a single day after formation and only about 20% last more than 1 month, so that the available time to assess the dam stability usually is not enough for reliable in-depth analysis and only techniques allowing rapid data collection and analysis are possible. Nevertheless, some consequences from landslide damming can be reduced with mitigation and prevention measures where the expected damming probability is high and the possible consequences catastrophic. Therefore, planning and prevention instruments, as risk and susceptibility mapping, are fundamental to reduce the consequences of natural hazard and increase the cost efficiency of environmental management. Unfortunately, this is undoubtedly a challenging topic, as the comprehension of the phenomena involves complex variables and, as a consequence, no method has been proposed so far to map landslide damming probability at the basin scale.
Landslides dams are often generated by the reactivation of ancient movements triggered in the past during different climatic and environmental conditions (Casagli and Ermini 1999;Canuti et al. 2004;Dikau and Schrott 1999;Borgatti and Soldati 2010;Crozier 2010). In such cases, they are now dormant and vegetated (Rosi et al. 2018) with the strength parameters at the sliding surface close to the residual ones and can be reactivated by natural causes (e.g., earthquake, rainfall, or snowmelt) as well as man-made activity. Therefore, all dormant landslides able to reach a river section along their path can potentially obstruct the stream and should be subject to investigation.
New landslides, instead, may potentially develop wherever suitable conditions are met within hillslopes. The spatial probability of occurrence is usually estimated by landslide susceptibility analysis, based on well-established methods and the Arno river basin is no exception with recent maps which are available (Catani et al. 2005;Catani et al. 2013;Convertino et al. 2013). However, the damming probability is highly dependent on landslide volume, a quantity that is difficult to predict with accuracy. Similarly to other natural hazards (such as meteorite impacts, earthquakes, floods or forest fires), landslides as well occur as stochastic process with a magnitude frequency that follows an exponential distribution decreasing with increasing magnitude Malamud and Turcotte 2006;Clauset et al. 2009). By determining the spatial autocorrelation properties of such empirical statistical distributions, Catani et al. (2016) computed the magnitudefrequency probability of landslide volume as a tool to predict the landslide magnitude probability as a random space variable. The authors propose a tool to map locally the exceedance probability of landslides magnitude and provide a method that can be used to predict the occurrence probability of new landslides with volumes bigger than a certain threshold, necessary to obstruct the valley floor in a river basin.
According to some researches (Swanson et al. 1986;Canuti et al. 1998;Ermini and Casagli 2003;Dal Sasso et al. 2014;Tacconi Stefanelli et al. 2016), landslide dam behavior can be forecasted through the computation of geomorphological indexes, composed by parameters characterizing the involved natural systems: the landslide (or the dam) and the river (or the lake, if present). Geomorphological indexes are a powerful classification and prediction tool but, being mostly empirical, depend on extensive studies and measurement efforts. In most of the cases, such indexes need parameters that are not always available and easy to obtain, like landslide velocity (Swanson et al. 1986;Dal Sasso et al. 2014).
The recently proposed Morphological Obstruction Index (MOI) (Tacconi ) is a bivariate index that requires only simple morphometrical parameters which are easily obtained from common digital elevation models. This tool shows a proven capability in assessing the probability of formation and the stability of landslide dams with respect to other popular indexes (Swanson et al. 1986;Ermini and Casagli 2003). The MOI expression combines two of the most important parameters, the landslide volume V l (m 3 ) and the valley width W v (m): According to MOI, analyzed landslide dams can be classified within three evolutionary domains: formed, not formed, and of uncertain evolution. The limits of these regions are drawn by two straight lines, the "non-formation straight line" and the "formation straight line" (Fig. 1).
The equation of the former is expressed as follows: where V l ' is called the "non-formation volume" and is the minimum landslide volume (m 3 ) able to potentially dam a river a width Wv. Lower volumes do not completely obstruct the river. The expression of the latter is the upper limit for not formed dams and the inferior boundary of the formation domain and is expressed as follows: where V l ", called the "formation volume", is the minimum landslide volume (m 3 ) to have the river valley dammed, with a confidence interval of 99%. Two datasets were used to test the index reliability: an Italian inventory of 300 landslide dams ) and a database from a completely different environment in the Cordillera Blanca, Peru (Tacconi ) with 51 cases. The MOI was able to classify correctly about two-thirds of the database elements.
In this paper, we propose a simple semi-automatic methodology to verify the damming susceptibility at basin scale from existing and neo-formed landslides with geomorphological indexes. These results are achieved employing an available landslides inventory and the occurrence probability method proposed by Catani et al. (2016). The mapping methodology will be applied to the wellstudied basin of the Arno River in Italy, where the mass movements are quite common and a set of data needed for the methodology application were available.

Study area
The Arno River basin is placed in Central Italy and affects territory of the Tuscany Region (98.4%) and the Umbria Region (1.6%) involving 166 municipalities (Fig. 2). The Arno river basin has a surface of 9116 km 2 and, regarding the altitudes distribution of the basin, 55.3% is lower than 300 m a.s.l., 30.4% at altitudes between 300 and 600 m a.s.l., and 14.3% higher than 600 m a.s.l. In the distal part next to the sea, the floodplain is smoothly joined with a wide coastal plain. Flat surfaces of different extensions can be found extensively on morphological high. The areas with steep slopes are extended around the N-NE borders of the basin. The hilly and mountainous areas predominate within the basin with 78% of the total area.
The Arno River basin is placed in the south-eastern side of the Northern Apennines chain. The latter is a fold and thrust belt system made up by the juxtaposition of several tectonic units shifted toward E-NE (Boccaletti and Sani 1998). The orogenic phase was characterized by a compressive phase until the Middle-Upper Miocene. Since Upper Tortonian, the tectonic changes progressively from compressive to extensive, braking the Apennines in a system of structural highs (horst) and tectonic pits (graben) with NW-SE alignment. The latter phase resulted in the emplacement of Neogene sedimentary basins, mainly of marine (to the West) and fluviolacustrine (to the East) origin (Martini et al. 2001). The drainage of the Arno River is strongly conditioned by the chain structures and results in a prevalence of NW-SE trending streams.
The Arno River basin consists mainly of flysch and rocks with prevailing pelitic fraction along the reliefs, and cohesive and granular soils in the hilly basins. Igneous, metamorphic, and calcareous rocks outcrop in limited areas of the basin. These geological settings clearly affect the typology and occurrence of surface processes, primarily through the differences in the mechanical properties of the prevalent lithology.
Just as in all the Italian territory, in the study area, landslides and mass movements are very common. In particular along the main mountain chains, the Apennine, about 30,000 cases are recorded.

Materials and methods
The valley width can be considered a static variable in the MOI equation, since this parameter does not change significantly over decades within each river stretch. From this assumption, according to Eqs. (2) and (3), if we evaluate the average river with W v within each river stretch, two threshold landslide volumes V l ' and V l " (non-formation volume and formation volume) able to block a river can be calculated at each section. V l ' is the minimum volume of formation, below which a landslide does not produce complete obstruction with a 99% of confidence, and V l " is the minimum value above which the dam is definitely formed.
Landslide dams, as all kind of landslides, are often reactivations of ancient movements started in the past. Through an updated landslide database, it is possible to estimate the landslide volumes with some assumptions and simplifications. Landslides with volume bigger than V l ' and V l " for their river section are identified as potentially prone to block the river in that point. Therefore, a "map of the damming susceptibility" for reactivation of existing landslides can be produced.
The prediction of probability for new landslides, with volume bigger than V l ' and V l ", is a more challenging task as the volume is a difficult value to be computed. Several natural processes (earthquakes, meteorite impacts, floods, and forest fires) occur as stochastic events with non-normal magnitude frequency distribution, showing an exponential decrease of frequency with the increase of magnitude (Malamud and Turcotte 2006;Turcotte and Malamud 2004). Mass movements have similar behavior and generally seem to follow a power law magnitude-frequency distribution (MFD), at least for the medium-and large-size events (Guzzetti et al. 2002;Van Den Eeckhaut et al. 2007). Dimensions below a given rollover threshold do not follow this law, probably because they are underrepresented in historical inventories (missed by field surveys, hidden by anthropic activities, vegetation growth, or weathering). Above the rollover dimension, the MFD can be expressed by a power law, even if various sources of errors may produce deviation and noise (Guzzetti et al. 2002;Malamud et al. 2004). Scale-dependent fractal equations can show if the rollover is real, due to fractal behavior, or only and artifact due to mapping resolution or undersampling errors (Stark and Hovius 2001). MFD of landslides from robust historical inventories can be used in natural hazard and risk assessment as a solid basis for the magnitude prediction through a statistical approach. An MDF can be modeled by a power law scaling within the surveyed area (Van Den Eeckhaut et al. 2007). In order to have more information about the spatial variability of the MDF and to represent local pattern of mass movement magnitude in a large geographic region, Catani et al. (2016) computed a set of MFD parameters as random spatial variables in the Arno River basin and spatialized them with geostatistics. The authors found that in the Arno River basin, the landslide volumes are distributed according to a power law scaling for volumes bigger than a cutoff value v min . Starting with the hypothesis that a subset of the whole dataset would be equally representative in that area in statistical terms, Catani et al. (2016) verified that the power law exponent varies according to geographic position. The scaling exponent can be treated as a random spatial function with autocorrelation properties, both at local and regional scale. It is locally stationary only for areas with very homogeneous environmental characteristics (e.g., geology, geomorphology, hydrology, vegetation, and land use). Based on this finding, Catani et al. (2016) propose a simple method to map the power law exponent spatial distribution (Fig. 4, 4), in order to create maps of exceeding probability of landslide volume to be used in risk analysis.
The method is based on the application of the following equation: where P(≥v) is the occurrence probability of a landslide with volume equal or bigger than v, α is the power law exponent of the volume frequency distribution for landslides, and v min is the lower cutoff volume for which the volume frequency distribution follows the power law. For the Arno River basin, this value is equal Fig. 1 Schematic plot of the non-formation line and formation line to 10 4 m 3 . Tacconi  found that the same magnitude v min = 10 4 m 3 as the minimum landslide volume was able to cause a detectable consequence, like a partial damming, on a river dynamic.
Replacing v with the boundary volumes V l ' and V l " in Eq. (4), a simple methodology to get the landslide occurrence probability, with a volume greater than the non-formation value V l ' or the formation value V l ", may be easily obtained. Thus, two "maps of damming probability" of uncertain formation and formation for each stretch of river network can be finally obtained.
The proposed semi-automated procedure has been developed entirely in a GIS (geographic information system) environment. The methodology adopted to obtain the maps of damming susceptibility and damming probability is summarized in the diagram of Fig. 3 and can be briefly described in four main steps in Table 1.
The Arno River basin was selected as the test area for its size and favorable morphological characteristics and because a set of core data, essential for the application of the presented method, were available. These data, exploited as a basis, are ( Fig. 4): 1. An updated landslide database; 2. A digital terrain model; 3. The river network; 4. A map of local alpha distribution ).
The inventory of the Tuscany region was recently updated with the help of persistent scatter interferometry (PSI) (Rosi et al. 2018). This database is the collection of occurrences over a large temporal scale and is an "historical inventory." About 27,500 out of the 91,700 landslides in the new Tuscany database belong to the Arno River basin and have been selected for use in this work. The area of landslides ranges from 1 × 10 2 to 5 × 10 6 m 2 . Many of them (about 90%) are characterized by repeated reactivations and extremely to very slow velocities. Most of the mass movements (about 98%) can be attributed to rotational or planar slow-moving slides. Other less frequent movements, as topples or falls, are rare and located in Fig. 2 Location, elevation distribution, and main river network of the Arno River basin Original Paper areas with specific geological settings and are not taken into account in this work. The typology of landslide movement is conditioned by the morphological and geological settings of the basin. As shown in Fig. 2, the Apennine chain bound the eastern margin of the basin and the area is characterized by earth slides due to a dominance of flysch formations and medium-high slope angles. In the north of the basin, there is a predominance of flow-type landslides for the presence of metamorphic rocks and slope with high angles. In the middle part, shallow landslides are prevalent within cohesive and granular soils. A digital terrain model (DTM) with 10 m of resolution and a vector layer of the actual river network with a total length of about 150,000 km are also employed. A synthetic river network, obtained by using hydrological extraction models from a DTM, was discarded

Main steps
Step description I.Preliminary operations -Flat areas and short river stretch removal; -Division of the river network in stretches of 300 m long.

II.Parameters computation
-Creation of transects perpendicular to the river stretches; -Landscape classification in geomorphological units; -Conferring of Wv to the river stretches; -V l ' and V l " calculation for each river stretch.
III.Damming susceptibility mapping -Draining surfaces creation; -Landslides selection and volume calculation, V l ; -Comparison between V l with V l ' and V l ".
because rivers have suffered throughout history anthropogenic changes in several places, channeling very different paths from the natural hydraulic flow. For the determination of exceeding probability of volumes, a map of the local distribution of power law exponents has been used, as shown below, according to Catani et al. (2016).
Step I: preliminary operations To reduce the processing time and improve the visualization effect, it is recommended to remove a series of unnecessary data. The river blockage takes place almost exclusively in hilly or mountainous areas and preferentially along steep slopes (Costa and Schuster 1988;Fan et al. 2014;Tacconi Stefanelli et al. 2015. For this reason, sections that run in flat areas (with less than 5°slopes and below 150 m a.s.l.) have not been considered in the elaborations. Moreover, short river stretches of the zero and first order have small upstream watershed basins and the resulting ephemeral streams have a complex dynamics for which landslide damming occurs with different time and scales. Also, the consequences of damming would be indeed negligible. According to this consideration, river stretches of the zero and first order with length less than an arbitrary threshold of 20 m were not considered. Furthermore, to make maps easier to display and manage, the river network has been divided into river stretches 300 m long.
Many GIS software allow to easily carry these operations and to outline watersheds and drained areas from a DTM through algorithms of flow direction. Often, however, the main water flow resulted by the algorithm does not match the real one, mainly for the intense human activity that changes the rivers' path in artificial channels. This problem can be overtaken by using a real river network to force the DTM to locally assume minimum values according to the imposed river path.
Step II: W v , V l ', and V l " computation The classic way of landform measurement during field survey or through stereo aerial photos is too time-consuming (for regional scale analysis) and subjective. Also, according to the valley shape (Uor V-shape) and dimension, it is not always possible to identify a sharp change between the valley and the slope in the cross profile.
The analysis of digital terrain has been evolved in the last decades and different algorithms have been developed to automatically extract terrain features using commercial GIS software or stand-alone programs (Guth 1995;Romstad 2001;Tarboton 1997). Methods to automatically extract landform information at broad scales are becoming widely used in geomorphological and natural science studies with the increasing availability of remote sensing data and GIS applications (Drăguţ and Blaschke 2006;Walsh et al. 1998;Wang et al. 2010).
The valley width, such as any other landform, is not an easy parameter to identify and measure. Wood (1996Wood ( , 2009) realized "LandSerf" software, now integrated as a module into SAGA GIS software, designed to automatically classify landforms from digital models. The module derives land-surface parameters from DEMs (i.e., slope, aspect, and curvature), using a multiscale approach, that are used within image processing for pattern recognition and texture analysis. During this processing, the method allows the landscape classification, dividing it into homogeneous morphometric units (peaks, ridges, passes, channels, pits, and planes) (Fig. 5a). This terrain analysis technique follows a common concept in DEM-based landform mapping according to which each discrete landform type has a characteristic combination of elevation derivatives, such as a "morphometric signature" (e.g., Iwahashi and Pike 2007;Pike 1988). Using the module proposed by Wood (2009), it is possible to identify the polygons representing the channels' morphological unit, which can be used as an objective tool to define the valley floor limits over a broad spatial scale (Fig. 5b). The ability to discriminate different geomorphologic landforms is more effective in mountainous areas with strong elevation differences, than in flat areas where the differences between different land forms are less clear. Nevertheless, this does not affect our study because, as already stated, the landslide dams mostly occur in hilly and mountainous areas and therefore, the flat areas can be excluded from the analysis.
A further step is to associate a valley width value, W v , for each 300-m-long river stretch. To measure the distance between the two lateral valley floor boundaries, the river network is sampled creating 500-m-long lines (hereafter "transects"), perpendicular to the river stretches, outdistanced by 20 m (Fig. 6a).
Then, the created valley floor polygons can be used to "cut" the perpendicular transects 500 m long by using a simple cut command in any GIS software (Fig. 6b). The transects' length resulting from the cut is rounded to the lower nearest 10, in order to follow a prudential principle, since narrower valleys are easier to dam. The valley width, W v , of each river stretch is then assigned equal to the median value between the transects' length intersecting it. The median value was considered more significant rather than an averaged value, because the latter would have suffered most of abnormal extreme values due to the often-irregular valley geometry.
Therefore, with the W v value, for each river stretches, the two boundary landslide values of "non-formation volume" and "formation volume," V l ' and V l ", can be easily computed applying the equations of "non-formation" (Eq. (2)) and "formation" straight lines (Eq. (3)).
Step III: damming susceptibility mapping Thanks to an updated landslide polygons archive, it is possible to assess which landslide, if reactivated, is big enough to dam its own valley floor by using the two boundary volumes V l ' (below which a landslide definitely does not produce complete river blockages) and V l " (above which the river valley is certainly dammed).
It is reasonable to assume that a reactivated landslide will move downstream by gravity, following a path like a surface water flow. Draining directions within each slope are easily computed along the river network with the GIS software (Fig. 7). Each landslide can then be associated to the river stretch that it should reach if reactivated, according to the belonging draining surfaces.  Table 2 Comparison between landslide calculated volumes, V l , with the boundary volume of non-formation and formation, V l ' and V l " For the computation of the landslide volumes, two different procedures are followed, depending on the type of movement: one for rotational slides and another for the rest of the movements.
Volumes of rotational slides are calculated, according to a geometrical model assuming a semi ellipsoidal shape, using the equation proposed by WP/WLI (1990), as follows: where D r is the depth of the sliding surface (m), L r the maximum distance between the toe of the sliding surface and the crown of the landslide (m), and W r the maximum distance between the sides of the landslide perpendicular to L r (m). For the remaining landslide types which are not rotational slides, including shallow slides, solifluctions, and flows, a planar slow-moving sliding surface with a constant average depth of the landslide was assumed for simplification since the overall variability from an average value is not significant in areas with similar environmental characteristics (such as climate, geology, tectonics) Cruden and Varnes 1996;Ho et al. 2012;Segoni et al. 2012). The volume of the landslides is obtained with a simple equation, as follows: where the landslide surface, A (m 2 ), is computed automatically with GIS software. A constant average depth of 1.0 m was assumed for landslide depth Dr., compatible with the average soil thickness for shallow landslides in the study areas as empirically achieved through direct measurements by Catani et al. (2010) and Bicocchi et al. (2015Bicocchi et al. ( , 2016Bicocchi et al. ( , 2019. The constant thickness assumption is a simplification empirically achieved through direct measurements to manage a large number of landslides but it is possible only in well-studied areas. In other regions, a specific study on non-rotational landslides is requested to apply this assumption or to select a different approach (e.g., landslide depth related to the total involved area). Each landslide is then classified by assigning two dimensionless values with the simple scheme in Table 2: a value of 2 is assigned if the computed landslide volume, V l , is bigger than the boundary value, V l ' (or V l "), whereas a 0 is assigned if it is smaller. If the boundary value V l ' (or V l ") is bigger than V l but smaller than the V l values increased by 20% (V l × 1.2), then it is assigned a comparison value of 1. Following a cautionary principle, the V l value increased by 20% (V l × 1.2) is used as an arbitrary value to prevent any possible underestimation during parameter sampling and because an increase in landslide body size after the reactivation due to entrainment is possible. The landslide volume computed using the two aforementioned procedures is based on some approximations, since they use geometric simplifications, but it does still reflect the magnitude of the process and therefore is useful for the purposes of the study.
A classification of the damming susceptibility for every mapped landslide is assigned through the combination of the two comparison values in the intensity matrix of Fig. 8. The matrix divides the severity of the damming susceptibility into five classes of a qualitative scale, i.e., very low, low, moderate, high, and very high, colored with dark green, light green, yellow, orange, and red respectively. The gray squares, corresponding to the high V l " values (1 or 2) and lower V l ' values (0 or 1), are not a possible combination, because V l " is always bigger than V l ' according to their formulation.
Step IV: damming probability mapping Although the reactivation of mass movements represents the more common hazard for the damming of a river valley, the study concerning the forecasting of new landslides is not less important. Therefore, a method to get a map of damming probability is introduced by using the occurrence probability equation (Eq. (4)) of a landslide with a volume equal or greater than a v value proposed Fig. 8 Predisposition matrix used to the assignment of the damming predisposition intensity to the mapped landslides Table 3 River stretches susceptibility classification and related probability of uncertain formation (P (≥V l ')) and formation (P (≥V l ")) of a dam  Catani et al. (2016). To solve the probability Eq. (4), three values are needed for each section in the river network: two boundary volumes V l ' and V l ", computed in the " Step II: Wv, Vl', and Vl" computation" section, and the power law exponent of landslides volume frequency distribution, α, for each river stretches.
In the same way as for a reactivated landslide, a new landslide will move downslope by gravity, following the main flow directions. For this reason, using the alpha value distribution map in Fig. 4, 4 , each river stretches were associated with the corresponding alpha values (linked to the probability of landslide occurrence) available within its watershed area (Fig. 7). To highlight maximum possible hazards, the maximum alpha value, α max , was selected for each river stretches.
Replacing in Eq. (4) α with α max and V first with V l ', then with V l ", a new formulation of the occurrence probability equation is proposed, as follows: With this equation, the occurrence probability of a landslide with volume bigger than the boundary values V l ' and V l " can be computed for each river stretches in the Arno River basin. The former probability represents the "uncertain-formation probability" to have a landslide with the minimum volume to potentially dam the river in that point, whereas the latter is the "formation probability" Fig. 9 Map of damming susceptibility of landslides in the Arno River basin by reactivation Fig. 10 Damming predisposition of landslides in the Arno River basin by reactivation Original Paper to have a landslide with a volume able to obstruct, with 99% of confidence, the river stretch for sure.
The resulting probability of non-formation and formation of a landslide dam can be divided into five arbitrary classes as shown in Table 3: very low, low, moderate, high, and very high, colored dark green, light green, yellow, orange, and red respectively.

Results and discussions
The assessment of damming susceptibility on the Arno River basin landslide database is shown in the map of Fig. 9, where it can be observed that most exposed areas to damming potential are the Mt. Morello -Pratomagno and the Mandrioli -Alpe di Catenaia mountain ridges.
In the class distribution shown in Fig. 10, the most frequent class is the very low, with 94.40% of the whole database, followed by the moderate with 4.34% and the remaining percentage divided among low (0.78%), very high (0.47%), and high (0.02%) classes.
Concerning the damming probability caused by new landslides along the Arno River basin, the uncertain-formation and formation probability have a very different geographical distribution as shown in Figs. 11 and 12. The values within the two maps, reported in Fig. 13, have a different distribution as well.
The values of uncertain-formation probability have an almost normal distribution with very recurrent low, moderate, and very low classes, with 36.19%, 28.95%, and 23.57%, respectively, and lower high and very high classes with 9.54% and 1.74%, respectively (Fig. 13a). These two last results mean that about 10% of all the river stretches have more than 15% chances (see Table 3) to be obstructed by a new landslide, resulting in a not formed, formedunstable, or formed-stable dam. All the classes are widespread within the basin with a higher concentration of higher probability around Mt. Morello -Pratomagno and Mandrioli -Alpe di Catenaia mountain ridges (Fig. 11).
The formation probability displays a much clearer classes division and spatial distribution (Figs. 12 and 13b). Almost all (80.82%) the river stretches belong to the very low class, with a probability lower than 0.5% (see Table 3) that a landslide of neoformation has the minimum volume, V l ", able to surely obstruct the river stretch. The remaining 19.18% of river stretches are divided into low (12.74%), moderate (6.11%), high (0.30%), and very high For a validation of the results' reliability, a comparison between the mapping methods and the actually occurred landslide dams in the basin (collected by Tacconi ) has been carried out as shown in Fig. 14. In the map, all the landslides with high and very-high damming susceptibility are displayed as dark red polygons and the areas with higher density distribution of damming formation probability are shown through a statistical interpolation. In the areas with the higher values, we can find the higher concentration of censed landslide dams (8 out of 10) and landslides classified with higher damming susceptibility (68% of them). Few of these landslides and two known landslide dams are placed in areas with moderate/low probability values, due to local morphological and/or geological characteristics.

Conclusion
The costs of reconstruction and losses on both economic and lives due to the consequences of a river obstruction by a landslide are considerably higher compared with the costs of a proper spatial planning and maintenances of the river slopes. A tool capable to delimit the areas with higher risk could drastically reduce the costs, allowing to "punctually" focus maintenance works, monitoring, planning activities, and preventive measures.
The main aim of this research was to propose a useful and easy tool to predict which areas have a higher damming susceptibility from the two possible natural threats: the reactivation of existing landslides and the formation of new landslides. Two simple methodologies were developed applying the Morphological Obstruction Index, resulting in a useful forecasting and planning tool. These methodologies can assess, quickly and with a few data, the damming susceptibility of each river stretch, connected to existing landslides, and the probability of obstruction by new landslides along a river network. The proposed methodology was used on a test area, the Arno River basin, and resulted in practical and realistic maps of the spatial distribution of the susceptibility to obstruction along the river course which are well fitting with known data on past damming.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.