GIS-based predictive models for regional-scale landslide susceptibility assessment and risk mapping along road corridors

ABSTRACT In this paper, a study aimed to assess the landslide susceptibility at a regional scale for the wide provincial territory of Matera (Basilicata region, southern Italy) and the relative risk along the main road corridors distributed in this area is presented. A heuristic-bivariate statistical predictive model was performed to assess and map the landslide susceptibility in the study area by using a polynomial function of eight predisposing factors, weighted according to their influence on the instability process. The resulting susceptibility map was successively used for assessing the landslide risk along the provincial road network. The importance of these roads, representing the main network connecting the urban centres, derives from the absence of an efficient integrated transportation system through the entire regional territory. The landslide risk was evaluated through a matricial approach, which has allowed to define the risk levels (low, medium and high) along road stretches by overlapping the consequences and hazard maps, by combining their corresponding classes in a matrix and by associating to each combination a risk level. The resulting landslide risk map provides support information for decision-making and for identifying the priorities for the design of appropriate mitigation plans.


Introduction
Landslide risk is defined as the expected degree of loss (in terms of loss of life, people injured, damage to properties and disruption of economic activities) due to a landslide in a given area and for a defined period of time (Varnes 1984). For a given category of elements at risk, the specific risk can be quantified as the product of vulnerability, amount of elements at risk and the probability of occurrence of a specific hazard scenario with a given return period in a given area (Van Westen et al. 2006).
The first step in evaluating landslide risk is assessing susceptibility, which is the propensity of an area to generate landslides. The aim of landslide susceptibility mapping is to highlight the spatial distribution of potentially unstable areas. The assessment of landslide susceptibility and risk may be carried out by using several methods according to the scale of the study, the data availability and the aims of the analysis.
In general, the methods for assessing and mapping landslide susceptibility are based on the analysis of the relationships among existing landslides and various factors predisposing instability, weighted according to their influence on the instability process. Based on Soeters and Van Westen (1996), Van Westen (2000), Ayalew et al. (2005) and Fell et al. 2008, these methods, according to the study scale, can be distinguished in (1) inventory-based methods and qualitative heuristic analyses for very small scale (1:750,000-1:250,000), (2) semi-quantitative index-based methods and (3) quantitative models, i.e. bivariate and multivariate statistical methods and training and membership-based models (Artificial Neural Network, Fuzzy Logic, etc.) for small regional scale (1:100,000-1:50,000) and medium scale analyses (1:25,000-1:10,000) and (4) deterministic and probabilistic approaches for large scale (1: . For elaborating the territorial data needed to perform a given susceptibility predictive model, GIS (Geographic Information System) is a tool of paramount importance since it allows, because of its computational power, to manage data with a high degree of spatial variability (Van Westen 2000). Therefore, GIS offers the possibility for a dynamic, integrated and on-going management of the territory and its sudden modifications.
Assessing and mapping landslide susceptibility is an established practice, often aimed at comparing different predictive methods, which offers the possibility for a dynamic, integrated and on-going management of the territory and its sudden modifications (Devkota et al. 2013;Kavzoglu et al. 2014;Pellicani, Frattini & Spilotro 2014;Shahabi et al. 2014;Pham et al. 2016). On the contrary, quantifying, in mathematical terms, the landslide risk can be very complicated, due to several aspects, related to the complexity in assessing the temporal probability of a specific landslide event with given intensity (hazard) and the probability of damaging a given element at risk, i.e. vulnerability (Glade 2003;Uzielli et al. 2008;Pellicani, Van Westen, Spilotro 2014;Abdulwahid & Pradhan 2016). Nevertheless, for a specific type of landslide mechanism and element at risk, a number of procedures for assessing and mapping risk have been proposed in the literature. That is the case of rockfall risk along roads (Pantelidis 2011;Volkwein et al. 2011;Pellicani et al. 2016). In the recent literature, quantitative risk assessment procedures, attempted to estimate the risk in terms of annual probability of loss of life of occupants of vehicles impacted by falling boulders, have been developed (Corominas et al. 2005;Pantelidis 2011;Ferlisi et al. 2012;Mignelli et al. 2012;Stock et al. 2012;Budetta et al. 2016). The risk assessment in probabilistic terms requires the analysis of (1) the probability of occurrence of a rockfall event with a given magnitude, depending on the rockfall frequencymagnitude relationship and triggering event frequencies, (2) the probability that a boulder reaches the element at risk depending on the propagation process along the slope and mobility of element at risk, (3) vulnerability which depends on rockfall intensity and characteristics of exposed assets and (4) economic value of elements at risk, in terms of damaging of road and vehicles, interruption of traffic and loss of life (Corominas et al. 2005;Agliardi et al. 2009;Mavrouli & Corominas 2010;Ferlisi et al. 2012;Wang et al. 2014). Among the different types of landslides (Varnes 1984), rockfalls are generally characterized by small size, but by relatively high magnitude due to the high falling velocity of boulders and accordingly by greater impact and damaging potential on elements at risk, especially those mobiles (vehicles). All the remaining types of landslides, i.e. flows and slides, affecting road corridors are more widespread on the territory than rockfalls, as they affect hillslopes with small slope angles and are characterized by different lithologies (clay, sand, flysch, etc.), unlike the rockfalls affecting mainly steep and rocky slopes. For these reasons, landslides have a greater impact on the fixed elements at risk (road and traffic). Nevertheless, the assessment of landslide risk along road corridors is poorly treated in the literature.
In this paper, a GIS-based procedure for assessing and mapping landslide susceptibility in the provincial territory of Matera (Basilicata region, southern Italy) and the associated risk along the provincial road corridors is presented. Therefore, the first aim of this study is to produce a landslide susceptibility map at a regional scale (1:100,000) for the entire territory of Matera Province (3479 km 2 ), to be used as basic thematic layer for mapping the landslide risk along provincial road corridors located in the study area. The methods used for landslide susceptibility and risk assessment are, respectively, a polynomial heuristic-bivariate statistical model (Pellicani, Frattini & Spilotro 2014) and a qualitative matrix approach. First, a set of thematic maps related to the predisposing factors was prepared and correlated with a portion of the landslide inventory map in order to obtain, by using a bivariate procedure (Van Westen 1993), the weights representing their importance on the instability process. Finally, the factors were weighted and combined among them in a polynomial function. This procedure was validated by comparing the resulting landslide susceptibility map with the remaining portion of the landslide inventory. Subsequently, the assessment of landslide risk was carried out qualitatively by overlaying the consequences and hazard maps and by combining in a matrix, the relative classes. The consequences were derived by combining the vulnerability and exposure maps; while the hazard was evaluated in function of susceptibility and landslide intensity, depending on size and velocity of the landslide process.

Matera territory and road corridors
The study area is the provincial territory of Matera, in Basilicata region (southern Italy), covering an area of 3479 km 2 (Figure 1(a)). The flat and hilly areas are predominant, respectively about 63% and 30%, while the mountain ones occupy 7% of the whole region. From a geological point of view (Figure 1(b)), the study area is mainly characterized by a regressive, generally coarsening upward, sedimentary body consisting of sand and conglomerate levels unconformably overlying the Pliocene -Middle Pleistocene Subapennine Clay Formation, belonging to Bradanic trough  (Caputo et al. 2010). The eastern sector consists of the same formation overlying the calcareous structure of the Murgia Plateau, belonging to Apulian foreland. The active tectonics of the southwestern sector, in correspondence of Apennine Chain, has resulted in the presence of chaotic structural conditions (different types of joints and overlying structure) and alternate lithologies, clayey flysch and calcareous-marly flysch, even at high altitudes (over 1000 m).
The road corridors of the Matera Province extend in the study area for 1324 km, connecting 31 municipalities among them. This road network exerts an important role for the entire transportation system of the provincial territory, due to the lack of other types of effective infrastructures (railway, motorway, aerial, etc.). In the last years, about 44% (584 km) of roads were affected by damages. Concerning the geological variability in the study area, the most typical and recurrent instability phenomena, interacting with road corridors, affect mainly the following outcropping lithologies ( Figure 2): (1) clays, exposed along hillslopes or covered by sandy (variously cemented) levels, (2) clayey flysch and calcareous-marly flysch, (3) sand and conglomerate and (4) debris deposits.

Landslide inventories
Two landslide inventories were used in this study (Figure 3(a,b)). The first was produced by River Basin Authority of Basilicata within PAI project (Plan for the Hydrogeological Asset) (Figure 3(a)). The procedure for delimiting the unstable areas and editing the landslide inventory map is characterized by (1) the acquisition of the available historical documents and data on existing landslides, (2) checking and recognition of landslide processes by means of interpretation of 1:33,000 scale aerial-photos (1991 and 1997) and field surveys, (3) identification and classification of the type, intensity and state of activity of landslides, (4) definition of morphological features and locations of landslides and, finally (e), representation at an appropriate scale. The PAI landslide inventory map is constantly updated; in particular, the version used in this study is updated to 2015.
The second landslide inventory includes the landslides affecting the provincial road corridors (Figure 3(b)) and was exclusively produced for the aims of this study, through field surveys along the 1324 km of road. The landslide polygons were digitized in GIS on georeferenced orthophotos (relative to 2013) at a 1:5000 scale. About 1280 landslides were detected, mapped and classified according to movement type. Seven landslide types were recognized: falls/topples (1%), translational slides (1%), rotational landslides (26.7%), rapid flows (0.3%), complex movements (10%), slow flows (21%) and widespread shallow landslides (25%). The first six types are in accordance with the classification defined by Cruden and Varnes 1996, while the last category was introduced to indicate areas characterized by the coalescence of different landslides not individually classifiable or by shallow landslides with poorly defined boundaries (Pellicani & Spilotro 2015). For 15% of the landslides, the movement type was not determined. For each landslide, the following set of information surveyed in the field was collected in an alphanumeric database: morphometry, location on hillslope and relation to road, geological and land use features, classification of movement type, state of activity and groundwater presence, damages produced on road and existing consolidation works.

Predisposing factor thematic layers
The starting point for performing a GIS-based susceptibility predictive model is the analysis of the spatial relationships among existing landslides and various factors considered to have an influence on the occurrence of landslides. The selection of causal factors differs depending on the scale of analysis, the characteristics of the study area, the landslide type and the failure mechanisms (Guzzetti et al. 1999;Van Westen et al. 2008). In this study, for assessing and mapping the landslide susceptibility in Matera provincial area, the following predisposing factors were selected and mapped: lithology, land use, hydrological factors, i.e. drainage density, curve number and potential erodibility, and morphometric factors, i.e. elevation, slope angle and aspect. For each of these thematic maps, the incidence of landslides in their classes was evaluated.

Lithology
Geo-lithological map at 1:100,000 scale is obtained by digitizing the geological sheets of Geological Map of Italy. This map consists of 11 outcropping lithological units, which are shown in Figure 4(a). By analysing the spatial distribution of landslides in the different lithological classes, it is noted that the lithologies mostly affected by instability phenomena are represented by heterogeneous deposits (22%), debris (17%) and clay soils (14%). While in areas with sand or sandy conglomerate outcrops, the landslides are before originated by cracking of sand or conglomerate (variously cemented) and later interest (progressive failure) the underlying clayey layer, with very rapid final evolution. The products of collapse, in the presence of rainwater or water deriving from the altered underground drainage, evolve, generally with seasonal periodicity, as debris or mud flows. Landslides in the other lithological classes have a lower incidence, except for the flyschoid soils (12%). Indeed, the lithological diversity and the chaotic structure of these soils define specific situations for each site.

Land use
Land use map compiled at 1:100,000 scale was obtained by Corine Land Cover (2012). In this study area, the land use map is composed by 31 units. In order to reduce the data redundancy, the units belonging to the same land use category were grouped together. The resulting land use map consists of the following nine classes: urban area, crops, pasture, arable areas, forests, orchards and olive groves, mixed vegetation areas, bare and water bodies (Figure 4(b)). Among these land use categories, the most affected by instability are mixed vegetation areas (24%), pasture (18%) and bare (17%).

Morphometric factors
3.1.3.1. Elevation. Digital Terrain Model (DTM) with 20 £ 20 m grid-cells (Figure 4(c)) was generated through an interpolation algorithm from contour lines with a 5 m interval and elevation points, which were extracted from the Apulia Regional Technical Map (CTR 2008), at 1:5000 scale. From the DTM standard morphometric parameters, such as slope and aspect, and hydrological parameters were derived. Considering the unstable areas, the mean elevation ranges between 200 and 350 m a.s.l.

Slope angle.
A slope angle raster with 11 classes was derived from the DTM 20 £ 20 m (Figure 4(d)). The maximum slope angle in the area is 65 , with a mean value of about 9 . By evaluating the statistical distribution of slope angles in areas affected by landslides, the resulting histogram has revealed that hillslopes with mean slope angles lower than 4 are not affected by landslides, while most of the unstable areas (56%) have a slope angle between 8 and 16 .
3.1.3.3. Slope aspect. The aspect has a great influence on the slope stability. Indeed, it induces the sunlight and, as a consequence, the temperature and the moisture of the soil. In the study area the slopes are mostly positioned in the direction NE-SW (Figure 4(e)). By evaluating the statistical distribution of aspect in areas affected by landslides, the resulting histogram has revealed that slopes with S-SW direction, where sunlight induces more heat and less humidity, are largely affected by landslides. The other relevant slope directions are NE, E, SE and W, while N and NW are rarely affected by landslides.

Hydrological factors
3.1.4.1. Curve number. The runoff curve number is an empirical parameter, developed by the USDA Natural Resources Conservation Service (Soil Conservation Service 1972), used in hydrology for determining the approximate amount of direct runoff from a rainfall event in a particular area (Manzo et al. 2013). In general, the soil predisposition to generate runoff depends on the land cover typology, the soil saturation and permeability. According to USDA-SCS method, CN values are defined for combinations of land cover classes and hydrologic soil groups (A, B, C and D). These groups are classified based on the runoff and infiltration capacity of the soil. Soils in hydrologic group A have low runoff potential. Soils that have a moderate rate of infiltration when thoroughly wet are in hydrologic group B. Hydrologic group C soils that have a slow rate of infiltration rate when thoroughly wet. Soils in hydrologic group D have a high runoff potential.
To produce the CN map, 12 land cover classes were selected by merging the 31 categories defined in the Corine Land Cover. The 11 lithological classes were reclassified in terms of hydrologic soil groups, in relation to their permeability properties, as follows: Asand, gravel and debris, Bheterogeneous deposits and carbonatic rocks, Csandstone, clayey flysch and calcareous-marly flysch, Dclays. Finally, CN values, derived from the USDA-SCS classification, were associated to each combination of the land cover classes and the four hydrologic soil groups (Table 1). In general, the lower CN values are relative to permeable surfaces, which absorb the rainfall, while the higher CN values correspond to impermeable soils, where rainfalls turn into runoff. The CN map is shown in Figure 4(f).
3.1.4.2. Drainage density. Geomorphology-based hydrological models, obtained extracting drainage networks from DTM, allow to reconstruct the hydrological response of a basin (runoff formation) subsequent to a precipitation, based on the link between hydrological response and geomorphological features of the basin. From this model, the drainage density of basins was derived by dividing the length of the stream network comprised in each watershed for the basin surface (Figure 4(g)).
3.1.4.3. Potential erodibility. Rill erosion evolving into deeper forms may be generated in the areas with sand and sandy conglomerate deposits, with various degrees of damage to the roads. Areas with potential erodibility were identified where the stream network overlaps on sandy soils. In particular the potential erodibility map, shown in Figure 4(h), was produced by selecting the portions of the stream network falling on sandy outcrops and by extrapolating for them the slope angle values.

Susceptibility modelling and performance analysis
To assess the spatial distribution of landslide susceptibility in the provincial territory of Matera, a polynomial heuristic-bivariate statistical model was performed. In Pellicani, Frattini & Spilotro 2014, this method produced reliable outcomes comparable with those derived from two multivariate statistical models (discriminant analysis and logistic regression). It is easier to manage and enables the introduction of expert opinion into the analysis with respect with other methods more computationally demanding in processing large data-sets, especially considering the scale of the study. The proposed method is derived from a heuristic 'ranking and rating' index-based approach (Stevenson 1977;Amadesi & Vianello 1978;Romana 1985;Gupta & Joshi 1990;Anbalagan 1992;Gupta & Anbalagan 1997;Budetta et al. 2008), modified through a bivariate statistical procedure (Van Westen 1993) for assigning the weights to the factors influencing slope instability. The susceptibility was obtained as polynomial function of predisposing factors, according to the following expression: where X i are the causal factors and n are the weights, expressing the importance of each factor on the instability. Every factor was classified into a number of relevant classes, associated with a score, ranging from 0 to 1, expressing how much each class contributes to slope instability. These scores were selected heuristically on the basis of the expert knowledge and the incidence of landslides on each class, previously evaluated (see Section 3.1). The entire susceptibility modelling was implemented using ArcGis software on a 20 £ 20 m grid-cell. The weights were evaluated through a statistical bivariate procedure, correlating each thematic layer with a portion of existing landslides (training set) contained in the acquired landslide inventory map. The weight to be assigned to each factor was calculated in the following way: where D A is the 'areal density', defined as the ratio of the whole area involved in landslide to the total study area, and D class is where D i is the 'landslide density per class', defined as the ratio of the area involved in landslide for each class i of the causal factors to the total area in landsliding, s i is the score of the same class and n class is the number of classes. By overlapping each thematic raster, classified according to the scores and weighted, and implementing the susceptibility algorithm in each grid-cell, the local value of landslide susceptibility was obtained. In order to obtain a subdivision of the study area into zones with different levels of landslide susceptibility (low, medium and high), these values were ranked according to threshold values. The thresholds between low and medium classes and between medium and high classes were obtained through a statistical procedure. In particular, the map containing the local values of susceptibility was crossed with landslide inventory map and two georeferenced matrices, containing the susceptibility values of, respectively, stable and unstable points were derived. By assuming that the two data-sets are normally distributed, the correlation was established for stable and unstable points and the mean and standard deviation of these two distributions were used to define the thresholds between the susceptibility classes. The threshold between medium and high susceptibility classes was calculated through this formula: where m u and sd u are the mean and standard deviation of stable point distribution and m s and sd s are the mean and standard deviation of unstable point distribution.
The threshold between medium and low susceptibility classes was evaluated as follows: According to the values of about S H = 0.6 and S M = 0.4, the local values of landslide susceptibility were reclassified and the final landslide susceptibility map was obtained.
The above-described prediction procedure was carried out considering a portion of existing landslides (training set), therefore the model results were validated using the remaining part of landslides (test set) (Chung & Fabbri 2003Carrara et al. 2008).The training data-set was defined by selecting randomly the 65% of unstable area. The remaining 35% of cells were then used for model performance evaluation. Performance refers to comparing the model predictions with a real-world data-set, for assessing its accuracy or reliability (Forbes 1995;Begueria 2006). The performance analysis permits to establish the degree of confidence of the model, which is of great importance for transferring the results to the final users. For the present study, the efficiency and reliability of the model was evaluated by using the success-rate curve. It measures a goodness of fit of the model (Chung & Fabbri 2003). The success-rate curve shows the percentage of landslides falling within the area classified as unstable for increasing values of the selected susceptibility index (Carrara et al. 2008). In particular, the y-axis is normally considered as the number of landslides, or the cumulative percentage of landslide area correctly classified, while the x-axis is the cumulative percentage of area classified as unstable by the model (Frattini et al. 2010). Indeed, for a given number of instability phenomena correctly predicted by the model, the smaller the area classified as unstable the better is the model performance.

Risk assessment
In the absence of detailed input data, the assessment of the landslide risk, in quantitative terms, can be affected by large uncertainty (Wang et al. 2014). For this reason, the spatial distribution of landslide risk along the road corridors of Matera Province was assessed and mapped using a qualitative matrix approach ( Figure 5). It represents a variant of the most known hazard-consequence matrix approach (Chowdhury & Flentjie 2003; Australian Geomechanics Society (AGS) 2007), in which risk is obtained by combining a set of hazard categories with a set of consequence categories. In this study, risk input maps were generated separately, subdivided into classes and afterwards combined pairwise in two-dimensional tables or matrices. The first matrix assesses the landslide hazard, i.e. the probability of a landslide of certain magnitude to occur, and is a function of the return time, which depends on landslide type and intensity. Due to the lack of temporal data associated to landslide mechanisms, the hazard was evaluated qualitatively as a function of the landslide intensity. For a given type of landslide mechanism, the intensity changes according to the areal extent and velocity of landslide. The hazard matrix was obtained qualitatively through the following steps: (1) Assessing the velocity range of landslide phenomena, by associating a velocity category to each type of landslide mechanism collected in the acquired inventory map (Table 2). (2) Evaluating the areal extent of landslide phenomena and subdivision into four classes (Table 3).
(3) Determining the landslide intensity by combining in a matrix the velocities and the areas with relative classes (Table 4). The combination between the four classes of velocity and the four classes of area was classified into four intensity classes: null, low, medium and high (Table 5). A landslide intensity map was generated by reclassifying the landslides contained into the acquired inventory in terms of intensity classes.  (4) Combining the four classes of intensity with the four classes of susceptibility in a matrix and classifying these combinations in terms of five hazard classes: very low, low, medium and high. A column was added to the hazard matrix to consider the area not affected by existing landslides, in which future landslides could occur. For these areas, a very low hazard level was associated to no susceptibility, low hazard level to low susceptibility class, medium hazard level to medium susceptibility class and high hazard level to high susceptibility class. The landslide hazard map was obtained by overlying the susceptibility map and the landslide intensity map.
Consequences are generally defined as the outcome or potential outcome to an element at risk arising from the occurrence of a landslide of certain magnitude (Glade & Crozier 2005). Therefore, consequences are a function of the amount of elements at risk and the vulnerability of the affected elements (Puissant et al. 2004;Andersson-Sk€ old et al. 2014).
The vulnerability of provincial road corridors crossing the Matera territory was assessed in function of the presence of previous criticalities, collected by Provincial Road Department, which needed repairs. In particular, about 44% (584 km) of roads were affected, in the last years, by instability and different types of road repairs were carried out. As the vulnerability is intended as the degree of damage to a given element at risk caused by an instability phenomenon, the road repairs, representative of the type and severity of damages, were ranked in five classes (Table 6), to which a score ranging from 0 (no works) to 1 (severe damages and structural consolidation works) was assigned based on expert knowledge. As the road stretches affected by damages were repaired by performing more than one type of work, among those summarized in Table 6, the vulnerability values were derived by summing, for each road section, the several scores associated at the work categories and by normalizing the resulting values from 0 to 1. Five vulnerability classes were associated to the following score ranges: null (V 0 ) to zero value, low (V 1 ) to 0-0.25, medium (V 2 ) to 0.25-0.45, moderate (V 3 ) to 0.45-0.75 and high (V 4 ) to 0.75-1.   The exposure of the road corridors was evaluated in terms of amount of potential traffic. This choice was derived from the assumption that the economic value and the reconstruction costs of the entire road network are constant (exposure in monetary terms) and performing a probabilistic analysis of vehicles and passengers distribution along the road network is not possible, due to the lack of detailed input data in relation to the analysis scale. The traffic along each road section, connecting two or more towns, was qualitatively assessed as a function of population of each of them. The exposure values were normalized and subdivided into five classes: very low (E 0 = 0-0.1), low (E 1 = 0.1-0.25), medium (E 2 = 0.25-0.4), moderate (E 3 = 0.4-0.5) and high (E 4 = 0.5-1). The consequences map was obtained by overlying the vulnerability map and exposure map and by combining in a matrix ( Figure 5) the five vulnerability classes with the five exposure classes and ranking the combinations into five classes: insignificant, minor, medium, major and catastrophic. Finally, the risk matrix was produced by combining among them the hazard and consequence classes and by associating to each combination the risk classes, i.e. low, medium and high ( Figure 5). The risk map was obtained by overlaying the hazard and consequence maps and associating the risk classes according to the risk matrix.

Results and discussion
Based on the incidence of landslides in thematic layer classes, the scores to be associated to each class were heuristically selected (Table 7). For the lithological classes, the highest scores were assigned to clays (0.9), heterogeneous deposits (0.7) and debris (0.6), respectively; medium values to clayey flyschs (0.5) and calcareous-marly flyschs (0.4); lowest scores to remaining lithological classes. The land use classes more predisposing to instability were considered mixed vegetation areas (0.9), pasture (0.8) and bare (0.7). Regarding the morphometric factors, the frequency distribution of landslides has determined the highest scores for elevations ranging from 200 to 350 m a.s.l., for slope angle values comprised between 8 and 12 (0.9), which could be representative of areas already affected by landslides, and between 12 and 16 (1.0), indicative of potential zones for landsliding initiation, and, finally, for hillslopes orientated to South-West (0.9) and South (0.8). The highest scores were also assigned to drainage density classes ranging between 0.2 and 0.3 (0.9), 0.3 and 0.4 (1.0), 0.4 and 1 (0.8). The score values do not increase proportionally to the class values, as some graphical mistakes in the representation of the basins could be present. On the contrary, the higher the CN value, the higher the assigned score. Finally, for the potential erodibility classes, expressed in degree, the same scores associated to slope angle classes were used.
Through the statistical bivariate procedure, the landslide density per class, the areal density and, finally, the weights were calculated (Table 7). In particular, the weight values have highlighted the predisposing factors most strongly related to the occurrence of landslides, which are lithology, curve number and slope angle. By combining the thematic maps and relative weights according to the polynomial algorithm, the landslide susceptibility map at 1:100,000 scale of the entire provincial territory of Matera was obtained (Figure 6(a)). This map, a raster with a spatial resolution (cell size) of 20 m, was divided into four classes, the first composed by the values near zero (no susceptibility), Table 6. Typologies of road repair works and corresponding scores.

1.00
Water regulation (lateral and transverse drainages to roadway) 0.75 Works for reparation of embankments protecting artefacts, of hydraulic works and bumps, for restoring deformed road surface with mixed stabilized and/or asphalt, etc.

0.50
Interventions for road cleaning (cleaning of clogged drains, removal of mud deriving from hillslopes, etc.)

0.25
No works 0      (4) and (5). Based on the susceptibility zonation, 35% of the provincial territory of Matera is highly susceptible to landslides (S H = 0.6 Ä 1), 43% is affected by medium susceptibility (S M = 0.4 Ä 0.6), 15% is characterized by low susceptibility (S L = 0.1 Ä 0.4) and the 7% is identified as not susceptible (S 0 = 0 Ä 0.1). In order to evaluate the predisposing factor classes affected by high level of susceptibility, the thematic maps were overlaid to the susceptibility map. The highest susceptibility level was observed in correspondence of hillslope with average slope angle of 16 , at elevation comprised between 200 and 400 m a.s.l., exposed at South-East, South and South-West, covered by arable and mixed vegetation on clayey outcropping soils. The reliability of model results was evaluated through the success-rate curve, obtained by plotting the cumulative percentage of area classified unstable (or susceptible) by the model (on x-axis) with the cumulative percentage of landsliding area correctly predicted as unstable (on y-axis). In this validation procedure, the test set containing the 35% of unstable points of landslide inventory, randomly selected, was used as potential future unstable areas. To obtain the curve, both the test set landslides (observed phenomena) and susceptibility map (predicted phenomena) were reclassified into two dichotomous values: 1 for observed and predicted instability and 0 for observed and predicted stability. Forty different susceptibility maps were produced by varying the cut-off value of 0.025 (i.e. the threshold value that subdivides the unstable areas, susceptibility values higher than the threshold, from stable areas, susceptibility values lower than the threshold) and the percentage of unstable area was obtained. These maps were then combined with the binary raster of landslides and the percentage of landslide area within the unstable zone was calculated. The resulting curve, shown in Figure 6(b), highlights a good predictive capability of the model. Indeed, considering the most susceptible 15% of the area of the susceptibility map, the 60% of future landslides are located in this area. The landslide risk analysis was carried out exclusively along the road corridors. For this reason, the susceptibility map was cut for considering the hillslopes above and below the road in a buffer of about 800 m width (Figure 7(a)) and the total landslide inventory, composed of landslides mapped in both the inventories (Figure 3), was used to create the velocity and area maps. Before assessing and mapping the landslide risk, the spatial data contained in the input maps were made uniform, in order to obtain raster maps with the same spatial resolution, i.e. 20 £ 20 m. The landslide inventory map was transformed into a raster map, classified in terms of landslide mechanism type, and then reclassified in terms of landslide velocity and area to obtain the intensity map (Figure 7(b)). While, the vulnerability and exposure raster maps (Figure 7(c,d)) were obtained by rasterizing the road vectors and reclassified them, respectively, in terms of road repair work category and amount of traffic.
By combining the landslide velocity and area classes into a matrix, the intensity classes were obtained. The majority of existing landslides have medium (56%) and large (43.8%) sizes and are characterized by slow (47.2%) and moderate (46%) velocity. According to intensity matrix (Table  4), the 28.2% of landslides has low intensity, the 34.3% medium intensity and the 32.3% high intensity. By overlying the landslide intensity map on the susceptibility map and by combining the corresponding classes into the hazard matrix, the hazard zonation was obtained (Figure 8(a)). Based on the hazard matrix, the highest hazard level was assigned both to areas affected by instability, consisting of medium or high intensity landslides, and predicted as medium or high susceptible areas and to those not yet affected by landslide but characterized by a high susceptibility level. Medium hazard levels were associated to areas with a medium level of susceptibility affected by medium intensity landslides or free by landslides, areas with a low level of susceptibility but affected by high intensity landslides and, finally, areas with a high level of susceptibility but affected by low intensity landslides. By overlying the hazard map with roads in vector format, a map with the hazard classes distributed directly on the road linear vector was produced. Based on the hazard zonation, 57% of the road corridors are affected by high hazard level, 3% is characterized by medium hazard and 40% by low or very low hazard.
The vulnerability zonation along road corridors has revealed that about the 56% of roads is free of damage, 8% by low vulnerability, 16% by both medium and moderate vulnerability and 4% by high vulnerability while the exposure zonation has highlighted that 15% of roads is characterized by moderate (E 3 ) amount of traffic, 17% by both very low (E 0 ) and medium (E 2 ) exposure, 20% by high exposure (E 4 ) and the 31% by low exposure (E 1 ). By comparing the vulnerability and exposure maps (Figure 7(c,d)) it can be noted that road stretches affected by highest levels of vulnerability (i.e. V 4 and V 3 ) are conversely characterized by low amount of traffic. Indeed, 78% of roads free of damage (V 0 ) has moderate and high exposure levels. The distribution of consequence levels along roads (Figure 8(b)) was achieved from the overlapping of vulnerability and exposure maps among them and from the combination of the corresponding classes into the consequence matrix. In particular, the assignment of consequence levels to each combination of vulnerability and exposure classes was carried out by assuming that the influence of exposure was more prevalent than vulnerability, since the later does not derive from a probabilistic analysis of potential damage degree but from the evaluation of distribution of past damage on roads. Based on consequences matrix, 29% of roads is affected by insignificant consequences, 19% by minor, 36% by medium, 9% by major and 7% by catastrophic consequences.
Following the scheme of Figure 5, the hazard and consequence maps were overlaid, the corresponding classes were combined in the risk matrix and the risk classes (low, medium and high) were associated to each combination, up to obtain the final landslide risk map along 1324 km of Matera provincial roads (Figure 8(c)). Road sections affected by very low and low hazard and by insignificant and minor consequences, in their different combinations (H VL -C 0 , H VL -C 1 , H L -C 0 and H L -C 1 ), or by medium hazard and insignificant consequences (H M -C 0 ) or medium consequences and very low hazard (H VL -C 2 ) are characterized by low risk. Conversely, medium and high hazard levels with major and catastrophic consequences (H M -C 3 , H M -C 4 , H H -C 3 and H H -C 4 ), medium consequences with both medium and high hazard levels (H M -C 2 and H H -C 2 ) or catastrophic consequences with low hazard (H L -C 4 ) involve a high risk level. Finally, a medium level of risk corresponds to the following combinations: very low hazard with major and catastrophic consequences (H VL -C 3 , H VL -C 4 ), low hazard with medium and high consequences (H L -C 2 , H L -C 3 ), medium hazard with minor consequences (H M -C 1 ) and high hazard with insignificant and minor consequences (H H -C 0 , H H -C 1 ). Figure 8(d) shows the percentages of the overall road network and of the stretches crossed by landslide bodies included in each of the three risk classes. The landslides mapped in the surveyed inventory (Figure 3(b)) were used. The road sections subjected to low landslide risk are 22% (290 km) on the total length, while the remaining 48% (642 km) and 30% (392 km) of the examined road corridors are affected, respectively, by medium and high risk levels. The comparison between the risk map and the landslide inventory recognized along roads has also revealed that about the 10% (136 km) of roads are crossed by instability phenomena. In particular, on the total of instability phenomena, the 49.5% of landslides affected sections where the risk was evaluated high and 41% sections where the risk was assessed medium. The remaining 9.5% of landslides crossed road stretches where the risk was mapped as low.
The proposed procedure has allowed to obtain the zonation of landslide risk, even if with limited hazard and vulnerability data, through a matrix approach. Although the proposed methodology for the risk assessment did not consider the prediction of landslide run-out and temporal probability, the results obtained allow to rank the road stretches in terms of risk levels and, consequently, to identify the priorities for designing detailed field surveys and appropriate landslide risk mitigation plans.

Conclusions
A procedure for assessing and mapping landslide susceptibility in the Matera provincial area and the associated risk along the provincial road network was presented. The main aim of this study is to provide a landslide risk map as support to local governments for identifying the priorities for the design of appropriate mitigation plans.This study was carried out for the provincial roads as they represent the main network connecting the urban centres of the Matera Province and is affected for the 10% of its total length by instability phenomena.
The landslide susceptibility was assessed by implementing a heuristic-bivariate statistical method based on a polynomial function of eight predisposing factors (i.e. lithology, land use, elevation, slope angle, aspect, curve number, drainage density and potential erodibility), weighted according to their influence on the instability process. The result validation procedure, by using the success-rate curve, has shown that the method, although is based on a semi-quantitative combined analysis of predisposing factors, has provided a reliable assessment of landslide susceptibility.
The landslide risk was evaluated, in the absence of detailed input data due to the scale of analysis, by a qualitative matrix-based approach. The proposed procedure has allowed to assess and map landslide risk, based on only two input data, through the following three sequential steps: (1) producing landslide velocity and areal extension maps from landslide inventory map, according to the movement type, and combining them for obtaining the landslide intensity map; (2) combining, pairwise, intensity and susceptibility maps and vulnerability and exposure maps, these last produced by roads input map, for obtaining the hazard and consequences maps, respectively; (3)overlying the hazard and consequences maps and combining the corresponding classes in a matrix in order to obtain the final risk map of the roads, subdivided in low, medium and high risk levels.
The comparison between the risk map and the distribution of landslides crossing the road stretches has shown that the majority of instability phenomena (about 91%) are located on road sections classified at medium and high risk. This is a significant result, as it shows that it is possible to use, for a regional-scale landslide risk assessment, a simplified-qualitative procedure, based on a few data, easy to find and manage, and reliable in relation to the scale of analysis.