Applying Simulated Seismic Damage Scenarios in the Volcanic Region of Mount Etna (Sicily): A Case-Study From the MW 4.9, 2018 Earthquake

An application for a quick earthquake damage scenario assessment is here presented as a potential tool for planning prevention actions or managing seismic emergencies in the volcanic region of Mt. Etna (Italy). As case-study, we considered the December 26, 2018 earthquake that, with a magnitude MW 4.9, represents the largest event occurring in the area during the last 70 years. The QUEST working group (the INGV macroseismic team) carried out a detailed survey in the damage area, collecting data on the number of buildings in the different vulnerability classes and related damage, with the aim to assign intensity. The maximum intensity reached degree VIII EMS along a narrow strip extending for 5 km astride the Fiandaca fault. In this paper, we simulated the damage scenario in the most struck municipalities of the epicentral area by testing different methodological approaches proposed in the literature using the information of the ISTAT census data collected by the Italian Institute of Statistics. We evaluated the damage level of the residential buildings and we validated the results comparing with the real damage data recognized in the field. Our analysis highlighted the difficulty of applying methods calibrated for larger earthquakes in tectonic domains, to small magnitude events in volcanic zones, where some operating assumptions must be introduced. Despite this, the results confirm the potential of the simulations based on statistical damage assessment methods also in these peculiar conditions, opening the way to finalized plans of pre- and post-earthquake interventions.


INTRODUCTION
Risk scenarios in volcanic areas are mostly referring to damage or disruption caused by lava flows, tephra fallout, or pyroclastic flows, i.e., in general to the eruption effects. This despite damage caused by volcano seismicity, whether or not related to an eruption, is a critical issue in these areas and often represents an under-studied aspect of the risk. In Italy, the analyses carried out at Vesuvius and Campi Flegrei are an example of the proper way to face the problem (Working Group, 2013) and are at the basis of the emergency plans issued by the Italian Department of Civil Protection for the Neapolitan volcanic district. In the study by Zuccaro and De Gregorio (2019), for instance, the damage expected from pre-eruptive seismic activity is evaluated with a uniform seismic input on the whole area, but not considering the characteristics of the attenuation of ground shaking in a volcanic area. Luckily, the occurrence of strong volcano-tectonic earthquakes in the Neapolitan volcanic district is quite rare, and almost limited to the island of Ischia (De Natale et al., 2019;Selva et al., 2019). The situation is different at Mt. Etna, the largest active volcano in Europe, where the contribution of volcanotectonic seismicity plays a more important role and significantly enhances the level of risk, since the earthquakes are frequent and often produce heavy damage .
Here, earthquake damage scenarios have been developed in the framework of the recent EU projects UPStrat-MAFA (Sigbjörnsson et al., 2016) and KnowRISK 1 . In these analyses, both the building vulnerability and the impact on network systems have been also considered Meroni et al., 2016).
While the estimation of damage to residential buildings is a consolidated practice for "tectonic" earthquakes, its application on volcanoes requires ad-hoc approaches taking into account the characteristics of local seismicity, such as the small-moderate magnitude (M max ∼5.3), the shallow depth of sources (<3 km), the high values of peak ground motion (PGM) parameters and the low-frequency content in the near field, and the strong attenuation of intensity in very short distances (Azzaro et al., 2006(Azzaro et al., , 2017Langer et al., 2016;Iervolino, 2018;Tusa et al., 2020). The main feature of the earthquake seismic scenarios at Etna is that the involved areas are relatively small (∼30 km 2 ) and characterized by a large variability of effects in a few kilometers (QUEST Working Group, 2019).
The problem is complicated by the fact that, for very shallow earthquakes, the geometrical spreading of the seismic intensity reflects the geometry of a linear source, i.e., the shape of the shaking area around the causative fault is elongated and not circular (Azzaro et al., 2013). Considering the anisotropic attenuation in earthquake scenarios, as illustrated below, introduces operative problems due to the wide variability of the shaking values also inside the same municipality, difficult to tackle for the statistical significance of vulnerability building data in an inhabited area. In general, data on residential buildings are indeed provided at a municipal scale, unless particular cases where strategical buildings (hospitals, fire stations, schools, etc.) are individually considered. The ISTAT National census collects data on dwellings or building every 10 years and releases them in an aggregated form for each municipality. The accuracy level of census data is an important aspect to obtain reliable earthquake scenarios.
In this paper we discuss some critical points about the damage models in the literature, in the light of the recent analyses for seismic risk at a national scale (Dolce et al., 2020). Notwithstanding damage estimation models have recently been harmonized (i.e., Lagomarsino et al., 2019;Zuccaro et al., 2020) to be integrated into the same platform (IRMA-Italian Risk 1 https://knowriskproject.com/the-project/ Map, Dolce et al., 2020) and to obtain a comparison of results, in the present work, we use three damage models in their original form, based on macroseismic intensity as a input parameter and respecting the philosophy on which the models were created and validated. In this way we avoid resorting to intensity vs. acceleration conversion laws that are affected by significant uncertainties, and considering local amplification effects (where present), already included in the macroseismic parameter. The calibration of damage models on a local scale, as in the case of Mt. Etna, is a crucial point, characterized by the extremely rapid attenuation of intensity, and hence damage effects. An advantage in using macroseismic models is the immediate validation of the results with the data collected directly in terms of intensity.
Finally, the application to December 26, 2018, M W 4.9 Fleri earthquake, the strongest event recorded in the last 70 years at Mt. Etna, allows verifying how damage models set for different seismotectonic contexts may depict the damage level caused by a volcanic earthquakes. To this end, the validation process is based on real damage data collected in a detailed macroseismic survey throughout the damage area .

Macroseismic Survey and Intensity Data
Two days after the intense seismic swarm accompanying the December 24, 2018 eruption on the summit area of Mt. Etna (Alparone et al., 2020), a strong earthquake (M L 4.8, M W 4.9) hit, at 02:19 UTC, the lower southeastern flank of the volcano. Due to an extremely shallow hypocenter, located at a depth of ca. 1 km b.s.l., the event produced heavy damage in the area between the towns of Acireale and Zafferana (Figure 1), with more than 1,100 homeless, but luckily without causing victims (QUEST Working Group, 2019). As many other shallow shocks at Etna, the earthquake was accompanied by remarkable effects of surface faulting along the Fiandaca fault, the southernmost structure of the Timpe tectonic system (Civico et al., 2019;Cucci et al., 2019).
The severity of the effects prompted the QUEST Working Group-the INGV team devoted to the macroseismic survey-to undertake a detailed inspection in the localities of the epicentral area with the aim to assess intensity according to the European Macroseismic Scale EMS (Grünthal, 1998). Since the dense urbanization of the area and the rapid attenuation of seismic intensity moving out from the epicentral area, the survey was carried out, in some key zones, building by building in order to consider properly the variability of effects as well as the building vulnerability and the associated damage. As a result, an intensity map of 44 localities was produced by Azzaro et al. (2020), 24 of them reporting damage. Briefly, the intensity in the epicentral area reached degree 8 EMS-the most damaged zone is between the villages of Fleri and Pennisi along the Fiandaca fault-but the intensity distribution is strongly anisotropic, with a preferential  propagation along the strike of the causative fault and a very strong attenuation in the orthogonal direction (Figure 2).
In the following, we refer to this data for the analysis of the damage scenario produced by the 2018 earthquake, basing particularly on the detailed forms compiled for surveyed localities.

Modeling Intensity Scenarios
As a first step for calibrating the damage scenarios, we simulated the intensity distribution of the December 26, 2018 earthquake by using a probabilistic approach based on the Bayesian statistics (Rotondi and Zonno, 2004). Briefly, the method calculates the decay of the macroseismic intensity conditioned on the epicentral intensity I 0 of the earthquake and the epicenter-site distance; as a result, it provides the probabilistic distribution of the intensity expected (I exp ) at a given site (Zonno et al., 2009). The intensity to be assumed as a reference value is given by the mode (the most frequent value) of the smoothed binomial distribution, whereas the uncertainty can be calculated by setting other probability thresholds (50, 75%, etc.). Being a probabilistic estimation, it does   1919From 1919to 1945From 1946to 1960From 1961to 1971From 1972to 1981From 1982to 1991From 1992 After 2001 1 or 2 from 3 to 5 more than 6

Isolated Contiguous
Good Bad (*) indirect parameter deducted from ISTAT data. not need applying other methods to calculate uncertainties, as for example a Monte Carlo approach.
The above procedure has been adapted to the Etna area by Azzaro et al. (2013) in order to take into account the particular features of attenuation of the volcano-tectonic earthquakes. For example, it is possible to consider both isotropic and anisotropic attenuation models, useful to represent point or linear sources, respectively.
The evident asymmetry in the distribution of the observed intensities (Figure 2) suggested us using the anisotropic attenuation model, in which the preferential direction of propagation of the seismic energy (i.e., minimum attenuation) corresponds to the causative fault, while the maximum attenuation is orthogonal. The decay of intensity along these two trends produces significant differences in terms of calculated intensities, reaching two intensity degrees at equal distances from the epicenter. Concerning the fault model adopted for the simulation, we considered a linear source with a length of 4.5 km, that corresponds to the NW-SE trending segment of the Fiandaca fault (Civico et al., 2019) hosting the instrumental epicenter and the maximum intensity area (I = 8 EMS); the dip of the fault plane is considered vertical.
In Figure 3A the intensity scenario is represented on a grid with a resolution of 0.0025 • (ca. 250 m), where the values of expected intensity I exp represent the mode of the smoothed binomial distribution. To account for the uncertainty affecting the modeled scenario, we also calculated the intensities expected at probability 25, 50, and 75% (I (values commonly used in EMS to indicate, for example, intensity as "6-7, " where the data can be interpreted in the same way as "6" or "7"). In order to estimate the overall reliability of this synthetic scenario, we used the deterministic criterion of validation of the absolute discrepancy diff to compare observed (I) and the calculated (I exp ) intensities at a site: where N is the number of sites. For the whole scenario the result is 0.825, while for the study area (dashed box in Figure 3A) the value is 0.455, i.e., the difference between I and I exp is comparable to the uncertainty often associated with the intensity estimate.
Given the aims of this study and the applications based on the ISTAT census data, we finally plotted the same scenario at an urban scale referring the expected intensities to the census sections (a subdivision of the municipality) considered in the following chapters. In Figure 3B, I exp was calculated at the centroid of each polygon representing the sections. This approach reduces the resolution of the scenario, so that the extension of the areas characterized by a given value of I exp is a bit different from the ones calculated in Figure 3A. Anyway, the result is sufficiently good, with a diff value between synthetic intensities obtained in the centroid of the sections with the observed intensities located inside a polygon, equal to 0.525.

THE ISTAT DATA
To assess the building vulnerability at an urban scale, we use the information for residential housing provided by the ISTAT census data, homogeneously collected every 10 years on the entire national territory. These data allow a reliable estimate of the total number of buildings and their corresponding volume. The number of buildings is published directly by ISTAT (2011), for each municipality and for each census section. On the contrary, in 1991, the number of buildings for the census section has to be inferred from the number of dwellings, using the average value of the item "No. of Dwellings per Building" associated with each record.
Edifices are described by multiple characteristics: structural typology, date of construction (or renovation), number of floors, position in the block, state of repair and quality of maintenance ( Table 1). The last parameter is deduced indirectly from other ISTAT data such as the presence of efficient systems and the characteristics of installations 2 . This data allows a vulnerability classification of buildings when there is no other information collected specifically for the same purpose.
Until the 1991 census, the ISTAT data on residential buildings were provided at the resolution of census section in a disaggregated way, being so possible to correlate the relevant information (as for instance in Meroni et al., 2000).
From the 2001 census onwards, due to the introduction of more restrictive privacy rules, the ISTAT census information provides aggregated values only, reducing the vulnerability assessment to a rough estimation. There are no more disaggregated information at a census section scale neither at a municipal level, and few typological features on age, materials, and the other factors are available in a disaggregated form solely at a provincial level; in practice, they are not usable for detailed vulnerability analyses (Crowley et al., 2009).
Notwithstanding the most recent dataset (ISTAT, 2011) do not provides information suitable for setting the vulnerability model at a local scale, it well describe the recent urban development. Indeed, there are areas in the country with a scarce building development or even affected by the depopulation, for instance rural districts or mountain villages (e.g., in the Apennines). In these cases, data from the 1991 ISTAT census can still be considered representative. On the contrary, other areas of Italy are characterized by a considerable urban growth occurred in the last years, as in the case of the study area. The examined municipalities recorded a significant demographic growth (Figure 4) ranging from 11.4% in Acireale to 43.3% in Viagrande, resulting in a substantial increase in the number of buildings of approximately + 32% from 1991 up to 2011. Moreover, as often happens in case of a rapid urban development, the vulnerability level of the recent buildings is difficult to characterize because of the need to know the age of the buildings (i.e., year of construction) to be referred with the relevant technical rules.
The characterization at a municipal scale of buildings of the last 30 years is therefore a critical step in our analysis. For this reason, different methods for estimating vulnerability are examined in the following.

VULNERABILITY AND DAMAGE ASSESSMENT METHODS
From the methodological point of view, there is in Italy (but also in Europe) a tradition dating back 1980's of studies correlating damage to buildings vs. macroseismic intensity, based on data collected after strong earthquakes. This led structural engineers to derive robust fragility curves or vulnerability functions to be used for estimating local scenarios. Since then, this approach has been continuously improved both in the method and in the input data for calibration (e.g., Spence and Le Brun, 2006), with the increasing use of census data capable of providing better estimations. Several examples are represented by applications in European Risk projects such as in Spain, Portugal, France, Greece, Romania, Turkey, etc. (Lantada et al., 2010;Riedel et al., 2015;Kassaras et al., 2018;Mosoarca et al., 2019).
In general, the scarcity of instrumental ground motion data at a local scale as well as uncertainties associated with GMPEs, explains why in Italy/Europe the use of macroseismic intensity data is normally used for damage scenarios at an urban scale.  Frontiers in Earth Science | www.frontiersin.org Among the methods available in the literature, we adopt the Bernardini et al. (2008, hereinafter M_1) and Giovinazzi and Lagomarsino (2001, hereinafter M_2) which allow to consider the whole information contained in the 1991 ISTAT census by providing a very detailed and precise vulnerability classification of both masonry and reinforced concrete buildings. Indeed, the accuracy of the input data makes the difference with the numerous methods able to group the buildings into vulnerability classes (Kassem et al., 2020). Even though the input data is very accurate, these methods cannot still provide information about the edifices of the last 20 years and need to be integrated (Figure 5).
Another method proposed by Di Pasquale et al. (2005, M_3) is directly applicable to the most recent ISTAT data (2011 census). This empirical approach is based on the relationships between structural types and age classes of the buildings and was widely tested by the technicians of the Department of Civil Protection (Di Pasquale et al., 2005).
For M_1 and M_2 methods, in order to characterize the vulnerability of the edifices built after 1991, we considered the 2011 ISTAT data, available as disaggregated variables at a provincial level only, by adapting them to the municipal scale. Naming P the number of buildings in the 1991-2011 period at provincial level and M the same number for each municipality here considered-Zafferana Etnea, Santa Venerina, Acireale, Aci Sant'Antonio, Viagrande-we calculated the vulnerability distribution of P , and added the normalized ( P / M ) values to the vulnerability distribution of each municipality, previously calculated on the 1991 ISTAT data. We assumed that the recent constructions over the last 20 years P follow the vulnerability distribution assessed at the province level. More details on this method can be found in the application to the 2012 Emilia-Romagna earthquake (Meroni et al., 2017).
With the exception of method M_3, these two methods also consider the year of seismic classification, after which the adoption of more restrictive seismic standards is conceivable. They include an additional parameter, namely the date of seismic classification of the territory, which defines a lower vulnerability for the buildings constructed with earthquake resistant design. Although the study-area was classified since 1914, with a revision in 1962, such an early classification does not guarantee an adequate vulnerability level compared to the present building seismic code.

Method M_1
The approach proposed by Bernardini et al. (2008) defines a score for homogeneous groups of buildings, consistent with a vulnerability assessment (Meroni et al., 2000) calibrated on more than 28.000 detailed buildings vulnerability forms collected during the main seismic crises occurred in Italy from 1983 to 2000 (GNDT database-National Group for the Defense Against Earthquakes).
The 1991 census (ISTAT, 1991) was taken as a starting point, since it provides disaggregated data at a resolution of the census section.
The mean vulnerability index I v , for each group of buildings, is defined by the relation: where i = 1 ÷ 6 and j =1 ÷ 6 (see Table 2B). For a given k structural typology (k = 1 ÷ 5), the factors Delta_i and Delta_j (Table 2A) refer to the i ranges of the construction age (or total retrofitting) of the buildings and to the j typological factors, respectively ( Table 2B). The factors Manut and Classif account for the state of building maintenance and the year of seismic classification of the municipality.
The corresponding EMS vulnerability class is determined according to the range of the vulnerability score shown in Table 2C (Bernardini et al., 2008). Vulnerability classes range from A (the weakest buildings, having the highest indices) to F (the most resistant ones, with the lowest scores).
The method proposed by Bernardini et al. (2007) uses macroseismic fragility curves to evaluate damage to residential   buildings, through the definition of five damage classes (D1 ÷ D5) according to EMS classification. The EMS adopts qualitative ratings to evaluate the frequencies of buildings with different degrees of damage, for each vulnerability class and intensity. For instance, the intensity degree VII is reached when "Many buildings of vulnerability class A suffer damage of grade 3; a few of grade 4. Many buildings of vulnerability class B suffer damage of grade 2; a few of grade 3.A few buildings of vulnerability class C sustain damage of grade 2.A few buildings of vulnerability class D sustain damage of grade 1." (Grünthal, 1998). According to this approach it is possible to estimate a damage grade µ D following the equations: where µ D is the mean value of D, a random variable of damage grade; I is the intensity and I V is the vulnerability index.
The fragility curves for damage distribution P D > d | I are modeled according to a Beta distribution, with a probability density function given by: in which, is the Gamma function, p and q are the parameters of the Beta distribution, defined as a function of the average value µ D and the variance σ D 2 from the relations: Method M_2 Lagomarsino (2001, 2004) proposed a vulnerability assessment method based on EMS by identifying seven distinct categories of buildings. A mean vulnerability index V m (k) is given by the combination of the building age and the structural type: four classes are defined for masonry buildings (k = 1 ÷ 4) and three for reinforced concrete buildings (k = 5 ÷ 7), as shown in Tables 3A, 4A. FIGURE 7 | Damage distribution obtained through the methods M_1 (Bernardini et al., 2007), M_2 (Lagomarsino and Giovinazzi, 2006), and M_3 (Di Pasquale et al., 2005). Black line indicates the average number of buildings in each damage class. The overall vulnerability index V is calculated from the typological score V m (k) and considering appropriate behavior modifiers derived from other ISTAT information (e.g., number of floors, structural context, maintenance status, see Tables 3B,  4B): where k = 1 ÷ 7 refers to building structural typology. The Delta(i,k) for i = 1 ÷ 3 score represents the behavior modifiers for masonry (Table 3B) and reinforced concrete buildings (Table 4B). They can either increase or decrease the initial value of the mean vulnerability index V m (k).
The damage model proposed by Lagomarsino and Giovinazzi (2006) is similar to the one adopted by the M_1 method. It classifies the building stock according to the vulnerability definition of the EMS and provides damage distributions, conditioned by the level of intensity, for each degree of damage of the EMS.
The fragility curves are modeled according to a Beta distribution, with a probability density function with the value µ D calculated as: µ D is the mean damage grade of a random variable D, I is the intensity level and V is the vulnerability index (A ÷ F).
Applications of this method have been developed for Etna area  and Portugal (Zonno et al., 2010;Mota de Sá et al., 2016;Sousa and Campos Costa, 2016).

Method M_3
This method, adopted in the 2005 by the Department of Civil Protection for the assessments of seismic risk in  Italy (Di Pasquale et al., 2005), subdivides the building stock into four vulnerability classes (A, B, C1, and C2) by means of a correlation between type of construction and age. The age classes are those of the ISTAT census, and therefore the method is immediately applicable to the 2011 data. Table 5A lists the vulnerability class as a function of the horizontal and vertical structural elements (adapted from Braga et al., 1982). A correlation between vulnerability class and the age of masonry buildings has been obtained through a statistical study of a sample of about 50,000 dwellings after the 1980 Irpinia M W 6.8 earthquake ( Table 5B). The reinforced concrete buildings are classified into the C2 vulnerability class.
Damage scenarios are estimated by the Damage Probability Matrices (hereinafter DPM), that is a statistical correlation among macroseismic intensity, vulnerability class and EMS damage grade (D0 ÷ D5). According to this approach, damage grade k (k = 1 ÷ 5) is evaluated through the ground shaking (expressed in terms of intensity) and the frequency of buildings in the vulnerability class. The scenario is obtained by adding in discrete terms the number of buildings in each vulnerability class (A, B, C1, and C2) given the intensity degree. The number of buildings is weighted by the probability p k of the damage grade k given by the adopted DPM.
The DPM are described by a binomial distribution: where p k is the damage probability of level k (k = 1 ÷ 5), for a given intensity degree I.
The binomial distribution is defined by the binomial coefficient (or "average damage") d(I), ranging between 0 and 1. Technically, using the only average damage d(I) it is possible to describe the whole damage distribution for each vulnerability class and intensity grade. The values of p k , for each intensity and vulnerability class, have been found through an error minimization procedure on the observed frequencies of damage level (observed DPM), relative to the same sample of structures surveyed after the 1980 Irpinia earthquake (Di Orsini, 1997, 1998). The assumed values of average damage d(I) are summarized in Figure 6.

ANALYSIS OF RESULTS AND COMPARISON WITH SURVEYED DATA
As described before, the intensity attenuation is characterized by a rapid variability at very short distances, and needs to be represented in detail when matched with the vulnerability data to produce reliable damage scenarios. The use of a single intensity value at a municipal scale is a rough approximation, while the ISTAT census section scale is a preferable option. Our analysis followed three main steps: (i) calculating the macroseismic intensity in the census section centroids by means of the specific attenuation model (see section The December 26, 2018 Earthquake: Observed and Simulated Intensity Scenarios), (ii) assessing the building vulnerability distribution for each census section according to the models M_1 to M_3, (iii) evaluating the damage distribution by means of macroseismic fragility curves and DPM.
Following the procedure used by Dolce et al. (2020) for the recent Italian risk maps we examined the estimations obtained by the M_1, M_2, and M_3 methods in terms of damage values that are comparable because are defined through the same EMS scale. Figure 7 shows the damage distribution in five classes (D1 ÷ D5) for the residential buildings of the entire area. It is noteworthy that the methods M_1 and M_2 are more conservative than M_3. Largest differences exist in damage classes D2 and D3, where damage evaluated with M_3 is double than the ones calculated with M_1 and M_2. Furthermore, the M_2 method is the most conservative of all and it does not estimate any collapsed building (class D5).
In order to validate the methodological approaches here proposed, we compare the resulting damage assessments based on the ISTAT data with real data on individual buildings acquired during the macroseismic survey .
The main problem in applying the aforementioned methods is the different areal coverage to which data are referred, complicated by the fact that the study area is densely urbanized, in many places without a solution of continuity. In practice, while at the ISTAT scale the urban settlements are "viewed" through administrative boundaries (i.e., municipal limits and census sections), the extent of the macroseismic survey is focused on the "locality, " that is a territorial unit significant from the statistical point of view of the macroseismic intensity (Grünthal, 1998). In these terms, a locality such as a town or village, typically consists of an historical center and more recent outskirts around them. The localities investigated in the survey are those reported in the national geographic gazetteer used in the Italian Macroseismic DataBase (DBMI15; Locati et al., 2019), which grants the unequivocal association of a locality with a pair of geographical coordinates. In this way all the intensity data available in the DBMI15 for a given geographic reference, are really representative of the seismic history of that locality.
As a result, there is only a partial correspondence between the data acquired on individual buildings during the macroseismic survey and the aggregated ones provided by ISTAT. This situation is shown in Figure 8 where the surveyed areas (in blue) are overlapped on the census sections (in gray).
The macroseismic survey collected data out of 1,278 buildings in the localities damaged by the December 26, 2018 earthquake, distributed in four different municipalities . The focus on the historical centers has the effect of considering mostly the old buildings or the more vulnerable structures in general, whereas the new urbanized zones result "sampled" only for the nearest outskirts, while the residential areas with sparse and isolated buildings are not considered in the macroseismic survey practice. Figure 9A illustrates the case of Fleri, a locality of the municipality of Zafferana Etnea which has been assigned an intensity I = 8 EMS. According to the 2011 ISTAT data, there are 453 buildings in the two census sections (red polygons) corresponding to this locality. During the macroseismic survey carried out in the areas marked by the green polygons, 205 buildings were inspected and classified in terms of vulnerability and damage ( Figure 9B). The comparison between the simulated and surveyed damage distributions is shown in Figure 9C. Considering that no damage is expected in the recently urbanized areas with low vulnerability buildings, the surveyed damage distribution can be considered representative for the entire zone.
In this case, it is evident that the M_3 simulation provides values better approximating the observed data ( Figure 9C). Further analyses on other localities show greater differences between the estimated and observed values, especially for the higher damage classes (D3 ÷ D5).
To quantify the comparison of the simulations with surveyed data, we calculated the index of the dispersion of the results by evaluating a sort of error err j D , defined as a sum of residuals, expressed in terms of numbers of buildings.
The error err j D for each j-th method (M_1 ÷ M_3) is: where n j D is the numbers of estimated buildings for the damage grade D (D0 ÷ D5) for the j-th method. The n 0 D is the number of buildings surveyed during the macrosesmic campaign for the damage grades D.
The minimum value of sum of residuals err j D defines the method better approximating the real damage distribution.
For the final damage scenario, we considered the whole investigated area (gray polygons of Figure 8). In order to perform a comparison with the observed damage, we extrapolated the damage distribution of the 1,278 surveyed buildings to the 4,793 buildings associated to the relevant census sections, assuming that the damage in the investigated area is representative of the overall damage distribution.
The minimum value of err j D is obtained through the M_3 method (Di Pasquale et al., 2005). This method has a err j D values (307) half of the values of M_1 and M_2 (661 and 582, respectively) and therefore the M_3 method is the one better approximating the distribution observed by the survey.
The comparison of the damage distributions for the overall damage scenario is illustrated in Figure 10. The simulations provide damage distributions more conservative than the surveyed data: the damage degree from moderate to heavy damage (D3 ÷ D5) is not properly estimated. The largest differences are obtained for class D3, considering that discrepancies for classes D4 and D5 can be neglected because the sample is too small.
The presence of buildings with a level of vulnerability higher than the average level can explain this difference. Indeed, the survey highlights the presence of many reinforced concrete buildings built in the timespan 1970-90's that are particularly vulnerable. To simulate this condition, we forced the damage scenario by worsening the vulnerability of a part of the reinforced concrete buildings. In detail, 50% of the buildings of class C was classified as B. The downgrade of vulnerability is a practice suggested by the EMS guideline (Grünthal, 1998) in case of bad maintenance conditions or constructive defects, as in the case of reinforced concrete buildings without or with moderate level of resistant design. The new simulated distribution shows a minimum sum of residuals err j D equal to 231 (the minimum value of all simulations), a result confirming the robustness of the method as well as the need to calibrate the vulnerability assessment with local data.
The number of damaged buildings in each census section calculated with the M_3 method is illustrated in Figure 11. The first map shows the value of collapsed or very heavy damaged buildings (D5 + D4): the presence of victims (dead and injured) can be associated with this spatial distribution. The unusable buildings, calculated as the weighted sum of buildings with heavy and substantial damage (D4 + 60%D3) (Dolce et al., 2020), are illustrated in the central map of Figure 11; while the remaining part (D1 + D2 + 40%D3) represents the slightly damaged buildings that can be repaired but imply economic costs (Lucantoni et al., 2001;Di Pasquale et al., 2005).
In general, the spatial distributions of damaged buildings, illustrated in Figure 11, are mainly concentrated along the Fiandaca fault and damage decreases with distance from the fault, reproducing the attenuation effects of the seismic shaking.

CONCLUSIVE REMARKS
In this paper we compared the results of the earthquake damage scenarios obtained by theoretical vulnerability models and fragility curves with the real data collected through a macroseismic survey after a strong shock occurred at Etna volcano of the December 26, 2018 (M W 4.9). We used three models available in the literature, used in the Italian territory and in Europe, to test if suitable in the case of volcano tectonic earthquakes. The application of these models point out a number of issues that has to be solved, the main ones being: 1. the macroseismic data are collected according to the consolidated procedure for assigning the intensity, that is referred to a well-defined locality and cannot be extended to a territory or, worse, to single buildings; 2. the organization of the ISTAT data is critical to perform correct vulnerability and damage assessments, and strongly influences the quality of the final estimates: the recent and updated (ISTAT, 2011) data are available only in an aggregated form, preventing their immediate use in sophisticated damage models.
Furthermore, the perimeter of the census sections in the study area has considerably changed in the last decades, making difficult to track the changes determined by the recent urban growth. Since the localities of the macroseismic survey are unevenly located compared to the administrative boundaries or the census sections, the geographical match with the ISTAT data is not immediate.
The analysis of this case-study highlights the following considerations: 1. the seismic damage assessment methods M_1 (Bernardini et al., 2008) and M_2 (Giovinazzi and Lagomarsino, 2001) can be applied to disaggregated ISTAT data; in Italy, this typology of data (ISTAT, 1991) fixes the situation dating back 1991. In general, these models can be used for settlements where the urban development in the last 30 years has being scarce (for example in the Apennines, Central Italy), but not in areas with recent urban expansion, such as the slopes of Etna. In these cases, the use of approximations for the characterization of the recent urbanized areas is necessary, but may lead to estimation errors; 2. the simpler vulnerability model M_3 (Di Pasquale et al., 2005) can be used with more recent data (ISTAT, 2011). Despite the limited number of ISTAT parameters and a rougher haracterization of the building vulnerability, this method produces reliable estimates when calibrated with real data collected through the macroseismic survey; 3. the considered vulnerability models are calibrated on data mainly collected in Central Italy and do not fully adapt to some typologies of residential buildings in the Etna region. In general, bad quality of materials, construction errors, poor observance of the rules, subsequent structural changes may contribute to increase the effective level of vulnerability, and deserve a correct calibration.
Validation remains an open issue until more detailed data will be available as, for instance, disaggregated updated ISTAT data or the AeDes data (Agibilità e Danno in Emergenza Sismica-"Postearthquake damage and safety assessment and short-term countermeasures") collected, buildings by buildings, in areas struck by strong earthquakes in Italy.
The common used damage estimation methods are calibrated for tectonic earthquakes and validated with data of strong events often presenting seismic sequences with cumulative damage (see the glaring example of the 2016-2017 Central Italy earthquakes in Graziani et al., 2019). In the present work we prove that it is possible to generate reliable damage scenarios in volcanictectonic regions, despite some different characteristics as the small-moderate magnitude, the shallow depth, the high values of peak ground motion parameters and their fast attenuation, and the low frequency content in the near field.
The necessary conditions are that the analysis has to be carried out on a detailed scale (census sections level) and the actual characteristics of the residential buildings have to be accounted for. At the present stage of our investigations, the M_3 method (Di Pasquale et al., 2005) provides a damage scenario better reproducing the effects of the 2018 Etna earthquake.
This application may also contributes to plan measures of intervention for the improvement of the building vulnerability in a densely populated areas such as the Etna volcano particularly exposed to seismic risk.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
VP and FM assessed damage scenarios and compared the results with the collected data. RA and SD'A produced intensity shaking scenario, collected, and analyzed macroseismic survey data. All authors contributed to the article and approved the submitted version.