Three-dimensional simulations of rockfalls in Ischia, Southern Italy, and preliminary susceptibility zonation

Abstract Ischia Island is a volcano-tectonic horst in the Phlegrean Volcanic District, Italy. We investigated rockfalls in Ischia using STONE, a three-dimensional model for simulating trajectories for given detachment locations of blocks. We propose methodological advances regarding the use of high-resolution LiDAR elevation data, the localization of possible detachments sources, and the inclusion of scenario-based seismic shaking as a trigger for rockfalls. We demonstrated that raw LiDAR data are useful to distinguish areas covered by tall vegetation, allowing realistic simulation of trajectories. We found that the areas most susceptibile to rockfalls are located along the N, N-W and S-W steep flanks of Mt. Epomeo, the S and S-W coast, and the sides of some steep exposed hydrographic channels located in the southern sector of the island. A novel procedure for dynamic activation of sources depending on ground shaking, in the event of an earthquake, helped inferring a seismically-triggered source map and the corresponding rockfall trajectories, for a scenario with 475 y return time. Thus, we obtained preliminary rockfall suceptibility in Ischia both in a “static” (trigger-independent) scenario, and in a seismic shaking triggering scenario. They must not be considered a risk map, but a starting point for a detailed field analysis.


Introduction
Landslides are natural phenomena which shapes the Earth surface, and in populated areas they can endanger lives and cause huge economic damages. Rockfalls are among the most dangerous types of landslides for their rapidity and destructive potential. Research on landslides focuses on methods that allows spatial zonations of the probability of occurrence (susceptibility), while hazard assessment requires knowledge of magnitude and temporal frequency. Additional information about vulnerable elements exposed to hazard allows to estimate the risk associate to landslides, in general, and rockfalls in particular.
For example, Copons and Vilaplana (2008) discussed the exposure of human lives and land to rockfalls in Andorra; Youssef et al. (2015) investigated the relationship between rockfalls and urban expansion in Saudi Arabia; Alvioli et al. (2021) produced rockfall susceptibility maps along the whole Italian railway network; Dorren et al. (2022) discussed mitigation of rockfall risk in the Swiss Alps.
Here, we performed a preliminary susceptibility zonation of rockfall susceptibility in the Ischia Island, with a three-dimensional model and high-resolution topographic data. We further implemented a novel mechanism to link spatial susceptibility with a seismic forcing with well-defined return time, adding a temporal component to get one step closer to full hazard assessment; we further analyzed the possible impact of rockfalls on roads and built up areas, in a critical way.
Rockfalls are widespread geomorphological processes that represent one of the main hazards in mountain areas (Whalley 1984), with mobilized volumes ranging from less than 1 to 10 5 m 3 (Evans and Hungr 1993;Hungr et al. 1999;Ruiz-Carulla and Corominas 2020) and with a very rapid to extremely rapid movement (Broili 1973;Dorren 2003). To describe them, we adopted the model STONE (Guzzetti et al. 2002), suitable for events including the detachment of individual blocks from a steep slope or a rocky cliff and following a path consisting in falling, bouncing and/or rolling on the ground. Thus, we do not discuss rock avalanches and other phenomena, here.
The most relevant input of STONE is a map of possible rockfall sources. This is a common input to many similar models, namely: RockGIS, (Matas et al. 2017); Rockyfor3D, (Dorren et al. 2022) and references therein; 2D CRSP, (Jones et al. 2000); Hy-STONE (Agliardi and Crosta 2003). All of these models calculate geometrical rockfall trajectories for given starting points; and none of them contains an physical mechanism to determine if and when blocks actually detach from a specific point. We applied STONE because we have been working with the model for several years already, in different settings (Guzzetti et al. 2003(Guzzetti et al. , 2004Santangelo et al. 2019;Sarro et al. 2020;Santangelo et al. 2021;Alvioli et al. 2021).
Several methods exist to gather information about slope stability and to map landslides features, including rockfalls sources and trajectories followed by individual falling blocks in individual slopes. They may involve direct field observation (Heckmann et al. 2016;Rossi et al. 2021), laser scanning (Copons and Vilaplana 2008), photogrammetry based on photographs acquired on the ground (Matas et al. 2022) or by aerial vehicles (Santangelo et al. 2019;Buyer et al. 2020;Giordan et al. 2020), and infrared thermography (Baron et al. 2014;Loche et al. 2022). Such information is useful to infer the likely location of future rockfall sources, and represent a valuable input for the application of three-dimensional models for rockfall trajectories. Simulation of boulders trajectories represents an essential tool in hazard/risk analyses (Lan et al. 2007(Lan et al. , 2010Pellicani et al. 2016).
Here, we present results of an analysis based on a new application of LiDAR elevation data and use of the model STONE (Guzzetti et al. 2002) for a probabilistic assessment of rockfall trajectories. We investigated various steps of the procedure involved in preparing input data for running the computer code, assessment of the results of simulations and their impact on built-up areas and roads, and a scenariobased hazard assessment for seismic triggering of rockfall events to identify the most rockfall-susceptible areas of the island, to plan detailed field studies aimed at full rockfall risk assessment in Ischia.
This work is organized as follows: Section 2 describes the geological setting, the historical seismicity and the historical flank instability of Ischia, and introduces the three-dimensional rockfall simulations. Section 3 describes and lists the data used for this study, specifying which were already available and which were newly developed. Section 4 describes the different steps performed in this work: interpolation of raw LiDAR data; values of input parameters, also in relation to vegetation; rockfall source localization; application of a new probabilistic model for triggering seismically induced rockfall. Section 5 describes the main results for the assessment of trajectories, Section 6 discusses results and draws conclusions about the simulation of rockfalls in Ischia island.

Geographical and geological setting
The Ischia island is one of the three volcanic complexes located in the Phlegrean Volcanic District (Orsi et al. 1999) which last erupted in 1302 (Civetta et al. 1991). It is an active volcanic field at rest that covers an area of about 46 km 2 and it is characterized by a remarkable ground uplift caused by the resurgence of an ancient collapsed caldera (Acocella and Funiciello 1999), whose rim is facing the southwest and northwest sectors of the island. (a) Location of the study area; the background represents seismic hazard levels in Italy, represented by maximum ground acceleration with exceeding probability of 10% in 50 years (Stucchi et al. 2004), for illustrative purposes. (b) A shaded relief obtained from 2 Â 2 m elevation data of the island, colorized with shades of green as a function of the difference between a digital surface model (including vegetation) and a digital terrain model (excluding vegetation). Brown polygons are buildings and black lines are roads, both from the OpenStreetMap project.
Being a volcanic island, Ischia shows different processes, e.g. fumarolic activity, earthquakes, slope instabilities and volcanic climax eruptions. Volcanic edifices experience slope instability as consequence of different solicitations such as (i) eruption mechanism and depositional process, (ii) tectonic stress, (iii) extreme weather conditions and (iv) seismic shaking. All these events induce the mobilization of unstable fractured volcanic flanks, initiating rockfalls (Carracedo 1996;Hurlimann et al. 2000;Delcamp et al. 2017;Roberti et al. 2021). The island is also densely populated, and the national inventory of landslide phenomena (Trigila et al. 2010;ISPRA 2018) in the area show rockfall runout/deposition areas with substantial overlap with buildings and infrastructure ( Figure 1). Thus, assessment of possible locations for initiations of new falling blocks, and their possible trajectories downhill, is of paramount importance.
The regional extensional tectonics is supposed to be related to the opening of back-arc basins caused by an East-retreating subduction of the Apulo-Adriatic lithosphere (Doglioni et al. 1996). The resurgent dome generated an uplifted block, the Mount Epomeo horst, located in the central and western sector of the island; it is marked by a system of subvertical faults striking NW-SE and NE-SW, on the edges of the dome and away from it, and and N-S and W-E mainly located at the borders of the dome (Acocella and Funiciello 1999). The southeastern part of the island is characterized by highly dipping ENE-WSW normal faults. Volcanism at Ischia started over 150 ka B.P. (Cassignol and Gillot 1982) and continued, with very long period (centuries to millennia) of quiescence, until the last eruption occurred in 1302 A.D. (de Vita et al. 2006). The oldest exposed rocks belong to a partially eroded volcanic complex, which crops out in the south-eastern part of the island, covered by more recent deposits, composed of volcanic effusive and explosive rocks, mostly trachytes and phonolites (de Vita et al. 2006). The large Mount Epomeo Green Tuff caldera formed during the eruption that took place 55 ka ago. In addition to volcanic soils, debris flow deposits cover the southern slope of Mount Epomeo and in the northern and western sectors of the island (de Vita et al. 2006). Southern sector collapses of Ischia, probably occurred between between 8.6 and 5.7 ka ago as a consequence of the resurgence, generated three major debris flow deposit units (Tibaldi and Vezzoli 2004).

Historical and recent seismicity
Historical seismicity of Ischia is dominated by the earthquakes of 1881 (D'Auria et al. 2018) and 1883 (Carlino et al. 2010), which mainly affected Casamicciola Terme and the other municipalities of the island. Both events, especially the second one, strongly affected the towns on the island causing the death of 2,313 people (Cubellis et al. 2004), and triggered several ground and environmental effects: large landslides North-Western slopes of Mt. Epomeo, ground cracks and variations of the springs' flow rate, chemistry and temperature (Alessio et al. 1996). Since 1883, before 2017 the seismicity of Ischia was characterized by small and shallow events, most of which were detectable only in Casamicciola Terme (D'Auria et al. 2018) On August 2017 a seismic sequence characterized by a starting Mw 4.0 earthquake occurred 1-2 km below the town of Casamicciola Terme (De Novellis et al. 2018). Given that the highest geothermal activity is located in the southwest sector of the island while the main seismogenic portion of Ischia is represented by the Casamicciola Terme area, the recorded seismicity does not seem to be correlated with the geothermal activity (Chiodini et al. 2004). Then it is reasonable to speculate that seismogenesis at Ischia is probably related to the structural dynamics of the northern part of the island Paoletti et al. 2013).

Slope stability
The North, North-West and South-West facing steep flanks of Mount Epomeo are considered the gravitationally least stable slopes of the island (Della Seta et al. 2011). Mass movement deposits are widespread on the island and were generated by rockfalls, slides, toppling, debris flows and debris avalanches (Mele and del Prete 1998;Tibaldi and Vezzoli 2004;del Prete and Mele 2006). Slope movements on the island are correlated with vertical uplift and volcanism (Fusi et al. 1990;Alessio et al. 1996;Mele and del Prete 1998). This relationship is documented by displacement of volcanic and non-volcanic deposits, whose age and original stratigraphic position are in some cases well constrained (de Vita et al. 2006).
Old mass movements have been correlated to volcanic eruptions (Buchner 1986), whilst recent events (del Prete and Mele 2006) have been related to exceptionally heavy rains (Mele and del Prete 1998) or to seismicity (Guadagno and Mele 1992;Molin et al. 2003;Rapolla et al. 2010). Rockfalls mainly occur as an effect of the gravitational instability on the cliffs along the coastline or on the major scarps generated by the rapid resurgence. Their activation was favored by the widespread rock fracturing (favored in turn by hydrothermal weathering), which allowed even huge blocks to be detached (Della Seta et al. 2011).
Rockfall dynamics can be defined as a two steps process: the detachment of a rocky block from a slope (first step) and the subsequent downslope movement (second step). In general terms, landslide hazard analysis should determine both the likelihood of occurrence of failures, both in time and space. Due to its characteristic and to its rapidity, a rockfall hazard assessment should also consider the travel distance of the falling blocks, and their trajectories (Guzzetti et al. 2004).
More specifically, a comprehensive analysis of rockfalls hazard should include a wealth of specific data: accurate source areas locations in a rock cliff (Loye et al. 2009), the volume of rock mass that can be released (Mavrouli and Corominas 2020;Francioni et al. 2020;Hantz et al. 2021), the dissipation of energy during rolling and rebounds on the slope surface ) and impacts on trees (Dorren et al. 2006;Lundstrom et al. 2009;Lu et al. 2021), penetration depth in the soil during impact (Pichler et al. 2005;Lu et al. 2019) and, possibly, a model for fragmentation of the block after each impact (Ruiz-Carulla and Corominas 2020).
In this research we did not adopt a physical model to infer the timing of the rockfall and the detachment processes. The existing methods to analyze the source area, investigating the precise location of the possible detachment points, the travel distance of the falling blocks, the expected volume of the blocks, and the frequency and size of rockfalls, are only applicable on very small areas and are beyond the scope of this work. Expert mapping of possible source areas by photo and satellite interpretation is an adequately strong and effective approach which can be used to identify the source points on relatively large areas (Santangelo et al. 2019, in combination with morphometric analysis (Michoud et al. 2012;Alvioli et al. 2021).
Here, we are interested in obtaining a preliminary rockfall zonation of the island. We also estimated the likelihood of built-up locations and roads to be hit by a block falling from upslope, according to the proposed zonation; risk assessment is not the primary focus of this work, though. To evaluate the rockfall trajectories, the following information is needed: (i) the point of origin of the falling rock, (ii) the underlying topography, (iii) the friction and energy restitution coefficients of the cropping out rocks downhill, useful to run a model for the description of how the mass would reach a rest state.
Eventually, we performed an assessment of seismic-induced rockfalls within an earthquake scenario with a specific return time. We proposed a new model to link a scenario-based seismic ground shaking to a probability of triggering rockfalls. For this additional step we need additional input, with respect to the data listed above. The seismic-triggering mechanism is based on spatially distributed peak ground acceleration (PGA), and it will be explained in detail in the following.

Materials
We used materials which were either already available to us (a small scale geo-lithological map; an inventory of rockfalls occurred on the island; a vector map of the buildings and roads, and a map of ground shaking for a specific scenario), or newly developed maps (a digital elevation model, and a map of sources prepared in expert way). Here follows a short description of such data.
1. Raw LiDAR data provided by the Italian Ministry of the Environment.
Interpolated data to obtain digital terrain and surface models, to run the computer code of the model STONE and to obtain spatial information about vegetation height; Figure 1). 2. Vector layer containing the street network and the polygons of buildings, from OpenStreetMap; Figure 1). The map contains about 580 km of roads, and more than 16,000 buildings. 3. Inventory of landslide phenomena in Ischia, selected from the inventory of landslide phenomena in Italy known as IFFI (Trigila et al. 2010; ISPRA 2018). The  IFFI inventory contains 308 landslide polygons in Ischia island, of which 100 are identified as rockfalls and diffused rockfalls ( Figure 2). 4. Map of sample potential source area for rockfalls on Ischia, from expert interpretation of Google Earth TM images ( Figure 2). 5. Geo-lithological map of Ischia, scale 1:50,000 Rapolla et al. (2010). Figure 3 shows the map, with eight lithological classes (listed in Table 1). 6. Map of PGA for the major Casamicciola earthquake of 1833, obtained from Rapolla et al. (2010), and corresponding to a return time of 475 years ( Figure 4). We are aware of a recent work (Albano et al. 2018) proposing a model for a PGA map of the 2017 earthquake in Ischia, also in relationship with ground movements. We did not consider that specific event, here, because we deal with scenarios associated with specific return times; we will discuss this choice further in the Methods section. 7. A land cover map, obtained from ISPRA (2012), useful to compare inferred vegetation areas with actual land cover, though at lower resolution (10 m) than the DEM adopted here.

Methods
The different steps performed in this work can be summarized as: (1) interpolation of raw LiDAR data, to obtain digital elevation models; (2) assignment of values for input parameters, based on geo-lithological data and presence of vegetation; (3) localization of prospective sources for the initiation of rockfall trajectories and run of the computer program STONE, for different options of the input source map; (4) application of a new model for selective triggering of prospective sources, and corresponding trajectories, due to seismic shaking. We illustrate these point in the following four paragraphs.

Interpolation of raw LiDAR data
We used data from the Italian "Geoportale Nazionale," 1 managed by the Ministry of Environment and providing different kinds of spatial data. In particular, the archive contains an extensive LiDAR survey covering a substantial portion of Italy, with data  Figure 3 shows the corresponding geo-lithological map. We set STONE to perform random sampling of values of the parameters in a ± 10% range around the nominal values listed here.
stored at intermediate processing level. Each point contains the x, y coordinates of the point, its elevation, and a binary classification flag. For this research, we selected point clouds covering the Ischia island and we interpreted the two values of the flag either as referring to reflections of the laser impulses from the ground (first returns labeled with flag "1") or from the upper vegetation, or other elevated surfaces (second returns, flagged with "2"). The two sets of data contain, respectively, 39,896,156 and 25,320,537 points. We interpolated separately the two point clouds, using the module r.surf.rst specifically designed to perform surface interpolation from vector points map by splines, within GRASS GIS (Mit a sov a and Hofierka 1993). Interpolation resulted in two digital elevation models, with a square grid of 2 m Â 2 m -the nominal resolution of the digital elevation model provided "as-is" from the same data source as the point clouds. In the following, we will refer to the first one as digital terrain model (DTM) and to the second one as digital surface model (DSM).
Moreover, we interpreted the point-to-point difference between DSM and DTM as due to vegetation, and exploited this information to infer modifications of ground parameters relevant for the simulations with STONE. We partially took into account disturbances due to the presence of anthropic structures and buildings using additional land cover data, which we correlated with point-to-point DSM -DTM differences, as explained in the next section.

Numerical input parameters
In STONE, the behavior of individual blocks rolling and bouncing down a slope on different types of substrata is controlled by three main parameters. Friction controls the rolling stage, while normal and tangential restitution control bounces. We inferred the values of the three parameters by expert comparison of the geo-lithological formations existing in the island and values of the parameters reported in the bibliography, mainly by (Guzzetti et al. 2002) and Alvioli et al. (2021). The selected values are listed in Table 1.
The criteria adopted to match lithotypes with parameter values were as follows (Guzzetti et al. 2004). Massive and/or thickly bedded rocks (e.g. lavas, limestone, thickly bedded siltstone, claystone etc.) are characterized by very high values of the normal and tangential energy restitution coefficients and very small values of the dynamic rolling angle, as well as pyroclastic deposits and breccias and pumice; whereas thinly bedded limestone, clay and reworked materials are characterized by intermediate values of the normal and tangential restitution coefficients and also of the dynamic rolling angle. Empirical field data revealed that the alluvial deposits show low values of the normal and tangential energy restitution coefficients, and a very high value of the dynamic rolling friction angle (Guzzetti et al. 2002(Guzzetti et al. , 2004. The program STONE performs a random modification of individual values of the parameters, within configurable intervals and probability distribution functions. Moreover, we assumed that vegetation interferes with trajectories of falling blocks reducing their ability to proceed downhill (Dorren et al. 2006;Masuya et al. 2009;Dorren et al. 2022). Thus, we modified the parameters exploiting information about the presence of vegetation obtained using the available LIDAR data (Section 4.1). The map of vegetated area used in Figure 1) to colorize the shaded relief of the island was obtained subtracting the DEM generated using the" second impulse" from that obtained interpolating the" first impulse," and considering only points with a positive altitude difference of more than 5 m, which actually exists in the island.
We assumed that vegetation interferes with trajectories of falling blocks only modifying the parameters of cells identified as containing vegetation; we did not set up physical barriers. Thus, in the grid cells with vegetation, we enhanced the friction Table 2. Comparison between the DEM difference used in this work as a proxy for the presence of higher vegetation, and land cover obtained from a national land cover map, at 10 m resolution, obtained from the ISPRA website. Values are percentages of each land cover type within the pixels in which the DEM difference is larger than 1, 5 or 10 m; pixels with buildings and roads (cf. Figure 1) were excluded from the calculations. coefficients in Table 1 by 50% and reduced both the restitution coefficients by 50%. We considered as vegetation only the grid cells with a positive difference between the DSM and DTM, larger than 5 m. We decided to use grid cells with DEM differences larger than 5 m as follows. We used a land cover map (ISPRA 2012; 10 m resolution), and calculated the percentage of cells falling in each land cover class where the DEM difference was larger than a threshold. Using three different thresholds, namely 1 m, 5 m and 10 m ( Table 2). We observed that 5 m provided a reasonable balance between the accumulated percentages of classes which potentially affect rockfall trajectories (low values in "broad leaved," "needle leaved," "orchards" and "vineyards") and classes which surely correspond to false positives for the DEM difference interpreted as upper vegetation, namely "artificial surfaces." Thus, we selected 5 m as a minimum difference to modify terrain parameters in STONE.

Localization of rockfall sources
The most relevant input of the program STONE is a raster grid specifying the location of sources. For small areas, i.e. at slope or small catchment scale, sources can be mapped on the field and by photointerpretation. On larger areas, either procedures are impractical and sources must be inferred in some probabilistic way. We follow a similar approach as in Alvioli et al. (2021), where the analysis concerned the whole Italian railway system. Sources were localized by morphometric analysis on a few sample locations mapped by experts, and statistical generalization to the whole area of interest identifies locations of possible future sources.
Here, we used the results of the generalization procedure of Alvioli et al. (2021) as a starting point, and applied additional cuts based on slope angle and relief (elevation range) values. Additional cuts are based on mapping of possible sources of rockfalls on Google Earth TM images, shown as blue polygons in Figure 1. The expert mapping analysis consists in a detailed investigation of the location of potential sources of rockfalls through the photointerpretation of satellite images spanning a period of about 10 years (Guzzetti et al. 2004;Santangelo et al. 2019Santangelo et al. , 2021. The mapped polygons were used to calibrate the statistical procedure providing additional information on the local conditions that may induce detachment of blocks. We analyzed the distribution of values of slope and relief within the mapped polygons, first, and of the whole study area, then. Figure 5. The latter can be calculated in different ways, as virtually any morphometric quantities derived from a DTM (Trevisani and Cavalli 2016;Voros et al. 2021).
Generally speaking, relief is defined in each grid cell as the difference between the maximum and minimum values of elevation in a neighborhood of that cell; thus, it clearly depends on the size of the neighborhood. We investigated the values of relief within mapped sources, in relation with values in the whole study area, for different sizes of the neighborhood, using the module r.neighbors in GRASS GIS. It turned out that two specific sizes, namely square neighborhoods of size 15 and 25 cells, have distributions that helps generalizing the morphometric characteristics of mapped sources to the whole area. This is clearly shown in the ratio of histograms corresponding to values of slope, 15 Â 15 relief and 25 Â 25 relief within the sources and in the whole area, in Figure 5.
From the shape of the ratios, we concluded that the probability of any grid cell in the study area to be a source of rockfalls is much larger for increasing values of slope, 15 Â 15 and 25 Â 25 relief. Thus we assigned each grid cell a probability proportional to the ratio source/whole area of slope as shown in the figure, P(S). In addition, we heuristically defined thresholds for relief, i.e. we set the grid cells with joint values of 15 Â 15 relief below 35 m and 25 Â 25 relief below 45 m to null probability. We refer to the relief-constrained probability as P(S, R), where R ¼ ðR 15 , R 25 Þ: The probability, in turn, controls the number of trajectories initiating from each grid cell of the DTM, providing a probabilistic output for the model. The number of simulated trajectories per source cell varies from a minimum of 10 trajectories, at lower probabilities, to a maximum of 100, for probability equal to unity. The resulting raster source map is called SRC Prob , where SRC stand for "source" and where Prob stands for "probabilistic."

A new model for seismically induced rockfalls
To investigate the impact of earthquake-induced rockfalls, we illustrate here a new model to link ground shaking to the natural predisposition of slopes to be affected by rockfalls. The method was developed within a national project called FRA.SI -"Integrated multi-scale methodologies for the zonation of earthquake-induced hazard in Italy," jointly developed by CNR IRPI, IREA and IGAG (https://frasi-project.irpi. cnr.it, in Italian). Basic assumptions of the method are that seismic shaking is one of the possible triggers of rockfalls which typically affect slopes with specific topographic characteristics, or where rockfall occurred already (Govi 1977). We briefly summarize the rationale of the method, here below.
We consider the probabilistic source map described in Section 4.3 as the baseline map, in which the different grid cells with non-null probability can be actually activated (i.e. trajectories initiating from these cells are actually simulated within STONE) as a function of ground shaking due to an earthquake. We selected peak ground acceleration (PGA) as a scalar, distributed measure of the intensity of an earthquake event. We defined the overall probability of a grid cell with given values slope angle S, relief R ¼ ðR 15 , R 25 Þ, and PGA of representing a source of rockfalls as follows: where the first factor on the right hand side, P(S, R), is the relief-constrained probability as a function of S, R 15 and R 26 , defined in Section 4.3, and the second factor is a damping factor. The latter is a linear function of PGA, and assumes values between 0 and unity; thus, it sets to zero the overall probability P seismic ðS, R; PGAÞ where PGA is null, and leaves P(S, R) unmodified where PGA ¼ PGA max . The values PGA min and PGA max in Eq. (1) were determined from the minimum and maximum values of expected PGA calculated from return times of 475 y and 975 y, all over Italy, with the model of Mori et al. (2020); Falcone et al. (2021), and are 0.0 and 0.81, respectively (in units if g, the Earth's acceleration of gravity). The rationale for using these minimum and maximum values rely on the possibility of using the same method of defining P seismic ðS, R; PGAÞ in a parametric way, as a function of values of PGA contained in specific maps, corresponding to specific earthquake events, with return times within the range 475 y and 975 y. Thus, the method is applicable both using scenario-based PGA maps.
For the specific PGA map available for the study area of interest in this work, we obtained a raster source map called SRC seismic , using Eq. (1). We stress, here, that we used the PGA scenario of Rapolla et al. (2010) because that was explicitly associated with a specific return time. This is consistent with the meaning of the parameters PGA min and PGA max of Eq. (1), which were calibrated against national PGA maps for the same specific scenarios. Use of a different PGA maps, corresponding to as specific event as e.g. in Albano et al. (2018), would be inconsistent with such definition, as it would require an absolute calibration of the parameters PGA min and PGA max against a number of specific events for which rockfalls were observed and mapped. This is beyond the scope of this work; such kind of study is underway and will be presented elsewhere (Alvioli et al. 2022b).

Simulation of rockfall trajectories
The software STONE, adopted here, assumes point-like masses falling under the sole action of gravity and the constraints of topography, and it calculates trajectories dominated by ballistic dynamics during falling, and bouncing and/or rolling on the ground in a probabilistic way. The assessment of trajectories derives as consequence of several of steps, the first of which consists in the identification of perspective rockfall source areas. Identification of potential rockfall sources is a key step of numerical, physically based simulation of rockfall runout Section 4.3).
The software STONE was designed to assess rockfall hazard at the regional and local scales using thematic data (e.g. geological, geomorphological, land use maps etc.), available for large scale or obtained in field surveys (Guzzetti et al. 2002). In the model, a kinematic simulation is performed computing the trajectory at discrete time steps. Within each time step the boulder can be in a free falling state, in a rolling state, or in a bouncing state (Guzzetti et al. 2002). The trajectory of each boulder is computed from the digital topography and it depends on the starting point, and a set of coefficients used to simulate the loss of velocity at the impact points or during the rolling state, i.e. dynamic friction where the boulder rolls, normal and tangential restitution coefficients values at each impact point.
Attribution of ranges values for dynamic friction and normal and tangential restitution coefficients depends on the lithotype was inferred from a geo-lithological map ( Figure 3). We assigned values by comparison of the information available from the literature for each lithological class (Guzzetti et al. 2002;Alvioli et al. 2021). For each simulation the program allows a random variation in the coefficients (within 10% of the nominal values listed in Table 1, uniformly distributed) and in the initial direction of motion, resulting in an output with probabilistic meaning.
Furthermore, since the vegetation can be considered an important element that can influence and reduce the distance that a rolling boulder can reach, the friction between the falling block and the terrain was increased by 50% and the normal and tangential restitution coefficients were reduced by 50% in all such cells. Thus, we have run STONE with two identical sets of input data, except the set of numerical parameters describing the soil response, corresponding to considering or not the effect of vegetation.  without considering the presence of vegetation and contemplating the effect of vegetation. We decided to use these thresholds on the basis of the results of diagrams shown in Figure 5. The analysis of Figure 5 highlights that in cells mapped as sources  (a) impact on buildings; red squares represents building which potentially could be hit by falling boulders and yellow squares represent building identified as potential hits by by STONE but ruled out as "false positive" in an expert way. (b) roads; the color ramp represents intersection with a very high (purple), high (red), medium (orange) and low (green) number of trajectories. the probability of observed values of slope and relief higher than a certain limit could be used to locate sources through a probabilistic analysis. Therefore, we considered the "static" source probability as a function of slope according to the normalized ratio of Figure 5b; probability; then, we used the values of relief calculated in a 15 Â 15 and 25 Â 25 cells neighborhoods as thresholds to establish which of the identified sources have more probability of representing a source of falling boulders (Figure 5d and f). We note, here, that considering relief only provides a correction to the probabilities obtained using slope alone.

Results & discussion
As highlighted in Figure 7a and b the presence of vegetation reduces the reach of boulders trajectories, by construction. As explained in section 4.2 we considered as vegetation the grid cells with a positive difference between the DSM and DTM larger than 5 m. We decided to use the threshold of 5 m because, as shown in Table 2, it represents a good compromise between larger percentages of land cover classes which can influence rockfall trajectories ("broad leaved," "needle leaved," "orchards" and "vineyards") and small percentage in the "artificial surfaces" class.
We stress that LiDAR data provide high resolution DEMs suitable for physicallybased slope stability models, including STONE and comparable models (Matas et al. 2017;Dorren et al. 2022;Jones et al. 2000;Agliardi and Crosta 2003). The strategy was applied before, though with no distinction between DTM and DSM (Santangelo et al. 2019). High-resolution elevation data are useful both to run morphometric analysis and locate potential rockfall sources and to carry out a detailed analysis of boulders trajectories. In the present case, availability of both DTM and DSM data allowed to include the effect of vegetation.
As one can see in Figure 6a and b, the areas most susceptible to rockfall events are those located along the South and South-West coast of the island, characterized by steep cliffs, the North, North-West and South-West steep flanks of Mount Epomeo and the exposed sides of the hydrographic channels located in the South sector of the island. This is a qualitative confirmation of the findings by Della Seta et al. (2011); del Prete and Mele (2006); the polygons denoted by diffused rockfalls in the IFFI inventory also overlap with the runout calculated with STONE. Moreover, Ischia is a densely populated island (about 1,380 people per km, 2 vs. an average of about 200 people per km, 2 in Italy) and, as shown in Figure 6, it is characterized by the presence of potential rockfall sources on several places. Thus, built up areas ( Figure 8a) and roads (Figure 8b) could be affected by rapid phenomena as rockfalls. Figure 6a and b shows the result of the overlap of cell-by-cells trajectory count (Figure 6b) with roads and buildings on the island (Figure 1), to assess the possible impact of rockfalls on infrastructure. We stress that, at this stage, trajectory count only represents susceptibility maps and intersection with vulnerable infrastructure cannot be considered as a risk map, because magnitude and the temporal components are not considered up to now. Moreover, potential source areas were identified through a statistical approach, derived from a morphometric analysis of terrain within polygons mapped from Google Earth TM images. Therefore a field survey to ascertain actual potential sources and would be useful to validate the inferred source areas. The same goes for evidences on the field of actual deposition areas of rockfalls occurred in the past.
Rockfall sources were established on the basis of slope and relief values. Therefore, the trajectory count output of STONE represent susceptibility maps, in which false positives are often possible. On the other hand, to go beyond traditional susceptibility maps and give the best possible estimate of rockfall hazard, one can consider a specific trigger, with a specific probability in timeor return time, as we did in this work. We considered a seismic trigger for rockfalls, thus approximating an assessment seismically-induced rockfall hazardgiven that a magnitude component would still be missing. The rockfall trajectories simulated on the basis of the PGA map for the major Casamicciola earthquake of 1833, obtained from Rapolla et al. (2010), and corresponding to a return time of 475 years, show that the source map differs from the probabilistic source map described in Figure 6 first of all for the number of potential sources: less sources have been recognized on the whole island in the SRC seismic in comparison of that recognized in the baseline ("static") probabilistic source map ( Figure 6). Furthermore the overall trajectory count in the result is smaller than in the full simulation, as one can see comparing Figures 6 and 9. Almost all of the earthquake-induced rockfalls are located on the steep N-NW flanks of Mount Epomeo, on the SE and SW cliffs of the coastline.

Conclusions
In this study we performed a statistical analysis to infer the potential sources of rockfall on the Ischia island, their possible trajectories both in a "static" (trigger- Figure 9. Simulation of seismically-induced rockfall trajectories, using a seismic trigger corresponding to an event with return time 475 y and with a peak ground acceleration with the spatial pattern of Figure 4. The model adopted in this work singles out rockfall sources on the basis of slope, S, and values of PGA, as defined in Eq. (1) (Alvioli et al. 2022a(Alvioli et al. , 2022b. independent) scenario corresponding to a susceptibility assessment, and in a seismic shaking triggering scenario, going in the direction of seismically-induced rockfall hazard. We performed our analysis using the 3D model STONE, which needs input data of potential rockfall sources, geo-lithological characteristics of the studied area to assign terrain parameters, and a DEM describing the topography. We can draw the following conclusions: We improved an existing morphometric method to generalize expert-mapped rockfall sources, previously applied at national scale and at 10 m resolution . In this work, at 2 m grid resolution, we included relief calculated with two different moving windows sizes along with slope angle, to generalize mapped polygons to additional potential sources of rockfalls with a probabilistic method. Raw LiDAR data are useful to distinguish areas covered by tall vegetation, allowing a more realistic simulation of rockfalls trajectories and their traveling distance. This in principle can be reproduced in all of the areas where LiDAR points are available from the Geoportale Nazionale, a large portion of Italy. Preliminary assessment of rockfall susceptibility in Ischia showed that the areas with highest susceptibility are located along the N, N-W and S-W steep flanks of Mount Epomeo, the S and S-W coast of the island. Many of the susceptible areas overlap roads and built-up areas. We implemented a new method to model a seismic trigger for rockfalls, for physically-based rockfall simulations with models that are otherwise completely independent from a specific trigger, and from time. This effective method allows to approximate seismically-induced rockfall hazard.
Many of the methods presented here can be further improved, either with additional modeling and/or with data from field surveys. In the inclusion of the effects of vegetation, we simply modified the terrain parameters describing the static and dynamic response in STONE to effectively account for the presence of vegetation; future work may focus on calibrating the parameters, with proper field data. A LiDAR campaign was performed recently by a few of us, and results will be reported elsewhere. The survey will also be useful to locate actually unstable area, to further constrain the maps of sources proposed here, and to pinpoint deposition areas or individual blocks from past rockfalls, to constrain the parameters of the model.
Eventually, we stress that the method presented here is amenable for application both a small, medium and large scale. While this work is an example application on the small scale, the same approach was applied at national level (Alvioli et al. 2022a) using ground shaking maps corresponding to different return times (Mori et al. 2020;Falcone et al. 2021). Ground shaking maps corresponding to a return time of 475 y was actually used to calibrate the seismic trigger for "static" rockfall sources adopted here. The same method can be applied at intermediate scale using ground shaking maps corresponding to specific earthquake events. The final aim of the seismic-trigger method is the application on near-real time, with ground shaking maps obtained after an earthquake event. Note 1. http://www.pcn.minambiente.it/mattm/en/