Probabilistic modeling of future volcanic eruptions at Mount Etna

1 Introduction

[2] Mt Etna (Italy) is one of the most active and hazardous volcanoes in the world [Behncke et al., 2005]. It is well known for the persistent activity from the summit craters, frequent lava flow-forming eruptions from vents situated on the flanks of the volcano, and the large population settled on and around the sides of the volcano that is at risk [Guest and Murray, 1979; Bisson et al., 2009]. The ever-expanding use of areas near the volcano increases the potential impact of future eruptions of Mt Etna on the regional economy and on the health and safety of its citizens (Figure 1). Eruptions can neither be prevented nor interrupted, but actions can be taken to minimize damage from them [Ganci et al., 2012b]. Reduction of risk to life and property can be undertaken by avoiding threatened areas and by taking protective measures to reduce the effects when and where vulnerable areas cannot be avoided [Scifoni et al., 2010]. Etna is continuously monitored using seismology, deformation, gravity, magnetism, gas emission studies, and petrology by Istituto Nazionale di Geofisica e Vulcanologia (INGV) and its observatory in Catania [Napoli et al., 2008; Greco et al., 2010; Bonaccorso et al., 2011]. Detection of volcanic precursors can generally identify the locality of impending volcanic activity, even though it often does not determine exactly the typology or timing of an eruption or even its certainty. Hazard-zonation maps can then be used for risk-based decision making in land-use planning and emergency management, as well as addressing the more general problem of quantitative volcanic hazards assessment [Wadge et al., 1994; Newhall, 2000; Sparks, 2003; Magill et al., 2005; Behncke et al., 2005; Favalli et al., 2009; Bisson et al., 2009; Crisci et al., 2010; Cappello et al., 2011a; Vicari et al., 2011a, 2011b]. Thus, effective monitoring of Etna volcano, combined with the preparation of emergency plans to face future eruptions, can help reduce the risk to lives and property in an area densely urbanized, where about 0.8 million people live [Behncke et al., 2005].

[3] Our purpose is to give a proper statistical treatment of the records of volcanic eruptions at Mt Etna to allow key at-risk areas to be rapidly and appropriately identified. In terms of hazard, persistent summit activity does not represent a significant threat to the towns located on Etna's slopes, although important tourist facilities and infrastructures close to the central craters have been repeatedly destroyed. However, periodic flank eruptions pose a serious danger to the lives and the property of the local populace [Behncke et al., 2005], especially when the eruptive vents open up at very low altitudes (up to 200–300 m above sea level). In terms of volcanological data, the catalog of Etna flank eruptions is fully reliable only for the last ~400 years, while summit eruptions are sufficiently documented starting from the first decades of the 20th century [Mulargia et al., 1985; Behncke et al., 2005; Tanguy et al., 2007]. Therefore, the statistical analysis of volcanic activity of Mt Etna was conducted with the dual aim of (1) constructing a spatiotemporal probability map for vent opening of future flank eruptions and (2) forecasting the expected number of eruptive events at the summit craters.

[4] Different approaches have been used to evaluate the probability of vent opening at Etna volcano by analyzing the distribution of spatial location of flank eruptions and their temporal frequency.

[5] For the spatial distribution, the density of vents has been calculated to identify areas where the highest concentration is located, which correspond to the zones known as NE, S, and W Rifts [Guest and Murray, 1979; Acocella and Neri, 2003; Salvi et al., 2006; Crisci et al., 2010; Rongo et al., 2011; Neri et al., 2011; Cappello et al., 2012].

[6] For the temporal distribution of flank eruptions, Mulargia et al. [1985] concluded that Etna's lateral eruptions follow a stationary Poisson process with a constant intensity. Salvi et al. [2006] refuted this result, proving that flank eruptions at Etna do not follow a simple Poisson distribution but are more likely represented by a nonhomogeneous Poisson process with power law intensity. Recently, Bebbington [2007] applied the Hidden Markov Models (HMM) to volcanic occurrences to identify the activity state (hidden) of the volcano most consistent with the observations and to forecast the next eruptions. However, the stationary distribution assumed by HMM cannot properly model the increasing trend of flank eruptions observed over the last 40 years at Etna [Allard et al., 2006]. Finally, the temporal analysis performed by Smethurst et al. [2009] revealed that flank eruptions follow a nonhomogeneous Poisson process with a piecewise intensity increasing nearly linearly since the mid-1900s.

[7] In this paper, we carefully assess the distribution of future vent locations using a probabilistic modeling that combines both the temporal and spatial analysis of flank eruptions in the past four centuries. Spatiotemporal probability maps for future vent opening on the flanks of Mt Etna in the next few years/decades are presented and discussed. We assume that the temporal patterns extracted from the records of past eruptions are useful to predict the future behavior of a volcano; likewise, we consider the regions with a high density of volcanic structures to be the most probable emission zones of future lava flows. Clearly the resulting map does not represent a deterministic forecast of future eruptions but a probabilistic estimate originating from a statistical analysis of the flank eruptive activity of Etna since 1600. Moreover, we show the results of the statistical modeling of the persistent summit activity during the last 110 years, considering Etna's summit craters both separately or as a whole.

2 Etna Volcanic Eruptions: Types, Frequencies, and Structural Features

2.2 Persistent Summit Activity During the Last 110 Years

[15] The summit area of Mt Etna has frequently undergone major morphological changes due to its persistent eruptive activity, both effusive and explosive [Neri et al., 2008]. Four summit craters are active today: the Voragine, the NE crater, the Bocca Nuova, and the SE crater (Figure 2). Their eruptive activity consists of the following: (1) continuous degassing; (2) mild Strombolian explosions accompanied by low rate (~0.1–1 m³ s⁻¹) lava effusion from vents inside or on the flanks of the summit craters; and (3) high-rate paroxysmal (up to 300–400 m³ s⁻¹) eruptions, characterized by lava fountains, lava overflows, and tephra-rich eruption columns several kilometers high.

[16] This morphology has not always been as it is today. Just 100 years ago, there was only one summit crater: the central crater (CC; Guest [1973]). Indeed, during the 17th–18th centuries, the summit of the volcano had a single deep and wide crater depression, which was formed because of the collapse of the summit area triggered by the huge 1669 eruption. Later on, during the 19th–20th centuries, the volcanic activity resumed within this crater, causing its gradual filling.

[17] During the almost persistent activity of the central crater (CC), the NE crater (NEC) formed in 1911. At the beginning, it was a pit located ~3160 m above sea level (asl) on the northern slope of the CC cone [Guest, 1973] and characterized mainly by Strombolian activity. NEC became very active in the second half of the 20th century, producing numerous and occasionally violent and long eruptions (especially in 1947–1960, 1974–1986, and 1995–1998), up to becoming a major pyroclastic cone as well as the highest point of the volcano (3330 m asl in 2007; Neri et al., 2008).

[18] Starting from the middle of the 20th century, the CC changed its internal morphology frequently due to the intracrater volcanic activity alternating with subsidence and collapse phenomena. One of these events led to the creation of the Voragine (VOR; also known as “The Chasm”) in 1945. VOR formed in the northeastern portion of the fairly flat floor of CC [Guest, 1973]. Major eruptive phases occurred at CC in 1960–1964, which culminated in the opening of eruptive fissures, several paroxysmal events accompanied by huge tephra emissions, continuous lava overflows, and the growth of some pyroclastic intracrater cones [Behncke et al., 2004]. Four years later, in 1968, the Bocca Nuova (BN) was formed, still inside the CC, a few tens of meters west of the VOR [Behncke et al., 2003]. BN widened and deepened gradually, always within the CC, becoming even larger than the VOR. Later, in 1998, the Strombolian activity resumed at the crater floor, generating several lava fountaining episodes (in 1999) that filled the crater and triggered numerous lava overflows [Behncke et al., 2003; Calvari et al., 2003]. After the 2001 flank eruption, BN and VOR were involved in subsidence and sporadic phreatomagmatic explosions, and progressively they collapsed, becoming almost a single crater depression whose edges are, in part, coincident with those of the CC [Giammanco et al., 2007]. Therefore, it is evident that VOR and BN are internal, minor vents of the CC itself (see dashed black line in Figure 2).

[19] The SE crater (SEC) formed during the 1971 eruption. Since its creation, SEC has been the most active of the summit craters of the volcano. At first, it was a degassing pit located close to the southeast base of the CC cone, at ~3070 m asl [Behncke et al., 2006]. During the first 25 years of activity, SEC erupted quite frequently and built a cone about 100 m high. The cone has grown dramatically in the periods 1996–2001 and 2006–2012, through over 143 paroxysmal episodes and effusive eruptions. Today it has reached an elevation of ~3300 m, forming also a huge new pyroclastic cone on its eastern slope [Alparone et al., 2005; Behncke et al., 2006; Neri et al., 2008; Vicari et al., 2011b; Ganci et al., 2012a].

[20] The progressive increase in eruptive activity characterizing the Mt Etna summit area during the last 110 years can be readily seen in Table 1. What is most striking is the fact that the appearance of the SEC marks a real change in the eruptive frequency. Over the past 41 years, the SEC has produced a number of eruptive events that correspond to more than double the events produced jointly by the CC and NEC in 110 years. Never in the recent past has Mt Etna been so active; to find similar levels of activity, we would probably have to go back several hundred years, when the historical chronicles, however, were certainly not as accurate as today.

Table 1. Number and Type of Eruptions for Each Summit Cratera

Crater	Paroxysmal Episodes	Strombolian Activities	Lava Flows
Central crater (since 1900)	28	3	16
NE crater (since 1911)	39	10	28
SE crater (since 1971)	154	13	133

• ^a The central crater includes the activity of Voragine and Bocca Nuova. The last column (Lava Flows) indicates the number of paroxysmal and strombolian events that were also characterized by (1) rheomorphic lavas formed during highly fed explosive activity; (2) lava overflows emerging from one of the summit crater; and (3) lavas erupted from short fissures opened along the external slope of the summit cone (above ~3000 m asl).

3 Statistical Analysis of Flank Eruptions

[21] Over the last four centuries, the completeness and accuracy of historical data at Mt Etna provides a robust basis for conducting a statistical analysis to detect spatial and temporal patterns in the distribution of past events and hence to formulate the appropriate spatiotemporal probabilistic model for future flank eruptions.

[22] We analyzed the spatial pattern of eruptive fissures using the Clark and Evans [1954] aggregation index R that measures the ratio between the expected nearest-neighbor distances under complete spatial randomness (CSR) and the observed mean nearest-neighbor distances:

(1)

where d_min(i) is the distance between the ith eruptive fissure and its nearest neighbor, A is the area of the region, and N is the number of structures considered in the calculation. In order to consider the eruptive fissures in full, d_min is calculated as the “minimax distance,” that is, the minimum value of the maximum distances between each end point of the ith eruptive fissure and all the end points of the jth eruptive fissure [Cappello et al., 2012].

[23] R values less than 1 highlight a cluster tendency, since the observed mean distance between neighboring fissures is less than that expected in random patterns. R equal to 1 points out that casualness is present among data, while R greater than 1 indicates uniformity. The theoretical maximum of R is 2.149 occurring when data are maximally dispersed. With historical data at Etna, we obtained an R value equal to 0.027, which proves the nonhomogeneous spatial distribution of eruptive fissures.

[24] In addition to the spatial analysis, we examined the historical records during the last four centuries to detect the presence of trends or patterns characterizing the temporal distribution of eruption occurrences. We calculated the cumulative number of eruptive fissure systems linked to flank eruptions over time and the inter-event periods, i.e., the repose times between subsequent eruptions (Figure 3). Both distributions underline the temporal nonstationarity of eruption occurrences, since a growing trend is evident for the number of flank eruptions, while a mainly decreasing tendency stands out for repose times during the last 30 years. This highlights that volcanic eruptions occurred with a higher frequency in this last period than in the previous ~370 years of activity. To formally prove that the sequence of eruptive events is not a stationary process, we applied the Kolmogorov-Smirnov test statistic D_n [Mulargia et al., 1985], which is based on the maximum difference between the theoretical (uniform) and the empirical cumulative distribution functions. Since the historical records of Etna's flank eruptions provided D_n = 0.267, which is higher than the critical value at a 5% significance level (0.174 = 1.36/√n for a large number n of events), we can reject the null hypothesis of uniform distribution (Figure 3).

[25] In order to evaluate the possible shift during time of the location of the eruptive activity, we also calculated the centroids closer to the eruptive fractures in the last 400 years, grouping them every 50 years (Figure 4). The result indicates that the increasing trend in volcanic activity is accompanied by a migration of the centroids that became progressively closer to the summit of the volcano, on the high eastern side. This fact is testified by the rapid morphostructural evolution of the summit area during the last decades: the former central crater (the single summit crater existing up to early 1900) was modified by the growth of the NE crater (formed in 1911 on the northern side), the Bocca Nuova (formed in 1968 inside the central crater and close to the its western rim, see Figure 2), and the SE crater (formed in 1971 on the southeastern side), which represent the connections of the three main volcanic rifts in the summit area [Neri et al., 2011].

[26] The statistical analysis conducted on flank eruptions in the last four centuries demonstrates the spatial nonhomogeneity and the temporal nonstationarity of flank activity, hence suggesting that the most appropriate space-time probabilistic model for future flank eruptions at Etna volcano is a nonhomogenous Poisson process (NHPP) with a time-varying intensity function.

4 The Space-Time Probabilistic Modeling of Flank Eruptions

[27] Probabilistic modeling provides a powerful tool to estimate the timing and location of future eruptions. Typically, modeling exploits statistical analysis on historical information to calculate probability estimates over a particular area of interest and construct a spatiotemporal map furnishing detailed recurrence rates, i.e., events expected per unit area per unit time.

[28] Here the area of interest on Etna volcano is 1170 km² (rectangular dotted area of 32.5 × 36 km in Figure 1), over which we superimposed a regular grid of potential vents yielding a resolution of 0.5 km. For each potential vent (x, y), the probability of at least a new vent opening in the unit area ΔxΔy (0.5 km × 0.5 km) within time Δt is calculated with the NHPP model as follows:

(2)

with λ_xyt the space-time varying intensity function calculated as

(3)

where d_i is the distance between the point (x, y) and the ith eruptive fissure, h is the smoothing factor, t_i is the elapsed time since the eruption at the ith fissure, N is the total number of fissures, and k is a constant that assigns larger weights to more recent events with respect to older ones.

[29] The smoothing factor h (also called bandwidth or window width) strictly influences the result of equation 3, controlling how λ_xyt varies with distance from historical eruptive fissures. We assigned it a value of 1 km by applying the explicit version of the Least Squares Cross-Validation (LSCV) proposed by Worton [1995], which is based on minimizing the integrated square error between the true and the estimated distribution:

(4)

where N is the total number of structures considered in the calculation and d_ij is the minimax distance.

[30] In order to calculate the constant k, we conducted a retrospective search to quantify how the distribution of the oldest eruptive fissures has affected the formation of the most recent ones. The approach used, called back analysis [Cappello et al., 2011b], allowed us to empirically retrieve the best value of the parameter k using only the known past data. We divided all eruptive fissures according to the age of their formation in two parts and used the first one as a learning data set and the second one as a testing data set. The cutoff time we chose to divide the eruptive fissures is 1981, since it determines two meaningful sets from a statistical point of view: 50 fissures were used for learning and 11 for testing (i.e., the 1983, March 1985, December 1985, 1986–1987, 1989, 1991–1993, 2001, 2002–2003, 2004–2005, 2006, and 2008–2009 eruptions). Moreover, from a volcanological point of view, activity at Etna since 1981 shows a considerably increasing trend that further justified our choice [Behncke and Neri, 2003; Behncke et al., 2005; Allard et al., 2006; Salvi et al., 2006; Neri et al., 2008, 2009; Cappello et al., 2011b].

[31] The best fit for the constant k was retrieved by evaluating the mean square error of each probability forecast with the Brier score [Brier, 1950]:

(5)

where p_i is the forecast probability of each grid cell and o_i indicates whether an eruptive fissure has been formed in that cell after 1981 (o_i = 1 if at least a fissure has occurred, and o_i = 0 if not). The Brier score runs from 0 to 1, with smaller values indicating better forecasts. We found that the value for k that minimizes BS is −0.00107 and used it in equation 3 for all eruptive fissures.

[32] The function λ_xyt is then rescaled so that its integral across the whole area equals the temporal intensity calculated with the power intensity function [Smethurst et al., 2009; and reference therein]:

(6)

where T₁ is the year of the oldest eruptive fissure, t is the year of interest, and δ and θ are unknown positive parameters. The parameter δ determines the shape of the power function: δ > 1 (resp. δ < 1) corresponds to an increasing (resp. decreasing) intensity with time, while δ = 1 reduces equation 6 to a homogenous process. The best estimates of δ and θ (Table 2) are found by minimizing the square differences of residuals between the expected amount (computed as the integral of λ from T₁ to t) and the actual number of eruptions over time (estimated as the cumulative number of events observed up to t) on the entire catalog (from T₁ = 1610 to T₂ = 2012; see Figure 5a). In order to indicate the reliability of the estimates of δ and θ, we also listed the 95% confidence intervals (minimum and maximum estimates in Table 2). Using the obtained values of δ and θ, we extrapolated the curves of annual recurrence rates for different forecasting periods p (1, 10 or 50 years) and calculated the minimum, best, and maximum temporal intensities evaluating equation 6 for t = T₂ + p (Figure 5b). Obviously, the estimated values of δ > 1 (Table 2) are only representative of the increasing phase of the current eruptive cycle. Indeed, it is not physically plausible that the intensity will increase continuously in the next years due to the cyclic variations in output rate and frequency of Etna's eruptive activity [Wadge et al., 1975; Hughes et al., 1990; Behncke and Neri, 2003; Allard et al., 2006].

5 Spatiotemporal Probability Maps for Future Vent Opening

[33] The spatiotemporal probability hazard maps for future vent opening for the three different forecasting periods (1, 10 or 50 years) are reported in Figure 6. The temporal recurrence rates are those obtained with the best estimates of δ and θ (Table 2). The highest value reached by the probability map of vent opening for the coming year is fairly low, at 0.00178. This is bounded by the temporal recurrence rate (0.216), which predicts about one eruption every 5 years. As is evident in Figures 6b and 6c, the probability estimates gradually increase when the forecasting periods are 10 years (maximum value of 0.018) and 50 years (maximum value of 0.091). Nevertheless, the global trend of the probability distributions is kept, showing that the areas closer to the summit of the volcano, above 2500 m above sea level, are where it is very probable that a new eruption will take place. Dissimilarity in the probability estimates can be noticed also in the areas at lower altitude, between 2500 and 1800 m above sea level. Below 1000 m, the expectation of at least a new eruption is very low, except in the South Rift, where the probabilities slowly decrease until the altitude of ~600 m above sea level. This effect is caused by the presence of the 1669 eruptive fissure.

6 Statistical Modeling of Summit Eruptions

[34] The persistent summit activity documented at Etna from the first decades of the 20th century until today represents a good source for a statistical analysis. Thus, here we provide a more comprehensive and accurate study in forecasting volcanic eruptions, estimating also the space-time annual rate of eruptive events at the summit craters (Figure 2): the CC (comprising VOR and BN), the NEC, and the SEC (including the new cone recently formed on its eastern slope).

[35] The data set for the statistical analysis comprises the various types of eruptions occurring since 1900 at the summit craters of Etna, including paroxysmal episodes and Strombolian activities (see Table 1). We conducted two kinds of analysis on this data set. First we calculated the recurrence rate of the summit activity considering the whole catalog of eruptions. Next we divided the data set into three and estimated space-time forecasts for each summit crater (CC, NEC, and SEC).

[36] As for the recurrence rate of the summit activity as a whole (without distinguishing the activity of each summit crater), we estimated the best values for the parameters δ and θ (Table 3) minimizing the square differences of residuals between the expected and the actual number of eruptions over time (Figure 7a). These best estimates were then used in the power intensity function described by equation 6 to obtain the curves of recurrence rates with the expected number of eruptions for the next 1 and 10 years (Figure 7b). Unlike the analysis of flank eruptions, we decided not to use forecasting periods longer than 10 years since they might provide inconsistent results starting from the short time series of documented summit eruptions (~110 years). The best values of recurrence rates for the coming year (12.327 eruptions/year) and the next 10 years (17.030 eruptions/year) estimate that slightly more than one eruption per month is expected in the entire summit area. In the worst cases, more than two (25.530) and about three eruptions (35.268) per month could occur, respectively, in the next 1 and 10 years. These estimates are strongly influenced by the extraordinary activity that took place at SEC in 2000 (66 short-term lava fountains) and recently between 2011 and 2012 (25 midterm lava fountains).

[37] For the space-time forecasts, we divided the catalog in three data sets selecting the eruptions occurring at CC, NEC, and SEC. With these three data sets, we used the same mathematical modeling, calculating the best expected number of eruptions on the basis of the historical catalogs (Figures 8a–8c), from which we derived the best estimates for parameters δ and θ (Table 4). The recurrence rates were again estimated for the coming year and the next 10 years (Figures 8d–8f). The spatiotemporal analysis suggests that the expected number of eruptions at the CC (from 0.461 to 0.484) is comparable with the rates obtained for the flank eruptions (from 0.216 to 0.219). Unlike those of the NEC and of the SEC, the curve of the annual rate for the CC shows a decrease in the speed of growth. At the NEC, the expected number of eruptions reaches annual rates of 1.453 and 1.641, respectively, for the next 1 and 10 years. But the summit crater showing the most exceptional activity is undoubtedly at the SEC. On the basis of the historical eruptions since its formation in 1971, we expected a number of eruptions equal to 14.582 for the coming year, that is, more than one eruption per month, and 22.010 eruptions in the next 10 years, namely slightly less than two events per month.

7 Conclusive Remarks

[38] We have demonstrated a new way to construct the spatiotemporal probability map of vent opening over a specific time period. The spatiotemporal probability map furnishes detailed recurrence rates, estimating the number of expected events per unit area per unit time. In order to find the areas where it is more likely that a new eruption will occur, the best approach consists in analyzing the past eruptive history of a volcano to deduce the possible location of future eruptions.

[39] We based the construction of the spatiotemporal probability map of vent opening on an exhaustive statistical analysis and a NHPP model applied to the Etna flank eruptions of the last 400 years. The spatial distribution of past eruptive fissures is studied, as well as the spatial shift in the effusive eruptions, which appears evident in the last 30 years of Etna activity. An innovation of our approach consists in forecasting the location of new eruptive vents by weighting the most recent events more heavily than the older ones. If a new eruption occurs, the map can simply be updated by merging past information with the timing and location of the new vent. In order to make major changes to the map, a statistically significant number of eruptions would have to occur at different unexpected locations in the study area. The spatiotemporal probability maps of vent opening are quite different in content and meaning from the susceptibility map produced by Cappello et al. [2012], which estimates only the most likely emission zones. Here, we also considered the frequency of flank eruptions to evaluate the future timing and areas of Etna prone to vent opening for the coming year and the next 10 and 50 years. We chose three different forecasting periods to highlight how the probability estimates change over time. Since the catalog of flank eruptions at Etna is well studied and documented only for the last four centuries, forecasting periods longer than 50 years would provide results lacking in significance from a statistical point of view.

[40] Even though the focus of this study is Mt Etna, the approach used to assess the most exposed areas to be affected by new eventual eruptions is meant to be general and applicable to any volcanic region, if available data are significant from a statistical point of view. The spatiotemporal probability map is extremely important in volcanic hazards assessment because the degree of danger presented by a hazardous event strongly depends on the eruption vent location. Hence, the more rigorous the evaluation of the spatiotemporal probability map of vent opening, then the greater the accuracy of the hazard map.

[41] Finally, the results of the analysis of the persistent summit activity during the last 110 years highlight that the hazard rate for eruptive events is not constant with time and differs for each summit crater of Mt Etna. The variations in the recurrence rates at the summit craters indicate that in the last century the main eruptive activity was initially located at the central crater, then migrated to the NE crater and later to the SE crater, as the new summit craters were formed. Most likely, the SE crater will remain the most active crater in the summit area for the near future.