Defining high-detail hazard maps by a cellular automata approach: application to Mount Etna (Italy)

The individuation of areas that are more likely to be affected by new events in volcanic regions is of fundamental relevance for the mitigation of the possible consequences, both in terms of loss of human life and material properties. Here, we describe a methodology for defining flexible high-detail lava-hazard maps and a technique for the validation of the results obtained. The methodology relies on: (i) an accurate analysis of the past behavior of the volcano; (ii) a new version of the SCIARA model for lava-flow simulation (based on the macroscopic cellular automata paradigm); and (iii) high-performance parallel computing for increasing computational efficiency. The new release of the SCIARA model introduces a Bingham-like rheology as part of the minimization algorithm of the differences for the determination of outflows from a generic cell, and an improved approach to lava cooling. The method is here applied to Mount Etna, the most active volcano in Europe, and applications to landuse planning and hazard mitigation are presented.


Introduction
Many volcanic areas around the World are densely populated and urbanized.In Italy, Mount Etna is home to approximately one million people, even though it is the most active volcano in Europe [Behncke and Neri 2003a].More than half of the events that have occurred in the last four centuries have reported damage to human property in numerous towns on the flanks of the volcano.In particular, the eruptions in 1669 and 1928 destroyed entire villages, and in the first of these, this included portions of Catania, the main city in the region [Chester et al. 1999, Crisci et al. 2003].
Over the last few decades, the vulnerability of the Etnean area has increased exponentially due to continued, sometimes chaotic, urbanization [Dibben 2008], with the consequence that new eruptions might involve even greater risk.In recent effusive crises, countermeasures based on the construction of embankments or digging of channels were adopted to stop or deflect the lava flow [Barberi et al. 1993, Barberi et al. 2003].In some cases, even exceptional actions were necessary in the attempts to stop the lava-front advance, which involved, for instance, the use of explosives and/or the dumping of concrete blocks to divert the lava flow [e.g., Barberi et al. 2003].These actions are generally decided and completed while an eruption is in progress, without any adequate assessment of the outcome, and often exposing the people involved to great risk.In addition, Mount Etna erupts from vents both at its summit and on its flanks [Acocella andNeri 2003, Behncke andNeri 2003b], so that future events might take place almost anywhere on the volcano [Guest and Murray 1979].
Finally, some eruptions unfold very rapidly (e.g., in 1928 and 1981), which leaves little time for any coordinated response [Guest et al. 1987].It is therefore imperative to predict the evolution of future lava fields as early and swiftly as possible, to enable the appropriate civil-defense actions, such as evacuation and other protective measures, and to allow more rational land-use planning.
A modern approach to individuate areas that might be affected by future lava flows is the application of algorithms that allow numerical simulations of either the lava flows [Ishihara et al. 1990, Del Negro et al. 2008] or the simple paths that might be followed by future lava flows [Felpeto et al. 2001, Favalli et al. 2009].For instance, in 2001, the path of the lava flow from an eruption that threatened the town of Nicolosi on Mount Etna was correctly predicted using a lavaflow simulation model [Crisci et al. 2004], which provided useful information to the Civil Defense authorities.However, this modus operandi can be difficult to replicate, and a-priori knowledge of the degree of exposure of the areas surrounding a volcano is desirable, to allow both the realization of preventive countermeasures and more rational land-use planning.
The simulation code used in the present study is the latest release of the Simulation by Cellular Interactive Automata of the Rheology of Aetnean lava flows (SCIARA) model, which has been greatly improved with respect to previous versions [e.g.Crisci et al. 1982].Specifically, a Bingham-like rheology has been introduced for the first time as part of the minimization algorithm of the differences [Di Gregorio and Serra 1999], which is applied for computing lava outflows from the generic cell towards its neighbors.The hexagonal cellular space adopted in the previous releases of the model for the mitigation of the problem of the anisotropic flow direction has been replaced by a Mooreneighborhood-square space, which produces an even better solution for the anisotropic effect.Furthermore, improvements have been introduced to better account for lava cooling.
Based on this new release of the SCIARA computational model, this study illustrates a methodology for the definition of flexible, high-resolution lava-invasion hazard maps applied to Mount Etna.Specific applications for civil-defense purposes and land-use planning are also outlined.

The SCIARA computational model
In general, a computational model for simulating lava flows should be well calibrated and validated against several test cases in order to be fruitfully applied for land-use planning and civil-defense purposes.Furthermore, the model should be as computationally efficient as possible, as a great number of simulations might be required over large areas [D'Ambrosio et al. 2006].Several versions of the SCIARA model, which is based on the cellular automata computational paradigm [Von Neumann 1966], have been released since the first issue [e.g.Crisci et al. 1982, Di Gregorio andSerra 1999].Cellular automata models are intrinsically parallel, so their implementation on parallel computers is straightforward and computational time shortens almost proportionally to the number of available processors [D'Ambrosio and Spataro 2007].
In the current version of the SCIARA model, a Bingham-like rheology has been introduced for the first time, as part of the minimization algorithm of the differences, and this is applied to compute lava outflows from a generic cell.The hexagonal cellular space adopted in the previous releases [Crisci et al. 2004] of the model for mitigation of the problem of the anisotropic flow direction has been replaced by a Moore-neighborhood-square space, which produces better solutions for the anisotropic effect.Furthermore, the new version takes better account of lava cooling.The model was tested with good results on the 2001 and 2006 lava flows on Mount Etna (Italy) (Figure 1), and on an ideal surface (an octagonal-base pyramid; Figure 2b,c), to evaluate the magnitude of the anisotropic effect.After preliminary calibration, the new release has been demonstrated to be more accurate than the previous versions.Tests performed on an ideal surface have indicated that the new SCIARA model solves the typical anisotropic problem of deterministic cellular automata for fluids on ideal surfaces.We briefly outline main specifications of the SCIARA model in the following [see Spataro et al. 2010, for details].

The SCIARA model specifications
As stated above, the new SCIARA model introduces a rheology that was inspired by the Bingham model, and therefore the concepts of critical height and viscosity are explicitly considered [Park and Iversen 1984, Dragoni et al. 1986, Ishihara et al. 1990].In particular, lava can flow out if and only if its thickness overcomes a critical value (i.e. the critical height), so that the basal stress exceeds the yield strength.Moreover, the viscosity is accounted for in terms of the flow relaxation rate, as this represents the parameter for the distribution algorithm that influences the amount of lava that actually leaves the cell, according to a power law of the following kind: log r = a + bT where, T is the lava temperature, and the a and b coefficients are determined by solving this system: Similarly, the critical height mainly depends on the lava temperature, according to a power law of the following kind: the coefficients of which are obtained by solving this system: where, T sol and T vent are the lava temperatures at solidification and at the vents, respectively.
It is known that, in general, deterministic cellular automata for the simulation of macroscopic fluids have a strong dependence on the cell geometry and the directions of the cellular space.To solve the problem, different solutions have been proposed in the literature, such as the adoption of hexagonal cells [Avolio et al. 1998, D'Ambrosio et al. 2003, Crisci et al. 2004], or Monte Carlo approaches [Miyamoto andSasaki 1997, Vicari et al. 2007].The first of these solutions, however, does not perfectly solve the HAZARD MAPS WITH CELLULAR AUTOMATA Here, we approach the anisotropy problem by the introduction of a fictitious topographic alteration along the diagonal neighbors of a generic cell (Figure 2a).In a standard situation of non-altered heights, the cells along diagonals have a lower height with respect to the cells belonging to the von Neumann neighborhood, even in the case of constant slope.This is because the distance between the central cell and the diagonal neighbors is greater than the distance between the central cell and the von Neumann neighborhood cells (Figure 2a).This introduces a side-effect into the distribution algorithm, which operates on the basis of height differences.If the algorithm finds a greater height difference along the diagonals, it will favor them by producing greater outflow.To solve this problem, we consider the height of diagonal neighbors taken at the intersection between the diagonal line and the circle with a radius w centered in the center of the central cell (Figure 2a).Under the commonly assumed hypothesis of an inclined plane between adjacent cells [Vicari et al. 2007], this solution allow there to be constant differences in the levels in correspondence to constant slopes, and the distribution algorithm can work correctly even for the Moore neighborhood.
We tested how this adopted version of SCIARA accounts for the anisotropy problem on an ideal surface, which was represented by an octagonal-base pyramid with faces inclined by an angle a = 5˚.The pyramid is represented by a 10 m cell size digital elevation model (DEM) of 203 columns and 203 rows.By locating the lava source at the top of the structure, flows along both the diagonals and the orthogonal directions of the square cellular space are observed.Figure 2 shows the results of two test cases in which a constant effusion rate, of 1 m 3 s −1 , is emitted for a total of 6 days, and no temperature loss is considered.The first simulation (Figure 2b) is obtained by considering the actual topographic heights of the cells, while the second (Figure 2c) takes into account the topographic corrections.As can be seen, the anisotropic problem is quite significant in the first case, in which, as expected, the diagonal flows reach the base of the pyramid more rapidly with respect to those on the orthogonal directions.In the second case, in which the topographic alterations are considered, the problem is very small, and all of the flows reach the base of the pyramid at the same moment.

The method for lava-flows impact prediction
The methodology for defining high-detail hazard maps proposed here relies on the application of a lava-flow computational model (i.e., the SCIARA model) for the simulation of new events, and on a new criterion for the evaluation of the impact of the simulations performed, in terms of the spatial hazard.
A preliminary analysis of the past behavior of Mount

a b
Etna volcano was carried out to classify the documented effusive events.We used this analysis to set up a database of simulations of possible future lava flows for Mount Etna.A first evaluation of the local hazard is given by the number of simulations that cover each site.However, a similar choice might be misleading, because different events can have different probabilities of occurrence, and these probabilities should be adequately taken into account.In most cases, the probabilities can be inferred from a statistical analysis of past eruptions, which obtains a refined measure of the hazard posed by the lava flows.We show an application of this approach to Mount Etna.

Application to Mount Etna
We adopted a similar procedure to that applied using a previous version of the SCIARA model that was described in Crisci et al. [2010] and was applied to the eastern flank of Mount Etna.Here, we tackled this problem with an elaborate approach for numerical simulation of Etnean lava flows that was based on the results of a high number of simulations of flows that erupted from a grid of hypothetical vents that covered the entire Mount Etna volcano.
On the basis of the documented behavior of the volcano, the probability of a new vent opening on Mount Etna was obtained in the form of a probability density function map [Cappello et al. 2011a], which considers both spatial and temporal components (Figure 3a).The spatial component was estimated through a Gaussian kernel that was applied over the main volcanic structures at Mount Etna [Cappello et al. 2011b], while the temporal rate was evaluated using an approach based on the so-called 'reposetime method' [Ho et al. 1991].
A regular grid of possible future vents was considered as the source for the simulations (Figure 3b).In addition, all of the flank eruptions of Mount Etna since AD 1600 were classified according to their durations and lava volumes [Crisci et al. 2010].The trend of the lava discharge rate during the effusive events is assumed to be as shown in Figure 4, which is in agreement with the recent effusive behavior of Mt Etna [Behncke et al. 2005].
An overall probability of occurrence, p e , was thus defined for each scenario, by consideration of the product of the individual probabilities of its main parameters: (1) where, p s denotes the probability of eruption from a given location, p c the probability related to the membership class of the event, and p t the probability related to the effusion rate trend of the event.Note that 0 ≤ p e ≤ 1, as this results from the product of the independent probability values [Crisci et al. 2010].Formally, if a given DEM cell of co-ordinates x, y was affected by n x,y ≤ N simulations, its hazard was defined as the sum of the probabilities of occurrence of the lava flows involved, p e (i) (i = 1, 2, ..., n x,y ): (2) Figure 3b shows the 500-m-spaced grid of vents that was considered in this study, for a total of 4,290 source locations.A subset of event classes (    the historical events considered by Crisci et al. [2010] and was considered for each crater, thus resulting in a total of 25,740 different simulations.These simulations were performed on an 8-node Apple Xserve Xeon-based cluster, for a total of 80 cores, and they were performed in ca. 10 days.The lava-flow hazard map obtained by applying Equation ( 2) to the database of simulations is shown in Figure 5, and this represents the probability that future eruptions will affect the entire Etnean area.

Validation of the map
For the validation of the hazard map obtained, we proposed the following method [Crisci et al. 2010].We simulated 13 well-documented eruptions on the topography used for all of the simulations presented here (a recent topography).As input data for the selected 13 simulations (e.g.lava discharge rate, duration), we used measurements taken from the existing literature.Afterwards, we evaluated the overlap between the affected areas and the hazard classes RONGO ET AL.

573
Figure 5. High-detail hazard map of Mount Etna from all of the simulations and from the application of Equation ( 2), representing the probabilities of the areas to be affected by future eruptions.identified in the hazard map obtained.It is worth noting that a more reliable validation can be obtained when Etna produces new events.
Initially, the extent of the portion of the i-th hazard class (i = 1, 2, ..., r) that is expected to be affected on average by future lava flows was evaluated as: (3) where, m(C i ) is the measure of the area of the i-th class and where: (4) is the mean probability of the i-th class, and #C i is the number of DEM cells within the i-th class.Moreover, the mean extent of the area affected by the n = 13 simulations, in comparison with the i-th hazard class, was evaluated as: (5) where, S j is the area affected by the j-th simulation.
Finally, the following function was considered: (6) This represents the mean distance between the observed and expected values defined above, where r is the number of hazard classes considered, and: − it is 0 if the two series fit perfectly; − its value becomes greater as the discrepancy increases.
The graphs in Figure 6 show the comparisons between the expected and observed areas affected by the 13 historical lava flows considered, re-simulated on the recent topography [Crisci et al. 2010], with respect to the hazard map classified in 10, 40 and 160 classes of hazard (Figure 6a-c, respectively).Figure 6d shows the mean distance between the observed and expected affected areas, d, with respect to the number of classes, which shows the convergence of d as the number of classes increases.

Applications for civil defense and land-use planning
The availability of a large number of simulations of lava flows with different magnitudes and venting throughout the whole volcano allows the straightforward extraction of selected scenarios on demand.This would be especially relevant once premonitory signals, such as localized seismicity and ground deformation, indicate a possible site of imminent eruption.For this purpose, the SCIARA simulation model has been integrated in a geographic information system application that allows embankments, HAZARD MAPS WITH CELLULAR AUTOMATA ( ) channels digging, barriers, etc., to also be take into account.
In a first civil-defense oriented application, it is possible to identify all of the source areas of lava flows that can affect a given area of interest, such as a town or a major infrastructure.This is obtained by querying the simulation database, by selecting the lava flows that affect the area of interest, and by delimiting their sources [see also Favalli et al. 2009].For the present application, we have chosen the town of Nicolosi, an important historical and cultural site that has many administrative buildings and tourist facilities.Figure 7 shows the vents that might originate eruptions that would affect the urban area of Nicolosi, and the resulting hazard (Figure 7a).A second application [Crisci et al. 2010] considers the effects of human intervention during an eruption, to assess the decrease in the hazard in the critical areas (Figure 7b).
On Mount Etna, emissions of lava flows from extensive eruptive fissures are relatively frequent.We have tested our system with respect to a hypothetical lava-flow emission from a fissure system on the east-northeast flank of Mount Etna, not far from the 1928 eruption site.We considered the emission of lava flow along a fissure system of ca.7 km long, which consisted of three en-echelon segments.The eruptive system was approximated by a subset of vents of the simulation grid, and all of the lava flows that originated from those selected from the simulation database (without the need to perform new simulations) and the resulting map is shown in Figure 8d.
Figure 8 also refers to a further application regarding temporal hazard mapping, by showing the evolution of the involved area with time.This application might be of fundamental importance for the assessment, from a temporal point of view, of how a hazard in a specific area evolves with time (e.g.day by day). Figure 8a-d also shows the inundated areas relative to 1, 3, 6 and 9 days, respectively, with the relative probabilities of occurrences should the fissure system considered be activated.This application constitutes real-time hazard assessment, as these maps show the temporal evolution of the probability of being inundated by lava flows.Similar maps can be used by Civil Protection managers to mitigate the hazard of an imminent/new event without any assumptions as to duration and emission rate of the upcoming event.

Concluding remarks
The macroscopic SCIARA model has here been fruitfully applied to derive a detailed lava-flow hazard map of Mount Etna that can be used for land-use planning and civil-defense purposes.The database of lava flow simulations presented is embedded in a geographic information system environment that allows the retrieval of specific scenarios at any time, and that constitutes a fundamental tool for landuse planning and in quantifying the impact of imminent eruptions and the efficiency of protective measures.
We showed that the latest SCIARA release, which reintroduces a square tessellation of the cellular space instead of the previously adopted hexagonal one, solves the problem of the anisotropic flow direction.This result is particularly significant, as SCIARA is a deterministic model, and as all of the previously proposed solutions refer to probabilistic cellular automata simulation models.This code has been used to simulate the 2001 and 2006 lava flows for Mount Etna (Italy), with high degrees of overlap obtained between the actual and simulated events, and perfect fitting was obtained in terms of the run-out.
The fundamental problem of assessing the impact of future eruptions of Mount Etna lies mostly in the lack of certainty concerning their duration, effusion rates and vent locations.The approach presented in the present study tackles this problem from two directions, by attributing to any given location on Mount Etna a statistical likelihood to (i) give rise to new eruptive vents, and (ii) be impacted upon by lava flows of future eruptions.
As the results obtained are strongly related to the morphology of the study area, each new eruption will require the creation of an updated DEM that incorporates the morphostructural changes that are induced by the eruption [e.g.Tarquini and Favalli 2010].Re-simulation will be necessary only for those events that affect the modified area, and a new, up-dated hazard map can then be obtained by simply reprocessing the new set of simulations.This is a relatively rapid procedure, even on sequential computers.In contrast, if Mount Etna were to produce eruptions for which the characteristics imply a modification of the previously determined representative set of lava flows, a new overall simulation phase would be required to obtain a correct hazard scenario.
Indeed, an apparent limit of our methodology lies in the knowledge that the risk of lava-flow inundation from Mount Etna is extremely low along the low-lying and coastal areas.During the past 1,000 years, lava flows have reached the sea on at least three occasions, in ca.1030 and ca.1160 AD [Tanguy et al. 2007], and in 1669.This is because the development of the lava tubes and ephemeral vents has a determining role in allowing the lava flows to extend to greater distances [Calvari and Pinkerton 1998].Unfortunately, with our current knowledge, the timing and location of the development of such features is not predictable, and therefore this cannot be simulated in the context of this methodology.A solution to this problem would be to simulate, in a manner similar to the approach followed here, the lava emission from a dense network of hypothetical ephemeral vent locations on the lower slopes of Mount Etna.This would provide a complete database for the full area down to the coast, from which the relevant scenarios can be extracted once the positions of the lava tubes and ephemeral vents are known.In general, the overall approach presented here can be applied to other volcanoes where there is also the risk of lava-flow inundation.

Figure 1 .
Figure 1.SCIARA simulations for the 2001 (a) and 2006 (b) lava-flow events at Mount Etna.In both cases the new model reproduced the events with a good level of agreement (area in red).

Figure 2 .
Figure 2. (a) Reference schema for the cell height determination in the Moore neighborhood.The heights of the cells along the von Neumann neighborhood correspond to DEM values.Those along diagonals are taken at the intersection between the diagonal line and the circle with radius w, such that the distance with respect to the centre of the central cell is constant for each adjacent neighbor; (b, c) SCIARA simulations on an octagonal-based pyramid with faces inclined by an angle a = 5˚, performed to evaluate the problem of the anisotropic flow direction when the actual cell topographic elevations are considered (b), and when the topographic corrections along the diagonals are considered (c).Note: the square lattice in panels (b) and (c) is only indicative of the cellular space orientation and does not correspond to the actual cellular space in terms of the number of rows and columns.a b c

RONGO ET AL. 571 Figure 3 .
Figure 3. (a) Characterization of new vents forming in the study region, with the different probabilities of activation, considered in this study.(b) Grid of the hypothetical vents defined as the source for the simulations to be carried out in this study.

Figure 4 .
Figure 4. Example of the representative evolution of the effusion rate during a typical Etnean eruption, for an event lasting 90 days, with a peak effusion rate of ca. 10 m 3 s −1 located at ca. 1/3 of the overall duration, and with ca.1/3 of the maximum value reached at ca. 2/3 of the duration.Note that the trend was obtained by uniting randomly chosen points between two Gaussian curves (dotted/dashed lines).

Figure 6 .
Figure 6.(a-c) Comparison between the observed and expected mean extents of the areas affected by the 13 historical lava flows considered and simulated on the present-day topography, with respect to the lava-invasion hazard map with 10, 40 and 160 classes, respectively; (d) Trend of the mean distance between the observed and expected affected areas, d, with respect to the number of classes, showing considerable convergence for high-density classed maps.The fitting is satisfied starting from d = 0.25 km 2 , which is obtained in correspondence with the 10-classed hazard map.

Figure 7 .Figure 8 .
Figure 7. (a) Map showing the population centers for the densely urbanised Nicolosi on the south flank of Etna, together with the location of a set of vents from which lava flows originated and the hazard zonation based on simulations.The eastern sector of Nicolosi appears to be the most exposed, with a narrow branch of medium-high hazard passing straight through the centre of the town.(b) Map showing the same area and the scenario from lava flows intersecting a hypothetical 2-km-long and 20-m-tall earth barrier to protect the centre of Nicolosi, which are re-simulated on a modified topography that includes the presence of the barrier.The embankment lowers the probability of lava flows inundating the centre of Nicolosi to some extent, as the hazard decreases.a b

Table 1 .
Events classification for the durations/volumes considered for the simulations carried out, as The probabilities of occurrence for each class.