Validation of physics-based ground shaking scenarios for empirical fragility studies: the case of the 2009 L’Aquila earthquake

This paper explores and validates the use of ground shaking scenarios generated via 3D physics-based numerical simulations (PBS) for seismic fragility studies. The 2009 L’Aquila seismic event is selected as case-study application, given the availability of a comprehensive post-earthquake database, gathering observed seismic damages detected on several building typologies representative of the Italian built environment, and of a validated numerical model for the PBS of ground shaking scenarios. Empirical fragility curves are derived as a function of different seismic intensity measures, by taking advantage of an improved statistical technique, overcoming possible uncertainties in the resulting estimates entailed by data aggregation. PBS-based fragility functions are compared to the corresponding sets of curves relying on updated ShakeMaps. The predictive capability of the adopted simulation strategies is then verified in terms of seismic damage scenarios, by respectively coupling PBS- and ShakeMap-based fragility models with the corresponding ground shaking scenarios. Comparison of observed and predicted damage distributions highlights the suitability of PBS for region-specific seismic vulnerability and risk applications.


Introduction
With the ever-increasing computational power, 3D physics-based numerical simulations (PBS) are becoming a more and more appealing tool to provide realistic site-specific scenarios of earthquake ground motion, in alternative to the commonly used empirical ground motion prediction equations (GMPE), based on statistical regressions from regional or 1 3 worldwide records. More specifically, PBS are the key approach towards generation of urban and regional risk scenarios, as it is the case for the ShakeOut (Porter et al. 2011) andHaywired (USGS 2017a, b) experiments in California, as well as for the Scenario Earthquake Shaking Maps (available at https:// www.j-shis. bosai. go. jp) suitable for prefecture emergency plans in Japan, in order to pinpoint target areas and target facilities needing maintenance for earthquake disaster prevention. As pointed out by Hirata (2017), the predicted damage in the Kumamoto prefecture, prior to the earthquake sequence of 2016, was in good agreement with the observed one, although countermeasures were not yet made effective at the time the earthquake occurred.
In spite of their exceptional interest in the context of seismic risk evaluations at urban scale (see e.g., Smerzini and Pitilakis 2018;Stupazzini et al. 2021;Riaño et al. 2021), the applicability of PBS ground motion scenarios is questioned because of their high-frequency limitations, both because of computational limits, that prevent using excessively fine numerical meshes and solve short propagation wavelengths, and because of missing details in the geological model description. To solve these issues, hybrid approaches have been proposed to couple results of low-frequency (LF) PBS with high-frequency (HF) components from stochastic or empirical Green's functions based approaches (Irikura and Miyake 2011). However, since the hybrid approaches have the limitation of disregarding the correlation between LF and HF portions of motion, new approaches have been devised that, based on machine learning techniques (Paolucci et al. 2018;Gatti and Clouteau 2020;Okazaki et al. 2021), aim at solving this issue. Recently, Paolucci et al. (2021) have published the BB-SPEEDset, a dataset including several thousands of simulated nearsource broadband accelerograms from validated case studies, that proved to have statistical properties similar to the NEar-Source Strong-motion records dataset NESS2 (Sgobba et al. 2021).
Validation of PBS results for engineering applications has already found considerable attention (e.g., Galasso et al. 2013;Bijelić et al. 2018;Tsioulou and Galasso 2018;Petrone et al. 2021a, b). However, in view of the seismic risk applications mentioned previously, it is also crucial to verify whether a PBS scenario is suitable to provide not only a reliable prediction of the level of damage observed in an urban area during a historical earthquake, but also a proper basis for the calibration of empirical fragility curves when instrumental information on ground shaking is not sufficient to reliably correlate the observed level of damage to the estimated ground motion intensity.
As a matter of fact, the spatial distribution of ground motion intensity is typically inferred either by ShakeMaps (Worden et al. 2020;Wald et al. 2021), in case a sufficient number of records is available, or by empirical GMPEs (Erdik 2017), as it is often the case for historical earthquakes with no instrumental records. In all such cases, there is a large level of uncertainty when a ground motion level is associated with the specific site where an earthquake effect is observed. Besides, the estimated ground motion is typically available only through its peak values, without information on other parameters related to the time history itself, such as duration and frequency content, and spatial variability of ground motion is neglected. Instead, once validated, the PBS ground motion scenario may provide a complete picture of the variability of the ground motion waveforms, supporting the derivation of empirical fragility curves conditioned on a wider set of intensity measures (IM), including multi-component input.
Motivated by these considerations, the primary goal of this paper is to investigate and validate the use of PBS ground shaking scenarios as input description for the derivation of empirical fragility curves, with application to a well-known case study for Italy, namely, the Mw 6.2, Apr 6 2009, L'Aquila earthquake. This case study is selected because of the availability, on one side, of a database of post-earthquake damage surveys with an unprecedented level of detail (Dolce et al. 2019;Rosti et al. 2021a, b), and, on the other side, of a validated numerical model for the PBS of earthquake ground motion (Smerzini and Villani 2012;Evangelista et al. 2017;Di Michele et al. 2022). Within this primary objective, a secondary objective of the paper is to propose new empirical fragility curves for L'Aquila, which distinguish from already published results owing to: (i) the use of the up-to-date version of the ShakeMap for the L'Aquila earthquake (v4 of Michelini et al. 2020 with respect to v3 of Michelini et al. 2008); (ii) the use of an enhanced statistical approach, which uses the Bernoulli distribution to characterize the random component of the probability model, overcoming possible uncertainty in the resulting fragility estimates driven by data binning; (iii) the adoption of alternative seismic intensity measures for the definition of the explanatory variable of the fragility functions.
This paper is organised as follows. In Sects. 2 and 3, the case study of the L'Aquila earthquake is introduced, with emphasis on the presentation of the detailed damage database used for the empirical fragility study, as well as of the ground shaking scenarios obtained by PBS through the high-performance numerical code SPEED (Mazzieri et al. 2013), in comparison with those provided by the ShakeMaps approach . Then, Sect. 4 presents and compares two different sets of typological fragility curves, obtained on the basis of either the latest version of the ShakeMaps or the PBS of the L'Aquila earthquake. Finally, in Sect. 5, the total building damage observed during the L'Aquila earthquake is compared with the predictions obtained by coupling both the ShakeMaps and the PBS ground motion maps with the corresponding set of fragility curves.

Case study: the 2009 L'Aquila earthquake
On April 6, 2009, a Mw 6.2 earthquake hit the city of L'Aquila, one of the largest urban centers in the Abruzzo region (Central Italy) with about 70,000 inhabitants, causing 308 deaths and vast destruction in the town itself and surrounding areas. The earthquake represents one of the largest shocks in Italy since the beginning of the instrumental age, after the 1980 Mw 6.9 Irpinia, the 2016 Mw 6.5 Norcia, the 1976 Mw 6.4 Friuli earthquake. The earthquake was generated by the seismic rupture of the Paganica segment of Upper Aterno Valley-Paganica fault system consisting of four NW-SE trending normal fault segments (Galadini et al. 2018). The causative fault borders the eastern side of the Aterno River Valley, a geological depression of tectonic origin related to the extensional tectonic stresses acting in Central Italy, as depicted in Fig. 1. Note that the L'Aquila historical center is located on the surface projection of the causative fault inside the Aterno valley, at only about 2 km from the epicenter of the earthquake, suffering therefore from the largest ground shaking owing to the close proximity to the seismic source and to the site amplification effects of the valley sediments.
The post-earthquake macroseismic survey revealed a maximum intensity degree of IX-X in the Mercalli-Cancani-Sieberg (MCS) scale in the towns of Onna and Paganica, while other 14 towns and villages, including L'Aquila, reached an intensity degree between VIII and IX (Galli et al. 2009). The isoseismal lines of maximum intensity showed a pattern strongly elongated in the NW-SE direction, from the epicenter toward the southern part of the Aterno Valley. This pattern reflects the geometry and the dynamics of the fault rupture, with prevailing directivity effects in the direction SE of the epicenter (see discussion in Sect. 3).
After the earthquake, a comprehensive survey campaign, coordinated by the Italian Department of Civil Protection (DPC), was carried out with the aim of collecting damage and usability information on the affected building stock. This information was then conveyed in a damage database (Rosti et al. 2018(Rosti et al. , 2020a(Rosti et al. , 2021a(Rosti et al. , b, 2022, now available in the Observed Damage Database-Da.D.O. (Dolce et al. 2019). This is a web-GIS tool, operated by the Italian DPC, collecting and archiving data on the construction and structural characteristics, as well as on seismic damage, of ordinary buildings inspected after seismic crises of national importance, from 1976 onwards.
The final L'Aquila database, including only data regarding residential buildings in municipalities which were completely surveyed (i.e. with completeness ratio larger than 90%, see discussion in Rosti et al. 2021a, b), consists of about 37,000 buildings, which are shown in Fig. 1. For each inspected building, the database collects information on the geographical coordinates, the construction material (masonry or RC), the typological building classification according to Table 1, the number of storeys (see height classification in Table 2), the construction age, the Damage State (DS) according to the European Macroseismic Scale EMS (Grünthal et al. 1998). Furthermore, as clarified in the following section, the database has been enriched by associating to each damaged building a value of ground motion IM coming from both the ShakeMap and the PBS computations. Referring to Rosti et al. (2018) for further details on the criteria for the attribution of the damage grade, it is important to underline herein that a global damage level is assigned to each building based on the maximum damage observed on the most damaged building component.
For the objective of this study, which is to compare the effect of different approaches for the characterization of ground shaking in the empirical fragility analysis, a subset of the damage dataset, corresponding to the detailed study area in Fig. 1 encompassing the L'Aquila municipality, is considered. An area of limited extent is  Table 1 Typological building classification according to Rosti et al. (2021a, b) In the first column, the percentage of buildings belonging to each class within the detailed study area are reported into brackets Building typology studied mainly because of the computational limits associated with the generation of simulated broadband time histories over a dense grid of sites covering the entire area of the SPEED model. Furthermore, the detailed study area covers the region of larger engineering interest, because of the high density of severely damaged buildings. The considered dataset, hereafter referred to as detailed dataset, counts 7987 residential buildings, 4564 (57%) of which refer to masonry, whereas 3423 (43%) are RC buildings (Fig. 2a). Note that in Table 1 the percentage of buildings belonging to each typological class is specified. Figure 2 shows the spatial distribution and percentage of the different damage levels, from DS0 (no damage) to DS5 (collapse), for both masonry (b) and RC buildings (c). Figure 3 shows the frequency of buildings in each DS according to the different building typological and height classes: buildings are grouped into four main classes-(a) Irregular layout or poor-quality masonry, with flexible or rigid diaphragms, w/wo connecting devices (e.g. tie-rods and/or tie-beams); (b) Regular layout and good-quality masonry, with flexible or rigid diaphragms, w/wo connecting devices; (c) RC seismically designed, pre 1981; (d) RC seismically designed, post 1981-and the corresponding height classes are considered to compute the statistics. The histograms indicate, as expected, that passing from more vulnerable classes, such as irregular masonry buildings, to less vulnerable classes (RC post-1981), there is a decrease of the percentage of higher damage states (DS4-DS5). Furthermore, there is a positive correlation between the percentage of DS4-DS5 with the height of the building, for any building class.
In the municipalities more distant from the earthquake source, post-earthquake surveys tend to be limited to damaged buildings, as they are carried out only upon the owners' request. As a consequence, resulting damage distributions could be biased by an underestimated number of undamaged buildings. To counteract this issue and to suitably account for the negative evidence of damage in sites less affected by the ground shaking, the post-earthquake dataset was integrated by 197′528 non-inspected buildings, located in the non-surveyed and partially-surveyed (with completeness ratio lower than 10%) municipalities, which were reasonably assumed undamaged (e.g. Rosti et al. 2022). The number of undamaged buildings was retrieved from national building census (ISTAT 2001). Allocation of undamaged RC buildings to predefined building types was immediate, as national building census classifies buildings based on construction material, construction age and number of stories. Differently, mapping of undamaged masonry buildings, identified by census parameters, to predefined building typologies was performed via the empirical correlation between census and typological building attributes reported in Rosti et al. (2022). 1 3

3D physics-based simulation of L'Aquila earthquake ground motion
A 3D physics-based numerical approach, through the spectral element code SPEED (Mazzieri et al. 2013, http:// speed. mox. polimi. it/), is used to simulate the seismic wave propagation during the L'Aquila earthquake and, hence, to construct the ground shaking  Figure 4 shows an overview of the SPEED model of the L'Aquila earthquake to highlight its main inputs. The model extends over a volume of 58 km × 58 km × 20 km (see also box in Fig. 1) and it is discretized using an unstructured hexahedral mesh capable To simulate the earthquake rupture, a kinematic finite fault model is adopted using the slip distribution proposed by Ameri et al. (2012), as documented in Evangelista et al. (2017). The seismic velocity model (see inset in Fig. 4), calibrated based on the available seismic microzonation data, is constructed by combing a horizontally layered crustal model and a 3D model of the Aterno valley. With respect to Evangelista et al. (2017), where the outcropping bedrock layer consists of very hard rock with shear wave velocity V S = 1700 m/s, the updated crustal model features a first layer (see Crust 1 in Fig. 4) with a parabolic gradient of V S from a minimum value of 800 m/s up to a value of 1700 m/s at 1000 m depth from the topographical surface. The modified velocity profile implies a V S30 value of 900 m/s, which is more consistent with the average distribution of V S30 values for the rock formations of Central Italy (Forte et al. 2019). In Fig. 4, a representative cross-section of the velocity model (see yellow line A-A′ superimposed on the mesh), passing through the L'Aquila downtown, is also shown to highlight the complexity of the wave propagation model. Note that ground topography is also accounted for in the SPEED simulation.
To overcome the frequency limit of the numerical model, the ANN2BB technique proposed by Paolucci et al. (2018) and further improved in Paolucci et al. (2021), is used to enrich the PBS signals at high frequency. The essence of this technique is the use of an Artificial Neural Network (ANN), trained on a strong-motion recorded dataset, to predict short-period spectral ordinates starting from the long-period spectral ordinates of the SPEED signal. The procedure allows for the generation of broadband ground motions which preserve the features of the SPEED signals at low frequencies and, at high frequency, provide a realistic (i.e. consistent with records) distribution of several IMs , as well as of the spatial variability . The last feature is particularly relevant when generating ground shaking scenarios for seismic risk studies, because it ensures that realistic features of spatial correlation and coherency are maintained at both long and short periods.
Referring to Evangelista et al. (2017) for a thorough comparison of the simulated ground motions with the available recordings, which represents a seminal step for PBS validation, albeit beyond the scope of this work, we focus herein on the analysis of the physics-based ground shaking scenario, in comparison with the one obtained from standard empirical approaches. To this end, it is relevant to compare the PBS scenario with the ShakeMap available for L'Aquila earthquake, as the latter represents the stateof-the-art approach for computing the input motion to estimate damage scenarios for historical earthquakes and for deriving empirical fragility curves within the national seismic risk assessment for Italy (Dolce et al. 2021;Borzi et al. 2021;Faenza et al. 2020). Such comparison aims at shedding light on analogies and differences between the two approaches, which may have an impact on the estimation of the fragility curves presented in the following section. Figure 5 shows the map of Peak Ground Velocity (PGV) obtained from PBS (left) and from the up-to-date version (v4) of the ShakeMap according to Michelini et al. (2020) (right), in comparison with the values recorded at the available strong-motion stations. The maximum horizontal component (Hmax) is shown to enable a consistent comparison between the two approaches, since ShakeMaps are released only for the Hmax component. To understand the differences between PBS and ShakeMap, it is worth recalling that the former is derived by computing the peak values directly from the waveforms simulated by SPEED on an arbitrarily dense grid of receivers, while the latter is generated by combining, through suitable geospatial interpolation algorithms, the ground motion values recorded at the available stations with the GMPE of Bindi et al. (2011), where data are not available. The Bindi et al. (2011) is the ground motion model selected for shallow active crustal regions in Italy within the ShakeMap workflow according to Michelini et al. (2020). Note that, in the area in Fig. 5, only 6 stations (triangles with thick contour) are used to adjust the field given by the GMPE. On both maps, the MCS intensity field available in the Macroseismic Intensity Database (DBMI15, Locati et al. 2021) is superimposed. On the SPEED map, at two selected stations, namely AQK and AQV, the simulated velocity waveform is compared with the recorded one (NS and EW component, respectively). From the comparison of Fig. 5 the following remarks can be made: • A good agreement is found between simulated and recorded waveforms in terms of peak values, frequency content and duration, although, at the stations of the Aterno transect AQV-AQG-AQA (not reported herein for sake of brevity, see Evangelista et al. 2017), simulations tend to provide lower amplitudes than recorded, most likely because of inaccuracies of the source model; • The spatial distribution of simulated PGV is in satisfactory agreement with the MCS intensity field: the areas of maximum PGV, which are found on the surface projection of the fault and inside the alluvial basin, correspond approximately to the sites with the highest intensity degrees. Consistently with the observed macroseismic field, the SPEED map shows a pattern oriented in the NW-SE direction with higher shaking toward SE, mainly because of the coupling of rupture propagation effects with site effects in the Aterno basin; • On the other hand, the ShakeMap is excessively homogeneous, being based predominantly on median estimates from GMPEs, and it fails to predict the strong anisotropy of ground shaking in the SE direction. Instead, the largest PGV values of the ShakeMap are found in the North sector of the Aterno valley (where the relatively high peak values recorded at the AQV-AQG-AQA array constrain the ShakeMap), in disagreement with corresponding MCS intensities of relatively modest degree (not larger than VII). • At larger distances from the fault (say more than about 15-20 km), the PGV values from PBS tend to be smaller than the ShakeMap, most likely because of the simplifications introduced in the modelling of the geological setting in the regions surrounding the L'Aquila epicentral area. As a matter of fact, the 3D numerical model follows a detailed geological description of the Aterno valley (see details in Evangelista et al. 2017), whereas a simplified geological model of the surrounding regions is considered, consisting of a horizontally layered crustal model (see Fig. 4).
Focusing on the detailed study area of Fig. 2, relevant for the fragility analysis presented in the following sections, the comparison between PBS and ShakeMap in terms of both PGA-Hmax and PGV-Hmax appears as in Fig. 6. To derive these maps, about 560 threecomponent broadband ground motion time histories are generated through the ANN2BB procedure on a regularly 500 m spaced grid of receivers. Starting from the simulated waveforms, a wide set of IMs, including Peak Ground Acceleration (PGA), Peak Ground Velocity (PGV), Peak Ground Displacement (PGD), Spectral Acceleration (SA) at different periods from 0 to 3 s, Arias Intensity (I A ), Cumulative Absolute Velocity (CAV) and Housner Intensity (HI), is computed. Finally, by applying a spline spatial interpolation, it is possible to associate to each building of the detailed damage database the corresponding IM value from PBS. Note that the availability of the entire ground motion time history allows to derive a broader set of IMs: more specifically, not only the standard peak measures as used in the ShakeMap (namely, PGA, PGV and SA at 0.3 s, 1.0 s and 3.0 s), but also integral measures which may offer further insights into their correlation with seismic damage. Although in this work we consider only PGA, PGV and SA, to enable a direct comparison with the standard approach based on ShakeMap, future studies will aim at evaluating the efficiency of a larger set of IMs.
The maps of Fig. 6 highlight further the differences in the spatial variability of ground shaking estimated by the two approaches. As shown in Paolucci et al. (2018), the PBS enriched by the ANN2BB technique at high frequencies provides a realistic spatial correlation of the peak ground motion values, with a relevant small-scale spatial variability, which reflects the physical features of the source rupture (not homogeneous across the fault plane) and of local site response characteristics. Furthermore, from a qualitative point of view, the spatial distribution of PGA shows a lower correlation range, as expected for short period IMs ). On the other hand, the ShakeMap provides a smooth pattern with limited spatial variability, also because of the coarse spatial interpolation grid with 1 km step. To quantify these differences, Fig. 7 shows the histograms of the PGA-Hmax (left) and PGV-Hmax (right) values computed from both approaches at the building locations of the detailed damage database. It is clear that the ShakeMap values have a more limited variability than PBS: standard deviation (in log10 units) decreases from 0.14 to 0.07 for PGA and from 0.12 to 0.04 for PGV. It is noted that the left tail of the frequency distributions from PBS, with minimum values of around 0.75 m/s 2 for PGA and 0.1 m/s for PGV, is shifted to lower values with respect to the ShakeMap results. This can be explained by the near-source location of the study area (see Fig. 6), where the coupling of smallscale asperities of the fault slip distribution with the geological heterogeneities induced by the presence of the Aterno valley implies a broader range of variability of predicted values from PBS. Such variability of ground shaking cannot be captured by the ShakeMap approach, where neither the complexity of the alluvial basin nor the details of the source rupture are accounted for.

Statistical approach and fitting technique
The derivation of fragility functions requires an appropriate statistical model and fitting technique to approximate observational data as a function of the ground motion severity.
A statistical model is defined by a systematic and a random component. The systematic component, i.e. the fragility function, is the mean of the response variable (i.e. damage in this study) as a function of the explanatory variable (i.e. the selected ground motion IM). The random component instead represents the conditional probability distribution of the response variable given the explanatory one.
In line with literature studies (e.g. Rota et al. 2008;Del Gaudio et al. 2017;Ader et al. 2020), the cumulative lognormal distribution is adopted for describing the probability of reaching or exceeding a preselected damage level, as a function of the seismic intensity measure: where ds i is the median value of the selected intensity measure associated with damage level ds i whereas β denotes the logarithmic standard deviation.
Given the availability of building-by-building data, the Bernoulli distribution (e.g. Rossetto et al. 2014;Ioannou et al. 2015) is selected for characterizing the random component of the statistical model: where Y ij is the damage response given the intensity measure threshold im j , n j is equal to 1, p ij is the exceedance probability defined by Eq. (1). The term Y ij is modelled as a binary outcome taking the value of 1 if the selected damage level is exceeded (i.e. DS ≥ ds i ) and 0 otherwise (i.e. DS < ds i ). The adoption of the Bernoulli distribution, which is a special case of the binomial distribution, allows for avoiding aggregation of damage data, which may introduce uncertainty in the resulting fragility estimates (e.g. Rossetto et al. 2014;Ioannou et al. 2015).
Fragility functions are simultaneoulsy fitted via the maximum likelihood estimate (MLE) approach and a unique constant dispersion value (β) is assumed for all damage states to prevent intersecting fragility curves (e.g. Lallemant et al. 2015;Porter et al. 2021;Nguyen and Lallemant 2021;Rosti et al. 2021a, b).
For each building typology, optimal parameters of the fragility model (i.e. ds i and β) result from maximising the logarithm of the likelihood function: where nDS is the number of damage levels and N is the number of data points. (1) 1 3

Fragility curves for PGA estimated from PBS and ShakeMap
Taking advantage of the statistical modelling discussed in Sect. 4.1, empirical fragility curves are derived for a wide catalogue of building typologies. Considering its extensive use in seismic vulnerability and risk applications (e.g. da Porto et al. 2021; Dolce et al. 2021), PGA is first selected for representing the severity of the ground motion shaking. As previously commented, for consistency between the PBS-and ShakeMap-based approaches, the Hmax component is considered. It is recognized that, to increase the accuracy of the fragility curves avoiding unreliable large probabilities of damage for very low IMs, it is necessary to condition the fragility curves at very low levels of ground motion IMs, accounting for the negative evidence of damage in the regions less affected by the ground shaking (Rosti et al. 2021a, b). To this end, ground shaking estimates are needed also at undamaged buildings in the non-surveyed and partially-surveyed municipalities, which fall largely outside the numerical model at more than 50 km from L'Aquila downtown (see Fig. 1). To overcome the lack of PBS results in this area, a hybrid strategy is pursued for defining seismic input at different building locations of the damage database, combining the PBS at the buildings in the detailed study area and the ShakeMap at the undamaged buildings in the non-or partially-surveyed municipalities. Table 3 collects the parameters (i.e. median values and logarithmic standard deviation) of the typological fragility functions in terms of PGA, estimated on the basis of PBS, supplemented by ShakeMap estimates in the regions less affected by the ground shaking. Typological fragility curves fully relying on ShakeMap estimates are also derived (Table 3). Comparisons between sets of fragility curves resulting from PBS and those entirely based on the ShakeMap are shown in Figs. 8, 9 and 10.
The two sets of fragility curves (PBS vs. ShakeMap) turn out to be consistent, with some limited differences. Specifically, PBS-derived fragility curves tend to be less conservative for masonry buildings (especially for the MH class), while a reverse trend is found for RC buildings (especially for post 1981 class). Such evidence of consistent results of both approaches is a key to support the use of input motions from PBS not only when ShakeMaps are not available, as for historical earthquakes, but also for conditioning fragility curves on IMs that can be derived solely from the availability of the ground motion time series (e.g. CAV, HI, I A ).

Fragility curves for alternative seismic intensity measures
One of the unquestionable advantages of PBS is the possibility of generating ground motion time histories at predefined grid points, thus permitting to locally characterize the ground shaking by a broad catalog of IMs, including peak, spectral and integral intensity measures. To exploit this potentiality, empirical fragility curves are also derived for IMs alternative to PGA, namely, PGV and average spectral acceleration (SAavg), in line with existing studies adopting average spectral quantities based on the arithmetic (e.g. Bradley 2010) and geometric mean of spectral accelerations over a range of periods (e.g. Kazantzi and Vamvatsikos 2015;Kohragi et al. 2017). In this study, SAavg is computed as weighted average of spectral acceleration values evaluated at 0 s, 0.3 s and 1 s. Although SAavg could be evaluated by referring to a different selection of periods (e.g. logarithmically/linearly equally-spaced, number of considered periods and definition of the period range), consideration of SA(0), SA(0.3) and SA(1) is motivated by the period values for which ShakeMap spectral acceleration estimates are available. Each spectral acceleration value  (1)) was then weighted based on the width of the associated period range (i.e. weights equal to 0.15 for PGA, 0.5 for spectral acceleration at 0.3 s and 0.35 for spectral acceleration at 1 s, e.g. Bradley 2010). With respect to spectral intensity measures evaluated at the building fundamental period (e.g. Rossetto and Elnashai 2003), spectral ordinates averaged over a range of periods allow for accounting for the variety of structural configurations within a given building typology. Furthermore, the adoption of the average spectral acceleration accounts for period elongation associated with structural damage and higher modes contributions (e.g. Rosti et al. 2020a). Analogously to Sect. 4.2, the ground motion IMs, represented by PGV and SAavg, respectively, are characterized by PBS in the area close to the L'Aquila municipality and by the ShakeMap in the less affected municipalities.
Based on the statistical procedure outlined in Sect. 4.1, sets of fragility functions are derived for each building typology, as a function of PGV (Table 4) and SAavg (Table 5). In the tables, parameters of the typological fragility functions entirely resting on ShakeMap estimates are also reported for comparison. It is noted that the dispersion values (β in Tables 3, 4 and 5, for the different IMs under consideration) of the PBS-derived fragility curves are systematically larger than those associated with ShakeMap, especially for RC buildings and for PGA. Average increments of β vary between around 4% and 7%, for PGA, and between about 8% and 18% for both PGV and SAavg, considering masonry and RC building typologies, respectively. We interpret the higher standard deviations for PBS as due to the greater variability in the spatial distribution of the IMs from the physics-based scenario, as clear from Fig. 7. As discussed previously, the variability of the ShakeMap results is likely underestimated in the near-source region (see for example the PGV spatial distribution in Fig. 7, right), because it does not account for the heterogeneities of both the slip distribution along the fault as well as the geological complexity of the area. Figure 11 shows empirically-derived fragility curves of some selected building typologies, as a function of different IMs, estimated by PBS (solid lines) and ShakeMap (dashed lines). As observed for PGA, PBS-derived fragility curves of URM building  typologies are generally less conservative than the corresponding ShakeMap-based ones. Difference between sets of PBS-and ShakeMap-based fragility functions is more significant when PGV and SAavg are considered (Fig. 11). PBS-based fragility curves of RC buildings seismically designed after 1981 are slightly more conservative in the lower ground motion range than the corresponding ShakeMap-based ones. For a given IM, the overall trend of PBS-and ShakeMap-derived fragility functions is nevertheless very similar, suggesting a general consistency among resulting fragility curves.

Validation against observed damage
Ground shaking scenarios from PBS and ShakeMap are respectively coupled with the corresponding fragility models for PGA, to verify their predictive capability in the context of seismic risk assessment. Referring to the detailed study area, building-bybuilding seismic scenarios are simulated in terms of physical damage. The availability of building-by-building data allows for a punctual definition of both seismic input and damage, avoiding the uncertainty driven by aggregated models (e.g. Nievas et al. 2022). The PBS scenario in terms of PGA is combined with the fragility model resulting from the hybrid strategy (PBS + ShakeMap). The ShakeMap ground shaking scenario ) is instead coupled with fragility functions fully resting on the ShakeMap . Figure 12 shows the spatial distribution of damage to buildings sited in the detailed study area (Fig. 12a) and in the L'Aquila historical centre (Fig. 12b), resulting from the joint use of ground shaking scenarios and consistent fragility models. Results are displayed in terms of probability of exceedance of preselected damage levels (i.e. DS1, DS3 and DS4), evaluated at the building level. In the figure, bar plots represent the frequency distribution of the exceedance probability values of each damage level. Bar colours correspond to the exceedance probability ranges of the discretisation adopted in the legend. For each building, a synthetic representation of damage (e.g. Dolce et al. 2003;Lagomarsino and Giovinazzi 2006;Rosti et al 2020b) is also provided in terms of mean damage (μ D ), defined as weighted average of the expected probabilities of occurrence (p k ) of the different damage levels (DS k ): Figure 13 shows the spatial and frequency distribution of mean damage in the study area and in the L'Aquila historical centre, resulting from the adopted simulation strategies. Although there is a general agreement between the statistical distributions of mean damage (4) D = 5 ∑ k=0 kp k Fig. 13 Comparison of mean damage predictions obtained from the two alternative simulation procedures: study area and L'Aquila historical centre from the two approaches, some differences are found especially at the scale of the historical center, where PBS tends to provide higher percentages of μ D ≥ 2.5-3.0, with a more significant clustering of μ D = 4-4.5 in the northern sector.
Besides contributing to the development of empirical fragility models, post-earthquake damage data (e.g. Dolce et al. 2019) are a unique opportunity for testing and validating the adequacy of seismic vulnerability and risk models for regional-scale applications (e.g. Smerzini and Pitilakis 2018;da Porto et al. 2021;Riga et al. 2021). In this study, the accuracy of the adopted simulation procedures to reproduce the observed seismic damage is globally assessed in terms of damage distribution within the detailed study area (Fig. 14). In the figure, predicted global damage distributions for masonry, RC and all buildings are compared to the observed ones. Predicted damage distributions are obtained by coupling PBS and ShakeMap ground motion scenarios with the corresponding fragility models. In Fig. 15, damage predictions are directly plotted against observations, highlighting possible overestimation or underestimation of the adopted simulation strategies with respect to each DS. Results show that both the adopted simulation procedures generally well reproduce the observed seismic damage in the detailed study area. Referring to the damage distributions for all buildings, predictions turn out to be in very good agreement with observations for all levels of damage from DS2 to DS5, with differences limited to 11% for both PBS and ShakeMap, indicating the general fitness of both procedures. Nevertheless, both procedures tend to overestimate null damage (DS0, with differences between 33% and 23% for PBS and ShakeMap, respectively), and to underestimate slight damage (DS1, with differences between 25% and 17% for PBS and ShakeMap). Larger discrepancies for the DS0/DS1 damage levels are likely driven by the need of constraining the lower tails of the fragility curves in the lower ground motion range (by taking advantage of undamaged buildings in the more distant municipalities), that, as already discussed, is a key and unavoidable step for reliable fragility assessment. Comparison between observed and predicted damage distributions is also detailed for each building typology (Fig. 16). In line with previous considerations, both simulation procedures allow to well reproduce the frequency of occurrence of damage levels from DS2 to DS5, which provide higher contribution to consequences estimation (e.g. da . Slight damage (DS1) is instead slightly underestimated in favour of null damage. This last finding can be explained by the need of constraining the tails of the fragility functions in the lower ground motion range, allowing for reliable seismic risk estimates (e.g. Rosti et al. 2021aRosti et al. , b, 2022.

Conclusions
This work validates the use of 3D PBS as a tool to characterize the ground shaking scenario to be used as input for empirical fragility studies on past earthquakes. The selected case study is the 2009 L'Aquila earthquake, for which a comprehensive database of postearthquake damage observations is available, together with relatively detailed ShakeMaps, constrained by few accelerometric records, to be used as a verification benchmark. Besides, a 3D numerical model, validated on the recordings of the event, is available encompassing a finite-fault kinematic source model as well as a detailed source-to-site wave propagation model embedding the topography of the region and the Aterno valley (Evangelista et al. 2017).
To the authors' knowledge, this work represents the first attempt to test the capability of the 3D physics-based numerical approach, which has attracted considerable research efforts in the last years, to provide region-specific ground shaking scenarios, to be used as input for the calibration of empirical fragility curves, as opposed to standard approaches relying on ShakeMaps or solely GMPEs. As a first step of this validation exercise, a novel set of fragility curves is empiricallyderived by statistical processing the L'Aquila post-earthquake damage database for several masonry and RC building typologies representative of the Italian building stock. Such database is presently the most complete in the Italian framework and likely one of the most complete in the Mediterranean context, where the building typologies considered in this work are commonly found, although with some obvious differences in terms of construction materials and building practice, especially for masonry structures.
The use of an improved statistical technique, relying on the Bernoulli distribution for modelling the random component of the statistical model and assuming a constant dispersion value for all damage states, allows for overcoming possible uncertainties in the resulting fragility estimates driven by data aggregation and for ensuring the ordinal nature of damage.
Fragility curves are calibrated by characterizing the ground motion intensity at the buildings located within the L'Aquila municipality, i.e. the area severely affected by the earthquake, using the broadband shaking scenarios from the PBS. The ground motion IMs selected for the fragility analysis include PGA, which is the reference ground motion intensity measure used in the platform for national seismic risk assessment in Italy, namely IRMA , as well as less standard IMs, such as PGV and SAavg (weighted average spectral acceleration in the range 0-1 s).
The fragility curves obtained from the PBS results are found in good agreement, for each building typology, with the ones obtained, through the same statistical approach, using the latest version of the ShakeMap (v4), released by the National Institute of Geophysics and Volcanology of Italy, that represent the reference approach to provide input motions for empirical fragility studies.
Finally, the two sets of fragility models, from PBS and ShakeMap, are coupled with the corresponding ground shaking scenarios to check the consistency of the predicted damage levels building-by-building and of their spatial distribution in the L'Aquila municipality area, with respect to the observed ones. Results point out that both approaches reproduce well the observed distribution of damage among the different classes, especially for damage states from DS2 to DS5 (errors limited to about 10%), that are expected to contribute the most to the loss estimation. Instead, some discrepancies, with differences up to around 30%, are found for lower damage states (DS0 and DS1).
Besides having provided a novel validated set of empirical fragility curves for different common building typologies in Italy, this study sheds light on potential advantages of simulation-based seismic shaking scenarios for empirical fragility studies, particularly when ShakeMaps are poorly constrained, or not available at all, owing to the lack of strongmotion recordings, as for historical earthquakes. More specifically, this work highlights the following main advantages related to the application of PBS to the problem at hand: • PBS are suitable to provide ground shaking scenarios with a realistic spatial correlation, reflecting the specificity of the regional seismic wave propagation features (e.g. alluvial basins, topography) as well as of the seismic fault rupture. Figure 6 and Fig. 7 demonstrate the capability of the PBS approach in providing a ground shaking scenario that, on the one hand, is consistent with the observed macroseismic pattern owing to the coupling of rupture propagation effects with local site response, and, on the other hand, is characterized by a realistic spatial correlation structure (see Infantino et al. 2021 for a quantitative evaluation of the spatial correlation structure of PBS). The latter is typically neglected by standard approaches for ground motion characterization, although it may affect results in the lower and upper intensity measure range (Rosti et al. 2020a), as well as the dispersion of the fragility curves. • Contrary to standard empirical models, the entire ground motion time history is available from PBS at any building site. This means that a wider portfolio of ground motion IMs can be computed, including not only peak values, as conventionally adopted in ShakeMaps (PGA, PGA and SA), but also integral measures (CAV, HI, I A ), reflecting in greater detail the variability of seismic shaking in amplitude, frequency and duration. It is, in fact, recognized that non-conventional or vector-valued ground motion intensity measures may improve the correlation with damage observations (e.g. Gehl et al. 2013;Masi et al. 2020). The sensitivity of fragility curves with respect to other intensity measures will be the subject of future studies. • PBS are suitable to obtain realistic site-specific seismic damage/risk scenarios for future earthquakes, especially in the near-source region, where the complexity of the fault rupture as well as the lateral variability of geological conditions may make unreliable the estimates based on standard empirical approaches. Such scenarios might be implemented in seismic risk platforms such as IRMA and used for civil protection activities.