Rockfall susceptibility and network-ranked susceptibility along the Italian railway

Rockfalls pose a substantial threat to ground transportation, due to their rapidity, destructive potential and high probability of occurrence on steep topographies, found along roads and railways. Approaches for assessment of rockfall susceptibility range from purely phenomenological methods and statistical methods, suitable for modeling large areas, to purely deterministic ones, usually easier to use in local analyses. A common requirement is the need to locate potential detachment points, often found uphill on cliffs, and the subsequent assessment of the runout areas of rockfalls stemming from such points. Here, we apply a physically based model to calculate rockfall trajectories along the whole Italian railway network, within a corridor of total length of about 17,000 km and varying width. We propose a data-driven method for the location of rockfall source points based on expert mapping of potential source areas on sample representative locations. Using empirical distributions of gridded slope values in source areas mapped by experts, we derived probabilistic maps of rockfall sources in the proximity of the railway network, regardless of a particular trigger. Source areas act as starting points of simulated trajectories in the three-dimensional model STONE. The program provides a pixel-by-pixel trajectory count, covering 24,500 km2, the largest homogeneous application of the model to date. We classified the map into a vector susceptibility map of the segments of the railway, for which we provide segment-wise rockfall susceptibility. Eventually, we considered a graph representation of the network to classify the segments both on the basis of rockfall susceptibility and the role of each segment in the network, resulting in a network-ranked susceptibility. Both maps are useful for subsequent hazard assessment, and to prioritize safety improvements along the railway, at national scale.


Introduction
Rockfall susceptibility along a transportation corridor is a recurring topic in landslide research. Conclusions based on quantitative studies are relevant for practical purposes, providing key information for risk mitigation. They may allow planning of protection barriers in key spots to optimize available resources and, as in the case motivating this work, they provide a quick estimate of the possible critical locations along a transportation network after a triggering event.
A susceptibility map does not contain information about the magnitude of possible events, nor their expected frequency over a period of time. Aiming to such a complete description (i.e., a hazard map) would require knowledge of the processes responsible for the mobilization of individual blocks or rock masses at the local scale. The input to prepare a susceptibility map is less demanding, though it would still not be trivial to gather the required data homogeneously over a large area, as we aim to in this work. The minimum requirement is knowledge of the point of origin of the falling mass, the underlying topography, and a model for the description of how the mass would reach a rest state, downhill. The majority of existing rockfall studies made use of physically based models, for they give a more accurate description of the whole process, and only a few examples exist of statistical models for rockfall susceptibility (Frattini et al., 2008).
Phenomenological (easier) ways to relate rockfall sources to their arrival point exist, for example using an empirical shadow angle (i.e., the angle between the distal limit of the shadow and the top of the talus slope) (Evans and Hungr, 1993). Here we committed to using the model STONE by Guzzetti et al. (2002). The model requires a digital elevation model (DEM) of the area, predetermined source locations, and numerical parameters controlling the energy restitution at the location of impacts of boulders with the topography. We used the trajectory count per grid cell produced by the model to calculate rockfall susceptibility.
To this end, in this paper we (i) applied a novel, datadriven method to specify source areas for rockfall trajectories (Alvioli et al., 2020b), based on expert-mapping of source areas in sample locations and statistical generalization across the whole study area; (ii) performed trajectory count simulation in a probabilistic way, using the model STONE, which we eventually classified as a susceptibility map along the national railway; (iii) performed a network analysis, within graph theory, trying to infer the impact of possible events on the functionality of the network as a whole. A national rockfall inventory consisting of more then 70,000 rockfall polygons helped validating the results.
The paper is organized as follows. Section 2 gives some context for this study, both for physically based rockfall modeling and for network analysis within graph theory. Section 3 describes the data used in this study. Section 4 describes the core method, and it is split into four separate paragraphs describing data preparation, how we singled out potential rockfall source locations, physically based rockfall susceptibility and analysis of the railway network within graph theory; a short paragraph describes the computational demand of this work. Results and discussion are split in a similar way, in Sections 5 and 6, respectively. Conclusions are drawn in Section 7.

Physically based rockfall modeling
Relevant stages of rockfalls are detachment of a rocky mass from a slope (usually a cliff face), typically by sliding, toppling or falling, and subsequent movement down slope, through a process which can assume different dynamics. These phenomena can vary among a rather broad spectrum, with two limiting cases represented by independent falls of individual fragments with no interactions among each other, and flow of granular masses with or without interaction with the substrate (Bourrier et al., 2013). Phenomena which do not fall strictly under the two limiting cases contain elements of both, and are the most difficult to model.
In this work we did not adopt a physical model for the detachment process. To our knowledge, no model exists that would cover the multifaceted mechanism of initiation of a rockfall (Collins and Stock, 2016;Matasci et al., 2018). One possibility to avoid a detailed description of this step would be direct mapping through observation of the rocky slopes, for example by laser scanners (Jaboyedoff et al., 2012;Riquelme et al., 2014;Li et al., 2019), photogrammetry (Buyer et al., 2020), and similar approaches combined with remotely piloted aircrafts (Peréz-Rey et al., 2019;Santangelo et al., 2019;Giordan et al., 2020). Such measurements and algorithms may provide detailed information about the fractured status and/or discontinuities of rocks, block orientations and distribution of joints, giving insight to infer locations of possible detachment areas (Turanboy et al., 2018;Menegoni et al., 2021;Mavrouli and Corominas, 2020).
The methods outlined above to identify potential rockfall sources are only applicable on small areas, and probably on individual slopes. Expert-mapping of possible source areas by photo-interpretation (Santangelo et al., 2019 is another approach which can be very effective, and can be applied to relatively larger areas. In this work, we aim to a very large area assessment, so none of these methods could be applied. A promising method to this end could be a multivariate statistical regression technique as the one outlined by Rossi et al. (2021), which was not implemented here, since we selected a data-driven method (Alvioli et al., 2020b), described in the following sections.
Here we only deal with individual masses following independent trajectories. These can be described by ballistic dynamics, during free fall, and bouncing or rolling on the ground. In particular, the model STONE adopted here assumes point-like boulders under the only action of gravity and the constraints due to topography and geology. The computer program STONE models the motion of three possible states in a number of time steps in a threedimensional fashion along the trajectory, until the motion comes to an end due to friction exhausting the initial kinetic energy.
Types of movements of rocky masses other than ballistic trajectories interrupted by bouncing or rolling are not included in this work, as they require models suited to describe avalanches and granular fluxes, while STONE does not include interactions between boulders, neither it accounts for their fragmentation. Recently, Matas et al. (2020) and Ruiz-Carulla et al. (2020) studied in detail the fragmentation mechanism of large boulders by directly observing blocks debris after they impacted the slope. Data about the angular distribution of the fragments helped calibrating numerical parameters of the model RockGIS (Matas et al., 2017), which implements a fractal model for the mass of the fragments, which is absent in STONE.
To go beyond a spatial susceptibility assessment, one should consider the magnitude of possible rockfalls (Guzzetti et al., 2003(Guzzetti et al., , 2004, in a given area, and their temporal recurrence. Magnitude may vary in a wide range of volumes, from cubic decimeters to millions of cubic meters (Hungr et al., 1999). Melzner et al. (2020) also investigated frequency-size distributions, finding a dependence on the method used to map rockfalls and, similarly to Corominas et al. (2018), recommended caution when considering the large-size region of the distributions. We did not consider the magnitude of rockfalls, here, as the model STONE does not directly account for that.
We also did not account for a specific trigger: as a matter of fact, rockfalls may be triggered by different phenomena, mostly by intense rainfall and seismic activity, with their occurrence influenced by a variety of factors (Macciotta et al., 2015). Static factors are steep terrain gradient and structural setting, and dynamic factors are thermal deformations and or freezing-thawing cycles of the rocks. Studies exist about correlations between rockfalls and rainfall, for example by means of probabilistic rainfall thresholds , or studies of the relation between rockfalls and rain, freezing periods, and strong temperature variations (Delonca et al., 2014). Inclusion in our method of a trigger-or time-dependent description is also beyond the model STONE, thus it falls out of the scope of this paper and needs to be considered elsewhere.

Analysis of the railway network within graph theory
Graph theory is the mathematical representation of pairwise relations between a collection of objects, and a whole science of graphs exists with examples in a broad spectrum of theoretical and applied research (Chartrand et al., 1993). One of major field of application of graph theory is the representation of a physical network. Literature on network and/or graph science is extensive, and it deals with a variety of phenomena. Examples are biological and social networks, power grids, telecommunications and, obviously, the world wide web and transportation networks (Barthélemy, 2011). Studies related to transportation networks often involve urban settings, including interdependent and/or multi modal networks of different means of transportation, and communication networks also influence the structure of cities (Alvioli, 2020).
The effect of disruption of a number of vertices or edges was also investigated in the literature; Mattsson and Jenelius (2015) contains a review on the resilience of transportation networks. A few general results are that redundancy in a complex network is not sufficient for tolerance upon failure, while tolerance is much higher for networks exhibiting scale invariance (Albert et al., 2000) as for the world wide web (Cohen et al., 2000). On the other hand, scale-free networks are subject to disruption upon removal of a few of the relevant vertices. Features of networks relevant for their resilience are the presence of loops (Katifori et al., 2010) and vulnerable spots (Bassolas et al., 2020).
In this work, we use basic notions of graph theory, useful to make sense of the relevance of a particular edge of the graph within the entire graph. We refer to "network" as the collection of physical links (track segments) and of nodes (stations or intersections), and to "graph" as the equivalent representation of edges and vertices. The elements of the physical network are associated with specific geographical locations, while the associated graph is an abstraction that only retains the topology of the physical network and, optionally, a few other quantities.
We further refer to betweenness centrality of vertices (Freeman, 1978) and edges (Lu and Zhang, 2013), a measure of centrality of each such element in the graph. Betweenness centrality (in the following, betweenness) of a vertex is given by the number of shortest paths going through the vertex, out of the total set of shortest paths connecting any two pair of vertices in the graph (all-to-all routes). In a similar way, we can define betweenness of edges. Values of betweenness of vertices and edges depend on the weights of edges in the graph. As the graph considered here is associated to a physical network, we weight edges with the actual length of railway links. Centrality of a vertex or an edge in the network gives a measure of its relevance for functioning of the whole network.
Eventually, we devised a new (simple) centrality index intended to quantify the impact on the network when one edge is removed, which mimics an interruption due to an external event. The new index helps performing a combined classification of the railway segments, based simultaneously on centrality of the segments and their rockfall susceptibility. We call the new classification network-ranked susceptibility. A similar approach, namely calculation of network risk assessment combining exposure to link failures and identification of link importance, was proposed by Cats et al. (2016) for a multi-modal public transport network in the Netherlands.

Available Data
This work involved both material which was already available to us and newly developed maps. Among the existing material, this study required (i) a digital elevation model, needed for the whole processing chain; (ii) a nation-wide inventory of rockfalls; (iii) a small-scale lithological map; (iv) a vector map of the railway network. Newly developed maps were a nation-wide slope unit and sample rockfall source area maps. Here follows a short description of such pieces of information.
• Digital elevation model. The highest resolution DEM freely available for the whole of Italy is TINITALY (Tarquini et al., 2007). The grid is 10 m x 10 m, in the reference system with EPSG:32633. Figure  1(a) shows a shaded relief map based on the 25 mresolution EU-DEM; its use in this work is limited to this figure.
• Inventory of landslide phenomena in Italy, known as IFFI (Trigila et al., 2010;ISPRA, 2018). The IFFI inventory contains over 620,000 landslide polygons, of which 70,576 are classified as rockfalls, and 4,051 are located in the vicinities of the railway track and were used in this work. Rockfall polygons source and runout areas.
• Vector layer containing the railway network track. Table 2 shows the distribution in elevation of the various links of the network, and Fig. 1(c) shows the slope units selected from this work along the railway track. The figure also shows the track on plain areas, where slope units are undefined.
• A map of sample potential source area for rockfalls in Italy, from expert interpretation of Google Earth™ images. Section 4.2 contains a description of the rationale and the criteria used during expert mapping.

Methods
Performing nation-wide rockfall susceptibility assessment, along the railway track, required careful planning. The long chains of actions, from data selection to the preparation of the final map, is in fact prone to propagating early-stage mistakes in unpredictable ways. Selection of STONE as a rockfall runout modeling software (Guzzetti et al., 2002;Agliardi and Crosta, 2003) was dictated by the need to compromise between a well-defined relationship between sources and runout, and overall computational demand.
To prepare a susceptibility map based on the results of the model STONE we adopted the following steps: (i) preparation of the data, (ii) selection of possible source areas, which is probably the most important input of the model, (iii) execution of the program, (iv) collection and classification of the results analyzing the susceptibility results from the perspective of the properties of the network. The flowchart in Figure 2 illustrates the procedure; the following paragraphs describe in detail all the steps performed in the analysis.

Preparation of data
We adopted the topographic subdivision, previously defined by Guzzetti and Reichenbach (1994) and slightly modified for slope units delineation (Alvioli et al., 2020a). The original subdivision contained 30 sections, obtained by quantitative analysis and visual interpretation of morphometric variables obtained from a digital elevation model. In the modified version, we neglected topographic units much smaller than the surrounding ones, and split one  Table 1. After (Fongo, 2018;Bucci et al., 2021). (c) The spatial buffer along the railway network adopted in this work. The buffer consists of the subset of slope units from Alvioli et al. (2020a) intersecting a buffer of 1 km around the whole railway track, depicted in brown. Slope units are in blue. The total number of selected slope units is 32,200, out of the over 330,000 polygons, for a planimetric area of 25,400 km 2 (cf. into two different ones, ending up with 29 topographic units, shown in Fig. 1(a). Use of topographic units was dictated by the need of similar terrain types, allowing to use the same soil parameters in each unit, for numerical simulations.
Utilization of topographic domains also allowed us to define potential source areas of rockfalls using the statistical generalization proposed in this work separately within the 29 sections, corresponding to 29 different terrain types. A lithological map of Italy at the 1:100,000 scale provided the necessary insight for the determination of different sets of input parameters for the software STONE in different lithological classes. These parameters control the behaviour of boulders during rolling or bouncing on the ground, along the simulated trajectories. Their final values were inferred from previous studies. Table 1: Numerical values of the parameters used in STONE. Figure 1(b) shows the corresponding lithological map (Fongo, 2018;Bucci et al., 2021). STONE performs random sampling of values of the parameters in a ± 5% range around the nominal values listed here. A second, much finer subdivision of the territory was adopted here, consisting in groups of slope units intersecting the railway. Use of slope units was both a natural choice, given that slope units are mapping units well suited for landslide studies, and a solution to the technical problem of defining meaningful sub-areas for execution of the simulations in parallel. We favor SU as suitable mapping units for the description of landslide phenomena (Alvioli et al., 2016;Camilo et al., 2017;Schlögel et al., 2018;Bornaetxea et al., 2018;Tanyaş et al., 2019a,b;Jacobs et al., 2020;Amato et al., 2020;Chen et al., 2020;Li and Lan, 2020). More specifically, we adopted SUs instead of a geometric buffer around the track, because rockfall trajectories initiated in a given SU will be bounded within the SU, with reasonable confidence. The set of SU used here were specifically calculated and optimized over the whole of Italy (Alvioli et al., 2020a) for this work, and for a series of additional works regarding selection of source areas of debris flows  and soil slides susceptibility (yet unpublished), along the railway network. The subset of SUs selected for this work is shown in Fig.  1(c). The selected polygons (cf. point 2a) cover an overall area of 25,400 km 2 .

Class ID
We prepared all of the data splitting the study area into sub-areas and working within sub-areas in parallel, computing-wise. In particular, for each group of SUs singled out as described above, we prepared an independent simulation with STONE. This step requires selecting the relevant portion of DEM, of sources areas and of grids containing the numerical values of parameters needed by STONE (Table 1).

Figure 2:
A flowchart of the procedure followed in this work to obtain maps of rockfall susceptibility and network-ranked susceptibility. The different steps involved in the procedure are described in detail in Section 3 (Available data), Section 4 (Methods) and Section 5 (Results).

Identification of rockfall source areas
Identification of potential rockfall source areas is a key step of physically based simulation of rockfall runout, because they substantially influence the extent of runout. A detailed investigation of the location of potential sources of rockfalls requires expert analysis of the cliffs in the study area (Guzzetti et al., 2004;Santangelo et al., 2019Santangelo et al., , 2020, which is typically a time consuming and expensive procedure. This makes identification of potential sources a limiting factor for systematic rockfall studies over large areas. Visual mapping of all of the potential sources by photointerpretation over the whole of the railway network in Italy would be impossible, in a reasonable time. Moreover, once a potential source is located on a DEM, it is desirable to assign a probability for the likelihood of that location to evolve into a rockfall. A straightforward way of selecting source areas for STONE is to establish a slope threshold above which any grid cell acts as a potential rockfall source (Guzzetti et al., 2003). Other criteria are based on a combination of geomorphological analysis and mapped sources (Agliardi and Crosta, 2003), a combination of a slope threshold and geomorphological analysis , or a multivariate statistical analysis (Rossi et al., 2021). In this work, we adopted a method to both locate sources and assign a probability of failure in a homogeneous way over a large area, on the 10 m-resolution DEM of Italy (TINITALY, Tarquini et al. (2007)), according to the following steps: 1a. selection of a 1 km-wide buffer region around the railway track; 2a. selection of the subset of slope units (SU) intersecting the buffer, out of the whole set of about 330,000 SUs calculated for Italy by Alvioli et al. (2020a); the subset of selected SUs is the final working buffer used here; 3a. expert mapping of potential rockfall source areas within the selected SUs, in sample locations we considered representative of the conditions that could trigger rockfalls in the specific topographic unit (Guzzetti and Reichenbach, 1994) under investigation; 4a. development of a data-driven procedure to identify source areas in different grid cells, with respect to the areas mapped as in 3a within the same topographic unit, using statistical generalization; 5a. visual analysis of the source areas map obtained from point 4a and assimilation in the final source areas map of additional sources in locations chosen in expert way.
Selection of the SU overlapping with a 1-km buffer around the railway track (cf. point 1a) ensures that the simulations will include all of the potential sources that could develop into rockfalls and whose runout could concern the railway network (this work does not attempt to account for any anthropogenic mass displacements in the vicinities of the railway track).
Expert mapping was performed on Google Earth™ images. Image interpretation to detect and map potential rockfall source areas can be carried out using different data and tools. We chose Google Earth™ images because they provide simultaneously morphological (elevation) and optical information, resulting in 2.5D images. They are less accurate than stereoscopic optical images, but the advantage is completeness and free availability at the national scale. This is far better than using 2D images in GIS, which have the same image availability as Google Earth™ , but come with no morphological information. The vector format editing tool allows mapping directly on the interpreted images, and exporting the polygons.
To map potential rockfall source areas, the geomorphologists considered that they correspond to sub-vertical rocky slopes showing talus or rockfall deposit underneath. They usually appear not vegetated, even if in some locations tall vegetation can hide potential detachment areas. For each topographic unit, the geomorphologists mapped all of the potential rockfall source areas within each SU in a significant subset of SUs, which were considered representative of the main litho-structural and morphological settings.
We used the information of expert-mapped sources to calculate the probability of any grid cell to be a source as follows. For a given topographic unit, we selected SUs containing mapped polygons. We calculated histograms of the distribution of values of slope angle in the cells within the mapped polygons, and within the whole slope unit. For each bin of 2°width, we took the ratio of the two histograms: values of slope within the sources over values within the whole slope units. The values of the numerator are smaller than at the values of the denominator, in each bin, by construction (the number of cells belonging to source areas falling in each bin are less or equal to the number of cells, in the same bin, from the whole slope unit). Since we took care of delineating all of the possible rockfall sources within each of the considered slope units, we are confident that the ratio represents the probability for a cell with slope in the 2°bin range of being a rockfall source, in that specific slope unit. This way, we have a probability point for each 2°slope bin. Considering all of the sampled SUs in the given topographic unit, we end up with a collection of probability values as the green dots in Fig. 3, for TU 1.2 & 1.3, which we chose as an illustrative  Table 1. The red curve, a non-linear 90% quantile regression corresponding to Eq. (2), is the adopted model. The horizontal line represents the limit under which probability is set to null. example.
The purpose of the procedure is to obtain probabilityslope relations which can be extended to the entire set of selected SU, in each topographic unit. Green dots in Fig. 3 represent the probability values from expert mapping, and curves show three possible approaches to define a functional dependence between probability and slope. The black histogram is the simple average, in each slope bin, of the observed probability values. The blue curve is a non-linear fit of the values estimated from observed data, with a function of the following form: where S is the slope in degrees, and a, b are parameters of the fit. The band around the fit is the area between the minimum and maximum curves obtained using a ± σ a and b ± σ b (with σ a,b the standard deviations provided by least squares minimization), shown to illustrate goodness of fit. The red curve is a non-linear quantile regression with the 90 th percentile, with a function containing a single parameter, c, as follows: Thus, Eq.
(2) is the curve dividing the observed sample in a way that 10% of the green dots in Fig. 3 fall above the curve. We decided to use the method of quantile regression, Eq.
(2), to assign the probability of a grid cell to represent a source area (cf. point 4a). We believe that it maximizes the information extracted from observed data, with respect to the average and the fit of Eq. 1. Moreover, we used a function with a fixed exponent (the fourth power of S/90) because using a parametric exponent as in Eq.
(1) the resulting curves systematically underestimated the probability of source presence for large values of slope. The procedure provides an individual curve P QR (S) for each topographic unit, obtained form expert mapping of potential rockfall source areas and from nonlinear quantile regression of probability values observed in the same topographic unit.
In each topographic unit, cell-by-cell values of probability obtained from Eq. (2) were directly related with the number of simulated trajectories considering the cell as the starting point of a falling boulder. We imposed a lower threshold probability under which the cell is not considered as a source. The threshold is the minimum between the value of P QR providing a source map covering at least 80% of the sources mapped by visual interpretation, and P QR = 10%. We simulated for each source cell a number of trajectories corresponding to the numerical value of P QR in percent; for example, a cell with probability X% would trigger X simulated trajectories. This allows to exploit the possibility in STONE of simulating a different number of falling boulders based on the probability of a given cell to be a source area. This method increases the degree of randomization of the simulations (together with the random selection of the value of input parameters within a given percent around the mean value provided as an input, built-in in STONE).
Eventually, we evaluated the degree of match between model and ground truth using hit rate, defined as HR = T P/(T P + F N ), where T P stands for true positives and F N for false negatives. Given that T P + F N = P , the total attempts of making a prediction (we do not consider negatives, here), the index gives a measure of the percentage of times a prediction was a hit, out of the total attempts.

Simulation of rockfall trajectories with STONE
After the selection of source areas, we performed the following steps to calculate rockfall susceptibility within a polygonal buffer around the whole Italian railway network, using the program STONE and the 10 m-resolution DEM used to identify sources of rockfalls: 1b. attribution of dynamic friction, normal and tangential restitution coefficients values, necessary for execution of the program STONE. This step was performed by expert assignment of values of coefficients available from the literature to each lithological class; 2b. preparation of data grids necessary for execution of STONE in each topographic unit, further split in rectangular domains containing contiguous SU polygons, as selected at point 2a in Section 4.2; 3b. execution of actual simulations with STONE. We performed distributed (parallel) runs within the 921 subareas described in point 2b; 4b. collection on the various raster maps containing the trajectory count per grid cell ("counter"), in each of the sub-areas of each of the 29 topographic units; 5b. transformation (classification) of the counter map into a probabilistic map.
Numerical values of dynamic friction, tangential and normal restitution are necessary for running the program STONE (cf. point 1b). We selected them as follows. We considered a lithological map at 1:100,000 scale, obtained from a digital map in vector format. The original map, from Italian Geological Service (ISPRA (Tacchia, 2004;Battaglini et al., 2012)), contained more than 5,000 unique descriptions of geological formations (Fongo, 2018;Bucci et al., 2021). Using expert criteria, with the help of about 400 geological sheets (digital photographs from the same source) at higher scale, the original unique description were classified into 19 lithological categories. The numerical values of coefficients to run the program STONE were obtained using previous estimations for similar lithologies available in the literature. (Guzzetti et al., 2002(Guzzetti et al., , 2003(Guzzetti et al., , 2004Agliardi and Crosta, 2003;Santangelo et al., 2019).
The selected values are listed in Table 1.
The actions of points 2b, 3b, 4b were performed in 921 sub-areas singled out by considering separately groups of adjacent SUs: (i) a DEM of the sub-area under investigation; (ii) a grid containing the positions of the rockfall source areas, specifying the number of trajectories to be simulated; (iii) three grids, containing the values of dynamic friction and of normal and tangential restitution coefficients, for each grid cell of the DEM.
The software STONE produces as an output three grids containing: (a) a count of the number of trajectories that crossed the cell during the simulations; (b) the maximum height of the trajectories, measured from the local elevation; (c) the maximum velocity of the simulated boulders in each cell. The three grids can be used jointly, in principle, to estimate rockfall hazard in each cell. As a matter of fact, for given source areas, (a) can be interpreted as the relative probability of having a trajectory crossing each grid cell, while (b) and (c) may be jointly considered as a proxy of the magnitude of the expected rockfalls in each cell. Estimating temporal frequency of rockfall events, instead, requires additional information (Guzzetti et al., 2003(Guzzetti et al., , 2004. In this work, we only used output (a), as an estimate of the relative spatial probability of rockfall events.
To convert the trajectory count per cell into a probabilistic index (cf. point 5b) we used the following rationale. We considered the national inventory of landslide phenomena in Italy (IFFI -ISPRA (2018) and Trigila et al. (2010)). We extracted polygons labeled in the inventory as "falls" or as "extended areas containing falls". We calculated empirical cumulative distribution functions (ECDFs) of the trajectory count, limited to the grid cells underlying the  Fig. 1(a) and Table 1). We used an approach based on ECDF to map the output of STONE, consisting in a cellby-cell trajectory count (abscissa), into a relative probability map (ordinate).
polygons extracted from IFFI inventory. In this way, we considered as relevant only the cells known to actually contain a rockfall. We used the ECDF to map all of the cells in the trajectory count grid with non-null value onto the [0,1] interval, which allows a probabilistic interpretation required for a susceptibility map.
Using an ECDF, instead of a simple normalization of the trajectory count values, has two advantages. First, it allows to attribute less relevance to cells with a small number counts, given the huge variance of the counts in the output grids. Second, using ECDFs makes results independent from the absolute values of the number of simulated trajectories per source grid cell and makes only relative values relevant, consistently with our method of assigning these numbers (cf. Section 4.2). We performed the counter-to-relative probability conversion operation separately in each of the 29 topographic units shown in Fig.  1(a). Figure 4 shows sample ECDFs obtained in four topographic units.

4.4.
Analysis of the railway network within graph theory The railway track running through Italy can be viewed as a collection of links, i.e. portions of track running without intersections from a starting point to an arrival point, and of nodes, i.e. the intersections. Nodes can be connected to multiple links. This collection of links and nodes represents a geographical network.
The links-nodes description can be further assimilated with a graph, whose vertices are the network nodes and whose edges are the network links. We performed a simple preliminary analysis of the graph corresponding to the railway network, calculating betweenness of vertices and All of the graphs contains three quantitative measures: (i) color of the vertices, displayed with a green-to-red palette, represents betweenness; (ii) color of the edges, displayed with a yellow-to-red palette, represents the edge betweenness; (iii) thickness of each edge is proportional to the actual (spatial) distance between the two vertices (nodes) it connects. Information content is the same in the three representations.
edges. This analysis describes the topology of the graph, and it is a starting point for discussing rockfall susceptibility in relation to the railway network operativity. Figure  5 is a graphical representations of such analysis. Figure 5 is a graph representation of the geographical network, where vertices are now connected by straight lines -graph edges. In Fig. 5, the dots colored with a green-to-red palette represent betweenness of vertices (cf. Section 2.2), with red corresponding to higher values. The representation of edges is twofold. The color of the edges, quantified with a yellow-to-red palette, represents betweenness of the edges: more central edges (i.e. edges crossed by a larger number of all-to-all shortest paths) have larger betweenness. Thickness of the edges, instead, is proportional to the physical distance between the nodes they connect. Edge betweenness in the (weighted) graph, is known to be one of the most important variables in networks (Barthélemy, 2011), and we used it in the subsequent analysis.
To understand the impact of possible interruptions due to rockfalls on the railway network, from graph theory point of view, we performed two separate analyses. First, we considered the impact on the graph (and, correspondingly, on the real network) of removing edges (i.e. railway links, consisting of a variable number of 1-km segments) from the network, one by one. Second, we considered jointly the classification of the railway segments into susceptibility classes and their role in the equivalent graph as described above.
Removing one edge from the graph, all of the remaining edges (as associated to the residual network) change their values of betweenness by a small amount. For each edge in the graph, we calculated total variation betweenness from the original graph to the graph obtained removing the one edge, and assigned this value as an attribute to the removed edge in the original graph. We looped over the whole set of edges, ending up with a new ranking of edges according to the newly calculated attribute. We used the newly defined quantity as an additional index for a joint classification of railway segments. The joint classification is based both on rockfall susceptibility and on the variation of betweenness index, and we named it network-ranked susceptibility.

Computing implementation
Simulations over a large area with the software STONE has the limitation that the software calculates the userdefined number of trajectories in a serial way. Despite each calculation is relatively fast, the large number of trajectories needed for this study would make it impossible to run in a single instance. We adopted a simple strategy to overcome this problem, splitting the study areas in many sub-areas (here, groups of contiguous SUs) and running several instances of STONE in parallel (cf. Section 4.3).
Raster grids were prepared for each sub-area in automatic way within GRASS GIS (Neteler and Mitasova, 2007;Casagrande et al., 2007). We used a multi-core machine, equipped with 48 computing cores and 330 GB RAM memory, to perform in parallel the GIS operations necessary to prepare the grids (cf. point 2b, Section 4.3), execution of the program itself (cf. point 3b) and collection of the 921 sub-results (cf. point 4b). The total running time for each complete run was about two days, using a masterslave strategy, which best balanced workload.

Results
This section describes the results of the physically based model STONE (Guzzetti et al., 2002) along the national railway network in Italy. We performed simulations in areas encompassed by all of the slope units (Alvioli et al., 2016(Alvioli et al., , 2020a overlapping with a 1 km buffer delineated around the railway track.
Results are described in three different steps: (i) agreement between modeled and observed rockfall source areas mapped by experts, and with sub-areas of the polygons in the IFFI inventory (cf. Section 4.2); (ii) assessment of the trajectory count output of the program STONE, performed calculating the agreement between runout areas and polygons in the IFFI inventory; classification of railway segments into rockfall susceptibility classes (cf. Section 4.3); (iii) analysis of the impact of the network classification for the railway network within the newly defined network-ranked susceptibility (cf. Section 4.4).
We believe it interesting to show and discuss separately the results of source area modeling, point (i), since they represent a relevant input to the physically based model used in this work and, probably, to similar approaches. The procedure we developed for the purpose was calibrated on a sample of expert-mapped source areas. The output of the statistical generalization is a grid with with probabilistic interpretation; comparison between modeled probability and mapped source areas used for calibration, thus, is meaningful. To this end, we calculated the number of grid cells encompassed by polygons representing mapped source areas in which the statistical procedure assigned non-null probability, and the number of cells with values of probability larger than 0.8, to check if larger probabilities are found within the mapped polygons.
Relevance of point (ii) is straightforward, as rockfall susceptibility along the railway track is the primary goal of this work. Point (iii) gives insight into the mutual relationship between susceptibility classes and the role of each railway segment within the network as a whole.

Results of the procedure to identify source areas
The procedure described in Section 4.2 applies to each of the topographic units adopted in this work as homogeneous domains for simulation with STONE. Figure 6 shows the probability-slope dependence in all of the TUs. Once source areas and the associated probabilities were estimated, we visually checked that no evident errors existed along the entire railway track. We systematically checked the results of the statistical classification of source areas along the entire railroad network, through expert visual analysis. In a few selected locations, in which expert geomorphologists considered that the infrastructure was potentially subject to manifest risk, additional source areas were delineated. Table 3 lists, for each topographic unit: the total area, the area of slope units considered here (i.e., the extent of the study area, the number and total area of the expert-mapped polygons of potential source grid cells.
For illustrative purposes, due to the large extent of the study area, we show details of expert mapping and source selection procedure, and their comparison, for one specific topographic unit, namely unit 1.

& 1.3 (Central-Eastern
Alps & Carso) of Guzzetti and Reichenbach (1994). Figure 7(a) shows an overall view of unit 1.2 & 1.3, the railway track, the polygons mapped by experts in this unit, along with the portion of IFFI rockfalls overlapping with the set of SUs in this unit. Similar settings were found across the whole study area, about 25,000 km 2 .
In Fig 7(b), the detail shows a zoom of a few slope units in the north-east of the unit, in which both expert-mapped and IFFI polygons are present (small orange box in the top figure). IFFI polygons are split into main body (purple pattern), their portion below the 10 th elevation percentile ("IFFI Lo", cyan) and the portion above 90 th elevation percentile ("IFFI Hi", green). Percentiles help showing that most IFFI rockfall polygons actually contain both the source, runout and deposition zone. We infer that from the fact that the upper percentile lies on top of slope units and approximately corresponds with expert-mapped sources, while the lower percentile lies at the bottom of slopes and stretches down to the valley, for the larger polygons. These characteristics were shared by all of the IFFI rockfall polygons that we explicitly checked visually. We use upper and lower percentiles to have a general overview of the performance of both source selection and simulation results, in the following.
The output of the generalization is a 10 m x 10 m grid map aligned with the TINITALY DEM, whose value in each cell is the probability for a rockfall trajectory to originate from that specific location. An example probabilistic source map selected within topographic unit 1.2 & 1.3 is in Fig. 8(a) A comparison between modeled probability and expertmapped source areas represents an evaluation of the per- Figure 6: Each panel is as in Fig. 3; topographic unit 1.2 & 1.3 is shown in Fig. 3; 2.3a and 2.3b were grouped into one section, for what concerns modeling of source areas; section 3.1 does not contain mapped sources; sections 6.3 and 8.2 contain no intersections between railway and slope units; this figure shows results for the remaining 24 sections.
formance of the calibration step of the source identification procedure. Performance is quantified by hit rate of modeled probability versus ground truth, represented by mapped polygons. We considered the modeled probability in two different ways. We either selected all of the grid cells for which the statistical procedure assigned any nonnull probability, and cells with values of probability larger than 0.8, as an arbitrary threshold to distinguish values of "high" probability.
We devised a second comparison, using rockfall polygons in the IFFI inventory (Trigila et al., 2010), intended as a validation step of the procedure for source identification. We calculated, within each polygon in the inventory, the portion corresponding to the upper 90% in elevation. This was our best (arbitrary) guess to what could have been the source area of each catalogued rockfall. Results for the comparisons in each of the topographic units, in terms of HR, are listed in Table 3. We did not evaluate true negative rate, the counterpart of hit rate, or a full confusion matrix: true negatives are unknown, here, because expert mapping was performed only in selected location, and we did not know if IFFI is a complete inventory. Hit rate was calculated within the slope units selected as study area for this work. Table 3 shows HR corresponding to four combinations of HR calculated as: (i) modeled sources, with any probability value, against the expert-mapped source areas ("Mapped Total"); (ii) modeled sources, limited to cells with highest probability, against the expert-mapped source areas ("Mapped P > 80%"); (iii) modeled sources, with any probability value, against the upper 90% elevation of IFFI polygons ("IFFI Total"); (iv) modeled sources, with highest probability, against the upper 90% elevation of IFFI polygons ("IFFI P > 80%").  .1) are below 25%. The corresponding sections have small or very small total SU area, and many are located in the vicinities of plains or coasts, which are probably the most difficult to model on a regional or national level as they differ the most from the typical settings in which we mapped sources in expert way.
Combination (ii): the validity of this comparison is prob- Table 3: Numerical evaluation of the statistical generalization for the probability of grid cells to initiate a rockfall trajectory. We list the total area of each section (Section ID, as in Guzzetti and Reichenbach (1994)), the area covered by slope units selected for this work Alvioli et al. (2020a), the number of expert-mapped source polygons and their total area. HR is hit rate (or true positive rate), the ratio T P /P = T P /(T P +F N ). "Mapped" refers to the cells underlying the expert-mapped source areas (colored in green-yellow-orange-deep orange, from high to small HR), and "IFFI" (colored in light to deep blue from large to small HR) refers to cells in the 90 th elevation percentile of the rockfall polygons in the national inventory (Trigila et al., 2010). ably more difficult to understand on general grounds. Out of the 26 sections (three out of the original 29 were ruled out by absence of railway track or absence of slopes prone to rockfalls), four have HR larger than 20%, 11 sections show smaller than 10%, five of them smaller than 5%.
Combination (iii): HR results of the comparison between modeled sources and polygons labeled as "rockfall" or "extended areas containing rockfalls" in the IFFI national inventory. In this case, compared with the region above the 90 th elevation percentile, within each polygon. Results reveal a lower degree of match with respect to case (i): 20 sections assigned with HR less than 75%, of which 13 smaller than 50%, of which one smaller than 25%.
Combination (iv): as in case (ii), interpretation of the comparison is less obvious than the one with any value of probability for mapped sources, case (iii). Out of the 26 relevant sections, one has HR larger than 20%, 13 sections show HR smaller than 10%, five of them smaller than 5%.

Results for rockfall susceptibility within the model STONE
The final product of this work is a classification of the railway track in Italy, split into 1-km segments, with susceptibility values for the occurrence of rockfalls. To this end, the most relevant points are as follows (cf. enumerated list in Section 4.3). We prepared maps of trajectory count per cell using the model STONE. Then, we classified all of the cells encompassed by the selected SUs (cf. step 2a of the enumerated list in Section 4.2); we used the map to calculate a unique value of susceptibility per 1-km segment of the track.
We converted the trajectory count map, containing wildly varying values, into the [0,1] interval -consistent with a (relative) probability, i.e. susceptibility. We did so using ECDF functions calibrated in the subset of grid cells in which rockfalls are known to have occurred, reported in the national IFFI inventory (Trigila et al., 2010). The Table 4: Hit rate for the validation of rockfall runout areas, i.e. obtained from the comparison of probability of occurrence P QR from reclassification of trajectory count map by means of ECDFs as described in Section 4. The ID in the first column corresponds to topographic units of Fig. 1(a), Table 1 and Table 3. "Total" HR values correspond to values of any non-null probability value, and partial HR values to probability between the specified intervals. The largest value among the partial ones is in bold. This comparison involves only the 10 th elevation percentile of each polygon in the IFFI inventory.
Section ID HR HR HR HR HR inventory contains polygons mapped as "rockfalls" or as "extended areas containing rockfalls". We calculated specific ECDFs in each topographic unit shown in Fig. 1(a). Figure 4 shows a few sample ECDFs corresponding to the first four topographic units in Tables 3, 4 and 5. The remaining ECDFs are very similar, and curves would overlap with those in Fig. 4. Figure 8(b) shows the final, classified trajectory count (susceptibility) in a sample location. The raster map is colored with a green-to-red color ramp from low to high values of susceptibility. In the figure, we compare the susceptibility map with the IFFI polygons existing in the area (cf. Fig. 7). We performed validation of the classified output in the susceptibility map as follows. We selected all of the rockfall polygons from the IFFI inventory within the slope units buffer, and further selected the DEM grid cells under the polygons and falling in the lowest 10 th elevation percentile -also shown in Fig. 8. We are reasonably confident that these cells represent rockfall deposit areas.
Comparison of probability maps among different topo-graphic units with deposit areas inferred from IFFI polygons, and with the whole IFFI polygons, provided values of hit rate (HR = T P/(T P + F N )) listed in Table 4 and labeled as "Total" HR. Similarly, we defined "partial" HR values, limiting the calculations within the individual classes delimited by (0, 0.2, 0.4, 0.6, 0.8, 1) values, whose results are also listed in Table 4 in five additional columns. Table 5 has a similar structure as Table 4, but it lists results corresponding to comparison of modeled source areas with the whole extent of rockfall polygons in the IFFI inventory. Topographic unit 2.3b does not contain rockfalls recorded in the IFFI inventory relevant to the railway track, while sections 6.3 and 8.2 do not contain intersections of the track with slope units. Thus, we reported no values for such zones in Tables 4 and 5.
The next step was the classification of the railway network split into 1-km segments and classified using a simple strategy. For each segment, overlapping many susceptibility values, we selected the largest value. This is a conservative and reasonable choice, given that if a single  19% 19% point of the segment would be hit by a rockfall, the whole segment would be unusable. On the other hand, splitting the whole link connecting two nodes into small segments of the same length guarantees spatial homogeneity. Unfortunately, available data was insufficient to perform a validation of the result of the classification of segments.
To perform an additional assessment of the results, we performed graph analysis, and evaluated its implications for the operativity of the railway network.

Impact of rockfalls on the railway network: networkranked susceptibility
In Section 4.4 we described an analysis of the railway network in terms of graph theory. We defined a new index for each railway segment, obtained as the variation of total edge betweenness on the network if the edge corresponding to the segment is removed from the graph. Figure 9(b) shows the result of the analysis; edges are colored with a blue color palette, for increasing value of the attribute calculated by removing one edge at a time. The figure clearly highlights the most relevant edges, topology-wise, and one can easily spot the difference between this analysis and edge betweenness, Fig. 5.
The result of Fig. 9(b) is an intermediate step towards the definition of a combined quantity, defined considering simultaneously the variation of betweenness, thus strictly related to network properties, and classification of the railroad based on rockfall susceptibility. The combination was performed in the simplest possible way, by crossing the classes of both quantities to obtain a new classification, the last result of this work.
Classes based on rockfall susceptibility of railway links, split in 1-km segments, was defined in Section 4.3, and Fig. 9(a) shows the results. Variation of betweenness on removal of one graph edge was defined in this section, and it is shown in Fig. 9(b). The latter quantity was categorized into five categories, using natural breaks.
We refer to the joint classification of susceptibility and segment relevance as network-ranked susceptibility. Results for network-ranked susceptibility, expressed in terms of kilometers or railway in each class, are in Table 6. The ta- Figure 9: (a) Classification of the railway network, split into segments of 1 km length, with the procedure developed in this work. Segments are colored with a red-togreen palette, corresponding to low-to-high susceptibility for rockfalls. (b) A measure of relevance of railway links in the Italian railway network, alternative to betweenness, calculated in Section 4.4. Relevance is shown with a blue color palette; blue denotes more relevance. Relevance is high if total variation of vertex betweenness is high, when the corresponding edge is removed from the graph (cf. Fig.  5). (c) Network-ranked susceptibility, defined as a combination of segment-wise rockfall susceptibility, modeled in this work with STONE and shown in (a), and relevance of the corresponding railway links, show in (b). Table 6 lists a breakdown of the total railway length, in kilometers, within the seven classes shown in this figure in different colors.
ble clearly shows the classification obtained by crossing five classes in the two considered attributes, i.e. segment relevance (score) and rockfall susceptibility. The 25 classes are further colored (i.e., further classified) using seven colors, in order of increasing ranked susceptibility: light brown, yellow, red, magenta, blue, cyan, green. We defined the coloring scheme in a totally arbitrary way, except that: (i) the lowest rockfall susceptibility class, corresponding to the largest portion of railroad, was split into two classes of network-ranked susceptibility, and (ii) the remaining portion of the railroad was (almost) evenly distributed into additional five classes.
We applied the same coloring scheme of Table 6 to the graphical representation of the classification of segments, Fig. 9(c). Comparison with the results of classification based on the only susceptibility highlights manifest differences. In fact, while reddish segments are consistently located in hilly and mountainous areas, plain areas now con- tains the "very low" and "low" classes of network-ranked susceptibility, while they just had the lowest class of susceptibility; this was a totally arbitrary choice. We discuss the possibility of optimizing the new network-ranked susceptibility in the next section.

Discussion
As for the Methods and Results sections, we discuss separately the source areas identification procedure, rockfall susceptibility assessment and implications from network analysis combined with susceptibility.
We stress here that we used the national landslide inventory IFFI (Trigila et al., 2010;ISPRA, 2018) for validation in several ways. Nevertheless, we did not actually use the inventory for building or calibrating the rockfall model adopted here. In fact, we calculated hit rates of the statistical generalization of expert-mapped source areas and hit rates of runout areas against proper subsets and sub-areas of the polygons in the IFFI inventory.

Identification of rockfall source areas
An illustration of the comparison between modeled source probability, expert-mapping sources and the upper elevation percentile of IFFI polygons is in Fig. 8, for one particular region in topographic unit 1.2 & 1.3. We can see that in this region non-null pixels of the probabilistic source map have a wider extent than both expert-mapping and IFFI upper elevation percentile. This is by construction, because the statistical procedure described in Section 4.2 was devised to conservatively assign non-null probabilities to any grid cell, above a rather low threshold, based on their slope. This is consistent with the observation that expert-mapped sources actually contains low values of slopes (cf. Figs. 3 and 6). Figure 8(a) shows a nice agreement between dark-shaded probabilistic map and both IFFI upper elevation percentile and expertmapped area, though not all areas match perfectly.
Nevertheless, values of model match against the mapped and IFFI polygons, in Table 3 are often rather low. We attribute that to the following reasons. To build a model of source areas, we hypothesized a relationship between the probability of a grid cell of initiating a rockfall and local slope. This represents a compromise between an acceptable overall time needed for the procedure over a large study area and a realistic product, but it certainly does not embed all of the local terrain properties that influence the expert criteria applied for mapping potential source areas. Moreover, we made specific choices and assumptions, often arbitrary, in the modeling process.
Errors may also arise from the discrepancy between the DEM used in the analysis and the apparent resolution of Google Earth™ imagery used for expert mapping, especially in the locations with steepest relief, the ones in which we are mostly interested in. In the analysis, we used a DEM generated from a triangulated irregular network (TIN); visual inspection of a shaded relief generated form the DEM highlighted locations in which the triangulation used to prepare the DEM is manifest, which surely affects the slope map, nation-wide.

The model STONE and rockfall susceptibility
The output of the model STONE is a raster map: the classified trajectory count per grid cell. We used it in several ways, to get to the final results of this work: (i) we compared the classified output with the relevant subset of the IFFI inventory; (ii) we split the railway network into 1-km segments and classified them on the basis of rockfall susceptibility; (iii) we studied the properties of the railway network within graph theory, considering the network in relation with the susceptibility results, which helped defining a joint classification.
Tables 4 and 5 list results of the validation against IFFI polygons. Values for the total hit rate in this case are rather large, which means that most of the runout areas in the inventory were actually overlapping with at least one trajectory. Considering the "partial" HR in update of the methods presented in this work will be implemented. This could most probably be in the procedure for the selection of source areas (for example as suggested at the end of Section 6.2, or as in Rossi et al. (2021)), or in the strategy for classification of the trajectory count. Table 5 correspond to validation against the whole extent of polygons in the IFFI inventory, at variance with Table 4. Values of "Total" hit rate are rather similar to the case where we compared to the 10 th elevation percentile of each polygon, in a few cases a slightly higher. Interestingly, in this case, values of hit rates are almost always equally distributed among the five sub-classes of probability. One reason for that could be due to spatial inaccuracies of the polygons in the inventory.

Results in
The trajectory count obtained here, Fig. 8(b), is a raster map with resolution of 10 m. Nevertheless, we consider as a final product of this work the vector map of the railway network, split into 1-km segments, classified with five levels of susceptibility. The reason for that is twofold. First, the original aim of this work was to prepare a map which could be used for ranking, or prioritizing, safety countermeasures at the national level. Second, we acknowledge that we necessarily made assumptions and approximations to be able to work homogeneously at the national scale; for this reason, a raster map with 10 m resolution could contain local inaccuracies, which are mitigated by the classification of 1-km segments. This implies that application of the classified map is intended to locate the segments with highest susceptibility at the national scale; conclusions at smaller (local) scales can still be obtained within the model STONE, but should be supported by higher resolution data.

Network-ranked rockfall susceptibility
Network-ranked susceptibility for rockfalls, shown in Table 6 and Fig. 9(c), was defined in this work as a straightforward combination of five classes in edge ranking and rockfall susceptibility. We already noted that the proposed classes produce somewhat striking results, Fig.  9(c), in that segments with negligible rockfall susceptibility (cf . Fig 9(a)) do not fall in the lowest class of network-ranked susceptibility. We stress that the seven classes (colors) shown in Table 6 were chosen almost arbitrarily -hence the apparent mismatch.
We consider the classification proposed here as a methodological proposal, and we did not attempt a higher-level optimization of final classes, which we could actually perform with suitable validation data. Validation of susceptibility requires accurate data about observed rockfall occurrences, possibly including information about whether rockfalls did or did not cross the railway track. Calibration and validation of network-ranked susceptibility, in addition, would require information about the actual operation of the network.
Knowledge of actual traffic on the railroad would imply different network properties. In fact, betweenness of nodes and edges was calculated here considering all-to-all shortest paths. If we replaced them with actual routes we would obtain a different network and corresponding graph (Kurant and Thiran, 2006). As another example, knowledge of the number of trains running on each railroad link per day could be used to define relevance of edges in a different way, beyond the simple topology of the network considered here. Additional information on the average number of passengers per link/day could help in a similar way. On the other hand, knowledge of real time presence of trains on each link of the network, though not difficult to consider in technical terms, is not immediately relevant for this work. That would require inclusion of some dynamical mechanism for the actual detachment of rocks uphill form each railway link, along with real time monitoring of possible triggering events.

Conclusions
We described a comprehensive modeling chain for the assessment of rockfall susceptibility along the railway network in Italy. To cope with issues related with the large size of the study area, extensive mapping and computational challenges, we necessarily adopted assumptions and approximations. These were discussed in details throughout this work, in which we described the performance of the different steps of the procedure and we highlighted weak points and possible improvements. An additional point to stress here is that, due to the large size of the study area and unavailability of specific data, rockfalls on anthropic slopes were not considered. On the other hand, we highlighted the merits of this work as an effective procedure for rockfall susceptibility assessment over a very large and geomorphologically diverse study area.
We described three main steps: (i) localization of rockfall source areas on a digital landscape, with different probabilities, (ii) simulation of rockfall trajectories stemming from any possible source area, and (iii) assessment of the impact on the operativity on the railway network, using a joint classification of network segments taking into account both rockfall susceptibility and network properties. The relevant findings were as follows: • The procedure for localizing source areas is a statistical generalization of the distribution of slope angle values within expert-mapped polygons. In general, values of hit rate were seldom in the higher quartile, and mostly in the third and second quartile. Results for the validation were similar. The procedure is suitable for identification of source areas in large areas (here, about 25,000 km 2 ), though improvement of calibration performance would require including additional morphometric variables besides slope.
• Rockfall trajectories were simulated using the program STONE. Validation of the classified runout in terms of hit rate, were rather satisfactory, in that the majority of values were in the higher quartile of the probability distribution. We further cast the susceptibility map onto the railway network, by splitting the railway link in 1-km segments and assigning a class of susceptibility to each segment. Validation of the susceptibility values for segments of the railway requires additional data, not available to us.
• The analysis of possible rockfall impact on the railway considered the effect of disruption of a railway segment on the entire network. That was quantified by a new quantity which we called network-ranked susceptibility, considering the joint classification of rockfall susceptibility and graph properties. We consider this step as a methodological proposal, whose optimization and validation would require dedicated data and additional knowledge of the actual traffic over the railway network.
We stress that, in this work, the national landslide inventory IFFI (Trigila et al., 2010) helped validating results, but did not enter neither in the runout modeling, nor in the identification of sources areas, which required expert source areas mapping. This step of the procedure leaves room for improvement, in that we only used distributions of slope angle values to characterize expert-mapped polygons.
Both steps (identification of source areas and simulation of runout) can, in principle, be extended outside of the buffer considered in this work, because most of the information and methods developed here are still valid. The buffer itself, consisting of slope units overlapping with the railway, can be easily enlarged, given that slope units are available for the whole of Italy (Alvioli et al., 2020a).
Eventually, we stress that a natural evolution of the analysis performed on the topological properties of the railway network is a real-time monitoring (still static, no traffic involved) of the network as a function of the real-time evolution of possible triggering events. This was actually the original motivation of this work. The segment-wise susceptibility map developed here, along with analogous map obtained for debris flows susceptibility and soil slides susceptibility, are key inputs for the prototype early warning system that was funded by the National railway company (RFI) within the framework for the implementation of the SANF (an acronym for national early warning system for rainfall-induced landslides) system (Rossi et al., 2012;Guzzetti et al., 2020).