PROBABILISTIC HAZARD ASSESSMENT OF BALLISTIC BOMBS FROM PAROXYSMS

PROBABILISTIC HAZARD ASSESSMENT OF BALLISTIC BOMBS FROM PAROXYSMS AND MAJOR EXPLOSIONS AT STROMBOLI Andrea Bevilacqua(1), Patrizia Landi(1), Paola Del Carlo(1), Augusto Neri(1), Antonella Bertagnini(1), Massimo Pompilio(1), Marina Bisson(1), Alessio Di Roberto(1), Willy Aspinall(2), Daniele Andronico(3) (1) Istituto Nazionale di Geofisica e Vulcanologia, Pisa, Italy. (2) University of Bristol, School of Earth Sciences, Bristol, United Kingdom. (3) Istituto Nazionale di Geofisica e Vulcanologia, Catania, Italy.

Our target is the statistical modeling of the hazard of ballistic bombs during major explosions and paroxysms at Stromboli. Overview We target to quantify: • How much is the hourly probability to be hit by a bomb as a function of the spatial location. • How much is the uncertainty affecting these estimates, due to the main epistemic uncertainty sources. On an annual basis, we aim to compare the total cumulative hazard with illustrative lethal annual tolerable risk levels, i. e. 10 -3 for workers and 10 -4 for the general public. (Deligne et al. , 2018 - J Applied Volc) In the following we are detailing a structured technique in five steps: KEY IDEA We decompose the hazard assessment as the product of three factors: 1) Maps of the affected distance range 2) Estimates of the affected sectors angles (direction&width) 3) Maps of the affected region based on three different models. 4) Conditional hazard related to three different values of bombs/m 2 5) Hourly hazard estimates and equivalent time to 10 -3 hazard. In this talk we illustrate a preliminary application to the hazard during the paroxysms. Hazard = P x R x H The probability to be in a region affected hit by a bomb, if a paroxysm occurs. In particular, our temporal estimates rely on the statistics of an updated historical dataset of major explosions and paroxysms after 1875 (Bevilacqua et al. , 2019 - Giornate di Studio su Stromboli; Neri et al. , 2019 - AGU Fall Meeting 2019).

Logic scheme of the hazard model (1) MAP of AFFECTED DISTANCE RANGE Ballistic bombs range [max. distance pdf] (2) AFFECTED SECTOR (direction&width) Affected circular sector(s) sector width [angular measurements] Affected circular sector(s) sector direction [angular measurements] ⊗ Epistemic uncertainty angular buffer uniform pdf 90% CONFIDENCE INTERVAL of the ANGULAR PDF of AFFECTED SECTOR COLOR LEGEND Blue – model input data Black – statistical outputs Red – sources of epistemic uncertainty

Logic scheme of the hazard model (1) MAP of AFFECTED DISTANCE RANGE Ballistic bombs range [max. distance pdf] Affected circular sector(s) sector width [angular measurements] ⊗ (3) MAP of AFFECTED REGION Impact density on the ground [three scenarios] (2) AFFECTED SECTOR (direction&width) Spatial probability of affected region, conditional on a paroxysm PROBABILITY MAP, with 90% CONFIDENCE VALUES Affected circular sector(s) sector direction [angular measurements] ⊗ Epistemic uncertainty angular buffer uniform pdf 90% CONFIDENCE INTERVAL of the ANGULAR PDF of AFFECTED SECTOR COLOR LEGEND Blue – model input data Black – statistical outputs Red – sources of epistemic uncertainty

Logic scheme of the hazard model (1) MAP of AFFECTED DISTANCE RANGE Ballistic bombs range [max. distance pdf] Affected circular sector(s) sector width [angular measurements] ⊗ (3) MAP of AFFECTED REGION Impact density on the ground [three scenarios] ⊗ Spatial probability of affected region, conditional on a paroxysm PROBABILITY MAP, with 90% CONFIDENCE VALUES (4) CONDITIONAL HAZARD ESTIMATES Hazard of ballistic bombs conditional on a paroxysm HAZARD VALUES, with 90% CONFIDENCE (2) AFFECTED SECTOR (direction&width) Time rate of paroxysms [hazard functions] Affected circular sector(s) sector direction [angular measurements] ⊗ Epistemic uncertainty angular buffer uniform pdf 90% CONFIDENCE INTERVAL of the ANGULAR PDF of AFFECTED SECTOR COLOR LEGEND Blue – model input data Black – statistical outputs Red – sources of epistemic uncertainty

Logic scheme of the hazard model (1) MAP of AFFECTED DISTANCE RANGE Ballistic bombs range [max. distance pdf] Affected circular sector(s) sector width [angular measurements] ⊗ (3) MAP of AFFECTED REGION ⊗ Impact density on the ground [three scenarios] Spatial probability of affected region, conditional on a paroxysm PROBABILITY MAP, with 90% CONFIDENCE VALUES (4) CONDITIONAL HAZARD ESTIMATES Hazard of ballistic bombs conditional on a paroxysm HAZARD VALUES, with 90% CONFIDENCE ⊗ (2) AFFECTED SECTOR (direction&width) Time rate of paroxysms [hazard functions] (5) HOURLY HAZARD ESTIMATES Hazard rate of ballistic bombs HAZARD VALUES on SPECIFIC TIMES with 90% CONFIDENCE Affected circular sector(s) sector direction [angular measurements] ⊗ Epistemic uncertainty angular buffer uniform pdf 90% CONFIDENCE INTERVAL of the ANGULAR PDF of AFFECTED SECTOR COLOR LEGEND Blue – model input data Black – statistical outputs Red – sources of epistemic uncertainty

ER BS BL OC KS SP T AT M BO (1) Map of affected distance range of major explosions and paroxysms Figure. (a, b) maximum range reached by ballistic bombs (a) of major explosions in 1996, 1998, 2002, 2009, 2010, and (b) the paroxysms in 1930, 1959, 2003 e 2007. LIT HI C (c, d) probability maps of being within the range of ballistic bombs of (c) major explosions, and (d) paroxysms. UNCERTAINTY RANGE Reported values are probability percentages. This includes a ± 10% buffer. We define a uniform pdf over the range of observed maximum distance d of the bombs from the craters. TH IC BL OC KS We consider the max. distance reached by pluridecimentric ballistic bombs. We do not include smaller bombs and falllout. LI f 0(d) = Unif(270, 1100) m major explosions f 1(d) = Unif(600, 2600) m paroxysms UNCERTAINTY RANGE In the sequel we focus mostly on the paroxysms.

(2 a) Estimates of affected sectors (direction&width) MODEL #0 Figure. Probability density functions of the directions of the ballistic bombs of the paroxysms in: 1930, 1959, 2003 e 2007. Reported values are probabilities over degrees. Peaks of probability are observed in the directions SW and NNE. MODEL #1 MODEL #2 The pdf is obtained through a Monte Carlo simulation of random directions uniformly sampled inside the sectors affected by past paroxysms. Each sample can be seen as a single bomb. This produces a population of angle values: a Gaussian kernel density estimator provides the pdf f 2(·) over the circumference. MODEL #0 is uniform axysimmetric. These results are obtained with ~104 samples. MODEL #1 samples a number of bombs that is proportional with the width of the sector, which means a constant number of bombs/degree. MODEL #2 samples an equal number of bombs in each sector;

(2 b) Epistemic uncertainty (direction&width) MODEL #2 equal number of bombs MODEL #1 equal density of bombs 2007 2003 1959 1930 Ginostra Reported values are probabilities over degrees. We include the 5 th and 95 th percentiles and the mean values. 1959 Figure. Probability density functions of the directions of the ballistic bombs of the paroxysms in 1930, 1959, 2003 e 2007. 1930 San Vincenzo Ginostra San Vincenzo We quantify the effects of two sources of epistemic uncertainty affecting the sectors of past paroxysms: • Unif(-20°, +20°) samples the uncertainty on the direction (i. e. the bisector is rotated) • Unif(0°, +20°) samples the uncertainty on the half width (i. e. the sector is sym. enlarged on both sides). These values are preliminary and displayed as an example. These results are obtained with a hierarchical Monte Carlo of ~102 x 104 samples, i. e. the pdf estimation described in the previous slide is repeated for 102 independent perturbations of direction and width of the sectors affected by each past paroxysm. Probability peaks are reduced and the differences between the models become less prominent.

(3 a) Map of affected region by the bombs during the paroxysms [without epistemic uncertainty] MODEL #0 MODEL #1 (a) (b) We run a second Monte Carlo simulation sampling the pdf f 1(·) of distance range D of the bombs from the craters. Due to the Law of Large Numbers, in each direction θ the radial interval [0, D] is affected with a probability equal to the value of the angular pdf f 2(θ) times the ratio of total angular amplitude α of the paroxysm over 360°. This is equivalent to sample with f 1(·) ⊗ f 2(·) a distance and a direction for a population of bombs affecting a total angular amplitude α. This is performed for each of the three models Model #0, #1, and #2. MODEL #2 (c) Figure. Probability maps of ballistic bombs. (a, b, c) assume different models. Reported values are the probability percentages to be in a region affected by the bombs. Both maps assume an angular amplitude equal to the average of the historically observed angles: α = mean(75°, 85°, 160°, 210°) = 132. 5°. These results are obtained with ~103 samples of the radial distance. The effects of epistemic uncertainty are included in the sequel.

Unif(-10°, +10°) perturbs the direction, and Unif(0°, +10°) the half width of past sectors. constant bombs/degree constant number of bombs i. e. the pdf estimation and the 103 radial distance samples are repeated for 102 independent perturbations of direction and width of the sector affected by each past paroxysm. MODEL #2 These results are obtained with a hierarchical Monte Carlo of ~102 x [103 +104] samples, MODEL #1 (3 b) Epistemic uncertainty (map of affected region) - Example A, perturbations of 10° 5 th percentile 95 th percentile Mean values viewpoint COA helipad near the summit helipad in Ginostra 5 th percentile Mean values 95 th percentile viewpoint COA helipad near the summit helipad in Ginostra

Unif(-20°, +20°) perturbs the direction, and Unif(0°, +20°) the half width of past sectors. constant bombs/degree constant number of bombs i. e. the pdf estimation and the 103 radial distance samples are repeated for 102 independent perturbations of direction and width of the sector affected by each past paroxysm. MODEL #2 These results are obtained with a hierarchical Monte Carlo of ~102 x [103 +104] samples, MODEL #1 (3 c) Epistemic uncertainty (map of affected region) - Example B, perturbations of 20° 5 th percentile 95 th percentile Mean values viewpoint C. O. A. helipad near the summit helipad in Ginostra 5 th percentile Mean values 95 th percentile viewpoint C. O. A. helipad near the summit helipad in Ginostra

(4) Examples of conditional hazard estimates of the bombs during the paroxysms We choose four locations, corresponding to four probability levels to be in the affected region: (a) viewpoint at 400 m elevation a. s. l. on the Regions (a) and (b) trekking path from «Punta Labronzo» Pa ∈ [40%, 90%] (b) the helipad and the shelters near the summit (c) the helipad in Ginostra village, These results are obtained with the uncertainty example B. could be also affected by the bombs during a Pb ∈ [10%, 45%] major explosion. Pc ∈ [15%, 30%] (d) the volcanic observatory in San Vincenzo (C. O. A. ), Pd ∈ [ 5%, 10%] The hazard estimate is the product of Pi x Rj, i. e. the chance to be in a region affected and hit by a bomb. Fieldwork highlighted variable estimates of impact/meter 2 for the bombs of pluridecimetric diameter: - A value of R 1 ≅ 25% has been measured after the paroxysms in 2007 and 2019 at the elipad near the summit We adopt this value as a preliminary high level of impact density during a paroxysm. Pistolesi et al. , 2011 Andronico et al. , 2019 AGU Fall Meeting - A bomb every 5 -10 m, i. e. R 2 ≅ 4%, and R 3 ≅ 1%, have been observed after major explosions in 2009 and 2010. We adopt these values as preliminary medium and low levels of impact density during a paroxysm. Impacts R 1 ≅ 25% bombs/m 2 Impacts R 2 ≅ 4% bombs/m 2 Impacts R 3 ≅ 1% bombs/m 2 Region (a), Pa ∈ [40%, 90%] Pa x R 1 ∈ [10%, 22. 5%] Pa x R 2 ∈ [16‰, 36‰] Pa x R 3 ∈ [4. 0‰, 9. 0‰] Region (b), Pb ∈ [10%, 45%] Pb x R 1 ∈ [ 2. 5%, 11%] Pb x R 2 ∈ [4. 0‰, 18‰] Pb x R 3 ∈ [1. 0‰, 4. 5‰] Region (c), Pc ∈ [15%, 30%] Pc x R 1 ∈ [ 3. 7%, 7. 5%] Pc x R 2 ∈ [6. 0‰, 12‰] Pc x R 3 ∈ [1. 5‰, 3. 0‰] Region (d), Pd ∈ [ 5%, 10%] Pd x R 1 ∈ [ 1. 2%, 2. 4%] Pd x R 2 ∈ [2. 0‰, 4. 0‰] Pd x R 3 ∈ [ 0. 5‰, 1. 0‰] Andronico&Pistolesi, 2010 Gurioli et al. , 2013 These values range approximately between 5 x 10 -4 and 2 x 10 -1.

More details on the historical dataset and the statistical models: Bevilacqua et al. , 2019 [Giornate di Studio su Stromboli] Neri et al. , 2019 [AGU Fall Meeting 2019] (5 a) Hazard functions of major explosions or paroxysms (b) (a) H 1 H 2 In this talk we focus on the paroxysms hazard function Figure. Hazard functions of major explosions and paroxysms. (a) considers all the paroxysms post 1875. In (b) blue and violet bold lines are based on the data post 1985. Thin lines are based on the data post 1875, dashed lines on the interval 1875 -1916. Red lines mark levels 6 months after the last event.

(5 b) Hourly hazard estimates of the bombs during the paroxysms H 1=(3. 4 ± 0. 8) x 10 -5 H 2=( 1. 1 ± 0. 3) x 10 -5 Poisson models renewal models The hourly hazard estimate is the product Hi j k = Pi x Rj x Hk, i. e. the hourly chance to be in a region affected and hit by a bomb, if a paroxysm occurs. Impacts R 1 ≅ 25% bombs/m 2 Impacts R 2 ≅ 4% bombs/m 2 Impacts R 3 ≅ 1% bombs/m 2 Region (a), Pa ∈ [40%, 90%] Ha 1 1 ∈ [26, 94] x 10 -7 Ha 2 1 ∈ [4. 2, 15] x 10 -7 Ha 3 1 ∈ [ 10, 38] x 10 -8 Region (b), Pb ∈ [10%, 45%] Hb 1 1 ∈ [6. 5, 46] x 10 -7 Hb 2 1 ∈ [1. 0, 7. 6] x 10 -7 Hb 3 1 ∈ [2. 6, 19] x 10 -8 Region (c), Pc ∈ [15%, 30%] Hc 1 1 ∈ [9. 6, 31] x 10 -7 Hc 2 1 ∈ [1. 6, 5. 0] x 10 -7 Hc 3 1 ∈ [3. 9, 13] x 10 -8 Region (d), Pd ∈ [ 5%, 10%] Hd 1 1 ∈ [3. 1, 10] x 10 -7 Hd 2 1 ∈ [0. 5, 1. 7] x 10 -7 Hd 3 1 ∈ [1. 3, 4. 2] x 10 -8 Impacts R 1 ≅ 25% bombs/m 2 Impacts R 2 ≅ 4% bombs/m 2 Impacts R 3 ≅ 1% bombs/m 2 Region (a), Pa ∈ [40%, 90%] Ha 1 2 ∈ [8. 0, 31] x 10 -7 Ha 2 2 ∈ [ 13, 50] x 10 -8 Ha 3 2 ∈ [3. 2, 13] x 10 -8 Region (b), Pb ∈ [10%, 45%] Hb 1 2 ∈ [2. 0, 15] x 10 -7 Hb 2 2 ∈ [3. 2, 25] x 10 -8 Hb 3 2 ∈ [0. 8, 6. 3] x 10 -8 Region (c), Pc ∈ [15%, 30%] Hc 1 2 ∈ [3. 0, 10] x 10 -7 Hc 2 2 ∈ [4. 8, 17] x 10 -8 Hc 3 2 ∈ [1. 2, 4. 2] x 10 -8 Region (d), Pd ∈ [ 5%, 10%] Hd 1 2 ∈ [9. 6, 34] x 10 -7 Hd 2 2 ∈ [1. 6, 5. 6] x 10 -8 Hd 3 2 ∈ [0. 4, 1. 4] x 10 -8 These values range approximately between 4 x 10 -9 and 9 x 10 -6. Lethal annual tolerable risk levels are 10 -3 for workers and -4 10 for the general public This enables the calculation of the equivalent time to a cumulative hazard above a specified threshold, for example 1‰, i. e. 10 -3 Region (a), Pa ∈ [40%, 90%], impacts R 1, rate H 1 [110, 390] hours ≈ [4. 4, 16] days Region (b), Pb ∈ [10%, 45%], impacts R 1, rate H 1 [220, 1500] hours ≈ [9. 1, 64] days Region (c), Pc ∈ [15%, 30%], impacts R 2, rate H 2 Region (d), Pd ∈ [ 5%, 10%], impacts R 3, rate H 2 [5900, 21000] hours [71000, 250000] hours ≈ ≈ [245, 870] days [3000, 10000] days ≈ ≈ [8. 0, 29] months [8. 4, 28] years Equivalent times to 10 -4 are 10 times shorter.

• Angular probability peaks are in the directions of Ginostra Village and Punta Labronzo. We included a model for the epistemic uncertainty affecting this. The analysis of the events in 2019 further supports these preferential directions. Conclusions • We considered four different locations on the island at different elevations and in different directions from the craters. In these locations, the conditional hazard of the bombs during the paroxysms ranges between 5 x 10 -4 and 2 x 10 -1. • Illustrative hourly rates of paroxysms after 6 months from the last paroxysm range between 1 x 10 -5 and 4 x 10 -5 depending on the dataset and the statistical models adopted. • Hourly hazard estimates range between 4 x 10 -9 and 9 x 10 -6 depending on the affected region, the bombs/meter 2, and the hazard function model. • On an annual basis, we preliminarily compared the total cumulative hazard with illustrative lethal annual tolerable risk levels, i. e. 10 -3 for workers and 10 -4 for the general public. A total risk of 10 -3 is equivalent to: (a) [110, 390] hours on the N trekking path at 400 m elevation, (b) [220, 1500] hours at the helipad near the summit, both assuming 25% bombs/meter 2 and the renewal models. Regions (a) and (b) are also exposed to the bombs during a major explosion, a hazard not yet included here. (c) [8. 0, 29] months near the helipad in Ginostra village, assuming 4% bombs/meter 2, (d) [8. 4, 28] years at C. O. A. in San Vincenzo village, assuming 1% bombs/meter 2, both assuming Poisson models. The talk does not necessarily represent official views and policies of the Dipartimento della Protezione Civile. On an annual basis, this provides the max. total time that a worker should stay exposed at these levels of hazard. Equivalent times to 10 -4 are 10 times shorter. People spending several months in Ginostra village may exceed their annual tolerable risk level.

Future Work • We are further improving the robustness of the historical dataset, through a critical review of available literature of XIX and XX century. • We are going to detail a similar analysis for the major explosions, which are more common but affecting a less extended region, and generally have a smaller number of bombs/meter 2. • We are going to refine the estimates of bombs/meter 2 and the epistemic uncertainty affecting the direction and width of the sectors of past events, and also the maximum distance of the bombs. • We are going to include more information specifically related to the paroxysms of Summer 2019. • We are going to further explore the differences between renewal models and Poisson models in terms of total cumulative risk over a year-long time interval. The talk does not necessarily represent official views and policies of the Dipartimento della Protezione Civile.