Regional debris flow susceptibility analysis in mountainous peri-urban areas through morphometric and land cover indicators

A method for assessing regional debris flow susceptibility at the watershed scale, based on an index composed of a morphometric indicator and a land cover indicator, is proposed and applied in 106 peri-urban mountainous watersheds in Bogotá, Colombia. The indicator of debris flow susceptibility is obtained from readily available information common to most peri-urban mountainous areas and can be used to prioritise watersheds that can subsequently be subjected to detailed hazard analysis. Susceptibility is considered to increase with flashiness and the possibility of debris flows occurring. Morphological variables recognised in the literature to significantly influence flashiness and occurrence of debris flows are used to construct the morphometric indicator by applying principal component analysis. Subsequently, this indicator is compared with the results of debris flow propagation to assess its capacity in identifying the morphological conditions of a watershed that make it able to transport debris flows. Propagation of debris flows was carried out using the Modified Single Flow Direction algorithm, following identification of source areas by applying thresholds identified in the slope–area curve of the watersheds. Results show that the morphometric variables can be grouped into four indicators: size, shape, hypsometry and (potential) energy, with energy being the component that best explains the capability of a watershed to transport debris flows. However, the morphometric indicator was found to not sufficiently explain the records of past floods in the study area. Combining the morphometric indicator with land cover indicators improved the agreement and provided a more reliable assessment of debris flow susceptibility in the study area. The analysis shows that, even if morphometric parameters identify a high disposition to the occurrence of debris flow, improving land cover can reduce the susceptibility. However, if favourable morphometric conditions are present but deterioration of the land cover in the watershed takes place, then the susceptibility to debris flow events increases. The indicator of debris flow susceptibility is useful in the identification of flood type, which is a crucial step in flood risk assessment especially in mountainous environments, and it can be used as input for prioritisation of flood risk management strategies at regional level and for the prioritisation and identification of detailed flood hazard analysis. The indicator is regional in scope, and therefore it is not intended to constitute a detailed assessment but to highlight watersheds that could potentially be more susceptible to damaging floods than others in the same region.


Introduction
Appropriate recognition of hydrogeomorphic hazards in mountain areas is crucial for risk management, since it provides the basis for more detailed studies and for the development of appropriate risk management strategies (Wilford et al., 2004;Jakob and Weatherly, 2005;Welsh, 2007).Besides the identification of the flood potential, it is important to distinguish between debris-flow-and non-debris-flowdominated watersheds since these constitute very different hazards.
There are several definitions for hydrogeomorphic processes.Wilford et al. (2004) distinguishes among floods, debris floods and debris flows, with sediment concentrations of 20 and 47 % as upper limits for floods and debris floods respectively.Santangelo et al. (2012) and Costa (1988) Published by Copernicus Publications on behalf of the European Geosciences Union.
differentiate water floods as Newtonian, turbulent fluids with non-uniform concentration profiles and sediment concentrations of less than about 20 % by volume and shear strengths less than 10 N m −2 ; hyperconcentrated flows as having sediment concentrations ranging from 20 to 47 % by volume and shear strengths lower than about 40 N m −2 ; and debris flows as being non-Newtonian visco-plastic or dilatant fluids with laminar flow and uniform concentration profiles, with sediment concentrations ranging from 47 to 77 % by volume and shear strengths greater than about 40 N m −2 .On the other hand, FLO-2D Software (2006) uses the terms mudflow (non-homogeneous, non-Newtonian, transient flood events), and mud flood (sediment concentration from 20 to 40-45 % by volume).Despite the variety of definitions, the characteristics of debris flows imply different hazard conditions from those related to clear-water floods, with debris flows being potentially more destructive.The higher destructive capacity is related to a much faster flow and higher peak discharges than those of a conventional flood -as well as high erosive capacity with the ability to transport large boulders and debris in suspension and the generation of impact forces comparable to rock and snow avalanches (Welsh, 2007;Santangelo et al., 2012).With a lower sediment concentration, debris floods and hyperconcentrated flows as presented by Wilford et al. (2004) and Santangelo et al. (2012) are less hazardous, since they carry less of the large boulders responsible for impact damage in debris flows and flow velocities are usually lower.They are, however, considered more dangerous than clear-water floods of similar magnitude (Welsh, 2007).Previous research on the identification of flood potential and areas susceptible to debris flows used quantitative methodologies such as logistic regression and discriminant analysis in addition to GIS and remote sensing technologies (De Scally and Owens, 2004;De Scally et al., 2010;Wilford et al., 2004;Rowbotham et al., 2005;Chen and Yu, 2011;Santangelo et al., 2012;Crosta and Frattini, 2004;Griffiths et al., 2004;Kostaschuk, 1986;Patton and Baker, 1976;Bertrand et al., 2013).These studies focused on the identification of basins or fan parameters to classify them according to their dominant hydrogeomorphic processes.A conclusion from these studies is that drainage basin morphology is an important control of fan processes (Crosta and Frattini, 2004) and that there are significant differences in morphometric characteristics between basins where the dominant process is debris flows and those mainly dominated by fluvial processes (Welsh, 2007).Morphometric parameters such as the basin area, Melton ratio and watershed length have been identified by several authors as reliable predictors for differentiating between debris-flow-and non-debrisflow-dominated watersheds and their respective fans (Welsh, 2007).However, the results of the analyses seem to be highly dependent on the geographical area where the methodology is applied, and in many cases the identification of morphometric parameters requires a previous independent classification of the watersheds, normally entailing stratigraphic observations, detailed fieldwork, aerial photo analyses and calculations.
When historical data on the occurrence of flash floods and debris flows are not available, the recognition of hydrogeomorphological hazards can be carried out through fieldwork analysis applying methods such as the one proposed by Aulitzky (1982) based on hazard indicators, or through stratigraphic evidence in conjunction with age control (Jakob and Weatherly, 2005;Giraud, 2005).However, such fieldwork and detailed geological and geotechnical analysis at the regional scale require significant resources and time, and may not be practicable in the extensive peri-urban areas of cities in mountainous areas such as those in the Andean Cordillera.Furthermore, urbanisation processes in the peri-urban areas of these cities make geologic investigation difficult.Moreover the history of the watershed may not be a conclusive indicator of current hazard conditions, since anthropogenic intervention can play a significant role in the hazard dynamics.This calls for a more rapid yet reliable assessment of the watersheds, allowing a prioritisation of watersheds where a more detailed analysis based on field data is to be carried out.
This paper proposes a method for regional assessment of debris flow susceptibility under limited availability of data in urban environments, where flash floods occur as debris flows, hyperconcentrated flows or clear-water flows as defined by Costa (1988).The proposed index is based exclusively on information derived from digital elevation models (DEMs) and satellite images to overcome the limitation often found in the availability of previous geological work such as stratigraphic analysis and fieldwork for large areas.
The ability of morphometric variables to identify debrisflow-dominated basins was tested.Morphometric variables and land cover characteristics were considered as factors that influence flood hazard and were combined in an index that can be interpreted as the potential susceptibility to which watersheds are prone, including the spatial differentiation of the dominant type of hazard.A key aspect of the index is the discrimination between debris flow and clear-waterflood-dominated watersheds in order to understand the level of threat that floods in the watersheds pose, and to support prioritisation of watersheds to be subjected to further detailed study.
The study area is the mountainous area surrounding the city of Bogotá (Colombia), where an accelerated urban process has taken place during the last decades, forming a periurban area mostly characterised by illegal developments.To overcome the lack of historic records and the infeasibility of carrying out detailed geologic fieldwork for the identification of hydrogeomorphological processes that allow validation of the susceptibility index, results are compared with an independent method based on the propagation of debris flows using a digital elevation model as well as with the few available flood records in the area. 2 Methods and data

Study area
This research focuses on the mountainous watersheds surrounding the city of Bogotá, the capital and economic centre as well as the largest urban agglomeration of Colombia, with an estimated 7.4 million inhabitants.The city is located in the Andean region (see Fig. 1).Several creeks drain the steep mountains surrounding the city and cross the urban area to finally drain into the larger Bogotá River.In this analysis the watersheds that drain into the mainstream of the Tunjuelo River basin, one of the largest tributaries of the Bogotá River, as well as the watersheds in the eastern hills were considered.This includes 66 watersheds in the Tunjuelo River basin and 40 in the eastern hills of Bogotá (see Fig. 1a).These are characterised by a mountainous environment with areas ranging between 0.2 and 57 km 2 .The area is mainly formed by sandstones of Cretaceous and Palaeogene age.This sedimentary rock forms surrounding mountains up to 4000 m in altitude, thus reaching to some 1500 m above the level of the high plain of Bogotá (Torres et al., 2005).The mean annual precipitation varies from 600 to 1200 mm in a bimodal regime with rainy seasons of April-May and October-November (Bernal et al., 2007).
In the study area, flooding is controlled by climate and physiography.However, land use practices and anthropogenic influences have increased flood risk not only through soil and land cover deterioration but also through an intensive occupation of floodplains.In the southern mountains of Bogotá, which belong to the Tunjuelo River basin, the urban and industrial growth has been accelerated from the 1950s.Between 1951 and 1982, the lower basin of the Tunjuelo River was the most important area for urban development.It was settled by the poorest population of Bogotá (Osorio, 2007), and its growth has been characterised by informality and lack of planning.The most devastating floods in Bogotá have occurred in the lower Tunjuelo River basin, involving not only the mainstream but also the tributary creeks where flash floods have caused human losses (DPAE, 2003a, b, c).The watersheds located in the eastern hills of Bogotá have different characteristics since most of the area corresponds to protected forests.However, informal urbanisation takes place in some areas.Additionally, mining has been common both in the Tunjuelo River basin and in the eastern hills, causing the deterioration of the environmental conditions of the watersheds.
The reaches of the creeks in the urban areas have been subjected to significant intervention and occupation.Most of the creeks in the eastern hills drain into the storm water system through structures with low hydraulic capacity (less than the return period of 10 years) (Hidrotec, 1999).The streams that drain into the Tunjuelo River have been severely modified, mainly in the reaches near the confluence, albeit without a comprehensive flood management plan.The flood control structures in the study area have been constructed in the mainstream of the Tunjuelo River, including a dry dam in the middle basin, three retention basins in the lower basin and levees.Additionally, there are two reservoirs in the upper basin of the Tunjuelo River that supply water to Bogotá.Conversely, flood control works have not been constructed in any of the watersheds in the study area except for the Chiguaza watershed, where levees and channelisation works were constructed in the confluence with the Tunjuelo River in 2008.
For the purposes of this study the area formed by the Tunjuelo River basin and the eastern hills will be considered as comparable, due to the lithologic composition and homogeneous flood management policy in the city.
In order to test the performance of the proposed morphometric indicator, additional to the watersheds in the study area, a subwatershed of the Chiguaza Creek in the Tunjuelo River basin (see Fig. 1a) was analysed, as well as two areas external to the study area: La Negra Creek and La Chapa Creek (see Fig. 1b and c).These additional watersheds were included in the analysis given the availability of records and previous studies.
La Negra Creek is located 67 km to the northwest of Bogotá.Among several records of inundation events, the most critical occurred on 17 November 1988.The flow moved along La Negra Creek to its confluence with the Negro River affecting an important area of the municipality of Utica (UNAL and INGEOMINAS, 2007).According to UNAL and INGEOMINAS (2007) the characteristics of this event reveal a concentration of sediment of 40 % by volume, which corresponds to the upper limit of the mud flood category according to FLO-2D Software (2006) or to a debris flood according to Wilford et al. (2004).
La Chapa Creek was chosen due to the high frequency of debris flows.Chaparro (2005) notes that La Chapa Creek is prone to debris flows characterised by the mobilisation of granular material of varying size, ranging from boulders of several metres in diameter to sand, embedded in a liquid phase -formed by water, fine soils and air -sometimes accompanied by vegetal material.The most recent event that took place in this watershed was recorded on video, which allowed the type of flow that dominates the watershed to be confirmed.

Methodology
Variability in the level of hazard is reflected in the proposed susceptibility index, where high values represent a higher potential for debris flow and therefore an increased hazard condition.Moreover, flashier conditions, which result from unfavourable morphometric and land cover conditions, contribute to high values of the index, providing an indication of potential for flash flood danger in a large area.The proposed index to represent the level of debris flow susceptibility at regional scale is composed of a morphometric indicator and a land cover indicator.
The units of analysis correspond to the watersheds delineated up to the confluence with the Tunjuelo River in the case of the streams located in the Tunjuelo River basin, and to the confluence with the storm water system in the case of the streams located in the eastern hills.
In order to develop the susceptibility index and identify whether it is appropriate, a methodology that can be divided into three stages was followed.The first stage addresses the development of the morphometric indicator, the second stage corresponds to the development of the land cover indicator and the third stage is the development of the susceptibility index.Figure 2 shows the main steps that were carried out to obtain the susceptibility index for the study area.
For stage 1, a model to calculate a morphometric indicator was developed by using principal component analysis on morphometric parameters that have been identified in the literature as important descriptors of flood potential and debris flow discriminators.Due to the poor availability of historical records in the study area, which can limit the validation of the proposed indicator, three independent methods were used to assess the appropriateness of the morphometric indicator (methods i, ii and iii in Fig. 2).The first method identifies debris flow source areas using two criteria (a and b in Fig. 2) and propagates the flow on a DEM using two angles of reach (ratio between the elevation difference and length from the debris flow initiation point to the downstream extent of the debris flow runout) (Horton and Jaboyedoff, 2008;Kappes et al., 2011), in order to identify the capacity of watersheds to transport potential debris flows to their fans.The binary result of the propagation reaching or not reaching the fan was used to classify the watersheds.The distribution of the values of the morphometric indicator and its component indicators was analysed grouping the values according to the classification obtained from the propagation results.Furthermore, a contingency table and its proportion correct (PC, fraction of watersheds that were correctly identified by the morphometric indicator) were used to establish the correspondence between the morphometric indicator and the classification from the propagation modelling to assess the skill of the morphometric indicator to identify the potential capacity of the watersheds to propagate debris flows.
In order to compare the morphometric indicator with field data, method ii was used (see Fig. 2).A flow type classification of 11 watersheds was carried out on the basis of the available studies, reports and the flood record database managed by the municipality.The flood record database contains 55 flood events from 2001 to 2012.Due to the short period of record of the database, robust frequency analysis is not feasible.Moreover, flood records are less frequent in the eastern hills and non-existent in the upper Tunjuelo River basin, which may be due to the low density of population in this latter area.However, the data contained in the database normally describe affected people, type of flow and damage, and provide relevant recent historical information on the type of hydrogeomorphic processes that take place in the watersheds.Watersheds where reports, studies or flood records clearly identify the occurrence or imminent possibility of debris flows were classified as debris flow watersheds (D), watersheds where a significant sediment concentration was identified in the past floods were classified as hyperconcentrated flow watersheds (H) and watersheds where the available reports describe the occurrence of floods without description of sediment sources and sediment concentration were classified as clear-water-flow watersheds (C).
The correspondence between the morphometric indicator and the classification obtained from flood records, studies and reports in the study area was assessed through a contingency table.
The two contingency tables (morphometric indicator vs. propagation classification and morphometric indicator vs. flood record classification) allowed assessing the representativeness of the indicator in terms of debris flow threat level.
Additionally, the morphometric indicator was calculated for two external watersheds and a subwatershed of the Chiguaza Creek in the study area.This constitutes method iii in Fig. 2. Since information on the dominant processes of these watersheds is available, they were used to assess the applicability of the indicator outside the study area in the first case and to add a valuable information to the analysis of the study area in the second case.
A qualitative indicator of land cover was developed in stage 2, which was combined with the morphometric indicator through a classification matrix and assessed through contingency tables in stage 3 (see Fig. 2).
The main input for the methods is a 5 m resolution raster DEM.This was constructed using contours that in the periurban area are available at intervals of 1 m.The contours were processed to obtain a triangulated irregular network that was subsequently transformed into a raster through linear interpolation.The details of each stage of the process are described in the following subsections.

Development of the morphometric indicator
Morphometric parameters used in the literature (see Table 1) were extracted for each watershed from the digital elevation model of the study area using ArcGis, SAGA and R functions.Many of the variables as listed in Table 1 are closely correlated.To reduce the dimensionality of the data set, principal component analysis was applied.A reduction of the variables is achieved by transforming the original variables to a new set of variables, the principal components, which are uncorrelated and which are ordered according to the components that retain most of the variation present in the original set of variables 2.2.1.These transformed variables were subsequently used to obtain the morphometric indicator.
Differentiation of debris flow watersheds or fans from those dominated by clear-water floods has been carried out by several authors, finding that morphometric variables are very valuable as discriminators of processes in watersheds (De Scally and Owens, 2004;De Scally et al., 2010;Wilford et al., 2004;Rowbotham et al., 2005;Chen and Yu, 2011;Santangelo et al., 2012;Crosta and Frattini, 2004;Griffiths et al., 2004;Kostaschuk, 1986;Patton and Baker, 1976;Jackson et al., 1987;Bertrand et al., 2013).On the other hand, research on the relationships between watershed characteristics and peak flood and flashiness has contributed to identifying morphometric variables that can help to describe the characteristics of the hydrologic response of a watershed (Patton, 1988).Table 1 summarises the morphometric variables that have been identified in the literature as appropriate discriminators of processes and descriptors of the hydrologic response of watersheds and that were chosen for the analysis.The parameters that are most commonly found in the literature as important discriminators of hydrogeomorphic processes are the area, the slope and the Melton ratio.However, other parameters, such as those derived from the hypsometric curve and the average of the multiresolution index (MRI Gallant, 2003), have also been included given their importance in the description of the evolution and erosion processes of watersheds in the case of the former and the description of the erosion areas in the case of the latter.
The hypsometric curve and the hypsometric integral are non-dimensional measures of the proportion of the catchment above a given elevation (Willgoose and Hancock, 1998).The hypsometric curve describes the landmass distribution and thus the potential energy distribution within the basin above its base (Luo and Harlin, 2003).This curve can be seen as an exceedance distribution of normalised elevation, where the probability of exceedance is determined by the portion of the basin area that lies above the specified elevation (Huang and Niemann, 2008).The hypsometric integral is defined as the area below the hypsometric curve.Values near to 1 in the hypsometric integral indicate a state of youth and are typical of convex curves.Nevertheless, mature s-shaped hypsometric curves can present a great variety of shapes but have the same hypsometric integral value (Pérez-Peña et al., 2009).In order to analyse the hypsometric properties of the watersheds, the procedure described by Harlin (1978) was used: the hypsometric curve was treated as a cumulative distribution function.The second, third and fourth moments were derived about the centroids, yielding measures of skewness and kurtosis for the hypsometric curves, which are represented by a continuous third-order polynomial function.
The multiresolution valley bottom flatness index (MRI) is obtained through a classification algorithm applied at multiple scales by progressive generalisation of the DEM combined with progressive reduction of the slope class threshold.The results at different scales are then combined into a single index.The MRI utilises the flatness and lowness characteristics of valley bottoms.Flatness is measured by the inverse of the slope, and lowness is measured by a ranking of the elevation with respect to the surrounding area.The two measures, both scaled to the range 0 to 1, are combined by multiplication and can be interpreted as membership functions of fuzzy sets (Gallant, 2003).
From the principal component analysis of the morphometric variables, the factor loadings, which represent the proportion of the total unit variance of the indicator which is explained by the principal component, were used to construct the weights of the indicators (OECD, 2008).In order to develop an overall morphometric indicator, the individual indicators obtained from the principal components were combined using as weights the variability explained by each principal component.
The appropriateness of the morphometric indicator to capture the level of debris flow susceptibility was assessed through its comparison with the debris flow propagation capacity of the watersheds; with the classification of 11 watersheds from available detailed studies and historic information; and through the analysis of the indicator obtained in two watersheds outside the study area and one subwatershed of the study area where debris flows have been confirmed.For the first two analyses, contingency tables were used; for the third, direct comparison of the values of the indicator was carried out.
Several hypotheses have been formulated to explain mobilisation of debris flows, and this aspect represents an active research field.The triggering mechanisms and the causal relationships are, however, still partially unknown (Salvetti et al., 2008).Approaches for the identification of debris source areas include the use of credal networks (Antonucci et al., 2007), the use of indices for predisposition factors to assess debris flow initiation hazard (Bonnet-Staub, 2000), empirical relationships (Baumann and Wick, 2011;Horton and Jaboyedoff, 2008;Blahut et al., 2010), the Melton ruggedness number (Rengifo, 2012) and the use of the slope versus area diagram as a topographic signature of debrisflow-dominated channels (Santos and Duarte, 2006).Two of these approaches to identify potential debris flow initiation points will be used in this paper for method i in Fig. 2. The first approach is based on the analysis of the break in the slope versus drainage area relationship, while the second uses an empirically determined critical condition in this relationship (Horton and Jaboyedoff, 2008).In both cases, the debris flow propagation areas were obtained through a propagation algorithm by considering two angles of reach (ratio between the elevation difference H and length from the debris flow initiation point to the downstream extent of the debris flow runout L) (Horton and Jaboyedoff, 2008;Kappes et al., 2011).
Regarding the first method of identification of debris source areas, the slope-area diagram is the relationship between the slope at a point versus the area draining through that point.It quantifies the local topographic gradient as a function of drainage area.Several authors have found a change in the power law relationship (or a scaling break) in slope-area data from DEMs at the point that the valley slope ceases to change below a certain drainage area.This has been inferred to represent a transition to hillslope processes and has been interpreted as the topographic signature for debris flow valley incision (Stock and Dietrich, 2003;Montgomery and Foufoula-Georgiou, 1993;Seidl and Dietrich, 1993).The same conclusion was made by Tucker and Bras (1998), explaining that different processes have an impact on the slope-area relationship, suggesting the possibility that slope-area data may be used to discriminate between different geomorphic process regimes.
Two distinct regions of the slope-area diagram are typically observed.Small catchment areas are dominated by rainsplash, interrill erosion, soil creep or other erosive processes that tend to round or smooth the landscape.As the catchment area becomes larger, a break in gradient of the curve occurs.This is where slope decreases as catchment area increases.This region of the catchment is dominated by fluvial erosive processes that tend to incise the landscape (Hancock, 2005).
The slope-area curve was constructed for two regions of the study area, corresponding to the Tunjuelo River basin and to the eastern hills of Bogotá.The break in the slope-area diagram was obtained using segmented regression, in order to determine a threshold to differentiate two regions, one dominated by erosive processes and the other dominated by fluvial erosive processes.This threshold will be used as the topographic signature of debris flow.
In the case of the second method that was applied to identify debris flow sources, corresponding to the procedure proposed by Horton and Jaboyedoff (2008), this applies criteria based on area, slope, curvature, hydrology, lithology and land cover.The slope criterion to identify debris flow source areas is based on the relationship between slope and drainage area shown in Eqs. ( 1) and ( 2), where β lim is the threshold slope in degrees and S U A is the upstream area in square kilometres.These equations were built on observations made by Rickenmann and Zimmermann (1993).Horton and Jaboyedoff (2008) denominated these criteria as threshold for extreme events given that the 1987 events on which the threshold is based were considered as extraordinary, and this denomination allowed differentiation from other set of points used by Horton and Jaboyedoff (2008). (1) In the method by Horton and Jaboyedoff (2008) every point located above the limits defined by Eqs. ( 1) and ( 2) is considered as critical.In the application of the method in Argentina (Baumann and Wick, 2011), the equations were bounded between 15 and 40 • since in the observations made by Rickenmann and Zimmermann (1993) in Switzerland all the triggering areas slope angles were below 40 • and contributing areas inferior to 1 ha were not considered as potential sources.Thus, the parameters used for detection of triggering areas are slopes in the range of 15-40 • , contributing areas superior to 1 ha and plane curvatures inferior to −0.01/200 m −1 under the condition that the point is located above the limit defined by Eqs. ( 1) and (2).
Using the threshold obtained from the analysis of the slope-area curve together with the criteria of curvature and minimum drainage area as proposed by Horton and Jaboyedoff (2008), initiation points were identified in the study area, constituting the first method of identification of debris flow initiation points.As a second method, the threshold of extreme events was used as a criterion for slope and area, and in addition the minimum area and curvature were used as recommended by Horton and Jaboyedoff (2008).The Modified Single Flow Direction (MSF) model (Gruber et al., 2009) was used to identify the areas that potentially could be affected by debris flows for the two groups of initiation points.The MSF is based on the Single Flow Direction (D8) algorithm and other standard functionalities of ArcInfo/ArcGIS to account for flow spreading allowing the flow to divert from the steepest descent path by as much as 45 • on both sides.The only required inputs are the source areas and a DEM.For a detailed explanation on the MSF algorithm see Gruber et al. (2009).As a stopping condition the MSF algorithm uses the angle of reach.The trajectory component of the MSF model usually provides the potential maximum inundation zones of a mass-movement event.Thus, it indicates which areas are more or less likely to be affected.However, the runout distance should also be based on a maximum.A reasonable angle of reach (H /L ratio) has to be evaluated on the basis of empirical data for the type of mass movement that is being modelled.Several efforts have been made to develop relationships to estimate the angle of reach, mainly using the volume of the debris flow.The minimum angle of reach that has been observed is 6.5 • (ratio H /L = 0.11) (Prochaska et al., 2008), and the highest and more repetitive is 11 • (ratio H /L = 0.19) (Rickenmann, 1999;Huggel and Kääb, 2003;Rickenmann and Zimmermann, 1993;Kappes et al., 2011).The two angles were used to test the sensitivity of the results, but larger and more fluid debris flows may show lower H /L ratios and consequently a larger flow reach.
Watersheds where the propagation area reaches the mouth of the drainage area using a ratio H /L of 0.19 are classified as debris-flow-dominated and labelled "0.19 H /L". Watersheds where the propagation reaches the mouth for a ratio H /L of 0.11 will be considered debris-flow-dominated as well, albeit with a more fluid flow.These are labelled "0.11 H /L". In this classification no distinction between hyperconcentrated flows and clear-water floods is made.Therefore, watersheds where the propagation area does not extend to the mouth of the drainage area will be classified as clearwater-flood-dominated.
In order to assess the validity of the MSF algorithm, the debris propagation results were compared with the extent of a well-documented debris flow event that occurred in the study area.

Development of the land cover indicator
The land cover indicator was constructed by analysing the characteristic land cover of each watershed, which was obtained from the classification of a LANDSAT Thematic Mapper 5 (TM5) image taken in 2001.The LANDSAT image was classified using a supervised classification algorithm.The Table 1.Morphometric variables used in the analysis.Note that L corresponds to the length of the streams in a watershed; H max and H min correspond to the highest and lowest elevation in a watershed respectively.

Variable Relevance Reference
Area (A) Correlated with discharge; proportional to sediment storage in the catchment; wide basins collect a large amount of water, which can dilute the flood, reducing the probability of debris flow forming in a first-order basin reaching the alluvial fan.Correlated with other morphometric parameters.Crosta andFrattini (2004), Baker (1976), De Scally and Owens (2004), Gray (1961), Shreve (1974) Perimeter (P ) Base for several watershed shape indices.Zavoianu (1985) Drainage Circularity coefficient The more circular a watershed, the sharper its hydrograph; this means the flashiness increases and therefore the threat of flooding is higher.

De Matauco and Ibisate (2004)
Elongation ratio E = 2/(L wshd (A/π) 0.5 ) Floods travel less rapidly; have less erosion and transport potential; and have less suspended load in elongated watersheds.
An elongated shape favours a diminution of floods because tributaries flow into the mainstream at greater intervals of time and space.Zavoianu (1985) Watershed width W wshd = A/L wshd Related to the size of fans Weissmann et al. (2005) Length-to-width ratio (LW) Measure of elongation Zavoianu (1985) Melton ratio M = (H max − H min )/A 0.5 Frequently used to discriminate among hydrogeomorphologic processes.
Welsh and Davies ( 2010 The classification identified areas covered by forests, grass, paramo vegetation1 , urban soil and water.From the land cover composition of each watershed a qualitative condition was derived. The natural susceptibility of a catchment to debris flow hazards due to geological, morphological and climatic predispositions can be enhanced by human activities and the effects of land use changes (Koscielny, 2008).In order to include this influence in the susceptibility analysis, the percentage of vegetation cover, urban area and bare soil was used to qualify the state of the watersheds.
Vegetation cover has been recognised as one of the factors related to frequency of debris flows (Jakob, 1996).Forests reduce hydrogeomorphic hazards since they retain organic and inorganic material; contain the transport of mobilised material, reducing the extent of destruction; intercept precipitation; and the stems of trees reduce the areas disturbed by snow avalanches, rockfalls, floods, debris floods and debris flows (Sakals and Innes, 2006).Runoff can be increased by deforestation, soil properties degradation and impervious surfaces construction (Koscielny, 2008) as a result of urbanisation.Likewise, erosion processes and slope instabilities can occur (Koscielny, 2008).The percentage of bare soil represents areas prone to erosion and normally associated with quarries that can provide a supply of sediment.
According to Schueler (1995) stream degradation occurs at approximately 10-20 % total impervious area.The increase in frequency and severity of floods due to imperviousness produces an increase in stream cross-sectional area.This occurs as a response of the stream accommodating higher flows through widening of the stream banks, downcutting of the stream bed or both.The channel instability triggers streambank erosion and habitat degradation.With respect to flood magnitude, this can be increased significantly by percentages of impervious cover larger than 10 %.Hollis (1975) found that peak flows with recurrence intervals of 2 years increased by factors of 2, 3 and 5 with 10, 15 and 30 % impervious development respectively.A threshold of 15 % was used to consider a high condition of urbanisation of the watersheds and therefore a high degree of degradation.In order to consider the degree of degradation related to bare soil, normally related to quarries in the study area, a threshold of 10 % was used.

Development of a composite susceptibility index
The resulting indicators of land cover and morphometry were combined using a matrix that allows classification of the catchments into high, medium and low susceptibility.Figure 3 shows the initial matrix used for the analysis.The corners corresponding to poor land cover and high morphometric indicator and good land cover and low morphometric indicator (cells a and i) were assigned a high and low susceptibility respectively, since they correspond to the extreme conditions in the analysis.The cells from b to h in Fig. 3 were considered to potentially correspond to any category (low, medium or high priority), and all the possible combinations of the matrix were tested assessing the proportion correct of a contingency table comparing the obtained susceptibility index and the classification of flow type from the flood records, where debris flows were considered the most hazardous type of events.Potentially 2187 combinations can be obtained by assigning the three susceptibility categories to cells b to h in the matrix shown in Fig. 3.Even though some combinations of the categories are not consistent with a progressive increase of susceptibility level from the bottom right corner of the matrix to the top left corner, all of them were tested.Under this procedure, the resulting matrix corresponds to the best fit of the susceptibility index and the classification of flow from flood records.

Results
The results obtained for each stage of the process are presented in the following subsections.The first subsection presents the results of the estimation of the morphometric indicator for the study area.This subsection includes the development of the morphometric indicator model based on the principal component analysis and the assessment of the appropriateness of the morphometric indicator.The latter covers the classification of watersheds according to the debris flow propagation capacity and the comparison of the morphometric indicator, with the propagation of debris flows described using the MSF model, with the 11 watersheds with confirmed flow type in the study area and with three additional watersheds with confirmed flow type outside the study areas.The second subsection shows the results of the development of the land cover indicator, and finally the third subsection shows the results of the combination of the morphometric indicator and the land cover indicator to obtain a final susceptibility index.

Morphometric indicator model
The results of the principal component analysis applying a varimax rotation carried out on the morphometric variables are shown in Table 2. From the Scree tests carried out on the eigenvalues obtained from the principal component analysis, the amount of principal components to be used was found to be four.These first four principal components account for 85 % of the variance in the data.From the analysis four groups of variables could be identified related to the size (inversely proportional to area), shape (proportional to circularity), hypsometry (proportional to hypsometric integral) and potential energy (proportional to the Melton number).
Using the factor loadings obtained from the first four principal components and scaling them to unity, the following equations were obtained: P hypso = 0.27Hs + 0.23Hi + 0.23Hk + 0.22DHs The transformation of the variables in the analysis was made in such a way that the higher the value of the component, the higher the flashiness or debris flow susceptibility.From the variability explained by each principal component, the morphometric indicator would be P morp = 0.28P shape + 0.20P hypso + 0.22P energy + 0.30P size . (7)

Assessment of appropriateness of the morphometric indicator
The slope-area relationship for the two regions of the study area (Tunjuelo basin and the eastern hills) and the two For the Tunjuelo, eastern hills and La Negra watersheds the break in the slope-area diagram according to the segmented regression is located in a range of 0.11-0.17km 2 for slopes between 0.14 and 0.27.This is in the range of the values found by other authors for transition from debris flows to alluvial processes (Santos and Duarte, 2006;Montgomery and Foufoula-Georgiou, 1993;Seidl and Dietrich, 1993).The points that belong to La Chapa watershed do not allow the identification of a threshold.The drainage area of this watershed is only 7 km 2 , which makes the identification of the break difficult (see Fig. 4).Despite the significant scatter of the values, the slope-area points of La Chapa Creek are located above the points of the other watersheds (see Fig. 4).Regarding the comparison of the points with the threshold of extreme events defined by Horton and Jaboyedoff (2008)  to and above the threshold for areas from 0.02 to 10 km 2 .None of the other watersheds reach the threshold of extreme events, although the points of La Negra watershed are close to the threshold for areas between 2 and 10 km 2 (see Fig. 4).The points of the Tunjuelo River basin, in general, lie lower than the points of the other watersheds and are thus more distant from the threshold of extreme events.It can be observed that the segmented regression of the Tunjuelo River basin is located below the segmented regression of the watersheds located in the eastern hills, with a difference of approximately 0.05 m m −1 in slope (see Fig. 4).
The propagation for initiation points that meet the slope-area thresholds was calculated using the MFS algorithm.However, this appears to overestimate the number of debris-flow-dominated watersheds.JICA ( 2006) identified slope failure areas related to debris flow occurrence in four watersheds located in the centre of the study area using aerial photographs from 1997 to 2004 (see Fig. 5).The method applied by JICA ( 2006) identifies recent slope failures, old slope failures and mass movements related to potential debris flow initiation.In order to assess the initiation points obtained from the two approaches applied in this study, these were grouped into clusters where the distance between points is less than 50 m, in such a way that the clusters represent an area that produces the same propagation trajectory as the individual points.
The photointerpretation carried out by JICA ( 2006) resulted in 108 areas of failure.The slope-area threshold procedure correctly identified 82 % of these slope failure areas, with 107 clusters lying on the areas identified by JICA (2006).In contrast, the extreme event threshold correctly identified 65 % of the slope failure areas, with 103 clusters lying on the slope failure areas.Regarding the amount of initiation clusters identified by each criterion, the slope-area threshold resulted in 389 clusters, while the extreme events criteria identified 299 clusters.The slope-area threshold results in a false positive rate of 72 %, and the extreme event threshold in a false positive rate of 66 %.The visual comparison of the initiation points is shown in Fig. 5.For the case of the slope-area threshold the clusters are scattered covering the mountainous area of the watersheds, and even if they intersect the failure areas the clusters cover significant areas out of them, without showing a pattern associated with the past landslides.In the case of the initiation points from the extreme events threshold, these are not scattered in the upper watersheds but concentrated in areas from which 65 % correspond to past failures.Even if the false positive rates for both methods are high, the overestimated amount and distribution of initiation points in the case of the slope-area threshold procedure lead to unrealistic results when the propagation is applied, with propagation areas occupying most of the area of  the watersheds.Therefore, the propagation was recalculated using only the points above the curve of extreme events.
To assess the performance of the MSF algorithm, the propagation area was compared with the survey of the inundation extent of the debris flow occurred on 19 May 1994 in the upper basin of the Chiguaza Creek.This event affected 830 people and caused the death of 4 people (JICA, 2006).Figure 6 shows the results of the propagation for the Chiguaza Creek.The inset shows the extent of the debris flow using the MSF algorithm as well as the observed extent.The comparison of the modelled and observed runout distance of this event shows a good agreement, although deviations from the modelled propagation areas exist in the final part of the runout (see inset).The deviations occur at bridges, which agrees with the analysis of the event carried out by JICA ( 2006), which concluded that obstructions in crossings had significantly influenced the trajectory of the flow.Simplified models like MSF cannot take the influence of bridges on the propagation of the flow into account.However, independent of the trajectory, the model seems to represent fairly well the downstream extent of the flow, which is the main result needed for the analysis carried out in this study, since the distance between the simulated and observed downstream limit is only 60 m.
Once the results of the MSF algorithm were obtained in the study area, these were used to classify the watersheds according to their capacity to propagate debris flows with two angles of reach.Figure 7 shows the distribution of the morphometric indicator against the classification of the watersheds according to the angle of reach.A clear differentiation can be seen for watersheds classified as 0.19 H /L (able to propagate debris flows to their fans with an angle of reach of 0.19), with the lower quartile located above the interquartile ranges of the other two classifications.However, the differentiation between watersheds classified as 0.11 H /L (able to propagate debris flows to their fans with an angle of reach of 0.11) and clear-water watersheds (C) is less clear.Even if the lower and upper quartiles of the 0.11 H /L watersheds are higher, the median value is smaller than the C watersheds.From this result, a qualitative subdivision into categories was made on the basis of the indicator.Low values from 0 to 0.35 correspond to watersheds unable to propagate debris flows to their fans according to the MSF algorithm, medium values from 0.35 to 0.61 correspond to watersheds where a propagation is possible with a reach angle of 0.11 and high values from 0.61 to 1 correspond to watersheds that can propagate debris flows with an angle of reach of 0.19.
Figure 8 shows the results of the morphometric indicator and its comparison with the results from the debris flow propagation algorithm.Figure 8a shows the values of the morphometric indicator and its classification according to Fig. 7, in which the flood records and the observed type of flow are overlaid.The definition of the observed flow type was possible for 11 watersheds, where the flood records, reports and available studies provide enough information to classify the watersheds into clear flow, hyperconcentrated flow and debris flow according to the method explained in Sect.2.2. Figure 8b shows the resulting propagation areas for different angles of reach.The corresponding classification of the watersheds is shown in Fig. 8c depending on whether or not the lowest point of the watershed is reached by the propagation areas according to the angle-of-reach condition.The comparison of the spatial distribution of the morphometric indicator and the available flood records shows that the area of highest density of flood records is located in the centre of the study area, where the morphometric indicator ranges from 0 to 0.61 (see Fig. 8a).
A contingency table was constructed to assess the skill of the morphometric indicator to identify watersheds with the capacity to propagate debris flows to the fans according to the MSF model considering 0.11 H /L watersheds less dangerous than 0.19 H /L watersheds since the former are more fluid.The results are shown in Fig. 9a.When the three categories of the morphometric indicator are compared with the three flow classifications from the MSF model for all the watersheds in the study area, the proportion correct given as the fraction of the watersheds correctly identified is 0.56.When the contingency table is reduced to 2 × 2 dimensions -this is when only the identification of clear-water and debris flow watersheds is assessed considering low values of the indicator associated with clear-water flows and high and medium values associated with debris flows for angles of reach of 0.19 and 0.11 -the proportion correct reaches 0.75.The contingency table to assess the skill of the morphometric indicator to identify the observed flood types in the study area is shown in Fig. 9b.The 3 × 3 contingency table for the 11 watersheds for which flood type classification was possible results in a value of 0.36 for the proportion correct, while the 2 × 2 contingency table provides a proportion correct of 0.55.
The values of the morphometric indicator obtained for the three watersheds outside of those analysed for the development of the morphometric indicator correspond to 0.3 in the case of the subwatershed of the Chiguaza Creek, 0.43 in the case of La Chapa Creek and 0.14 in the case of La Negra Creek.

Land cover indicator
The three factors used to qualify the state of the watersheds (percentage of vegetation cover, percentage of urban area and percentage of bare soil) are shown in the ternary plot in Fig. 10a, where five areas were identified.As explained in Sect.2.2.2, limits for intensive degradation of the watershed were established taking into account the percentage of vegetation cover and bare soil cover (15 and 10 % respectively) an additional limit was introduced in the ternary plot of 50 % of vegetation cover that delimits area D in Fig. 10a, which represents watersheds with low urban use but high bare soil with low vegetation cover.Watersheds corresponding to zones C, D and E in Fig. 10a were grouped into watersheds in poor condition, watersheds in zone B correspond to fair condition and watersheds in zone A to good condition.The position of the dots in the ternary plot represents the conditions of the watersheds of the study area.Most of the dots are located in zone A. However, highly urbanised watersheds with poor vegetation cover and bare soil can be identified.The spatial distribution of these watersheds can be observed in Fig. 10b where a critical area can be localised in the lower part of the Tunjuelo River basin.

Combination of indicators to obtain a final susceptibility index
From the tests on the possible combination matrices defined by the structure showed in Fig. 3, the highest proportion correct that was obtained was 0.75 considering the three susceptibility classifications and the three types of flow obtained for the 11 watersheds where information was enough to carry out the classification.The debris flows were assigned the most dangerous condition.A proportion correct of 0.91 was obtained when only distinction between clear-water flows and debris flows was considered.The optimum matrix is shown in Fig. 11a and b shows the contingency matrices.
Figure 12 shows the resulting classification of the watersheds applying the matrix shown in Fig. 11a.In this, observed occurrence of floods was superimposed on the susceptibility classification where each dot represents a recorded flood.The spatial distribution of the flood events clearly concentrates in the watersheds located in the lower basin of the Tunjuelo River where there is a cluster of watersheds classified as medium and high susceptibility.watersheds have a low indicator.This is mainly due to the size indicator, which in comparison with the size of the watersheds in the areas assigns low values, with the lowest being the indicator of La Negra Creek, which has an area of 68.4 km 2 (the largest area in the analysis).It is important to take into account that the composite morphometric indicator involves not only the capacity of the watershed to propagate debris flows but also the flashiness; this means that watersheds with the characteristics to propagate debris flows are not necessarily the flashiest.
From the results of the size indicator shown in Fig. 13b, it can be observed that in general watersheds classified as 0.19 H /L exhibit high values of the indicator.However, the size indicator does not discriminate between processes.This could be due to the scale of the analysis, since all the analysed watersheds can be considered small.In the principal component analysis, the size indicator has the highest weight in the total morphometric indicator (see Eq. 7).Several studies have shown that drainage area is correlated with other morphometric parameters; for example, De Scally and Owens (2004) suggest that drainage area acts as a surrogate for the channel gradient, and Gray (1961) and Shreve (1974) showed the correlation between the length of the mainstream and drainage area.Similar to the findings of other authors (Mesa, 1987;Gray, 1961;Shreve, 1974) a high determination coefficient was found between the logarithm of the stream length and the logarithm of drainage area (R 2 = 0.92).The same behaviour is exhibited by the logarithm of length of the watershed (R 2 = 0.90), logarithm of watershed width (R 2 = 0.95) and logarithm of perimeter of the watershed (R 2 = 0.96).The empirical relationship between length of the mainstream (longest stream) and the area is known as Hack's law (Hack, 1957).The exponent of the power law may vary slightly from region to region, but it is generally accepted to be slightly below 0.6 (Rigon et al., 1996).In the study area this exponent corresponds to 0.59.Several authors have tried to explain the relationship between mainstream length and basin area (Mantilla et al., 2000).The conclusion reached by Willemin (2000) indicates that there is some aspect of the evolution of fluvial systems not yet understood that somehow takes into account three geometric components (basin elongation, basin convexity and stream sinuosity), none of which is particularly well correlated with basin area, and produces a robust relationship between mainstream length and basin area.This conclusion is coherent with the findings of this study, where elongation does not show a strong correlation with basin area or a trend to more elongated basins with increasing size of the watershed, as it is not in the same principal component (see Table 2).
The energy indicator, which provides a measure of the potential energy, is composed of the relief ratio, the mean watershed slope, the stream slope, the Melton number and the mean of the MRI.As suggested by Gallant (2003), the MRI can lead to identify similarities and differences between catchments, which in this analysis correspond to the energy of the watersheds.High Melton numbers have been previously used as an effective discriminator of debrisflow-dominated watersheds.However, the threshold for the Melton number varies significantly depending on the region, ranging from 0.5 (Welsh and Davies, 2010) to 0.75 (De Scally and Owens, 2004).Despite the variability of its components, the energy indicator clearly distinguishes 0.19 H /L watersheds.Some superposition of values occurs, but the interquartile range of 0.19 H /L is separated from the interquartile ranges of the other two classifications (see Fig. 7).In terms of energy it is more difficult to distinguish between 0.11 H /L watersheds and C watersheds.However, the mean and the first and third quartiles of the energy indicator for 0.11 H /L are higher than in the case of C watersheds, but with a wider range of superposition.The high values of the energy indicator for the subwatershed of Chiguaza, La Chapa and La Negra creeks is consistent with the processes that take place in the watersheds.
Regarding the hypsometric indicator, since the hypsometric integral decreases as mass is removed from the watershed, it follows that an inverse relationship between hypsometric skewness and the hypsometric integral exists (Harlin, 1984).This condition was found in the study area with a determination coefficient of 0.71.The same behaviour is exhibited by the density skewness (R 2 = 0.82) and hypsometric kurtosis (R 2 = 0.45), where small values are characteristic of large integral values and small skewness.The density kurtosis shows no correlation with the hypsometric integral; this is reflected in the low correlation of this parameter with the corresponding principal component in the analysis (see Table 2).Headward erosion that starts at the lower reaches would represent a higher possibility of debris flow affecting the urbanised fans of the watersheds; therefore this increase in susceptibility would be represented by high hypsometric integrals, low hypsometric skewness and negative density skewness.Furthermore, according to Cohen et al. (2008) higher hypsometric integral values (greater than 0.5) represent catchments dominated by diffusive erosion processes (concave down hypsometric curve) while lower values (less than 0.5) represent fluvially dominated catchments (concave up hypsometric curve).Therefore the hypsometric integral is linked to erosion processes, landform curvature and landscape morphology.
The boxplots of the hypsometric indicator (see Fig. 13b) do not show a differentiation according to the classification of watersheds based on the capacity to propagate debris flows.The superposition of the values of the hypsometric indicator obtained for Chiguaza, La Chapa and Negra creeks shows interesting results, mainly for the case of La Chapa Creek where the value can be classified as high in comparison with the other watersheds.Linking this result with the slope-area plot (Fig. 4) where no fluvially dominated area was identified for La Chapa Creek, it can be inferred that there is a dominance of diffusive processes (characteristic of debris flows) in this watershed that is captured by the morphometric indicators.Therefore, the hypsometric indicator may contribute to explaining the dominance of processes that supply sediment in the watershed.The availability of sediment is one of the determining factors for the occurrence of debris flows.However, its assessment requires extensive fieldwork and detailed sediment source analysis.This assessment is not replaced by the hypsometric indicator, but for the scale of the analysis this indicator is considered to significantly contribute in the susceptibility recognition.
For the case of the shape indicator, drainage density was found to be correlated with the principal component related to the shape of the watersheds, which confirms its relation with the physiographic characteristics of the watersheds (Gregory and Walling, 1968).The boxplots of Fig. 13b show that the capacity to transport debris flows is independent of the shape indicator.However, the values of the indicator for the three watersheds used as external test areas (Chiguaza, La Chapa and Negra) are in the range of high values; particularly Chiguaza and Negra creeks show a very high value.These two watersheds are very similar in terms of shape, hypsometry and energy.
High values of the indicator involve small-area, highenergy watersheds with shapes that contribute to flashiness and hypsometric characteristics that imply erosive processes.

Debris flow propagation
The analysis of the slope-area curves shows that on average the slope in La Chapa watershed is higher for a given drainage area than for the other watersheds considered.If the same drainage area, e.g. 1 km 2 , is considered for the three watersheds with segmented regression fit shown in Fig. 4 namely Tunjuelo River basin, the eastern hills and La Negra Creek -the slope values from the slope-area curves are 0.1, 0.15 and 0.16, respectively, which means that on average for this drainage area the local slope in the Tunjuelo River basin is milder than in the eastern hills, with the latter being slightly milder than the local slope in La Negra Creek.In the case of La Chapa watershed the value of slope for a drainage area of 1 km 2 is 0.4.This result is important given that La Chapa Creek has a confirmed debris flow dominance, followed by La Negra Creek where concentrations in the transition from hyperconcentrated flows and debris flows have been identified.High values of the morphometric indicator are concentrated in the watersheds located in the northeast of the study area.This behaviour is in agreement with the characteristics of the slope-area diagram shown in Fig. 4, where on average the watersheds in the eastern hills have higher local slope for a given area than in the Tunjuelo Basin watersheds.This condition reflects a difference in energy between the two areas that is captured by the morphometric indicator.
The differences in the threshold of extreme events and the location of the slope-area points belonging to each watershed imply that applying the threshold of extreme events reduces the amount of initiation points in comparison with the use of the threshold obtained from the slope-area relationship for dominance of debris flow processes.The comparison of initiation points obtained from the slope-area and from the threshold of extreme events with the failure areas obtained from photointerpretation shows that the slope-area and initiation points seem to overestimate the amount of initiation points in the study area.It is important to highlight that the points correspond to values of local slope averaged in a range of area; therefore, even in the case of the Tunjuelo River basin individual points that meet the extreme event criteria can be identified.However, the amount is less than in the case of the other watersheds.
The results of the MSF algorithm using the threshold of extreme events show that, independently of the classification of flow type of the watershed based on the type of flow at the mouth, other types of flow can occur in other areas of the watersheds, as is the case with the Chiguaza Creek where the extent of the propagation was compared with a field survey providing a good correspondence between the two.The classification of the type of flows at the mouth of this watershed is clear-water flow.However in upper areas where the supply of sediment is high, the morphometric conditions favour debris flows and the land cover is characterised by areas with bare soil.

Land cover indicator, composite susceptibility index and comparison of results
Even if the morphometric indicator provides insight into the expected behaviour and dominant processes of the watersheds reflecting the propagation capacity of the watersheds with a proportion correct of 0.56, it does not fully explain the distribution, characteristics and occurrence of the flood events in the study area.The proportion correct of the contingency matrix comparing the classification obtained from the morphometric indicator and the flow type from flood records yields a value of only 0.36.When the land cover indicator was included in the analysis on the basis that the land cover can exert a positive influence in the case of vegetated surfaces, but also can enhance the susceptibility conditions when urban and bare soil areas are significant, the proportion correct of the contingency matrix comparing the resulting susceptibility indicator and the flow type from flood records increased to 0.75.
It is important to consider that the mountains of Bogotá, mainly in the south of the city and in some localised areas of the east, have been subjected to illegal urbanisation processes.The processes involved in informal settlement entail the construction of houses in the creeks, in some cases not only in the protection buffers but also in the channels.Furthermore, urbanisation requires river crossings that in many cases are not technically designed and constitute dangerous obstructions to the flow as presented in Sect.3.1.2.Another important aspect to consider is the accumulation of waste material in the channels, which during flood events is transported by the flow, and obstructions are common in highly urbanised watersheds in the study area.
The inclusion of the land cover influence in the analysis helps to explain not only the highly deteriorated conditions of some of the watersheds located in the south of the city where floods are frequent, but also the lower occurrence of flood events in some watersheds in the east of the city where the presence of forests and protected areas has contributed to preserving the natural conditions of the watersheds.This suggests the importance of taking land cover into account when assessing the susceptibility to different types of flash floods in peri-urban areas of cities in mountainous areas.

Conclusions
A susceptibility indicator composed of a morphometric indicator and a land cover indicator was used to classify the flash flood susceptibility of 106 watersheds located in the mountainous peri-urban areas of Bogotá (Colombia).Morphological variables recognised in the literature to have a significant influence in flashiness and occurrence of debris flows were used to construct the morphometric indicator.Subsequently, this indicator was compared with the results of simplified debris flow propagation techniques and with the flood type classification carried out in 11 watersheds of the study area, and it was assessed in three additional watersheds to those analysed in the development of the morphometric indicator.These comparisons were made in order to assess the appropriateness of the morphometric indicator.A susceptibility index for each of the catchments was subsequently obtained through the combination of the morphometric indicator and a land cover indicator.An important consideration during the analysis is that watersheds that are prone to debris flows are more dangerous than other flashy watersheds.
The derived susceptibility index is not absolute, but relative, and is useful in applications at regional scales for preliminary assessment and prioritisation of more detailed studies.A limitation of the method is that it does not take sediment availability into account, which is a determining factor for debris flow occurrence.Even if some morphometric indicators could be related to erosion and sediment availability, this factor should be assessed through other techniques.
The morphological variables that were identified to enhance debris flow hazard were analysed through principal component analysis, finding that the 20 variables could be summarised in four component indicators related to size, shape, hypsometry and energy of the watersheds.Size of the watersheds is the component that has the highest weight in the development of the final morphometric indicator.This result is in agreement with previous research that identifies this parameter as relevant in the identification of hazard.
The use of the slope-area curve to identify debris flow source areas showed an overestimation of potential sources when compared with other methods using empirical thresholds.However, it provides valuable information on the processes occurring in a watershed.The slope-area diagram obtained regionally can provide insight in the susceptibility at morphometric level when curves are compared between watersheds in different areas.In the case of the study area, the comparison of the slope-area curves of the Tunjuelo basin and the eastern hills watersheds allowed us to conclude that the latter exhibit on average a higher slope for a given area, which is reflected in the energy indicator that is linked to the capacity to transport debris flows.
The energy indicator was shown to distinguish watersheds with the capacity to transport debris flows to their fans.This indicator involves parameters previously successfully used to identify debris-flow-dominated watersheds.While the prevalence of debris flows in a watershed should be confirmed using detailed information on geology and geotechnics, this parameter can be taken as an initial assessment and used for prioritisation of where to focus such detailed studies.
The use of size, shape and hypsometry indicators in addition to the energy indicator contributes to including valuable information in the analysis to integrally assess the watersheds.Size includes information regarding flashiness as well as shape.Hypsometry was found to be a promising indicator regarding the geomorphic evolution of the watershed and the erosion.
Despite the ability of the morphometric indicator to identify the capability to transport debris flows, it was found not to be sufficient to explain the records of past floods in the study area.The land cover indicator was included, with the objective of involving in the analysis not only the benefit of vegetated areas but also the enhancement of hazard conditions produced by urbanisation and soil deterioration.The indicator produced by the combination of the morpho-metric indicator and the land cover indicator improved the agreement between the results of the classification and the records of past floods in the area.This implies that, even if morphometric parameters show a high disposition for debris flow, land cover can compensate and reduce the susceptibility.However, if favourable morphometric conditions are present but deterioration of the watershed takes place, the danger increases.

Figure 1 .
Figure 1.Location of the study areas.

Figure 2 .
Figure 2. Schematic representation of the methodology.
Correlated with base flow, peak flood discharge and flood potential.Baker (1976),Patton and Baker (1976) Watershed length (L wshd )Has been used to differentiate between watersheds prone to debris flows and debris floods in combination with the Melton ratioWilford et al. (2004) Watershed mean slope (S) Related to flashiness of the watershed.Used to discriminate between debris flow and clear-water-flood-dominated watersheds.Al-Rawas and Valeo (2010), De Matauco and Ibisate (2004) Mainstream slope (StrS) Used to discriminate between processes in watersheds.Welsh (2007) Relief ratio R Ra = (H max − H min )/L wshd Used to describe debris flow travel distance and event magnitude.Chen and Yu (2011) Shape factor SF = A/L wshd Related to flow peak and debris flow occurrence Chen and Yu (2011), Al-Rawas and Valeo (2010), Wan et al.(2008)Mainstream length (L str ) Used to discriminate between processes in watersheds.Chen et al. (2010) ),Sodnik and Miko (2006),Saczuk (1998),Rowbotham et al. (2005),Wilford et al. (2004) Hypsometric integral(HI)   Linked to the stage of geomorphic development of the basin; indicator of the erosional stage; related to several geometric and hydrological properties such as flood plain area and potential surface storage; the hypsometric curve has been used to establish empirical correlations between the hypsometric parameters of a watershed and its observed time to peak.Used to differentiate between processes in the watershed.Pérez-Peña et al. (2009), Harlin (1978), Luo and Harlin (2003),Willgoose and Hancock (1998),Hurtrez et al. (1999) Hypsometric skewness (Hs) Reflects the amount of headward erosion attained by streams; high values are characteristic of headward development of the mainstream and its tributaries, representing the amount of headward erosion in the upper reach of a basin.changes are concentrated and whether accelerated forms of erosion, like mass wasting, are more probable in the basin's upper reaches.When density skewness equals 0, equal amount of change is occurring, or has occurred, in the upper and lower reaches of the watershed.Harlin (1984) Density kurtosis (DHk) Relates to the mid-basin slope.Harlin (1984) Average of the multiresolution index -MRI mean (MRIm)Discriminates between depositional regions and erosional regions.different spectral wavelengths were extracted from the LANDSAT image for training samples, with known land cover obtained from the inspection of a highresolution Google image.The reflectance data of the training samples were used in a recursive partitioning algorithm from which a classification tree is obtained and applied to all pixels of the LANDSAT multiband image to establish separability of the classes based on the spectral signatures.

Figure 3 .
Figure 3. Matrix of classification of susceptibility.

Figure 4 .
Figure 4. Slope-area diagram for the study area and comparative areas.This figure shows the log slope versus log area for each pixel in the watershed areas.To increase readability the value of the slope is averaged in bins of 0.2 log of the drainage area.The black line corresponds to the curve of extreme events given by Eqs.(1) and (2).

Figure 5 .
Figure 5.Comparison of failure areas detected by JICA (2006) and initiation points identified through the slope-area and the extreme event threshold.

Figure 6 .
Figure 6.Affected area in the Chiguaza Creek on 19 May 1994 compared with propagation areas obtained from the MSF model.

Figure 9 .
Figure 9. Contingency table to compare the watershed classification according to debris flow propagation capacity from the MSF model and the morphometric indicator, and the flood type classification from available information and the morphometric indicator.

Figure 10 .
Figure 10.(a) Ternary plot for classification of watersheds according to land cover.The description of the zones of the plot is as follows: (A) low percentage bare soil, low percentage of urban soil and high percentage of vegetated areas; (B) high percentage of bare soil, low percentage of urban soil and high percentage of vegetated areas; (C) low percentage of bare soil, high percentage of urban land and low percentage of vegetated areas; (D) high percentage of bare soil, low percentage or urban soil and low percentage of vegetated land; (D) high percentage of bare soil, high percentage of urban area and low percentage of vegetated cover.(b) Classification of watersheds according to land cover.

Figure 13 Figure 12 .
Figure13shows the boxplots of the composite morphometric indicator and the individual indicators for size, energy, hypsometry and shape.The indicators were grouped according Figure 13.(a) Composite morphometric indicator, (b) indicators based on morphometry.Note that 0.19 H /L and 0.11 H /L correspond to watersheds that can propagate debris flows to their fans considering angles of reach of 0.19 and 0.11, respectively.

Table 2 .
Principal components and corresponding variables.The symbol column shows the abbreviation used in the formulas, and the loading column corresponds to the correlation of each variable with the corresponding principal component.Variables belonging to the PC1 were log-transformed, and variables with the symbol * were transformed as 1 − (value − minimum input value)/(maximum input value − minimum input value).
, the slope-area points of La Chapa watershed are located close