Quantitative analysis of the risk to road networks exposed to slow-moving landslides: a case study in the Campania region (southern Italy)

This paper shows the results of a study aimed at quantitatively estimating—in terms of direct (repair) costs, at large scale (1:5000)—the slow-moving landslide risk to a road network assumed as undamaged as well as the consequences to the same network in damaged conditions. The newly conceived methodological approaches address some challenging tasks concerning (i) the hazard analysis, which is expressed in terms of probability of occurrence of slow-moving landslides with a given intensity level that, in turn, is established based on empirical fragility curves, and (ii) the consequence analysis, which brings to the generation of time-dependent vulnerability curves. Their applicability is successfully tested in a case study in the Campania region (southern Italy) for which both very high-resolution DInSAR data and information gathered from in situ surveys on the severity of damage sustained by the selected road sections are available. Benefits associated with the use of the obtained results in informed decision-making processes are finally discussed.


Introduction
Road infrastructure plays a key role in the economic development of a society. For this reason, ensuring its functionality and safety conditions over time is one of the most important and, at the same time, demanding tasks that central and local authorities are asked to undertake. Indeed, owing to their typical linear extent, roads often develop within different geological contexts, each of them prone to given landslide types that originate risks to traveling persons and to roads themselves, the latter being associated with socio-economic impacts (including indirect costs) such as prolongation of travel time and missed trips (Argyroudis et al. 2019;Hackl et al. 2018). Addressing this issue involves carrying out activities aimed at predicting and (eventually) preventing the above risks, taking into account operational and economic constraints (Fell et al. 2005;Winter 2019). As far as the forecasting activities are concerned, they may consist of qualitative or quantitative risk analyses.
The qualitative risk analysis "uses word form, descriptive, or numeric rating scales to describe the magnitude of potential consequences and the likelihood that those consequences will occur" (Fell et al. 2008). Implemented methods make use of scoring systems and ranking matrices to be applied to a large number of slopes. The obtained outcomes (i) facilitate decision-makers to compare (i.e., evaluate and rank) the estimated "relative" risks and (ii) foster the prioritization of slopes requiring follow-up actions (i.e., study, repair, or maintenance), the cost of the latter being preliminarily assessed (Fell et al. 2005;Pantelidis 2011;Wong 2005). Applications of scoring systems tailored for roads mainly deal with cut slopes prone to first-time landslides, such as rock falls (Budetta 2004;Budetta and Nappi 2013;Bunce et al. 1997;Ferlisi et al. 2012;Li et al. 2009;Vishal et al. 2017); anyway, scoring systems are also conceived and applied to roads threatened by soil cuttings/embankments (Liang et al. 2006;Lowell and Morin 2000;ODOT 2001;Wong 1998) and natural slopes (Escario et al. 1997). Recently, Pellicani et al. (2017) tested the applicability of ranking matrices to the main road network of a large area of southern Italy affected by different landslide types, including existing ones (e.g., active landslides slowly moving on buried sliding surfaces).
The quantitative risk analysis (QRA) is "based on numerical values of the probability, vulnerability, and consequences, and resulting in a numerical value of the risk" (Fell et al. 2008); accordingly, QRA allows one to estimate "absolute" risks in terms of probability of a given level of loss, accounting for the uncertainties associated with input parameters (Fell et al. 2005;Ho 2004;Macciotta et al. 2016). Formal methods aimed at analyzing the factors (hazard and consequences) concurring to the risk definition are well-described in the scientific literature Fell et al. 2005;Ho et al. 2000;Wong 2005;Wong et al. 1997;among others); on the other hand, QRA requires input data more accurate (both in quality and quantity) than those adopted in qualitative risk analyses, along with high-resolution digital elevation/terrain models and detailed information on exposed elements (van Westen et al. 2008). The QRA outcomes can be presented in terms of annual risk, e.g., the probability that a particular person-for instance, the most exposed one-may lose his/her life or the expected repair (direct) costs (e.g., €/annum). The scientific literature provides several examples referring to first-time landslide (i.e., rock fall or debris flow) risk to people on moving/stopped vehicles along roads (Budetta 2002;Budetta et al. 2016;Bunce et al. 1997;Hungr et al. 1999;Ferlisi et al. 2012;Lentini et al. 2019;Mavrouli et al. 2019;Roberds 2005;Unterrader et al. 2018;Winter 2018;Winter and Wong 2020;Wong and Winter 2018) whereas very limited in number are the contributions specifically oriented to the quantitative estimation of the risk to roads affected by existing landslides.
For a study area in Belgium and based on focus and semistructured interviews with involved stakeholders, Vranken et al. (2013) quantitatively estimated the costs of measures (in €/annum) to repair or prevent the landslide-induced damage. Lu et al. (2014), updating the work of Catani et al. (2005), presented a GIS-based landslide risk zoning map that makes use of remote sensing data to provide-on a pixel basis-the expected losses (in €) to the built environment (including roads) in slowmoving landslide-affected areas of the Arno river Basin (Italy) considering five temporal predictions (2, 5, 10, 20, and 30 years). Peng et al. (2015) performed a risk analysis in the Three Gorges area (China) taking into account the repair costs for roads; the results are expressed in 10 3 ¥ per map cell over a period of 10 years. Mavrouli et al. (2019) presented a procedure aimed at quantitatively estimating (as multiples of a unit cost equalling € 1000) the annual risk to given stretches of a road network in Spain exposed to different dangers, including slow-moving landslides with permanent or episodic activity.
In all encountered cases, the expected losses to or repair costs of exposed roads are computed (i) implicitly assuming that roads themselves are undamaged (as initial condition) and (ii) considering a priori established time intervals. In reality, roads exposed to slow-moving landslides are often already damaged before carrying out the risk analysis; furthermore, the road damage severity (and vulnerability itself) is time-dependent since it may increase as cumulative displacements of interacting slow-moving landslide bodies progressively increase.
These issues are addressed in this paper whose main novelty relies on the joint use of remote sensing and road damage data for the quantitative estimation-at large scale (1:5000)-of slow-moving landslide hazard and related consequences (vulnerability × reconstruction cost) to either an undamaged or a damaged road network exposed to slow-moving landslides. In this regard, the processing of synthetic aperture radar (SAR) images via differential interferometric algorithms (DInSAR) represents a well-established cost-effective non-invasive technique capable of providing displacement time series of affected areas (Antronico et al. 2013;Bianchini et al. 2012;Cascini et al. 2010;Colesanti and Wasowski 2006;Crosetto et al. 2018;Herrera et al. 2013;Tofani et al. 2013;Wasowski and Bovenga 2014), in turn useful to characterize the monitored slow-moving landslides from both geometric and kinematic points of view Calvello et al. 2017;Cascini et al. 2013;Castaldo et al. 2015;Cigna et al. 2013;Di Maio et al. 2018;Frattini et al. 2018;Gullà et al. 2017;Journault et al. 2018;Raspini et al. 2013Raspini et al. , 2017Rosi et al. 2018). More recently, thanks to the increased availability of very high-resolution sensor datasets (i.e., COSMO-SkyMed and TerraSAR-X), the DInSAR data started to be used within procedures aimed at analyzing the consequences induced by slow-moving landslides on both buildings (Bianchini et al. 2015;Ferlisi et al. 2015Ferlisi et al. , 2019aLu et al. 2014;Nicodemo et al. 2017Nicodemo et al. , 2020Peduto et al. 2016Peduto et al. , 2017Peduto et al. , 2018Peduto et al. , 2019 and road networks (Infante et al. 2018(Infante et al. , 2019Nappo et al. 2019;North et al. 2017;Bovenga 2014, 2015). In the latter case, road damage data can be profitably collected by filling ad hoc predisposed fact-sheets during field surveys, provided that a system for classifying the damage severity has been established, as operated in the present work for a road network in the Campania region (southern Italy).

The proposed methodologies
Quantitative risk analysis The proposed methodology for the quantitative analysis of the slow-moving landslide risk to an undamaged road network ( Fig. 1) involves carrying out sequential activities according to

Original Paper
Landslides 18 & (2021) the general framework provided by Fell et al. (2005). In particular, these activities are essentially aimed at (i) characterizing the slow-moving landslides affecting the study area, (ii) analyzing their probability of occurrence, (iii) predicting the consequences to the exposed road network, and (iv) quantitatively estimating the risk (in terms of repair costs to be expected over a given period of time). To these aims, spatial data have to be managed in a geographic information system (GIS) environment .
As shown in Fig. 1, the slow-moving landslides inventoried in an official map are first selected and categorized according to their type and state of activity (Cruden and Varnes 1996); the extent of affected areas is also retrieved. Exposed (or at risk) stretches of road are then identified by overlaying the landslide inventory map with the graph of road network and each of them is later surveyed in order to detect the severity level of visual damage (if any) to be assigned to road sections (within stretches). The latter are linear elements, traced orthogonally to the road centerline at 1:5000 scale, which indicate where the landslide-induced damage-on average-concentrates. Indeed, a given exposed stretch may have more than one section exhibiting either the same or different damage severity levels.
The damage severity is classified based on a system, adapted from the one first proposed by Mansour et al. (2011) and later updated by Mavrouli et al. (2019), including four levels: -D 0 (negligible): road pavement deformation and cracks are absent or rarely visible (Fig. 2a); -D 1 (from very low to low): deformation and cracks locally affect the pavement of road without effects on its functionality (Fig. 2b); -D 2 (from moderate to severe): deformation and cracks substantially affect the road pavement, partly or entirely involving the traffic lanes and/or the roadside, so that the reduction of speed limits is required (Fig. 2c); -D 3 (very severe): deformation and cracks definitively compromise the road pavement continuity, partly or entirely involving the traffic lanes and/or the roadside, so that traffic restrictions are required (e.g., yield on the oncoming traffic) (Fig. 2d).
Given a road section, a buffer symmetrically disposed with respect to the road centerline is introduced in order to associate a recorded damage severity level with a certain value of a representative landslide intensity measure (IM). The buffer width is 40 m in the direction orthogonal to the road centerline, according to the ground resolution of radar sensors and with the general intent to collect as much DInSAR data as possible, whereas the buffer length (along the road centerline) is variable according to the number of damaged road sections and their relative position. In this regard, some indicative examples are shown in Fig. 3.
As far as the IM is concerned, examples in literature dealing with slow-moving landslides analyzed at large scale suggest using a representative velocity along the steepest slope direction (v slope ). In particular, this suggestion is quite in agreement with Mansour et al. (2011) who propose adopting a generic annual displacement rate of the landslide body, without specifying the main direction of movement; on the other hand, Picarelli (2011) emphasizes the role of cumulative displacement (over a given period of time) in the onset and development of damage to exposed facilities. Mavrouli et al. (2019) basically adhere to both proposals assuming as IM a A generic IM value can be retrieved on the basis of DInSAR data, such as those obtained by way of the persistent scatterer interferometry (PSI). The latter is a technique that involves carrying out a multi-interferogram analysis of multi-temporal SAR images in order to extract long-term high phase stability benchmarks of coherent PSI point targets called persistent scatterers (PS) (Costantini et al. 2008;Ferretti et al. 2001). Accordingly, a v slope_k vector (and related modulus) derives from projecting the velocity vector pertaining to the k-th PS within a reference area (e.g., a landslide-affected area on the whole or a part of it, such as the buffer considered in this work) from the sensor line of sight (v LOS_k ) to the steepest slope direction; this implies assuming for each PS a prevalent translational movement (Cascini et al. 2010(Cascini et al. , 2013Cigna et al. 2013;Vecchiotti et al. 2017).
Considering that the LOS projection along the steepest slope direction can be biased by errors (Cascini et al. 2010;Colesanti and Wasowski 2006), each PS must be distinguished according to the own scaling factor, namely the constant value by which the modulus of v LOS_k must be multiplied in order to obtain the modulus of v slope_k . This issue was addressed by Cascini et al. (2013) who observed that, for data acquired on single satellite orbits (either ascending or descending), a scaling factor equalling 3.3 represents an acceptable (upper) threshold to select the most reliable projected PS velocity values. In other words, only those PS whose scaling factor is less than 3.3 are "projectable" and, then, can be used for the following quantitative analyses; the remaining PS have to be discarded and appointed as "not projectable" (see also Plank et al. 2012 andHerrera et al. 2013).
In this work, for those buffers covered by at least one projectable PS in either ascending or descending orbit, the v slope modulus to be associated with a given buffer is computed as the root mean square PS velocity along the steepest slope direction according to the equations (Cascini et al. 2013): in which the weight values are established based on the PS coherence (i.e., the higher the PS coherence, the higher the weight value).
In the Eqs. (1a) and (1b), k refers again to the k-th PS within the buffer; N is the total number of PS within the buffer; w ck is the weight computed on coherence of the velocity of the k-th PS within the buffer; w cN is the sum of w ck ; C max and C min are the maximum and the minimum coherence values of the used dataset, respectively; C k is the coherence value of the k-th PS within the buffer; ε min is a given small number used in order to not discard the k-th PS with C k = C min . In the analyses, ε min value was fixed equal to 0.2 thus assigning a weight of 20% to the smallest coherence value.
The availability of both IM values and (related) road damage data allows for the generation of empirical fragility curves. To this aim, the frequency of occurrence of each level of damage severity is first calculated for different classes of IM values. Then, using methods commonly adopted in different engineering fields (Ferlisi et al. 2019b;Fotopoulou and Pitilakis 2013;Mavrouli et al. 2014;Negulescu and Foerster 2010;Negulescu et al. 2014;Peduto et al. 2017Peduto et al. , 2018Peduto et al. , 2019Pitilakis and Fotopoulou 2015;Saeidi et al. 2009Saeidi et al. , 2012Shinozuka et al. 2000;Zhang and Ng 2005), the probability P() for a road section (randomly selected from a homogenous sample of road sections) to reach or exceed a certain damage severity level (D i ) for a given value of the selected IM parameter is calculated using a cumulative log-normal distribution function: The fragility parameters (median IM i and dispersion β i ) of the standard normal cumulative distribution function Φ[_] are computed using the maximum likelihood estimation method (Shinozuka et al. 2003). This method is based on a statistical approach and requires only binary information (damage or no Fig. 3 Examples showing the criteria adopted to define the length of buffer/s in the case of a one road section (the buffer length L coincides with the stretch length); b two road sections located on the boundary of the slow-moving landslide-affected area (the buffer lengths are equal in value, i.e., L 1 = L 2 ); and c three road sections, two located on the boundary and one within the slow-moving landslide-affected area (L 1 ≠ L 3 ; L 2 = L 1 + L 3 ) damage). In particular, the likelihood function (L) can be expressed as (Shinozuka et al. 2000): In Eq. (3), P(Damage ≥D i IM j j ) represents the probability of reaching or exceeding a certain D i for a given value of IM j (j = 1, …, Y, being Y the total number of IM j values), whereas x j is the realization of the Bernoulli random variable X j , whose value equals 1 or 0 depending on whether or not the road section sustains the considered damage severity level under the IM value equal to IM j . The two parameters IM i and β i can be computed by solving the following equations to maximize ln(L) and hence L: Once the fragility curves are generated, vulnerability curves-relating the selected IM parameter with the average damage (μ D ) expected to the road sections-can be derived by fitting the μ D (IM j ) data obtained as (adapted from Pitilakis and Fotopoulou 2015): where P i (IM j ) is the discrete probability (to be retrieved on the basis of the fragility curves for a given value of IM j ) associated with a damage severity level (D i ) whose numerical index equals d i (taken for this application as 1, 2, and 3 for D 1 , D 2 , and D 3 , respectively). According to Lagomarsino and Giovinazzi (2006) the tangent hyperbolic function can be used as regression model, so that: In Eq. (6), c 1 , c 2 , c, 3 , and c 4 are four coefficients that must be determined for the considered sample of road sections.
The fragility curves also allow one to retrieve the intensity thresholds (T i , i = 1, 2, 3) which represent values of IM in correspondence of a probability of reaching or exceeding a given damage severity level equalling the 5% (Zhang and Ng 2005). This choice can be motivated by the need to limit IM to values compatible with the effectiveness of risk mitigation measures designed according to a performance-based approach that involves reducing the evolution of slope displacements within an estimated time interval, by limiting their magnitude to a target value (Galli and di Prisco 2013). Based on these thresholds, four different ranges of the IM or intensity levels-associated with a nominal scale-are established (Table 1). Then, by carrying out the probabilistic analysis of the IM values used for the generation of fragility curves, the probability of occurrence of slow-moving landslides with a given intensity level (P (SML),d , with d = 0, ..., 3) is retrieved.
The next step involves generating the μ D -RRC curve based on an official price list coeval to the year for which the most recent road damage data are available. The RRC, representing the relative repair cost (ranging from 0 to 1), is given by the ratio between the repair cost and the reconstruction cost (per unit road length, e.g., 1 m).
Finally, the slow-moving landslide risk to the road network (R (RN) ) is estimated on a quantitative basis. This risk corresponds to the repair costs [in €] of a road network identical to the one under study from a starting time t 0 (in which the road network is assumed as undamaged) to a reference time t* (equalling the average time required for the full development of damage-i.e., from D 0 to D 3 -in at least one road section within the observation period). It ranges between a minimum (R (RN),min ) and a maximum (R (RN),max ) value according to the following equations: where L b is the length (in meters) of the b-th buffer covered by PSI data (being A, B, C, D the total number of buffers associated with landslides of either negligible, or low, or moderate, or high intensity level), L DInSAR is the overall length (in meters) of the buffers covered by PSI data, L rs is the overall length (in meters) of the entire sample of exposed road stretches, L rn is the overall length (in meters) of the entire road network. Based on the μ D -RRC curve, the values of RRC d,min and RRC d,max are associated with the minimum (μ D_d,min ) and the maximum (μ D_d,max ) values of μ D for a given intensity level, respectively (d = 0, …, 3). For instance, referring to the intensity level "moderate" and making use of the vulnerability curve, the value of μ D_2,min relates to the intensity threshold T 2 whereas the value of μ D_2,max to T 3 . Furthermore, it must be noticed that RRC 0,min has to be associated with the costs to be incurred in order to guarantee the ordinary maintenance of the road network, whereas RRC 3,max is equal to 1. UC* is the unit cost (UC) of reconstruction to be referred to the time t* (in years) by using the compound interest formula: wherein UC 0 is the UC at t 0 and r is the nominal annual interest rate.

Quantitative consequence analysis
To estimate the repair costs of the road network under consideration (which is already damaged) in a time t greater than t*, one can proceed as shown in Fig. 4.
In particular: (a) the "actual" relative repair cost (RRC actual ) of the road network is first calculated considering the entire sample (X) of road sections with recorded damage by way of the equation: wherein RRC d is the relative repair cost for μ D = d (d = 0, …, 3) and P, Q, R, and S are the total number of road sections (each one associated with a buffer having a length L b ) whose recorded damage equals D 0 , D 1 , D 2 , and D 3 , respectively; (b) then, by using the μ D -RRC curve, the value of μ D, actual corresponding to the calculated value of RRC actual is estimated; (c) this datum (i.e., μ D, actual ), in turn, allows one to retrieve-based on the vulnerability curve-the value of IM that, on average, is representative of the road sections interacting with slow-moving landslides (v slope, actual ); (d) assuming that the landslide bodies are moving according to a v slope that is keeping constant over time, in a time interval Δt = tt* (wherein t > t*, having conservatively hypothesized that the damage actually exhibited by the road network developed in t* years) the average cumulative displacement exhibited by the landslide bodies along the steepest slope direction will be equal to Δs = v slope, actual × Δt; (e) since s = v slope, actual × t* (where s is the average cumulative displacement exhibited by the landslide bodies in a time t*, starting from t 0 ), in order to have at t* an average cumulative displacement equalling s + Δs, the average velocity must be equal to v 1 = (s + Δs)/t*; (f) based on the vulnerability curve, the knowledge of v 1 in value allows one to estimate the average damage expected in a time t (μ D, actual + Δt ) considering an average velocity equal to v 1 ; (g) finally, the value of μ D, actual + Δt can be introduced in the μ D -RRC curve to obtain the value of RRC actual + Δt from which, by subtracting the value of RRC actual , the increment in value (from t* to t) of RRC can be retrieved.
On the other hand, the hypothesis that the landslides are moving with a constant v slope allows for the point-by-point generation of a vulnerability curve referring to the time t based on the available one which refers to t*. To this aim, it is possible to proceed as described above (sub-points d), (e), and (f)) by arbitrarily choosing different v slope values (more values are chosen, more accurate will be the estimate).

Case study: Road network and available dataset
The analyzed road network develops within a territory (Fig. 5) extending for about 1600 km 2 and located within the national park of "Cilento, Vallo di Diano, and Alburni" (south-western part of the Campania region, southern Italy).
From a geological point of view (Fig. 5a), owing to the long-time and complex lithogenetic and orogenetic history, several lithostratigraphic units in form of nappes and/or irregular sequences can be distinguished in the territory under consideration. As deeply discussed by Santangelo et al. (2005), such units can be first categorized into internal units (mainly constituted by marly calcarenites, calcilutites, clay, often siliceous, sandy clays, sandstones, and conglomerates) and external units (including carbonate and terrigenous sediments), according to the original position before the tectonic deformation. One of the most widespread group in the study area is represented by the neogenic synorogenic units of the Miocene age (made up of clays, sandstones, and conglomerates with wild-flysch facies) deposited in a basin formed on top of advancing thrusts. Finally, the quaternary postorogenic units include continental and marine sediments whose deposition took place after the final emersion of the area in the Late Pliocene-Early Pleistocene. The territory is studded with small towns, some of them representing renowned touristic spots especially in the summer season; globally, it counts 59 municipalities (Fig. 5b).
As for the road network, it mainly consists of urban roads and suburban secondary roads composed by single carriageways with two lanes (one per each traveling direction). The graph of the analyzed network is shown in Fig. 5c, whereas the length of the composing (either State or former State) roads are synthesized in Table 2.
For the territory shown in Fig. 5a, a landslide inventory map at 1:5000 scale (Fig. 6a) is available as a result of the activities carried out in 2012 by the former "Sinistra Sele" River Basin Authority within the Hydrogeological Setting Plan -Landslide Risk excerpt (Italian Law 365/2000). It can be noticed that slow-moving landslides (summing up to 14,843, representing the 83% out of the total inventoried slope instabilities) can be distinguished among rotational/translational slides, lateral spreads, (earth) flows, deep-seated gravitational slope deformations (DGSD), and creep phenomena (Cruden and Varnes 1996;Hungr et al. 2014). The extent of affected areas is shown in Fig. 6b, whereas Fig. 6c shows the number and state of activity of different inventoried slowmoving landslide types.
As for DInSAR data, the available dataset is provided by the Italian "Ministry of the Environment and Protection of the Territory and the Sea." It results from the PSI processing of images acquired by COSMO-SkyMed (very high-resolution) radar sensors on both ascending (42 images from May 2011 to March 2014) and descending (42 images from October 2011 to December 2013) orbits. The spatial distribution over the study area of PS is shown in Fig. 7a for the ascending orbit and in Fig. 7b for the descending orbit, along with the corresponding values of the average annual velocity (in mm/year) computed along the LOS of the radar sensors.
As for road damage data, they were collected by filling-in ad hoc predisposed fact-sheets during field surveys. These fact-sheets are comprised of different sections that allow gathering information on the location of the exposed road stretch, the geological context, the slow-moving landslide type and its state of activity, the recorded damage (with explanatory photos) and its severity level, the PSI data in terms of displacement time series along the steepest slope direction (Fig. 8).
Further damage data were collected by managing the images available in the Google Street View archive for all the considered road sections (with one image per section at least) dating back to September 2008 (since the damage surveys ended in July 2018, this implies that the observation period spans 9.9 years). An example is shown in Fig. 8 with reference to a road section in the municipal territory of Pollica along the former State road 267 damaged by a rotational slide. In this case, Google Street View provides two images that allow us to detect an increase in the damage severity (from D 1 to D 3 ) within an 8-year period. The same road section is

Results
According to the procedure shown in Fig. 1, by overlapping the landslide inventory map (Fig. 6a) with the graph of the road network (Fig. 5c), 549 road stretches were first identified and later surveyed (from March 2018 to July 2018) in order to detect the related damage and classify its severity. Among these stretches, 102 crossing urban centers were discarded because their damage was not straightforwardly attributable to landslide movements or had been caused by tilting/failure of the structures (e.g., walls, sheet piles) retaining the road embankment. As far as the remaining 447 road stretches are concerned, a database including X (486) road sections whose damage severity covers all the considered levels (D 0 = 321, D 1 = 58, D 2 = 79, D 3 = 28) was generated (Fig. 9a). The damage mainly occurs in correspondence of road sections located on stretches intersecting the landslide body or their head, whereas few cases refer to stretches in correspondence of the landslide foot, mainly for damage severity levels ranging from D 1 to D 3 (Fig. 9b). Overall, referring to the position of the stretches with respect to the landslide-affected areas, the damage survey led to the following percentages of stretches that exhibit a damage severity exceeding or not exceeding the D 0 level: head (38.7% exceeding, 61.3% not exceeding); body (38.8% exceeding, 61.2% not exceeding); foot (9.0% exceeding, 91.0% not exceeding).
A subset x (318) of X was not considered for analysis purposes on the basis of well-established criteria (e.g., road sections not covered by DInSAR data or for which DInSAR data interpretation was not reliable since, for instance, the scaling factor of the projection operation from LOS to the steepest slope direction exceeded the fixed threshold).  The sample Y = X − x (168) includes road sections whose damage severity level (D 0 = 67, D 1 = 29, D 2 = 48, D 3 = 24) can be associated with a value of IM. To this aim, given a road section (and the related buffer), the v LOS_k velocity pertaining to each PS associated with the COSMO-SkyMed image dataset was first projected along the steepest slope direction according to Cascini et al. (2010Cascini et al. ( , 2013; then, the v slope value was obtained by using the Eqs. (1a) and (1b). Figure 10a shows that as v slope values, on average, increase the level of damage severity increases as well. This aspect is even clearer in Fig. 10b where it can be noticed that the road sections with (i) D 1 damage severity have v slope values mostly lower than 8 mm/year, (ii) D 2 damage severity have v slope values ranging from 4 to 24 mm/year, and (iii) D 3 damage severity are associated with v slope values ranging from 8 up to more than 32 mm/year.
Referring to the Y sample, the empirical fragility curves (Fig. 10c) and the vulnerability curve (Fig. 10d) were generated. As for the fragility curves, it is worth stressing that the one associated with the D 1 damage severity represents the probability that "at least a damage from very low to low" will be sustained by an exposed road section, arbitrarily chosen from the Y sample, Fig. 9 Numerical distribution of the surveyed road sections a according to the recorded damage severity level and b taking into account the position of corresponding stretches with respect to the slow-moving-landslide-affected areas Fig. 10 a Damage severity level recorded to the surveyed road sections vs. v slope , b class frequency of recorded damage severity levels according to v slope , c empirical fragility curves, and d empirical vulnerability curve when subjected to a given v slope value; the same meaning applies to other fragility curves (Shinozuka et al. 2000). Figure 10c shows that up to v slope values of about 10 mm/year the probabilities of reaching or exceeding D 1 or D 2 are both over 90% whereas D 3 does not exceed 10%. The latter reaches the highest probabilities for v slope values of about 30 mm/year. Table 3 shows the values of the fragility parameters. Based on the vulnerability curve (Fig. 10d), it can be observed that v slope values of about 5, 10, and 30 mm/year can be associated with numerical indices of the average damage expected to the road sections (μ D ) equalling 1, 2, and 3 respectively. The computed values of the coefficients of the vulnerability curve are summarized in Table 4.
The fragility curves allowed retrieving the threshold values (T i ) of IM, here assumed as the values of IM in correspondence of 5% probability of reaching/exceeding a given damage severity level; accordingly, we obtained: T 1 = 1.1 mm/year; T 2 = 3.5 mm/year; T 3 = 8.9 mm/year. Therefore, four different intensity levels were defined (see also table in Fig. 11). Then, the probabilistic analysis of the IM values used for the generation of the fragility curves was carried out. The distribution law that best fits the available data (Fig. 11a) is the negative exponential one (Fig. 11b). As shown in Fig. 11c, the values of P (SML),d (d = 0, ..., 3) pertaining to each intensity level can be easily obtained on the basis of the cumulative distribution function of the random variable IM.
The next step dealt with the generation of the curve relating the expected average damage (μ D ) with the relative repair cost (RRC). To this aim, we used the price list provided by the "Ente Nazionale per le Strade" (ANAS 2018a, b) which allowed us to retrieve the repair costs (to be associated with μ D equalling 0, 1 and 2), and the reconstruction cost (to be associated with μ D equalling 3) per unit road length (1 m). These costs are synthesized in Table 5 along with the corresponding RRC values, whereas Fig. 12 shows the obtained μ D -RRC curve. For the sake of simplicity, this curve was assumed as time-independent; this means that all cost items listed in the third column of Table 5 increase (or decrease) by the same rate with time.
Based on the μ D -RRC curve, RRC d,min and RRC d,max (d = 0, …, 3) values were estimated provided that corresponding μ D_d,min and μ D_d,max (d = 0, …, 3) values had been previously determined by way of the vulnerability curve, for a given intensity level. The obtained values are synthesized in Table 6.
Focusing on the 24 road sections recording a D 3 severity level based on the damage survey, the availability of antecedent Google Street View images allowed us to retrieve the reference time t* (equalling 7.4 years) by averaging the time required by each of the above road sections to move from D 0 (if any) to D 3 . In this regard, 10 road sections out of 24 exhibited the full development of damage within the observation period. Then, the UC* value (equalling € 254.94) was computed according to Eq. (8) by applying a nominal annual interest rate whose percentage value equals 1.56% on the basis of data provided by the Italian Institute for Statistics (ISTAT 2019) about the change in construction costs of roads over time (from 2005 to 2017); on the other hand, UC 0 was posed equal to the unit cost at the time when the damage survey was carried out, which can be computed on the basis of the price list of ANAS (2018a, b) so obtaining UC 0 = € 227.35 (see also  Table 5).
Finally, according to Eqs. (7a) and (7b), the risk was estimated; in particular, it ranges between a minimum (R (RN),min ) and a maximum (R (RN),max ) value equalling respectively: To carry out the quantitative consequence analysis (QCA) with reference to the road network under consideration (which is already damaged) in a time t greater than t*, we proceeded based on the framework shown in Fig. 4.
The "actual" relative repair cost (RRC actual ) was calculated by using Eq. (9), taking into account the entire sample (X) of road sections. In particular, the obtained value of RRC actual (equal to 0.178) was estimated considering the total length of buffers distinguished according to the recorded damage severity level (Table 7).
Based on the curve shown in Fig. 12, a μ D, actual equalling 1.02 was obtained in correspondence of the computed value of RRC actual ; then, the vulnerability curve (Fig. 10d) allowed us to retrieve a v slope, actual = 4.9 mm/year associated with μ D, actual .
Assuming a forecasting time interval Δt = tt* = 5 years, for example, the average cumulative displacement of the slow-moving landslide bodies expected during this Δt equalled Δs = v slope, actual × Δt = 24.5 mm; on the other hand, s = v slope,actual × t* = 36.2 mm.
Once the parameters s and Δs were known, it was stated that, in order to have at the time t* an average cumulative displacement s + Δs, the average v slope must be equal to v 1 = (s + Δs)/t* = 8.2 mm/year. In turn, v 1 allowed us to estimatebased on the vulnerability curvethe value of the expected average damage at the time t = t* + Δt = 12.4 years; in particular, μ D, actual + Δt = 1.75.
The latter value was introduced again into the curve of Fig. 12 to retrieve the value of RRC actual + Δt = 0.417 from which, by subtracting the RRC actual value, the increment of RRC from t* to t (ΔRRC) was finally obtained: Based on the above presented analytical steps and taking account of the methodological details provided in the sub-section entitled "Quantitative consequence analysis," μ D was estimated for different time periods greater than t*. In this regard, Fig. 13 shows the vulnerability curves obtained with reference to four time periods, including the one associated with t*. These curves allow highlighting how the vulnerability of the road network increases as cumulative displacements of interacting slow-moving landslides (and the induced road damage) increase. For instance, considering a v slope value equal to v slope, actual (4.9 mm/year), the μ D value progressively moves from 1.02 (at t* = 7.4 years) to 1.75 (at t* + 5 years) or to 2.30 (at t* + 10 years) or to 2.64 (at t* + 15 years).

Discussion and conclusions
This paper showed the results of a study aimed at quantitatively estimating-at large scale (1:5000)-the risk (QRA)/direct consequences (QCA) to an undamaged/damaged road network exposed Fig. 11 a Class frequency, b probability density function (PDF), and c cumulative distribution function (CDF) of v slope . The probabilities of occurrence of slow-moving landslides of a given intensity level are summarized in the included table Thickness of different layers refers to secondary suburban roads according to Domenichini et al. (1993); the road width was posed equal to 7 m to slow-moving landslides. With reference to a case study in the Campania region (southern Italy), the analyses benefited from the availability of very high-resolution DInSAR data and information gathered from in situ surveys on the severity of damage to selected road sections. The proposed methodological approaches allowed us to successfully address some open issues in the scientific literature concerning the hazard analysis (e.g., in terms of probability of occurrence of slow-moving landslides with a given intensity level, in turn established based on empirical fragility curves) and the consequence analysis (e.g., in terms of time-dependence of the vulnerability curves).
The QRA results were expressed in terms of repair costs (in €) ranging between a minimum and a maximum value, to be referred to the average time required for the full development of damage in at least one road section. The obtained range of risk values takes into account the high degree of uncertainty that, at large scale, affects the QRA. In this regard, uncertainties can be associated with the errors inherent to slow-moving landslide mapping, the subjective nature of the road damage severity assessment (which is based on the judgment of the expert carrying out the in situ surveys or interpreting the Google Street View images), the lack of knowledge on precise values of cumulative displacements required for road sections interacting with slow-moving landslides to move from a damage severity level to the next one, the unavailability of DInSAR data in both ascending and descending orbits, the variability of mechanical characteristics of materials forming the roads from one exposed section to another. Anyway, the QRA results can be profitably used for risk management purposes. For instance, the potentially involved decision-makers could decide that the risk-to which the road network under study is exposed-has to be mitigated by reducing the hazard of the DInSAR-covered road sections associated with a high-intensity level, independently of the recorded damage severity, by way of slope stabilization works. In particular, the latter should be designed on the basis of a performance-based approach that involves immediately limiting IM to values not exceeding the T 1 threshold (1.1 mm/year), which then remain constant over a certain time period (e.g., 30 years, corresponding to the nominal life of interventions). Should all damaged road sections have been repaired, the use of Eqs. (7a) and (7b) allows one to obtain maximum and minimum residual risk values respectively equalling € 753,147 and € 233,255 (keeping fixed the intensity threshold values and posing P (SML),3 = 0). Accordingly, the implementation of the above risk mitigation strategy could allow for a reduction of risk (from preto post-interventions) equalling 67.8% (for the maximum value) and 76.2% (for the minimum value).
On the other hand, if the decision about the risk mitigation is postponed, the methodological approach adopted for QCA can turn out to be useful. For instance, looking at Fig. 13, the potentially involved decision-makers could decide to wait 5 years before implementing the interventions (or, similarly, 5 years after 2018, when damage surveys were carried out). This would imply that the reconstruction cost of all exposed road stretches (whose total length equals 59,137 m) moves from € 13,445,016 (when UC = € 227.35) to € 14,526,736 (when UC = € 245.65, taking account of the nominal annual interest rate). Since μ D increases from 1.02 to 1.75 in 5 years-and, accordingly, the RRC from 0.178 to 0.417-the repair cost of the road network would change from € 2,392,234 to € 6,058,039 with a percentage increase equalling the 153.2%.
Of course, the above decisions on whether and when a risk mitigation strategy has to be implemented requires carrying out a wider study including the quantification of the indirect (social and economic) consequences associated with the slow-moving landslide-induced road damage. This could allow also for the vehicular traffic control before and during the implementation of interventions. However, estimating indirect consequences is beyond the scope of this paper.
Further refinements of the work done deal with carrying out QRA/QCA at a detailed scale (> 1:5000), provided that more accurate data (both in quality and quantity) are collected for single   slow-moving landslides and affected road stretches. This could lead to the prioritization of risk mitigation measures and their proper choice/design. Finally, it is worth mentioning the potential exportability of both the proposed methodological approaches and the obtained results (i.e., fragility/vulnerability curves) that, once further validated, could stand as a reference for QRA/QCA concerning other similar road networks in similar geo-environmental contexts that are widespread in southern Italy.

Acknowledgments
The processed DInSAR data were provided by the Italian Ministry of the Environment and Protection of Land and Sea within the nationwide Italian "Piano Straordinario di Telerilevamento". The authors wish to thank the editor and two anonymous reviewers who helped improve the paper in its revised version.

Funding Information
Open access funding provided by Università degli Studi di Salerno within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/ 4.0/.