Modelling Cumulative Seismic Damage at the Urban Scale: Calibration and Validation Across Italian Earthquake Sequences

1. Introduction

The analysis of seismic damage scenario in urban areas is an essential tool for risk mitigation, as it investigates the hazard, the vulnerability of the exposed area and the expected damage, enabling the development of effective prevention policies, preparedness strategies and emergency response. A key challenge in this context lies in accounting for seismic sequences, where structures are subjected to repeated shaking which progressively weakens their seismic performance, leading to damage states that result from the accumulation rather than from a single event. This problem has been, and continues to be, researched from multiple perspectives across its numerous components. These include, among others, the development of models to forecast space-time clustering of earthquakes in seismic sequences [1,2]; the formulation of methodologies to analyse the structural features that influence the building’s seismic vulnerability (the Vulnerability Index Method among others [3]), and the external factors that might increase the vulnerability with time [4,5]; as well as models that calculate the expected damage on a structure as a function of its vulnerability and of the seismic input intensity (Damage Probability Matrices [6], vulnerability and fragility curves [7,8,9,10,11], Machine Learning algorithms [12]); and finally modelling the damage accumulation process. On this topic, in physics-based earthquake engineering, various numerical models have analysed the behaviour of different construction typologies subject to repeated shaking [13,14,15,16,17]. Different tools (software) have also been developed to simulate damage scenario ([18] and references therein) and, only more recently, a first software for state-dependent seismic damage scenario was proposed [19], highlighting the growing need to address this aspect in seismic risk assessment.

The numerical simulation of potential future scenarios is inherently linked to the study of observed damage distributions from past seismic events. In fact, in historical seismology considerable effort has been devoted to analysing the effects documented in historical sources in order to differentiate the contribution of individual earthquakes in historical sequences [20]. Recently, a methodology was also proposed to investigate the effect of cumulative seismic intensity on source parameter estimation [21,22]. These aspects are particularly relevant, as the seismic parameters collected in historical catalogues (CFTI5Med [23,24], CPTI [25,26]) form the basis for seismic hazard assessment.

Data from modern and recent seismic damage scenarios are collected through surveys carried out by the Italian Department of Civil Protection together with groups of experts from research centres (e.g. the INGV QUEST [27]). A significant portion of damage data observed in Italy has been made available through the Da.D.O. platform developed by Eucentre [28,29], with the aim to provide researchers with homogenous and standardized information for analysis. In fact, over time different operational inspection tools have been adopted, until the introduction of the AeDES form [30,31]. The datasets in Da.D.O. comprise extensive data on the damage data recorded on large building stocks at the end of major Italian seismic sequence from Friuli 1976 to Mugello 2019 events. However, the information generally refers to the final damage state of the buildings at the end of the sequence, without detailing the eventual intermediate evolution of damage during the sequence.

The first case in Italy where damage was recorded after subsequent earthquakes was during the Central Italy 2016-2017 sequence [32]. These observations have prompted several studies on their implication for damage scenario simulation and macroseismic assessment, as they showed how the large variability in building stock properties, pre-existing conditions, and external factors make both damage surveys and macroseismic intensity evaluations particularly complicated.

In fact, modelling a scenario taking into account all these observations as well as the structural behaviour of each building affected by a sequence would be too computationally expensive. In [33,34] a simplified methodology has been proposed to simulate the effect of seismic sequences over large urban areas considering the accumulation of damage and the evolution of vulnerability. In their work they simulated a seismic sequence by means of an agent-based model, similar to the L’Aquila 2009 sequence, and applied it to the areas of Avola and Catania, Sicily. In [35], building from this methodology, we proposed a different vulnerability evolution function and leveraged observed data on seismic sequences and associated damage scenarios to calibrate the methodology parameters.

In this paper, in Section 2 we first briefly illustrate the methodology proposed in [35] and the novelties proposed in this study. Section 3 presents the results obtained for the calibration study on the 2009 L’Aquila sequence. In Section 4 and Section 5 we show the extension of the results for the Garfagnana-Lunigiana 2013 scenario and the Central Italy 2016-2017 case studies. Finally, in Section 6 we discuss the differences and similarities observed in these three case studies highlighting potentialities and limitations of the methodology.

2. Methodology: Simulations of Seismic Damage Accumulation Scenario

In this study we further investigate an optimized version of the methodology applied in [35] for the simulation of seismic damage accumulation scenarios in urban areas. In the previous study, building on the model originally proposed by [33,34], we simulated numerically a cumulative damage process on an urban environment due to a seismic sequence, according to the following assumptions: (a) each earthquake of a seismic sequence induces in buildings an accumulation of damage which progressively increases their seismic vulnerability; (b) due to the increased vulnerability, the building can reach after subsequent shocks, both with smaller or larger magnitudes, a higher damage level than expected; (c) the damage observed at the end of a sequence is equal to the sum of the damage levels reached at each earthquake of the sequence.

In [35] we applied the methodology to the L’Aquila 2009 seismic sequence. To simulate the sequence, we imported in the model the instrumentally derived macroseismic intensity maps (“ShakeMaps”) generated by the INGV [36,37]. We chose to select, out of all the recorded events, only the eighteen with magnitude greater than Mw 4.0 occurred from the Mw 6.1 mainshock on April 6 to the Mw 5.0 on April 9. As the building stock we used the data from the L’Aquila 2009 survey collected in the Da.D.O. [28,29], both to characterize the initial conditions (position with respect to each ShakeMap, and the initial vulnerability class) and to compare the observed damage levels distribution with the results of our simulations. We focused on masonry buildings as they are the most representative of the observed dataset. We recall that the damage is classified on a discrete categorical scale from D0 (no damage), to D1 (light damage), D2 (moderate damage), D3 (severe damage), D4 (very severe damage) and D5 (collapse); while the initial vulnerability of masonry buildings is categorized into classes from A (most vulnerable) to C1 (least vulnerable), based on the structural features of the vertical and horizontal components and the presence of chains/tie rods/ ring beams, according to Table 1. Additionally, in this paper the vulnerability classes are further partitioned in subclasses based on the combinations of the structural features. As the fragility model to assess the damage after each event and to assign the initial numerical vulnerability value, we employed the following (Eq.1) recently developed model [10], hereinafter “L21”.

$${\mu }_D=\{\begin{array}{c}2.5\left[1+\tanh \left({\displaystyle \frac{I+3.45V-11.7}{0.9+2.8V}}\right)\right]forV\ge 0.32\\ 2.5\left[1+\tanh \left({\displaystyle \frac{I+6.25V-12.6}{1.8}}\right)\right]forV<0.32\end{array}$$

(1)

Schematically, for each earthquake i of a sequence of N events, the damage level ${DL}_i$ is calculated as the most probable damage from the L21 model, as a function of the intensity $I_i$ of the event and of the vulnerability $V_{i-1}$ of the building, that is the updated vulnerability after the previous event: ${DL}_i=f\left(V_{i-1},I_i\right)$,$i=1,\dots N$. The total damage level ${TDL}_i$ is calculated as the sum of all the damage levels up the i-th event included: ${TDL}_i={\displaystyle \sum _{j=1}^i}{DL}_j$. Lastly, the vulnerability is updated as a function of the initial vulnerability $V_0$ and of the total damage level reached after the i-th event: $V_i={g(V}_0,{TDL}_i)$. For the first event, of course, the damage level is calculated as a function of the initial vulnerability ${DL}_1=f\left(V_0,I_1\right)$ and ${TDL}_1=$ ${DL}_1$, assuming that all buildings start from an undamaged state.

As the vulnerability update rule, in the previous study we proposed an exponential function (Eq. 2), hereinafter “Exp”, and a new class of vulnerability A-, with range [1.15,1.25], which all buildings might reach, as they become more vulnerable than their initial state due to the accumulation of damage.

$$\mathrm{br}-\mathrm{to}-\mathrm{break}{V}_i^{EXP}=V_0+(\max \left(A^{-}\right)-V_0){\displaystyle \frac{e^{k({TDL}_i-3)}-e^{-3k}}{e^{2k}-e^{-3k}}},\mathrm{k}=1/V_0$$

(2)

After the last event of the sequence, the final total damage level assigned to each building is converted into a categorical value, from D0 to D5. Moreover, to consistently compare the results of our simulations with the observed data, we considered for each building a global observed damage grade, D_Global, also formulated in [10], which represents the weighted average of the damage levels observed at five structural elements (vertical and horizontal structure, stairs, roofs and infills), according to the survey data available in the Da.D.O. [28,29] .

As the metric of error between the simulated and observed damage distributions, we use the Root Mean Squared Error (RMSE) as follows:

$$RMSE=\sqrt{{\displaystyle \frac{1}{6}}{\displaystyle \sum _{i=0}^5}{\left(F_i^{obs}-F_i^{sim}\right)}^2}$$

(3)

where $F_i^{obs}$ and $F_i^{sim}$ are, respectively, the observed and simulated frequencies for each final damage level, from D0 to D5.

To account for factors not included in the model, we introduced some randomness in the assignment of the initial value of vulnerability $V_0$ to each building, within the ranges in Table 1, and by randomly selecting a fixed percentage of buildings to be damaged for each event. To achieve the latter, each damage scenario was simulated for a number of runs (N.seeds) and, for each earthquake, a percentage p of buildings was randomly chosen for each intensity value of the ShakeMap associated with the event. Therefore, depending on the number of earthquakes in the sequence, on the number of unique values in each earthquake’s ShakeMap and on the number of buildings, a total averaged number of buildings, T, was affected by the sequence in our simulations. All the results were averaged over the number of runs N. seeds =100. The choice to introduce the percentage parameter p was also motivated by the results obtained from the simulation of the mainshock-only scenario, where we observed an overestimation of the light and moderate damage levels D1, D2 and D3, and an underestimation of buildings in D0, D4 and D5. Aiming at finding the value of p corresponding to the minimum Root Mean Square Error (RMSE) between the observed final damage distribution, in terms of D_Global, and the simulated distribution obtained at the end of the sequence, we varied p from 10% to 100% with a step size of 10%. A finer search would have been demanding for memory usage and computational time. For the L’Aquila 2009 case study, we found a minimum value of RMSE (2.52) at p=20%.

Although the preliminary results showed a good agreement between the simulated and observed distribution, we acknowledged that the model presented some limitations that needed to be further investigated in order to improve the methodology’s robustness and generalizability. Among those already discussed in the previous study, in this paper we present the results obtained with an optimized version of the Python script, in terms of memory usage and computational time, that allowed us to further analyse the role of the percentage p of buildings, the number of earthquakes of the sequence, their magnitudes and the number of buildings. Firstly, in order to further investigate the effect of the control parameter p, the optimized version of the code allowed us to refine the minimum RMSE value search for continuous values of p using the minimize_scalar function of the SciPy library (in the bounded method since p varies from 10% to 100%). Secondly, we attempted analysing the relation of p with respect to the intensity. We therefore replaced the constant parameter p with a function p(I), such that the fraction of buildings randomly sampled for each event depends on the value of intensity associated with their position, and tested different functional forms: linear, quadratic, power law and sigmoidal. Lower RMSE values, consistent to those obtained for constant p values, were found using the linear function p(I) = mI. In this work, therefore, we will show directly the results for the optimized values of m. The methodology has been firstly applied to the L’Aquila 2009 sequence, as calibration study. With the aim of investigating also the contribution of foreshocks and to ensure consistency when applying the model to different sequences, we decided to include all the events of magnitude greater than Mw 4.0, occurring within the time window characterized by seismicity rates exceeding the background level. The L’Aquila input sequence, thus, in this work, consists of 5 additional earthquakes with respect to the one in the previous study [35]. Following the analysis of the results for the calibration case, we applied the methodology to two additional seismic sequences, characterized by different number of buildings, earthquakes and values of magnitude, observing similar outcomes.

3. New Results on the L’Aquila 2009 Calibration Case Study

As detailed above, for the L’Aquila 2009 scenario, we included in the sequence selected for our simulations the Mw 4.0 foreshock on March 30 2009 and four aftershocks occurred in the months of April, June and July 2009. In Figure 1 we show the sequence used as seismic input for the L’Aquila 2009 scenario in this work, which consists of a Mw 4.0 foreshock (EN 1), a Mw 6.1 mainshock (EN 2), and of an aftershock sequence of four events of Mw ≥ 5.0 (ENs 12,14,16,19), four events with 4.5 ≤ Mw <5.0 (ENs 3,8,13,20), and 13 events with 4.0 ≤ Mw < 4.5.

For the details on the seismotectonics and for the building stock properties we refer to [35] and references therein.

Firstly, we simulated the L’Aquila scenario varying the value of p from 10% to 100% with a step of 10% and found a minimum RMSE values (2.67) at p=20%. As shown in Figure 2, the results are consistent with the outcome of the simulations run for the 18-event scenario in [35]; however, the minimum RMSE value is larger. Secondly, by applying the optimization procedure, we found a minimum RMSE (0.71) at p=16.06% (red dot in Figure 2). The significant decrease in the RMSE value was reflected in an improved agreement between the observed data and the distribution of final damage levels simulated for p=16.06% and, in particular, with the percentage of buildings in the damage levels D0 and D5.

We then replaced the constant parameter p with the function p (I) = mI, such that the fraction of buildings randomly sampled for each event depends on the value of intensity associated with their position. In this case, the minimum RMSE between simulated and observed distributions is searched as a function of the value m. For the L’Aquila 2009 scenario, we found the minimum RMSE value 0.73 for m=2.37. Compared to the p=16.06% results, no significant difference could be noticed in the distribution of the simulated final damages with respect to the observed data. In Figure 3 we show the averaged (over 100 runs) distribution of the final damage levels obtained for m=2.37.

In this study we also focused on analysing the evolution of the cumulative damage process in our simulations. We analysed the temporal evolution of the averaged (over 100 runs) percentage of buildings in each damage level (D0-D5), after each event of the sequence. As shown in Figure 4, for m=2.37, after the mainshock (EN 2) approximately 20% of the buildings transitions from the undamaged state (D0) to light (D1) and moderate (D2 and D3) damage states, with a slight increase also observed in D4. During the sequence, the percentage of buildings in D1 and D2 increases with a steady upward trend. Notably, the percentage of buildings in D3 exhibits a decrease as the percentages of buildings in higher damage levels (D4 and D5), while remaining limited in absolute terms, display a gradual increase. After the Mw 5.4 aftershock (EN 14), we observe a steep increase in the percentages of buildings in all damage states. An analogous trend is observed for p=16.06%, where only small differences in the percentages can be noted. Although the simulations indicate that the mainshock has the most significant impact on the overall damage evolution, it does not appear to directly generate the highest damage states (D4 and D5). Instead, these levels seem to result from the progressive deterioration of buildings which had reached moderate damage states after the mainshock, due to the cumulative effects of the subsequent events. On the one hand, this outcome might suggest a limitation of our model as it appears to underestimate the expected contribution of the mainshock to severe damage. On the other hand, the observed evolution is consistent with the assumption of cumulative damage: subsequent events, while individually less intense, cumulatively cause the buildings to exceed the thresholds for the highest damage states.

The results obtained in the L’Aquila 2009 case show that, within the investigated range of intensities and for the considered building stock, the introduction of an intensity-dependent parameter has led to an improvement with respect to the preliminary results presented in the previous study [35]; however, on the other hand, the outcomes remain substantially equivalent to those obtained using the optimized constant p value. In fact, given the maximum and minimum intensities from the ShakeMaps associated to the L’Aquila 2009 sequence, I_min = 1.0 and I_max = 8.8, the corresponding percentage of buildings associated with these intensity values for m=2.37 are 2.37% and 20.86%, respectively. Since the spatial distribution of surveyed buildings is significantly concentrated in the epicentral areas of each ShakeMap, where the highest intensity values prevail, the adoption of either a constant or variable p parameter does not produce substantial differences in the results, as in the end the percentage of buildings is in both cases approximately 20%.

4. Case Study 1: the Garfagnana-Lunigiana 2013 Sequence

In order to evaluate the generalizability of the methodology and to analyse potential dependencies on the characteristics of the calibration case study, such as the number of events with Mw ≥ 4.0, the magnitudes of the mainshock and of the strongest foreshocks and aftershocks, we initially chose to apply the model to a sequence characterized by both a smaller number of events and lower magnitudes. This sequence would additionally provide an opportunity to test whether the observed damage distribution, collected in the Da.D.O., actually represents the effect of the mainshock alone.

Among the datasets available in the Da.D.O., the sequences matching these criteria were the Pollino 1998 (a Mw 5.6 mainshock and no other event with Mw ≥ 4.0 [39]), the Emilia 20031 (a Mw 4.7 mainshock followed by an aftershock with Mw ≥ 4.0 [25,26,40]) and the Garfagnana-Lunigiana 2013 (a Mw 5.1 mainshock followed by three aftershocks with Mw≥ 4.0). The latter was eventually chosen as it enabled us to simulate the damage accumulation process on more than two events, and, although one of the main objectives of the Da.D.O. was to provide homogeneous observed damage data, the macroseismic survey for this sequence adopted the same compilation form (AeDES) used in the L’Aquila 2009 sequence.

4.1. The Garfagnana-Lunigiana 2013 Sequence

The sequence shown in Figure 5 consisted of a Mw 5.12 mainshock on June 21 2013 followed by three aftershocks with magnitudes greater than Mw 4.0. The INGV recorded over 2400 earthquakes from June 15 to the end of the following month, elongating from the epicentre of the mainshock towards NE (data from [38]). The mainshock was actually the strongest event recorded in Italy in 2013 [25,26]. Although the mainshock magnitude was moderate, the shaking was felt over a broader area in northern and central Italy. The epicentral area is localized in the northern part of the Tuscany region, across the provinces of Pistoia, Massa-Carrara and Lucca, and it is known historically as the Garfagnana and Lunigiana.

The investigation of the focal mechanism of the 2013 June mainshock, a normal dip-slip faulting on a low (40°-50°) NNW-dipping fault plane with a slight right-lateral strike-slip component [41,42,43], provided insights on the complex seismotectonics of this area of the Northern Apennines, at the northern margin of the Apuan Alps. This region is characterized by an extensional regime related to the opening of intramontane basins and the rollback of the Adriatic plate with normal faults trending NW-SE, generally dipping SW or NW, and often organized in segmented systems connected by transfer zones. Seismicity in this area includes events of both deep (>35 km) and shallow (≤ 35 km) hypocentral depths [37]. The zone between the extensional basins (grabens) of Lunigiana and Garfagnana represents the NW prolongation of the Etrurian Fault System [44]. Most of seismic history occurred within this so-called transfer zone - “a land of earthquakes” (“terra di terremoti” in Italian), where at least one earthquake exceeding the damage threshold (Is > 6 MCS) has occurred approximately every thirty years from 1976 to 2013 [45].

The September 7 1920 (Mw 6.5 and Imax 10 MCS) is the strongest event occurred in the Northern Apennines [23,24] and was decisive for the development and application of seismic prevention policies. The area was first classified as an area of application (Royal Decree-Law n. 1315 of September 23, 1920 [47]) of the regulations in force at that time3: no-build zones were established, maximum building height was restricted, and construction methods and materials were regulated. However, [50] highlights how centralized reconstruction processes are not always well conceived: in this case, the imposed use of reinforced concrete disregarded local building practices. These areas were later classified as seismic zone 2 with the new regulations (Royal Decree n. 431 of March 13, 1927 [51]). Moreover, after the October 1995 event, the regional law of the Tuscany region (n. 56 of July 30, 1997 [52]) was promulgated to promote interventions aimed at risk mitigation and prevention in the municipalities seismically classified. As will be discussed in the following section, the presence of a documented recent seismic history, together with the implementation of risk mitigation measures, may have been effective in preventing more severe damage during the 2013 sequence [45].

4.2. The Macroseismic Survey and the Observed Damage

The QUEST macroseismic group surveyed 27 localities with effects of macroseismic intensity greater than 5 MCS-EMS [45]. The survey focused on external damage and on residential buildings. Most of the localities consist of old medieval historical centres of stone buildings (vulnerability classes A and B), partially abandoned and in a poor state of maintenance, and of new constructions in reinforced concrete and reinforced masonry built after 1920 and 1995. Highest damage levels were reported in the oldest building, while on the new constructions light damage was observed. In the report, there is a note on the effects of the Mw 4.4 aftershock of 23^rd of June which caused a damage grade D4 (roof collapse) in a type A building in the Metra locality, in the municipality of Minucciano.

The QUEST macroseismic survey report [45] is the reference study for the Garfagnana-Lunigiana 2013 dataset in Da.D.O. [28,29]. The original dataset in Da.D.O. consists of 3258 residential buildings, localized in 20 municipalities, in the provinces of Lucca, Pistoia and Massa-Carrara, and it is predominantly composed by masonry constructions (Figure 6). Calculating the Completeness Ratio (CR) [8], compared to the ISTAT 2011 census [46], in this dataset only one municipality, Casola in Lunigiana, was surveyed for more than half its buildings (CR ≈ 74%).

Following the procedure adapted in the L’Aquila 2009 scenario, we applied a preprocessing to the Da.D.O. Garfagnana-Lunigiana 2013 dataset, selecting only masonry buildings for which data on both observed damage and structural characteristics defining the vulnerability class for our analysis were available. After preprocessing, the dataset used in this study consists of 2544 masonry buildings. In Figure 7 we show the spatial distribution of the operative dataset (squares coloured according to the vulnerability class) with respect to the earthquake sequence (stars).

In Figure 8, Figure 9, Figure 10 and Figure 11, we present a descriptive statistical overview of the features of the dataset selected for the Garfagnana 2013 scenario. Differently from the L’Aquila dataset, in this case the data on the year of construction are largely incomplete. As shown in Figure 8, among the masonry buildings considered in this study, approximately 61% were built before 1919, while for the remaining 39% this information is not compiled. Moreover, information on the reconstruction periods is also largely missing. Therefore, in this case, a descriptive analysis of the structural features defining the vulnerability class (i.e. vertical structure, horizontal structure, and presence of chains) across construction or reconstruction periods would not be meaningful. Instead, the discussion focuses on their overall distribution and on their relationship with observed damage. We note that, given the smaller size of this dataset, we chose to retain as many buildings as possible during the preprocessing stage, prioritizing those with the essential information required for this study (i.e. vulnerability features and observed damage), over the exclusion of records lacking additional structural attributes. Figure 9 illustrates the distributions of the three main features defining the initial vulnerability class of the buildings: (a) the quality of the masonry of the vertical structure, (b) the horizontal elements, and (c) the presence of chains. In panel (d) we show the distribution of the subclasses defined according to Table 1. This building stock is dominated by unreinforced masonry walls of low construction quality, mostly lacking seismic tying elements (chains/tie rods or ring beams); the horizontal structural systems are heterogeneous and the most common typologies are floors supported by beams with semi-rigid and deformable slabs, indicating a prevalence of systems with low to moderate in-plane stiffness, which limits seismic load redistribution capacity.

In Figure 10, we illustrate the distribution of the geometrical features of this dataset. Compared to the L’Aquila (see [35] for details), we observe similarity in the overall predominant building profile, although in this case the number of storeys and the average floor heights are slightly higher.

Finally, in Figure 11 we show the distribution of the observed damage levels (D0 – D5) across vulnerability classes. In diagonally hatched bars the damage levels are referred only to the vertical structure, which is the damage assigned according to Da.D.O.; in dotted hatched bars, we used the global damage formula proposed in [10] which averages the damage levels observed to the following structural elements: vertical and horizontal structure, stairs, roofs and infills. A comparison between the observed damage distributions and the EMS-98 damage probability distribution (see in [53]) is noteworthy: the vertical damage distribution is more consistent with the probability distribution associated with intensity I = 7 EMS, whereas the global damage distribution is closer to intensity I = 6 EMS.

4.3. Results

In this section, we present the results obtained for the Garfagnana-Lunigiana 2013 case study, comparing the two formulations of the procedure: first assuming a constant value of p, then considering p as a variable parameter linearly dependent on seismic intensity, p(I).

Following the procedure detailed and applied previously to the calibration case study of L’Aquila 2009, we first simulated the Garfagnana-Lunigiana 2013 scenario varying the parameter p from 10% to 100% with a step of 10% and found a minimum RMSE value of 0.72 for p=20%. Secondly, we applied the optimization process and found that the new minimum RMSE of 0.34 is reached at p=19.06%. In Figure 12 we show the values of RMSE against the p parameter, ranging from 10% to 100% with a step of 10% in black, and the optimized minimum value in red. We observe that the trend is similar to the one observed for the L’Aquila 2009 scenario (Figure 2).

Finally, the optimization procedure, aimed at finding the minimum RMSE between simulated and observed distributions as a function of the value m, found the minimum RMSE value 0.32 for m=3.01. Given the maximum and minimum intensities from the ShakeMaps associated to the Garfagnana-Lunigiana 2013 sequence, I_min = 3.6 and I_max= 7.2, the percentage of buildings associated with these intensity values for m=3.01 are 10.84% and 21.67%, respectively. In Figure 13 we show the distribution of final damage levels simulated for m=3.01. Compared to the L’Aquila 2009 results, in the Garfagnana-Lunigiana 2013 case, the differences between the outcomes obtained for p=20%, p=19.06% and m=3.01, are almost imperceptible, both in the distribution of the final damage levels and in the temporal evolution of the damage. Between p=20% and p=19.06%, we could only notice a slight increase of the percentage of A buildings in D0 associated with a decrease in damage levels D1 and D2.

Analysing the temporal evolution of the averaged (over 100 runs) percentages of buildings in each damage level after each event, as shown in Figure 14, for m=3.01, we observe that, after the mainshock (the first event) less than 20% of the buildings change from D0 to levels D1 and D2. The second event (Mw 4.0) increases the percentage of buildings from D0 to D1 and causes (less than) one building to change from D2 to D3.

The third event (Mw 4.4) cause an increase across the damaged levels up to D4, and finally the fourth event (Mw 4.5) causes (less than) one collapse (D5). We recall that our working assumption is that the observed damage distribution used for comparison corresponds to the final state at the end of the seismic sequence. With regard to the temporal evolution and progressive accumulation of damage, in this case there is evidence in the macroseismic survey of only one building collapse, which occurred after the third event [45]. In comparison, our model captures a worsening of the damage after the third event, while simulates a potential collapse only after the fourth event. Furthermore, the simulations indicate that the mainshock primarily caused slight to moderate damage, whereas the second event is responsible for higher damage grades.

5. Case Study 2: the Central Italy 2016-2017 Sequence

As an additional case study, we chose to apply the methodology on a more complex sequence, with a larger number of events and higher values of magnitude, as well as a more extensive building stock. Among the datasets in Da.D.O. [28,29], the Irpinia 1980 and the Central Italy 2016-2017 were the largest in terms of number of buildings. The latter was selected for this study due to the complexity of the seismic sequence and, similarly to the two previous case studies, because of the adoption of the AeDES form.

5.1. The Central Italy 2016-2017 Sequence

The Central Italy 2016-2017 sequence spanned across two years from August 2016 to December 2017 with a Mw 6.5 mainshock on October 30, 2016, preceded by a foreshock sequence with a Mw 6.0 (August 24), two Mw 5.4 (August 24 and October 26), a Mw 5.9 (October 26) earthquakes, and followed by a long aftershocks sequence with four events of magnitude greater than Mw 5.0. In Figure 15, we show the 72 earthquakes with Mw ≥ 4.0 recorded from August 24, 2016, to December 3, 2017. The August 24 Mw 6.0 caused extensive destruction, leading to 299 fatalities. The Mw 6.5 mainshock is the largest recorded in Italy after the Mw 6.9 Irpinia 1980. Moreover, the Central Italy sequence produced extensive surface ruptures, respectively after the strongest events on August 24 (Mw 6.0) [54] and October 30 (Mw 6.5) [54].

The Central Italy 2016-2017 seismic sequence is located between the source areas of the Umbria-Marche 1997-1998 and the L’Aquila 2009 sequences. The seismic activity in this area is associated with the extensional tectonic regime of the northern and central Apennines, characterized by NW-SE trending normal faults systems, as previously detailed. The 2016-2107 sequence activated a complex network of multiple fault segments extending over 80 km [56,57]. Following the Mw 6.0 shock of August 24, with an epicentre located near to the town of Accumoli, seismicity migrated northwest towards the city of Norcia, culminating in the Mw 6.5 mainshock, and subsequently extended further toward southeast to the Campotosto area in January 2017. These small historical centres, mostly located on steep morphology, are characterized by masonry buildings, which typically feature irregular stone walls with low mechanical properties (e.g. double-leaf walls with weak connection and rubble infill), and upper floors and roofs originally of timber. These localities have been historically affected by destructive events: in 1627, the I₀=7-8 MCS and Me 5.3 event caused severe damage in Accumoli; in 1639, the I₀=9-10 MCS and Me 6.2 event devasted Amatrice; and in 1703 the I₀=11 MCS and Me 6.7 event which destroyed Norcia. Moreover, seismic design regulations have been in place since the early 20^th century (1915 for Amatrice, 1927 for Accumoli, 1984 for Arquata). However, during the 2016-2017 sequence, some buildings exhibited a notable vulnerability. As discussed in [58], the comparison between Norcia and Amatrice following the 2016-2017 sequence led to the concept of “forgotten vulnerability”, whereby the long time elapsed since the last destructive earthquake resulted in a loss of seismic memory, reduced maintenance, retrofitting and spread of construction malpractices (e.g. replacement of the historical timber floors with heavy reinforced concrete slabs without proper strengthening of the masonry walls).

5.2. The Macroseismic Survey and the Observed Data

During the Central Italy sequence, for the first time, observed damage data were collected through survey after each damaging event [32,59]. However, as mentioned in the introduction, even on modern times, with tools and preparation, assessing on site the evolution of damage on large urban areas remains complicated.

First phase: Amatrice, August 2016

The first rapid macroseismic survey was carried in MCS scale [60]. Contextually, the localities were evaluated also in the EMS scale [61,62]. While the MCS survey aimed at rapidly assessing the observed damage scenario for civil protection, the macroseismic assessment in EMS scale allowed to evaluate properly the effect of the vulnerability in the damage distribution. During this first phase, 291 localities were surveyed jointly by DPC-IGAG and INGV-QUEST in MCS scale, 150 in EMS scale4. In the epicentral area the intensity reached I₀= 10 EMS: in the localities of Amatrice and Arquata del Tronto, most buildings with vulnerability class A suffered damage grade D5 (collapse), many reached D3 and D4 levels, while in buildings of reinforced concrete heavy damage on not structural components was observed (except some collapse in Amatrice). The damage distribution elongated toward NNW-SSE in the near field with I=8 EMS. In this phase, it is noted in [61] a potential effect of damage accumulation in San Pellegrino di Norcia, due to Mw 5.4 on August 24. Moreover, the damage distribution was also affected by site effects, as detailed in the previous section.

Second phase: Norcia, October 2016

The October 20 Mw 6.5 mainshock, preceded by the Mw 5.4 and Mw 5.9 on October 26, worsened the damage scenario in the localities affected by the August shocks and damaged the areas northward of Norcia. In the updated QUEST report [63], increases in intensity up to three degrees of the EMS scale were observed. The municipality of Accumoli, for example, was assigned an intensity I=8 EMS after the August 24 shocks (several partial collapses and limited total destruction), while the October 30 event led to 80% of buildings with vulnerability A and 30% of buildings with vulnerability B to near-total destruction, reaching I=10 EMS.

Third phase: Campotosto, January 2017

Finally, in January 2017, four events with Mw ≥ 5.0, caused additional damage southward in the area of Campotosto (from 5 EMS after the October 30 mainshock to 8 EMS) [64,65]5.

The building stock

The dataset in Da.D.O. consists of 89268 buildings, localized across four Italian regions (Marche, Umbria, Lazio and Abruzzo), in 481 municipalities. In Figure 16 the construction typology composition of the Da.D.O. dataset is shown, highlighting that approximately 75% of the buildings are masonry structures.

Following the approach used in the two previous case studies, the Da.D.O. Central Italy 2016-2017 dataset was processed, selecting only masonry buildings with complete data for our analysis. After preprocessing, the dataset used in this study consists of 65164 masonry buildings. In Figure 17 we show the spatial distribution of the operative dataset (squares coloured according to the vulnerability class) with respect to the earthquake sequence (stars).

A descriptive statistical overview of the features of the Central Italy 2016-2017 dataset is illustrated in Figure 18, Figure 19, Figure 20 and Figure 21. Similarly to the L’Aquila 2009 dataset (see [35]), the information on the period of construction and eventual reconstruction of the building stock is well documented. As shown in Figure 18, more than half the masonry buildings were built before 1919; among these, approximately 60% have documented a reconstruction period. Therefore, in this case, the descriptive analysis of the structural features defining the vulnerability class (i.e. vertical structure, horizontal structure, and presence of chains) is carried out across construction periods, as well as with respect to their overall distribution and their relationship with observed damage.

Figure 19 illustrates the distributions of the three main features defining the initial vulnerability class of the buildings: (a) the quality of the masonry of the vertical structure, (b) the presence of chains, and (c) the horizontal elements. In panel (d) we show the distribution of the subclasses defined according to Table 1. This building stock shares comparable structural characteristics with the previous datasets, with a higher percentage of good quality masonry buildings.

The distribution of the geometrical features of this dataset is presented in Figure 20. Compared to the L’Aquila 2009 (see [35] for details) and Garfagnana-Lunigiana 2013, we observe an overall similarity in the most common building profile, although in this case we note a higher percentage of buildings with a small number of storeys (1-2) and with an average floor height in range 2.5-3.5.

Finally, in Figure 21 we show the distribution of the observed damage levels (D0 – D5) across vulnerability classes. In diagonally hatched bars the damage levels are referred only to the vertical structure; in dotted hatched bars, to the global damage. Similarly to the previous scenarios, we observe that when transitioning from vertical to global damage, the distribution of the former reflects a higher intensity level (I=7 EMS), whereas the latter corresponds to a I= 6 EMS.

5.3. Results

In this section, we present and discuss the results of the simulations for the Central Italy 2016-2017 scenario. Consistently with the approach detailed in the calibration and first case study, we first consider p as a constant parameter; then as a linearly dependent function of the intensity.

The results obtained from the 10% step analysis and from the optimization differ even significantly less than the results for both the L’Aquila 2009 and the Garfagnana-Lunigiana 2013 case studies. As illustrated in Figure 22, the 10% step analysis identifies a minimum RMSE (3.53) at p=20%, consistently with the previous cases, while the optimized minimum RMSE (3.51) at p=20.65%.

Finally, the minimum RMSE value of 3.74 was found for m=3.05. Given the maximum and minimum intensities from the ShakeMaps associated to the Central Italy 2016-2017 sequence, I_min = 0.0 and I_max= 9.8, the percentage of buildings associated with these intensity values for m=3.05 are 0% and 29.89%, respectively. In Figure 23 we show the distribution obtained for m=3.05, and in Figure 24 the temporal evolution of the percentages of buildings across damage levels after each event.

Notably, in the Central Italy 2016-2017 scenario, the minimum RMSE values are higher than those obtained in the calibration and first case study. In fact, we observe an underestimation of the percentages of buildings with light damage levels D1 and D2, and an overestimation of the percentage of buildings in D0, and in the moderate to severe damage levels, in particular in D5, similarly to the L’Aquila 2009 case. The analysis of the temporal evolution of the damage percentages provides insights into the potential underlying causes of this larger error. As illustrated in Figure 24, after EN 1, that is the Mw 6.0 foreshock of August 24, 2016, approximately 20% of the buildings experience damage, transitioning from D0 to levels D1, D2, D3, and D4. Although at the scale of the figure is not noticeable, after EN 1, the percentage of building in D4 increases of 0.5%. After EN 1, we observe an increase in the damaged states (and a decrease in the undamaged state) until the Mw 5.4 foreshock (EN 4), where the percentages reach a first plateau. A second steep increase in the damaged levels starts after EN 18 (another Mw 5.4 foreshock) and reaches a second small plateau after the EN 21 (Mw 5.9). Finally, the third steep increase occurs after EN 28 (the Mw 6.5 mainshock), after which the percentages stabilize with minor upward trend and a slight increase after the ENs 55 (Mw 5.1) and 56 (Mw 5.5).

6. Discussion

The aim of this work was to investigate the influence of several parameters in the model proposed in [35], in particular the parameter p, the number of earthquakes, their magnitudes, and the number of buildings. To achieve this, we first improved the implementation of the model, reducing computational costs in terms of time and memory usage; then we further studied the calibration study; finally, we applied the model to two additional sequences. In this section we discuss and compare the outcomes of the three case studies analysed in this work, also with respect to the preliminary results presented in the previous study.

In the L’Aquila 2009 scenario, the effect of the number of earthquakes can be noted, as we showed in the results obtained by changing the number of events in the input sequence, for the same configuration of the model.

We recall that the decision to introduce four additional recorded events to the L’Aquila 2009 18-events sequence was motivated by the necessity to define a criterion to apply in the selection of the earthquakes across the three case studies. In the L’Aquila2009 case, we observed a slightly lower RMSE values in the 18-events, with respect to the 23-event scenario. This might suggest the need to test alternative magnitude thresholds to Mw 4.0, to analyse whether they contribute significantly to the damage accumulation.

The Garfagnana-Lunigiana 2013 simulations, on the other hand, were affected by the lower number of events and buildings. In fact, although the RMSE values in this scenario were the lowest with respect to the other two case studies, in this case we observed the largest standard deviation in the stability plots (not shown in this paper). As shown in Table 2, the average (over the number of runs) total percentage of buildings affected at least once during a sequence, T, is slightly more than 50% of the building dataset, while for the L’Aquila 2009 and the Central Italy 2016-2017 is close to 100%. Moreover, in the Garfagnana-Lunigiana 2013 scenario, we observe an underestimation of the buildings in damage levels D4 and D5. Although also in this case, as discussed in the previous study, the ShakeMap available in the INGV Archive for this sequence are supposedly overestimated (see discussion in [35]), we noted that the magnitude of the mainshock indicated in the metadata of the associated ShakeMap is Mw 5.1, while other catalogues evaluated a Mw 5.3.

The simulations of the Central Italy 2016-2017 scenario were the ones affected by the larger RMSE values. This was likely due to the large number of events with low to moderate magnitude with respect to the magnitude of the strongest events, which may have dominated the observed scenario, leading to saturation effects. Therefore, in our simulations, including the contribution of the smaller events might have led to an overestimation of their cumulative effect. Moreover, as detailed in Section 4.2, according to the QUEST reports, a portion of the buildings in the Central Italy 2016-2017 dataset, particularly those affected by the final phase of the sequence, located in the area of Campotosto, was still undergoing reconstruction after the L’Aquila 2009 event. In our simulations, however, consistently applying the same procedure to all case studies, these buildings were assumed to be in an undamaged state before the sequence. Under this assumption, one would expect an underestimation of the damage in our simulations, compared to the observed data. Instead, the results show the opposite behaviour, with an overestimation in the highest damage levels (D4 and D5), likely linked to the limitations previously listed.

Overall, the refinement of the code and the implementation of a finer search for the optimal constant value p, minimizing the RMSE, have resulted in an improvement, showing consistency with the preliminary results and achieving significantly lower RMSE values for L’Aquila 2009 and Garfagnana-Lunigiana 2013 case studies, whereas more limited improvements are observed for the Central Italy 2016-2017 case study, as discussed above.

Finally, the aim of introducing p as a function of the intensity of the earthquake was to investigate whether this parameter could be related to physical variables of the model. The simplest assumption, also from an implementation perspective, was to consider a linear dependence on the intensity. We recall that, in our simulations, p represents the fraction of buildings randomly sampled – at each event and intensity level - to be affected by the sequence (not necessarily damaged, as damage depends on the vulnerability of the building and the intensity of the area of ShakeMap it is positioned in). Our simulations show that the value of m, which corresponds to the minimum RMSE between simulated and observed data, is consistently around 3 in all three case studies. To be noted, we tested different bounds for m, including the maximum and minimum intensities for each case study, and in all cases the minimum RMSE was found at around 3 (the results presented in this paper refer to bounds [1,10]). Thus, we could choose this value as a fixed parameter of our entire procedure of damage estimation, in order to manage the modelling of the cumulative seismic damage at urban scale for geographical areas with assigned seismic hazard where we do not have any damage dataset available for comparison. In this case, we could plan to simulate synthetic (but realistic) sequences of seismic events and to evaluate the effects on existing buildings through the methodology presented in this paper and calibrated on the three considered case studies. But this, of course, will be subject of future works.

6. Conclusions

The assessment of damage accumulation in urban areas caused by repetitive shocks is a challenging topic, actively investigated for its implication on seismic risk mitigation. In this paper, we studied a simulation methodology able to reproduce the damage distribution observed in buildings after real recorded seismic sequences. Starting from the preliminary results obtained for the L’Aquila 2009 sequence, we further applied an optimized procedure on the Garfagnana-Lunigiana 2013 and Central Italy 2016-2017 sequences. The new results show that, once appropriately calibrated, our methodology can consistently simulate damage scenarios in agreement with observed data, taking into account the damage accumulation process and the associated building seismic vulnerability evolution. This study has also highlighted the effect of the number of earthquakes, their magnitudes and the number of buildings on the simulated scenarios. Future research might focus on analysing the effects of simulated seismic sequences on geographical areas with significant seismic hazard for which no damage datasets are available. The analyses could be developed at either urban or municipality level, and could represent a very useful tool to support recovery and safety plans for vulnerable areas or emergency management activities of Civil Protection departments.