Lava flow hazard assessment
The methodology, developed to achieve a reliable and comprehensive assessment of lava flow hazards at Mt. Etna, is founded on the MAGFLOW numerical model for simulating lava flow paths16,17,18. Numerical simulations are a powerful tool to explore various eruption scenarios, as they can be used to estimate the extent of the lava field, the time required for the flow to reach a particular point and the resulting morphological changes19,20,21. The numerical simulations of lava flow paths are based on our knowledge of Etna's past eruptions, derived from the integration of historical and geological data and by adopting a high-resolution updated digital elevation model (DEM). Accounting for a more comprehensive range of possible scenarios, through the combination of field data, numerical modeling and probabilistic analysis, gives a better understanding of effusive eruptions and their effects. However, such a combination is rarely straightforward and requires an accurate evaluation of diverse stages. For this reason, we developed a methodology based on four different tasks (Methods): (i) assessment of the spatiotemporal probability of future vent opening, (ii) estimation of the occurrence probability associated with classes of expected eruptions, (iii) simulation of a large number of eruptive scenarios with the MAGFLOW model, and (iv) computation of the long-term probability that a lava flow will inundate a certain area. We applied this four-stage methodology to estimate lava flow hazard at Mt. Etna both for flank eruptions and summit eruptive activity.
Hazard assessment for flank eruptions
The historical records of flank eruptions at Etna date back over 2000 years, but are continuous and reliable only from the 17th century22. The eruptive behavior of Etna in the last 400 years is undoubtedly irregular, with important fluctuations in frequency, type of eruptions and lava discharge rates, both in the long (centuries) and short (decades) terms6,23,24,25. During the past century, major flank eruptions causing significant damage to crops and disruption in towns and villages have occurred in 1910, 1923, 1928, 1971, 1979, 1981, 1983, 1991–1993, 2001, and 2002–20035,8,22. An evident increase in the lava output and frequency of eruptions began in 1971. Since then, 17 eruptions have occurred mostly along the main rift zones (Fig. 1), showing that the volcano has been in a strongly active period during the last 400 years8.
In flank eruptions at Mt. Etna, lavas erupt from newly formed fissures and vents, hence, the potential spatiotemporal distribution of new vent openings must be estimated as part of the analysis. In order to identify the most probable emission zones of future lava flows, we analyzed the spatial location of past eruptive vents, as well as the eruption frequency within a time window. Recently, we conducted a detailed study22 of the main structural features of flank eruptions at Etna, including the outcropping and buried eruptive fissures active in the past 2000 years, the dykes in the Valle del Bove depression7,26 (Fig. 1a), and the main faults that can potentially be used as pathways for intruding magma27 and/or influence the surface stress field of the volcano. It is noteworthy that shallow feeder dykes guided by faults are not very rare on Etna, such as in the 1928 eruption27, while in other volcanoes worldwide the use of faults as channels is uncommon28. All these geological and structural features were used to construct the spatial probability map of vent opening at Etna27. Thanks to the completeness and accuracy of historical data over the last four centuries, we demonstrate spatial non-homogeneity and temporal non-stationarity of flank eruptions on Etna, showing that effusive events follow a non-homogenous Poisson process with space-time varying intensities8. Here, we assess the distribution of future vent locations using a probabilistic modeling (Methods) based on the locations of the pre-1600 eruptive fissures, dykes and faults (for which an exact date cannot be established), and on the spatial positioning and temporal sequencing of the post-1600 flank eruptions (that are accurately documented and well-dated). We calculate the recurrence rates (events expected per unit area per unit time) and produce a spatiotemporal probability map of new vent opening for the next 50 years (Fig. 2a). The distribution of probabilities is strongly non-homogeneous. The highest probability of new eruptions occurs in the areas closer to the summit of the volcano, at elevation higher than 2500 m. Diversity in the probability estimates also occurs in lower elevation areas, between 2500 and 1800 m. Below 1000 m, the expectation of new vent opening is very low, except in the South Rift, where probabilities slowly decrease to an elevation of ∼600 m.
To estimate the probability of occurrence of eruptions, a complete revision of location and opening dynamics of all flank eruptions occurring in the past 400 years was carried out. For each lava flow field, we collected the main quantitative volcanological data, concerning the eruption start and end, the eruptive style, the area covered by flows, the lava volume emitted, and the fissure producing the lava flow. From an accurate analysis of the distribution of flow durations and lava volumes, we selected a time barrier at 30 days and two divisions for the total volume emitted, at 30 and 100 × 106 m3. Thus we obtained six eruptive classes, five of which are populated (Fig. 2b). This classification of expected eruptions permitted calculation of local occurrence probabilities for each class, and the deduction of essential information for the lava flow simulations (Methods).
The eruptive scenarios were simulated using MAGFLOW29,30, which is based on a physical model for the thermal and rheological evolution of the flowing lava. To determine lava flow emplacement, MAGFLOW requires several input parameters: a digital representation of the topography, the chemical composition of the lava, an estimate of the effusion rate, and the location of the eruptive vent. As digital representation of the topography, a 10-m resolution DEM of Etna updated to 2007 was utilized. For the chemical composition of the lava, the typical properties of Etna's basaltic rocks31,32 were used (Methods). The effusion rate, e.g. the discharge rate controlling how a lava body grows, was derived from the six eruptive classes obtained from the characterization of expected eruptions. Since the numerical simulations require flow duration and lava volume to be defined, short and long time periods were set at 30 and 90 days respectively, while 30, 100 and 200 × 106 m3 were fixed as values for total volume of lava emitted. These values of duration/volume produced six possible bell-shaped curves16 representing the flux rate as a function of time (Fig. 2c). Locations of eruptive vents are the nodes of a regular grid with boundaries shown in Fig. 1. For every potential vent, six simulations were executed, each one representative of a particular eruptive class.
Finally, assessment of the long-term hazard related to lava flows at Mt. Etna rests on combining the spatiotemporal probability of future opening of new eruptive vents, occurrence probabilities for each class of expected eruption, and the overlapping of a large number (28,908) of numerical simulations (Methods). The hazard map for the next 50 years is shown in Fig. 3. The highest hazard level is reached in the well-delimited zone inside the Valle del Bove, a 5 × 7 km morphological depression on the eastern flank of the volcano able to capture lava flows emitted from the eruptive vents below the summit craters, and along the upper portions of the South and North-East Rifts. Only eruptions exceptional in duration and lava output could form lava flows capable of travelling across the entire valley and seriously threatening the towns of Zafferana Etnea and Santa Venerina in the eastern sector. Other areas more likely to be threatened by lava flows are those on the southern flank, including densely populated areas near the towns of Nicolosi, Pedara, Trecastagni, and Ragalna. Towns located at lower elevations could be equally threatened by lava flow inundation, though the hazard progressively diminishes at greater distances from volcano summit. On the northern and western sectors of the volcano, the main population centers exposed to lava inundation are Linguaglossa, downslope of the North-East Rift, and Bronte, located along the likely trajectory of lava flows emitted from the West Rift.
Validation of the lava flow hazard map is an additional task in our four-stage methodology, giving a quantitative evaluation of its reliability. As historical eruptions provide the only available information, we propose a retrospective validation using more recent past eruptions to indicate potential future events33. In particular, we calculate a fitting score to measure the overlap of the probability distribution of lava flow hazards at Etna up to the year 1981 against the historical lava flow paths that occurred after 1981 (Fig. 1). The good agreement between the observed and expected inundation areas obtained for the hazard map up to 1981 confirms the consistency of our results33.
Hazard assessment for the summit eruptive activity
The summit area of Mt. Etna has frequently undergone major morphological changes due to its persistent eruptive activity, both effusive and explosive (Fig. 1a). A single, large Central Crater (CC) existed until the beginning of the 20th century, but soon afterwards the summit area was affected by repeated subsidence and collapse phenomena alternating with intracrater volcanic activity and the birth of new summit vents34. The North-East Crater (NEC; formed in 1911) opened a few hundred meters north of the CC rim. The Voragine (VOR; 1945) and the Bocca Nuova (BN; 1968) opened inside the CC. Finally, the South-East Crater (SEC; 1971) is the youngest and currently the most active summit vent, as also testified by the growth of a large new cone on its southeast flank, identified as the New South-East Crater1 (NSEC; 2007). The progressive increase in summit eruptive activity during the past 110 years is confirmed by the results of the statistical analysis8, indicating that the hazard from eruptive events is not constant with time and differs for each summit crater of Mt. Etna. A clear increase in the eruptive frequency is evident from the middle of the last century and particularly from 1971, when the SEC was formed.
The probability that a lava flow will inundate a certain area was estimated assuming an effusive eruption originated from a vent within the summit area of Mt. Etna. This was done using the methodology adopted for flank eruptions, including the development of the spatiotemporal probability map of future vent opening, the occurrence probabilities associated with the classes of expected eruptions, and numerical simulation of a number of eruptive scenarios.
Firstly, we calculated a spatiotemporal probability map of future vent opening for the next 10 years using the spatial location of the existing four summit craters (CC, NEC, SEC and NSEC) and the frequency of eruptions since 1900 (Methods). This produced the spatiotemporal probability map for future vent opening shown in Fig. 4a.
Subsequently, the classes of expected summit eruptions were identified. We used flow duration/lava volume distributions for paroxysmal events and Strombolian activity at the summit craters1,2, to distinguish ten classes. For the paroxysmal events, lava flows and lava fountains occurring at VOR and SEC since 1998 were considered. We established three different duration intervals (≤4, 4–8 and >8 hours), and emitted lava volumes (≤1, 1–2 and >2 × 106 m3), resulting in nine possible eruptive classes1 (Classes I–IX; Fig. 4b). The last class of expected eruptions (Class X) was identified on the basis of Strombolian activity from 1955 to 1996. Classification of expected eruptions enables probabilities of occurrence to be calculated (Methods), as well as deriving the ten effusion rate trends associated with each class for the MAGFLOW simulations. For Classes I–IX of paroxysmal events, we used the bell-shaped curves shown in Fig. 4c. Since the Strombolian episodes are often characterized by a continuous mild effusive activity, we assigned to the Class X a constant rate equal to the Mean Output Rate (MOR ≈ 0.5 m3 s−1; see Fig. 4b) for the entire duration (150 days) of the eruption.
To evaluate areas more likely to be affected by lava flow paths we computed ten simulations, corresponding to the expected classes of eruptions (Classes I–X), for each of the potential vents of the grid defined around Etna's summit craters (Fig. 4a). The topographic base was the 10-m resolution DEM of Etna updated to 200734, while the chemical composition of the lava was defined using the properties of products erupted from Etna's paroxysmal events35.
The hazard map for summit eruptions was obtained by combining the spatiotemporal probability of future vent opening, the probability of occurrence for each eruptive class and a large number of simulations (5,950) computed with the MAGFLOW model. The resulting hazard map (Fig. 5) estimates probable areas of inundation by lava flows issuing from vents within the summit area over the next 10 years. The inundated area measures about 51 km2, while the maximum distance of the simulated lava flow paths is reached at the end of Valle del Bove (∼700 m elevation), and in the direction of Linguaglossa, ∼1200 m elevation. Close to the summit area, the northern tourist facilities (Piano Provenzana) located above 1800 m are exposed to lava inundation, although no simulated flows erupted from the summit vents inundated the southern tourist facilities (Rifugio Sapienza) at ∼1900 m. New pyroclastic cones formed on the southern flank during the 2001 and 2002–2003 eruptions, created an effective topographic barrier to these lava flows, diverting lava west or east of the cones.
 The flanks of several basaltic volcanoes are highly populated, and the towns, cities and associated infrastructures may be impacted by lava flows during effusive eruptions. Lava flow hazards on the flanks of such volcanoes (e.g., Mauna Loa, Mount Etna) are of great concern for the protection of civilians because of the frequent eruptions and the growing vulnerable population [e.g., Trusdell, 1995; Rowland et al., 2005]. It is difficult to assess lava hazard through computer simulations of lava flow [e.g., Crisci et al., 2004] because of the evolving topography during an eruption: as new lava flows are emplaced, they produce new topographic features that can influence the paths of subsequent lava flows. Lava flows are unconfined, multiphase and multicomponent flows in which temperature, rheology, and effusion rates all vary in time and space. Given the high complexity of these processes, in most real case applications it would be impractical to perform a direct numerical simulation of the complete set of flow conservation laws. Because the real-time forecast of lava flow paths is essential during volcanic crises, a number of “simplified” models have been developed; these models are mainly based on the assumption that lava follows the steepest gradient downhill and spreads laterally according to its critical thickness and mass conservation [e.g., Wadge et al., 1994]. Other models are based on a probabilistic definition of lava flow spreading [e.g., Favalli et al., 2005]; they generally involve low CPU execution times and do not require accurate specification of the physical properties of lava (e.g., flow rate, temperature, cristallinity, viscosity, yield strength). Other models are based on cellular automata or neural network approaches [e.g., Crisci et al., 2004] or on various simplifications of the governing physical equations [e.g., Dragoni et al., 1986; Harris and Rowland, 2001; Costa and Macedonio, 2005].
 Here we use a modified version of DOWNFLOW [Favalli et al., 2005]. The previous versions of the DOWNFLOW code [Favalli et al., 2005, 2006] predicted which areas would be inundated by lava flows erupting from a single vent. The new DOWNFLOW version adopts probability distributions for vent opening and lava flow length distributions to create maps of lava inundation probability.
 We applied this new model to lava flow hazard assessment on Mount Etna volcano, located on the East coast of Sicily, Italy. Mount Etna is characterized by numerous (near yearly) eruptions accompanied by effusive and moderate explosive activity [Behncke and Neri, 2003; Neri and Acocella, 2006; Allard et al., 2006]. Eruptions from both summit and flank vents at variable effusion rates have produced numerous composite lava fields [e.g., Harris et al., 2000; Branca and Del Carlo, 2004; Calvari and Pinkerton, 1998; Romano and Sturiale, 1982]. Much of the erupted lava has flowed into a deep scar resulting from flank instability on the eastern flank of the volcano (Valle del Bove) [Borgia et al., 1992].
 Existing hazard zonation maps of Etna have been compiled considering the main volcanological parameters of all eruptions since A.D. 1600 (e.g., including vent location, age, type of eruption, volume of lava flow and tephra) and the areas invaded by lava flows since 700 B.C.; they assume that the characteristics of future eruptions will be similar to those of past eruptions [Guest and Murray, 1979; Andronico and Lodato, 2005; Behncke et al., 2005]. We propose a different hazard assessment approach based not only on the frequency of coverage of past lava flows, but also on the numerical simulation of possible lava flow paths, taking into account the actual topography over which lavas flow. Updated high-resolution digital elevation models (DEMs) of Mount Etna have recently been derived from LiDAR (Light Detection and Ranging) data [Mazzarini et al., 2005, 2007]. Using vent and lava flow length distributions based on the statistical analysis of past eruptions, DOWNFLOW performs simulations of all possible future lava flow paths and produces detailed quantitative maps of the probability of inundation by lava flows. This is possible because (1) DOWNFLOW accounts for the most relevant emplacement characteristics of real lava flows, such as filling of topographic depressions, overriding of topographic obstacles and lava flow spreading, without requiring detailed input parameters; (2) the code is able to compute invasion areas from each vent in a very short computational time; and (3) simulation outputs require relatively little disk space. Lava flow paths from all possible vents can therefore easily be simulated and stored, results can be postprocessed to produce quantitative hazard maps for different scenarios, and lava catchment areas can be identified rigorously by considering all the vents from which lava flows can invade a target area.
2. Mount Etna Digital Elevation Model
 Airborne altimetric LiDAR data were used to generate a high-resolution (2 m) DEM of most of Mount Etna from data acquired during flights on 29 September (0800–1100 local time) and 30 September (0800–1030 local time) 2005. LiDAR data have already proved useful in analyzing the morphology of rapidly evolving features of the landscape, for example active volcanoes [e.g., Mazzarini et al., 2005, 2007; Ventura and Vilardo, 2007], as well as in lava flow modeling [e.g., Harris et al., 2007].
 An Airborne Laser Terrain Mapper 3033 (Optech®) laser altimeter (see Table 1 for instrument characteristics) was used to survey an area of 616 km2. Measurements were made at a frequency of 33 kHz, resulting in a mean ground point density of 1 point per 2.4 m2 (generally one point every 1.5 × 1.5 m2). The LiDAR survey consists of more than 2.57 × 108 data points. These points (echoes) were distributed along 34 NNE-SSW trending strips covering part of Etna's western flanks and the majority of its eastern flanks. The resulting DEM, geocoded to a UTM-WGS 84 projection, has an elevation accuracy of ±0.4 m and a horizontal accuracy of ±1.5 m.
|Operating scan frequency||33–40 kHz|
|Operating altitude||80–3500 m|
|Maximum scan angle||±40°|
|Horizontal accuracy||1/2000 × altitude at 1σ|
|Elevation accuracy||±15 cm at 1.2 km and ±35 cm at 3000 m at 1σ|
|Operating wavelength||1.064 μm|
 At the LiDAR operating wavelength (λ = 1.064 μm), atmospheric absorption is minimum; however, during an eruption, ash and gases can greatly hamper data acquisition. Likewise, persistent degassing such as that which feeds a constant plume at Etna can cause unwanted atmospheric effects, especially near the summit craters. Our September 2005 LiDAR survey (with flight altitudes ranging from 1500 to 3500 m) occurred during a period of low volcanic activity characterized by clear atmospheric conditions. A 10 m step DEM was derived from the LiDAR data. Outside the region covered by the LiDAR survey, the LiDAR DEM was merged with a 10 m step DEM of all of Italy created by Tarquini et al. .
3. Lava Flow Simulation
3.1. DOWNFLOW Probabilistic Code
 DOWNFLOW [Favalli et al., 2005] computes maximum slope paths on a stochastically perturbed topography. It is based on an assessment of the steepest descent path, which is derived multiple times (generally thousands of times) on a randomly perturbed topography in an attempt to quantify all possible flow paths. For each emission point, only two model parameters must be input. The first parameter (Δh) represents the maximum vertical perturbation (either positive or negative), which is applied to the DEM during each iteration. The second parameter (n) represents the number of steepest descent path calculations performed for any single run from a user-defined emission point [Favalli et al., 2005, 2006].
 Favalli et al.  performed a dimensional analysis of steady state lava flows on an inclined plane; results indicated that lava spreading is controlled by the competing viscous and gravitational/self-gravitational (hydrostatic) forces. These forces can be characterized by a vertical scale length H that approximates the size of obstacles which can be overridden by the flow. Stochastic variations in ground elevation by a characteristic vertical height equal to 2Δh ≅ H employed in the computation of maximum slope paths therefore account for first-order variations in lava flow spreading. An appropriate value of Δh can be obtained through the best fit of simulated and actual lava flows; such a value relates to the characteristic height of obstacles capable of diverging the flow. In the case of Mount Etna, the best fit to lava flow maps for the 1991–1993, 2001, and 2002 flank eruptions is given by Δh = 3 m and n = 10,000 [Favalli et al., 2005]. For example, the simulated inundation area of the 2001 flow covers 90% of the actual flow. We thus assume that these Δh and n values are appropriate for typical flank eruptions at Etna. Comparison between real and simulated flow fields at Etna and Nyiragongo indicates that the DOWNFLOW approach yields realistic emplacement areas [Favalli et al., 2005, 2006].
 We applied DOWNFLOW to Mount Etna, for which a high-resolution LIDAR-derived DEM is available. We checked the sensitivity of DOWNFLOW against the DEM resolution in the case of the southern 2001 lava flow at Mount Etna [Favalli et al., 2005]. The 1998 DEM of Tarquini et al. , with a cell size of 10 m, was used as a topographic base. The DEM root mean square (RMS) vertical error in the study area is 1.43 m. By resampling the original DEM, we generated a series of new DEMs with different cell sizes (up to 120 m) and vertical errors (up to 4.23 m), as reported in Table 2. We then performed simulations with DOWNFLOW for each newly generated DEM, using the 2001 vent as the starting point, with Δh = 3 m and n = 10,000. To compare the simulations with the actual flow, for each simulation we calculated the parameter μ
where AS is the area covered by our simulation and AR is the area covered by the actual 2001 flow; μ is 1 when these two areas coincide. In the case of the 10-m DEM we obtained a μ value of 0.52, which indicates that there is good agreement between the simulated flow and the actual flow. As shown in Table 2, the agreement between the simulated flow and the actual 2001 lava flow is almost constant up to grid cells of 90 m (or up to RMS vertical errors of 3.27 m), whereas μ values gently decrease for higher cell dimensions. The probabilistic nature of DOWNFLOW is such that intrinsic errors in the adopted DEM are overcome, since random perturbations of the topography already add random “errors” to the topography.
 Recently DOWNFLOW was applied on DEMs of coarser resolution. In particular, DOWNFLOW has successfully reproduced the 1977 and 2002 lava flows at Nyiragongo Volcano using a 90-m pixel DEM derived from the Shuttle Radar Topographic Mission (SRTM; http://www2.jpl.nasa.gov/srtm) as a topographic base [Favalli et al., 2008]. This indicates that DOWNFLOW can be used with the SRTM digital elevation model for most volcanoes in the world where high-resolution data are not yet available.
 However, the model allows neither a description of the time evolution of lava flows and of the effects of effusion rates, nor the definition of flow lengths. Accordingly, the length distribution for the simulated lava flow must be input.
3.2. Probability Distribution of Vents
 Although eruptions at Mount Etna occur from both the summit and flanks, we focus on the hazard due to flank eruptions, which have proved to be the most dangerous in the past [Chester et al., 1985; Mulargia et al., 1985; Crisci et al., 2003; Andronico and Lodato, 2005; Behncke et al., 2005] and often have the highest effusion rates [Harris et al., 2000; Wadge, 1981]. According to Guest and Murray  and Chester et al. , certain flank areas of Mount Etna are more likely to produce eruptions than others; these areas correspond to zones with a high density of cones and eruptive fractures, i.e., the NE, S and W rift zones [e.g., Kieffer, 1975].
 The probability of a vent opening in a certain location is usually based on the assumption that future vents will form in zones with the highest density of old vents [Guest and Murray, 1979; Behncke et al., 2005]. The location of more than 400 vents was thus extracted from LiDAR-based DEMs of Mount Etna and published maps [Mazzarini and Armienti, 2001]. From this discrete set of points we estimated the vent opening probability distribution through a kernel smoothing technique [Bowman and Azzalini, 2003] using a Gaussian kernel with a bandwidth of 1 km. In Figure 1 we plot the smoothing kernel over the vent spacing distribution. The distribution of vent spacing (vent separation) was calculated by computing the distance between each vent and the nearest neighbor vent. Results show that there is a very low probability that a vent will open less than 80 m from another vent. The probability distribution for vent opening was thus assessed using the vent density map and the 80 m proximity threshold (Figure 2). This vent opening distribution does not, of course, take into account the possibility of vents or vent clusters opening in areas where vents are currently not exposed at the surface. This possibility is taken into account, as described in the following sections, by using a uniform vent distribution.
3.3. Lava Flow Length Distribution
 Following Favalli et al. , data from Guest , Behncke et al. , Allard et al. , and Neri and Acocella  were plotted on a flow length versus vent elevation diagram (Figure 3). In the last two decades eruptions have occurred above 2000 m and have been characterized by a maximum runout distance of 8.75 km [Behncke et al., 2005].
 According to Guest  and Lopes and Guest , at Mount Etna there is an inverse correlation between vent elevation (hv) and the length of the erupted lava flow (L): the lower the elevation the longer the lava flow. The points in the hv-L space for the lava flows at Mount Etna (Figure 3) are bounded by a straight line defined by the following equation [Favalli et al., 2005]:
where hv (m) is the vent elevation and Ll (m) is the maximum runout distance of the lava flow front.
 For our purposes, we need to know the probability PLij that a lava flow originating from a given location i reaches a site j. We thus assume PLij = 0 if the distance along the flow path from vent i to site j (Lij) is greater than Ll, as defined in equation (1); PLij = 1 when Lij ≤ 0.5Ll, and it increases linearly from 0 to 1 when Ll > Lij > 0.5Ll. Equation (1) gives the maximum runout distance for lava flows erupting at a given elevation. It actually overestimates maximum lava flow lengths for vents at altitudes <2000 m. For example, using equation (1), the maximum length of lava flows erupted at an elevation of 1500 m is 16,700 m. To quantify the impact of length overestimation on the final probability maps, let us assume that the real length of lava flows erupted at 1500 m is 14,000 m instead of 16,700 m. This would produce a variation of only 2% in the probabilities of inundation in the final hazard map.
4. Maps of the Probability of Lava Inundation
 Maps of the probability of inundation by lava flows were compiled for Mount Etna using DOWNFLOW [Favalli et al., 2005]. The maps represent the probability that a given location will be inundated by lava flows in the case of a flank eruption. The inundation probability maps are based on (1) the probability distribution of future vents (PV) and (2) the probability length distribution of future lava flows (PL).
 The probability of inundation (Pj) by lava flows at a generic site j is defined as
where the sum is over all the vents i, PVi is the probability that a vent forms at location i, PSij is the probability that site j is inundated by a lava flow erupted from a vent at location i