Empirical fragility assessment using conditional GMPE-based ground shaking fields: application to damage data for 2016 Amatrice Earthquake

Recent earthquakes have exposed the vulnerability of existing buildings; this is demonstrated by damage incurred after moderate-to-high magnitude earthquakes. This stresses the need to exploit available data from different sources to develop reliable seismic risk components. As far as it regards empirical fragility assessment, accurate estimation of ground-shaking at the location of buildings of interest is as crucial as the accurate evaluation of observed damage for these buildings. This implies that explicit consideration of the uncertainties in the prediction of ground shaking leads to more robust empirical fragility curves. In such context, the simulation-based methods can be employed to provide fragility estimates that integrate over the space of plausible ground-shaking fields. These ground-shaking fields are generated according to the joint probability distribution of ground-shaking at the location of the buildings of interest considering the spatial correlation structure in the ground motion prediction residuals and updated based on the registered ground shaking data and observed damage. As an alternative to the embedded coefficients in the ground motion prediction equations accounting for subsoil categories, stratigraphic coefficients can be applied directly to the ground motion fields at the engineering bedrock level. Empirical fragility curves obtained using the observed damage in the aftermath of Amatrice Earthquake for residential masonry buildings show that explicit consideration of the uncertainty in the prediction of ground-shaking significantly affects the results.


Introduction
Accurate assessment of seismic risk for buildings at a territorial scale depends to a large extent on the availability of reliable and accurate fragility curves. The process of fragility estimation needs to make the most out of background information, expertise, available data and modelling/analysis capacities (Shinozuka et al. 2000;Jaiswal et al. 2011). The seismic fragility curves can be distinguished, based on the type of data used for deriving them, into four categories (Rossetto and Elnashai 2003); namely, analytic (e.g., Porter et al. 2007;Baker 2015;, empirical (e.g., Shinozuka et al. 2000;Jaiswal et al. 2011), based on expert opinion (e.g., Mosleh and Apostolakis 1986;Porter et al. 2007), and hybrid (e.g., Singhal and Kiremidjian 1998). Several advantages and shortcomings can be associated to the empirical seismic fragility assessment. The main advantage is that, being based on observed damage, the empirical fragility curves are-potentially-able to provide a realistic picture of post-earthquake damage. For instance, they can consider the possible soil-structure interactions, the cumulative damage due to aftershocks (or back-toback seismic events), the degradation in structural behaviour due to aging, and the effect of flexible diaphragm and secondary members (e.g., infills); just to name a few. The shortcomings are mainly associated to the difficulties in creating a homogenous (e.g., in terms of building class characterization, soil conditions, …, etc) spatial database of observed damage, inaccurate estimation of seismic ground shaking, site effects and spatial correlations in recorded ground motions and observed damage (e.g., Crowley et al. 2008;Colombi et al. 2008;Rossetto et al. 2014;Ioannou et al. 2015;Lallemant et al. 2015), the errors associated to the assignments of the observed damage state (e.g., Colombi et al. 2008;Rossetto et al. 2014) and finally a potential bias that could be introduced by underestimating the number of undamaged buildings (post-earthquake surveys tend to identify the damaged buildings) and hence the total number of buildings. In fact, the use of an established and standardized damage scale, in which different damage data should be converted, is one essential requirement for the derivation of empirical fragility curves. This permits to update the existing fragility curves as new data is gathered and to compare the vulnerability of different building classes. The European Macroseismic Scale (EMS-98, Grünthal 1998) is an example of a standard damage scale, very often used to report the empirical fragility curves in Europe.

A brief history of empirical fragility assessment for buildings in Italy
Several works have focused on empirical assessment of the vulnerability of Italian building stock (mainly masonry and reinforced concrete). Braga et al. (1982) used the MSK-76 damage scale (Medvedev-Sponheuer-Karník) to derive damage probability matrices (DPM) showing the probability of being in each of the five damage levels of MSK-76 for three building vulnerability classes and a 12-level damage-dependent seismic intensity scale (MSK-76). Sabetta et al. (1998) derived both the vulnerability (i.e., mean normalized damage versus seismic intensity relations, see Calvi et al. 2006) and the fragility curves by adopting PGA and Mercalli Cancani Sieberg (MCS) macroseismic scale as measures of seismic intensity. This study adopted the five damage levels of MKS-76 scale as measure of damage for three building vulnerability classes. Orsini (1999) has derived continuous vulnerability curves expressed as cumulative distribution of the seismic intensity level (a parameter-less scale of intensity, PSI, Spence et al. 1992) related to the excursion of a given damage level. Dolce et al. (2003) have derived damage probability matrices and scenario-based mean damage maps based on MSK-76 scale considering the site effects for four building vulnerability classes for town of Potenza in Southern Italy. Di Pasquale et al. (2005) have reported the damage probability matrices for four vulnerability classes as a function of MCS intensity. Lagomarsino and Giovinazzi (2006) have derived vulnerability curves reporting the mean damage ratio versus seismic intensity expressed in the EMS-98 scale for masonry and RC buildings. Zuccaro and Cacace (2009) mapped the Italian building inventory in terms of damage probability matrices based on EMS-98 scale. Zuccaro and Cacace (2015), instead, proposed a methodology to reduce the variability in the vulnerability classification of EMS-98 through the application of vulnerability modifiers. These modifiers expressed the shift in vulnerability (expressed in terms of a synthetic vulnerability index) for various "sub-classes" identified within a given vulnerability class with respect to the overall vulnerability index for that class. Rota et al. (2008) derived empirical fragility curves based on damage survey data for 150,000 buildings from past Italian earthquakes, occurred in the period spanning from Irpinia (1980) to Molise (2002) earthquakes. These empirical fragility curves are reported using EMS-98 as the damage scale and PGA as the ground shaking intensity. The seismic intensity is calculated as a single mean estimate evaluated for the corresponding municipality based on Sabetta and Pugliese (1996) ground motion prediction equation (GMPE) and assuming rock as the soil type. Liel and Lynch (2012) employed post-earthquake investigations following the L'Aquila earthquake to quantify the seismic fragility of RC buildings. The seismic intensity in this study is characterized by PGA derived from the ShakeMap (Michelini et al. 2008) for L'Aquila's main event. De Luca et al. (2015) and Del Gaudio et al. (2016) derived empirical fragility curves based on a database of post-earthquake building surveys conducted on 131 RC buildings located in the town of Pettino after the 2009 L'Aquila Earthquake. The EMS-98 compatible damage levels were deducted from the standard Italian survey sheet AeDES (Baggio et al. 2007) and seismic intensity was expressed in terms of PGA contours of the L'Aquila main event ShakeMap. Del Gaudio et al. (2017) have derived empirical fragility curves based on the large database of post L'Aquila 2009 Earthquake damage surveys conducted by the Italian Civil Protection (Dolce and Goretti 2015) on RC buildings. De Luca et al. (2018) have used the observed damage for infilled RC structures in the aftermath of the 24th August 2016 Amatrice Earthquake for Bayesian updating of existing empirical fragility curves for central Italy. Rosti et al. (2018) have used the database of damage data collected after the 2009 L'Aquila (Italy) earthquake, to derive damage probability matrices for several building types representative of the Italian building stock. Ioannou et al. (2020) have employed a Bayesian framework to study the importance of considering the uncertainty in the ground motion intensity on the shape of empirical fragility curves obtained based on post-disaster data aggregated at the municipality level for 1980 Irpinia Earthquake.

The focus: empirical fragility estimation based on conditional GMPE-based ground-shaking fields
One important aspect that emerges from the available rich literature on empirical vulnerability and fragility assessment is that accurate estimation of the ground shaking intensity is crucial (as much as accurate damage estimation) towards accurate and reliable empirical fragility assessment (e.g., Hsieh et al. 2013;Yazgan 2015;Lallemant et al. 2015;Ioannou et al 2015). Gaussian random ground shaking fields, whose first-and second-moment statistics are described by the ground motion prediction equations (GMPE) and the compatible spatial correlation models (i.e., models describing the intra-event variability in the residuals of the same GMPE), are used quite often in the literature (e.g., Park et al. 2007;Goda and Hong 2008;Jayaram and Baker 2009;Goda and Atkinson 2010;Goda 2011;Sokolov and Wenzel 2011, Weatherill et al. 2015. Based on a procedure used originally for generation of ground shaking fields in the frequency domain (Vanmarcke et al. 1993), Park et al. (2007) and also Crowley et al (2008) proposed an updating procedure that yields updated or "conditional" Gaussian random fields based on available recordings for various stations (a.k.a., correlated ground motion fields with "anchor points"). The conditional Gaussian random ground shaking field has been used mostly for portfolio loss assessment (e.g., Park et al. 2007;Crowley et al. 2008;Miano et al. 2015Miano et al. , 2016; however, it has not been used extensively for empirical fragility assessment. Straub and Der Kiureghian (2008) discuss the differences between empirical fragility assessment based on correlated observations (such as ground shaking for a given event) and that based on independent observations. Hsieh et al. 2013 use the 444 free-field ground motion recordings of the Chi-Chi Taiwan Earthquake to generate the map of ground shaking. They use a kriging model to interpolate the ground-shaking for the points with no recording stations nearby. Lallemant et al. (2015) recommend updating the ground-shaking field based on the recorded ground shaking data. They perform the updating by selecting the ground shaking field realizations that match the recorded ground shaking data (two stations in the example). Ioannou et al. (2015) perform updating using a kriging model in the case where the ground motion recordings are available at the location of a subset of the surveyed buildings. Yazgan (2015) uses the conditional ground shaking fields in a Bayesian framework to update the parameters of prescribed fragility models based on the likelihood of observing damage in the location of the surveyed buildings for a given class. It is interesting to note that the procedure employed in this work for updating the statistical properties of the Gaussian random fields does not rely on the assumption that recording stations need to be close to the buildings' location, as discussed later on. In fact, in the case of Amatrice Earthquake, most of the ground motion recording stations are located far field (this is shown later in the paper). Nevertheless, the updating procedure can include all the stations. The updating procedure can be carried out automatically and with no significant computational effort (the procedure can be implemented using Matlab's toolboxes).
This work presents a procedure for the generation of conditional GMPE-based ground shaking fields for derivation of empirical fragility curves. Section 2 presents the methodology employed in this work for deriving empirical fragility curves. Section 3 presents the application of the methodology in deriving empirical fragility curves based on damage data for masonry buildings available right after the Amatrice Earthquake of 24th of August 2016. Section 4 is dedicated to the final conclusions drawn from the results outlined in Sect. 3.

Basic assumptions for empirical fragility assessment for a class of buildings
The concept of fragility curve is applicable to Poisson-type limit state excursion (see e.g., Der Kiureghian 2005; Jalayer and Ebrahimian 2020). In the context of seismic risk, the uncertainty related to earthquake occurrence can be classified as leading to Poisson-type limit state excursion. Using the damage survey data from spatially-distributed surveyed buildings belonging to the same class for fragility modelling implies that the building-to-building variability within the same class is inevitably going to be considered. In other words, part of the dispersion captured by the resulting fragility curve is due to the differences between the surveyed buildings, within the same class, and not just due to the uncertainty in the ground motion representation. It can be argued that using a single class fragility curve is equivalent to the assumption that each of the surveyed buildings is replaced by a nominal average building for that class. One consequence of such assumption is that the building-to-building variability is going to be considered in "disguise" and as a contribution to record-to-record variability related to the ground motion. On the other hand, the ground motion intensity level at each building location (estimated in this work based on the simulation procedure discussed in the following paragraphs) cannot be independent as they correspond to the same earthquake event and they can exhibit various levels of correlation based on the mutual distance between the buildings. In other words, the ground-shaking levels registered at the location of each surveyed building can show significant correlation for close-by buildings. Excluding such correlation may lead to unrealistic and un-conservative fragility curves (Straub and Der Kiureghian 2008). Unconservative because the fragility curve that is derived based on the assumption of independent intensity values is going to rely on a larger information content compared to the fragility curve derived with explicit consideration of possible correlations between the intensity values-of course if such correlation is significant (this is the case when the buildings are close by). The spatial correlation in the ground motion is modelled herein by considering the spatial correlation structure in the residuals of the ground motion prediction equation (GMPE). Although assigning a damage grade to buildings based on rapid visual survey techniques (e.g., visual inspections, use of satellite imagery) is certainly subjected to errors (e.g., human/device error of judgment and/or an inevitable degree of subjectivity), the uncertainty in assigning the damage grades to surveyed buildings and the fragility model parameter uncertainties are not considered here.

Empirical fragility assessment for a class of buildings
The empirical fragility is calculated within an updated robust reliability assessment framework (Papadimitriou et al. 2001;Beck and Au 2002;Jalayer et al. 2010). This framework is also known in literature as "predictive fragility" (e.g., Sasani and Kiureghian 2001;Straub and Der Kiureghian 2008). The concept of robust reliability is used for calculating "robust" empirical fragility curves. The term robust in this context implies that the uncertainty in the evaluation of the ground shaking is considered (see Miano et al. 2018;Jalayer and Ebrahimian 2020 for examples of robust analytical fragility assessment). Let the survey dataset D Damage define the set of N damage states D i (e.g., the EMS-98 damage grade, i can vary between 0 and 5) for the buildings surveyed. Let PGA = {PGA i , i= 1:N} denote the vector of ground shaking values (expressed in terms of peak ground acceleration) estimated at the position of each surveyed building. The robust empirical fragility is calculated by considering the uncertainty in generating a plausible field of PGA values conditioned on both the vector of PGA values registered at various accelerometric stations, denoted by D PGA , and the observed damage pattern D Damage for the portfolio of the buildings surveyed. The updated empirical fragility conditioned on available data D Damage and D PGA can be written as follows: (1) where D i denotes a given damage grade; pga = x is a given level of ground-shaking; D Damage (the damage survey data) is the survey data acquired for the portfolio of buildings (for all the building classes); f(PGA|D PGA , D Damage ) denotes the joint probability density function (PDF) for the ground-shaking field; Ω(PGA) denotes the domain of all plausible PGA fields; the joint PDF represents the joint probability of observing the ground shaking levels for building survey locations f(PGA|D PGA , D Damage ) given the damage survey data D Damage and the registered PGA data D PGA , where PGA = [PGA i , i= 1:N] (see Sect. 2.3 for the derivation of f(PGA|D PGA , D Damage )). In general, it is quite complicated to calculate the integral in Eq. (1) analytically. It is usually estimated using standard Monte Carlo Simulation (MCS). In the adopted MCS approach, N sim realizations of the ground shaking field PGA are generated based on the joint PDF f(PGA|D PGA , D Damage ). The simulation procedure yields estimators of the updated empirical fragility and its standard deviation (see also Jalayer et al. 2016: where P(D > D i |x,D PGA ,D Damage ) is the updated empirical fragility curve (for a given building class) and σ[P(D > D i |x,D PGA ,D Damage )] is the standard deviation (error in the estimation) of the updated empirical fragility curve based on empirical damage data D Damage and registered PGA data D PGA .

Generation of GMPE-based ground shaking fields
The joint probability density function f(PGA) for the vector of PGA = [PGA i , i= 1:N] values at the location of each building of interest for a given earthquake scenario can be evaluated by employing a ground motion prediction equation (GMPE). A full probabilistic representation (based on the first two moments) of GMPE can be expressed in terms of multi-variate Normal distribution which is identified by its expected value vector M and covariance matrix Σ. That is, once the first two moments are given, several realizations of the ground shaking field can be generated. Herein, the ground motion prediction model proposed by Bindi et al. (2011, also known as ITA10) for the peak ground acceleration (PGA, the geometric mean of two horizontal components) as the intensity measure is used. A recent work by Michelini et al. (2019) ranks ITA10 quite favourably as the GMPE to use for Italian earthquakes from shallow active crustal zones. It is worth noting that there are more recent pan European BND14 (Bindi et al. 2014a, b) and Italian GMPE's (ITA18, Lanzano et al. 2019) available. The latter includes the recent destructive Italian earthquakes such as Emilia Sequence 2012 and the Central Italy Sequence 2016. To authors' knowledge, no well-established spatial correlation models have been calibrated so far to the residuals of these more recent GMPE's. For instance, a very recent Italian correlation law  seems to have been calibrated to the residuals of a recent regional GMPE for Northern Italy , Emilia 2012 recorded ground motions) and not to ITA18. Nevertheless, the methodological procedure presented in this work is quite versatile and can use more recent GMPE's and their related correlation structures as they become available. However, it is not imperative for this procedure that the adopted GMPE has the registrations of the earthquake of interest (in this case Amatrice Earthquake) embedded. In fact, the updating procedure (a sort of statistical base-line correction) aims to adjust potential biases with respect to a GMPE that does not include the registrations of the earthquake of interest.
The functional form of the adopted GMPE model is the following: where E[log 10 PGA] is the expected value (first moment) for the (base 10) logarithm of peak ground acceleration (PGA, in cm/s 2 ); e 1 is a constant term, F D (R jb ,M), F M (M), F s and F sof represent the distance function, the magnitude scaling, the site amplification and the style of faulting correction, respectively. M is the moment magnitude, R jb is the Joyner-Boore distance in km (or the epicentral distance when the fault geometry is unknown-generally when M<5.5). The proposed equation for the distance function The covariance matrix, Σ, is defined as the sum of two inter-event and intra-event components: where σ intra represents the intra-event variability and σ inter represents the inter-event variability (both parameters are tabulated in Bindi et al. 2011); e is the all ones matrix and R is the matrix of correlation coefficients. R is composed of unit diagonal terms and off-diagonals equal to ρ jk , j≠k (both varying from 1 to N; where N is the total number of buildings surveyed). The covariance matrix is obtained according to the following formulation of ρ jk (Esposito and Iervolino 2012): where h jk represents the distance between sites j and k and b(T) is a tabulated coefficient equal to 10.8 km for PGA. It is to note that the above expression for the correlation coefficient has been calibrated to the residuals of the Bindi et al. (2011) GMPE adopted herein and therefore is reasonably consistent for the purpose used in this study.

Considering stratigraphic and topographic factors
Two alternatives are considered for considering the stratigraphic site effects; namely, (a) the coefficients imbedded in the GMPE (here Bindi et al. 2011, ITA10); (b) application of stratigraphic amplification factors to ground shaking at bedrock. Landolfi et al. (2011) and later Tropeano et al. (2018) propose stratigraphic coefficients that consider non-linear soil column propagation effects. Landolfi et al. (2011) have merged empirical, semi-empirical and analytical datasets to compute non-linear relationships quantifying stratigraphic amplification for different classes of subsoil profiles. That is, site effects are evaluated through the stratigraphic amplification factor, which is directly multiplied by the reference (i.e., propagated to bedrock) peak ground acceleration from the GMPE by Bindi et al. (2011) to obtain the peak acceleration at surface (where the subscribes B, C, D represent the soil class): It is also possible to apply topographic amplification/deamplification factors (S T ) to the GMPE. S T depends on the shape of slopes, because irregular surface geometry affects the focusing, defocusing, diffraction and scattering of seismic waves. This can lead to a change in amplitude, frequency and duration of ground motion compared to flat ground conditions (Sanchez-Sesma 1990;Paolucci 2002). Following Garcìa-Rodriguèz et al. (2008), a geometrical parameter more suitable for small scale studies seemed to be the slope curvature, which can be obtained from the digital elevation model (DEM) of the area (at 20 × 20 m resolution for the case-study presented here). This index permits to mark the concave and the convex features of a landscape, with negative and positive values respectively, accounting for attenuation in valleys and the seismic waves focusing on ridges. The effectiveness of this parameter was also validated by the numerical study of Torgoev et al. (2013) and adopted in seismic slope stability analyses by Silvestri et al. (2016) and reported in Table 1.
It should be noted that the above-mentioned stratigraphic and topographic factors are going to be applied herein in a deterministic manner. It can be shown that this changes the mean vector M by a constant (i.e., the inner product of a vector of constant factors and the median in the arithmetic scale) and leaves the covariance matrix Σ unaltered. (9)

Updating the generated ground shaking fields based on registered PGA data
One interesting feature of the method adopted for generating the ground shaking field realizations is that it can be updated based on the registered values. It was assumed in Sect. 2.3.1 that the PGA values at the location of each surveyed building are distributed as a joint multi-variate lognormal distribution. One of the specific characteristics of a joint Normal distribution for a vector of variables is that any given partition of the vector conditioned on the remaining components of the vector is still going to be a joint Normal distribution. With specific reference to the case of the vector of log 10 PGA values denoted as data D PGA , let the vector of mean values M and the covariance matrix Σ be partitioned as follows (Park et al. 2007): where M 1 is the mean (of the base 10 logarithm) vector of PGA = [PGA i , i= 1:N cl ] values according to the adopted GMPE; M 2 is the mean vector of calculated log 10 PGA at the stations within the area of interest (according to the adopted GMPE); Σ 11 is the covariance matrix for the calculated (from the GMPE) log 10 PGA for the surveyed buildings of class CL; Σ 12 = Σ 21 is the cross-covariance matrix for the log 10 PGA values calculated (from the GMPE) at the location of the surveyed buildings and those calculated at the location of the stations; Σ 22 is the covariance matrix for the log 10 PGA values calculated at the stations. As described briefly above, the conditional distribution of the calculated log 10 PGA values given the registered log 10 PGA values at the stations is a joint Normal distribution with mean vector M 1|2 and covariance matrix Σ 11|22 : where D PGA is the vector of the registered log 10 PGA values for the stations. Note that M 1|2 and Σ 11|22 represent the first two moments of the updated joint lognormal pdf f(PGA|D PGA ).

Updating of the generated ground shaking fields based on observed damage
Once the first two moments of the updated ground shaking field corresponding to the lognormal PDF f(PGA|D PGA ) are obtained, various plausible realizations of the groundshaking random field can be generated. As the formulation in Eq. 1 indicates, the groundshaking field is also conditioned on the observed damage survey data for the portfolio of surveyed buildings as reflected by the joint PDF f(PGA|D Damage ,D PGA ). A rejection sampling approach is adopted to reflect the conditioning of the random field joint PDF on the damage data. That is, the random field realizations that lead to non-monotonically increasing fragility functions are discarded. This is achieved by checking whether the gradient of the fragility curves' slope is negative or not for all the classes and for all the damage levels considered. As a result, those generated random fields (for a given portfolio of buildings) which do not lead to meaningful fragility curves for all the classes and all the considered damage levels are discarded. This essentially means that the random field realizations that are not compatible with the observed damage pattern are not plausible. It is to note that, this updating procedure assumes that the damage survey is carried out with no error.

The empirical fragility curves
This section describes how the fragility curve P(D > D i |x,PGA k ,D Damage ) can be obtained for a given realization of vector of ground shaking values at the location of surveyed buildings PGA k generated based on pdf f(PGA|D PGA , D Damage ) and given the damage survey D Damage . The empirical fragility curves P(D > D i |x,PGA k ,D Damage ) (k= 1:N sim ) have been calculated according to a logistic regression probability model (Wiesberg 1985, for instances of application in fragility assessment see e.g. Charvet et al. 2014;Lallemant et al. 2015;Yazdi et al. 2016;De Risi et al. 2017). This type of regression is suitable for cases where the dependent variable is binary (i.e., either 0 or 1). Thus, it is especially suitable for estimating the probability of exceeding a damage state D i . That is, for damage state D i , each surveyed building can either: a) exceed the designated damage state denoted as 1 or b) NOT exceed the designated damage state, denoted as 0.
Denoting probability of exceeding damage state D i as a function of the ground-shaking level PGA= x j as π i (x j ), the likelihood of having r i buildings that exceed damage state D i for the N CL buildings surveyed for class CL can be expressed as follows: N cl is the number of surveyed buildings for building class CL; π i (x) is expressed through the following relation: It is to note that D cl is the vector of observed damage grades for the buildings in class CL. Note that D cl is a subset of damage data D Damage for the entire portfolio of buildings surveyed. Figure 1 shows an example of empirical fragility curve obtained by employing logistic regression. As it can be observed from Fig. 1, the fragility is expressed as the conditional probability of exceeding damage level D5 for the class of buildings CL= MBC4 (Masonry buildings with tie rods or ring beams with number of stories > 2, described later in Sect. 3.1.2). The points correspond to the N MBC4 buildings surveyed for the class MBC4. The buildings whose damage state is estimated to be greater than or equal to D5 are assigned the value of 1 and rest of the surveyed buildings are assigned the value of 0.

The study area
Between August and October 2016, Central Italy was stricken by three damaging earthquakes. The first Mw 6.0 event occurred on August 24th at 01:36 UTC close to Accumoli village (referred to as Amatrice Earthquake); it was followed by a long seismic sequence, which 2 months later produced a Mw 5.9 aftershock on October 26th at 19:18 UTC at 3 km West of Visso and a Mw 6.5 event on October 30th at 06:40 UTC, 6 km North of Norcia (see Ebrahimian and Jalayer 2017 for more details about the Central Italy seismic sequence; see Sextos et al. (2018) for more details about observed damage).
This paper focuses on the 3 municipalities and 4 fractions, namely, Amatrice, Accumoli, Arquata del Tronto, Tina (Accumuli), Illica (Accumuli), and Pescara del Tronto (Arquata del Tronto),Trisungo (Arquata del Tronto), affected by the 24th August 2016 event as shown in Fig. 2. The figure shows a simplified geological map of the area of interest, overlaid with the PGA contour map provided by INGV (Italian Institute of Geophysics and Volcanology) official ShakeMaps and the fault projection.
According to the PGA distribution provided by INGV, ground acceleration values as high as 0.5 g affected several small towns in the vicinity of the epicentre. Among these towns, Amatrice is the one that has had the most widespread destruction and the highest number of fatalities. The two hamlets of Tino and Illica have been considered as examples of the stratigraphic and the topographic amplification of ground motion. The masonry buildings located in the above-mentioned towns were strongly damaged by the August 24th event and its aftershocks.

Damage survey data D Damage : satellite imagery (Copernicus) and visual survey
The European Macroseismic Scale EMS 1998 (Grünthal 1998)  generation of "damage grading" maps, made possible by comparing pre-and post-event satellite images. With specific reference to damage data observed for Amatrice Earthquake, Masi et al. (2017) indicate that satellite imagery is more reliable for damage grades higher than or equal to Grade 4 and less reliable for damage grades lower than Grade 4. Therefore, although it is quite efficient for performing damage survey for vast areas, its use for a complete damage classification should be ideally accompanied by other means of damage surveying such as field surveys and areal photos. Visual Survey: It is to note that the Copernicus-EMS maps provide the damage grades for the buildings. However, these maps do not provide information about the sub-division of buildings into different classes. Hence, with the objective of complementing the damage grading maps provided by Copernicus, the same portfolio of buildings is also subjected to visual survey. The visual survey was based on photography available from field trips done at the end of September 2016 (courtesy of G. Forte and A. Santo), videos provided by drone (Feliziani et al. 2017, taken days after the Amatrice Earthquake) for areas that were difficult to access and google street view (https ://www.googl e.it/stree tview , last access 01/06/2017, used only for verifying the building classes and not the damage grades). The visual survey done building by building is arguably one of the most reliable means for assessing the occurred damage. However, the visual survey carried out herein is limited by the number of available photos/videos.

Structural characterization of the buildings
The portfolio of surveyed buildings is limited to residential masonry buildings. Further breaking down was based on the recommendations in Rota et al. (2008). In Rota et al. (2008), the masonry buildings are classified according to four parameters: number of floors, the presence of tie rods or tie beams, the type of horizontal structure (flexible or rigid floors), regular or irregular masonry layout. However, due to limitations posed by visual survey, the breaking down into more detailed classifications within the portfolio of residential masonry buildings in this work is limited to two factors: the number of storeys and the presence of tie rods or ring beams. As highlighted in a study conducted by Sorrentino et al. (2019), around 1/5th (20%) of the portfolio of masonry buildings in Amatrice and Accumuli presented either ring beams or tie rods. This percentage is confirmed for Amatrice also in a pre-earthquake survey (Fumagalli et al. 2017). In fact, the percentage of buildings that present either tie rods or rings beams according to visual survey in the present study for Amatrice and Accumuli is 28% (for all the seven towns considered this percentage is around 24%). The pre-and post-earthquake surveys conducted in the literature (Fumagalli et al. 2017;Sorrentino et al. 2019) indicate that the two municipalities of Amatrice and Accumuli are for the most part consisted of masonry buildings with following characteristics: up to 5 storeys (with a prevalence of 3 storey buildings); multi-layered masonry walls in the ground floor; units approximately dressed or undressed; medium to low mortar quality. It seems that, for the most part, the surveyed buildings could be classified as irregular (due to insufficient connections between the horizontal and vertical structures, between the masonry layers and the low to medium strength mortar used). Therefore, it has been assumed that the surveyed portfolio consists of irregular masonry buildings. As a result, four distinct classes of masonry buildings have been defined: MBC1 can be associated to the union of classes IMA2 (Masonry-irregular layout-flexible floors-without tie rods and ring beams-1 or 2 stories) and IMA4 (Masonry-irregular layout-rigid floors-without tie rods and ring beams-1 or 2 stories) of Rota et al. (2008); MBC2 to the union of classes IMA6 (Masonry-irregular layout-flexible floors-without tie rods and ring beams-more than 2 stories) and IMA8 (Masonry-irregular layout-rigid floorswithout tie rods and ring beams-more than 2 stories); MBS3 to the union of classes IMA1 (Masonry-irregular layout-flexible floors-with tie rods or ring beams-1 or 2 stories) and IMA3 (Masonry-irregular layout-rigid floors-with tie rods or ring beams-1 or 2 stories); MBC4 to union of classes IMA5 (Masonry-irregular layout-flexible floors-with tie rods or ring beams-more than 2 stories) and IMA7 (Masonry-irregular layout-rigid floors-with tie rods or ring beams-more than 2 stories). Figure 3 reports the numeric break-down of the portfolio of buildings into four classes, classified by the town to which they belong.
It is to note (Fig. 3) that the three main municipalities Amatrice, Accumuli and Arquata del Tronto show similar proportionality between the different masonry buildings classes considered. This confirms a rather homogenous structural characterization of the buildings in the largest centres of the studied area. Tables 2 and 3 report the number of buildings for each class and for each damage grade obtained based on Copernicus EMS and the visual survey, respectively. Figure 4 shows the damage grades provided by Copernicus for the surveyed buildings for the towns of (a) Amatrice, (b) Accumoli, (c) Arquata del Tronto, (d) Pescara del Tronto, (e) Illica, (f) Tino, and (g) Trisungo (the buildings belonging to damage levels D0 and D1 are grouped with label D0 + D1 in Fig. 4). It is worth noting that damage levels D0 and D1 could not be distinguished using the survey techniques adopted in this work. Consequently, fragility curves for D1 are not reported hereafter. It is to note that the surveyed building stock does not include all the masonry buildings of the seven towns under investigation. The buildings included in the survey are those for which both    Figure S1). In the work of Fiorentino et al. 2018, it is stated that 48.6% of the masonry buildings of the historical centre of Amatrice had collapsed; while 19.2% experienced partial collapse. Associating collapse to D5 and partial collapse to D4, the collapse and partial collapse percentages are 37.9% and 20.7% based on visual survey and 40.0% and 22.1% based on Copernicus. Figure S2 (in the electronic supplementary file) shows the histograms of number of cases > Di (i= 2:5) for visual survey and Copernicus for all four classes of masonry buildings considered. An overall fine agreement can be observed (especially for D3, D4 and D5).

Site effects for the surveyed buildings
Post-earthquake field recognition identified, as the most affected area, the valley of Tronto river. This valley is host to several municipalities and hamlets. The local geological and geomorphological setting can be sketched as shown in Fig. 5. Indeed, the villages occupy either the valley close to the river (e.g., Trisungo) or the cliffs overlooking it; with the latter being located usually at the top of small ridges and ancient erosional terraces (e.g. Amatrice, Accumoli, Arquata del Tronto, Fig. 5a) or located on the slopes (e.g. Pescara del Tronto, Illica, Tino, Fig. 5b). The villages located on cliff-type morphology (Fig. 5a) are almost always bordered by steep slopes (25°-35°) with heights varying from 20 to 80 m. For these areas, buildings labelled with D4 and D5 damage grades are widespread and mainly localized near the steep escarpments and in the narrower part of the ridges as shown in Figs. 4b (Accumuli), (1) Amatrice; (2) Accumoli; (3) Pescara del Tronto 4c (Arquata del Tronto) and Fig. 5 part 1 (Amatrice). These buildings were affected by seismic waves' focalization due to topographic shape effects (e.g., Grelle et al. 2018). These topographic effects were not present in lowland areas of the valley which suffered less damage, relatively speaking (see Trisungo, Fig. 4g). On the other hand, other towns shown in Figs. 4a (Amatrice), 4d (Pescara del Tronto), 4e (Illica), 4f (Tino), suffered widespread damage due to both topographic and stratigraphic effects (see Fig. 5 part 3). Some of these towns lie on slopes characterized by few meters of soft soils resting on a stiffer material (see Fig. 5b, Accumuli); where stiff arenaceous formation of the Laga Flysch is buried by few meters of weathered deposits and colluvium mainly made of silty sands. In the case of Pescara del Tronto (Fig. 4d), the hamlet lies on debris and travertine sands resting above a limestone bedrock. Figures S3 (in the electronic supplementary file) shows the overall damage break-down with respect to geologic units for MBC 1-2 and MBC 3-4, respectively. To account for stratigraphic amplification effects at a small scale, Eurocode 8 seismic soil classes are attributed to lithologic units. This is done according to the suggestions provided by Forte et al. (2017) for a comparable geo-lithological setting and is reported in Table 4. Both adopted approaches (a) and (b) for considering the stratigraphic amplification are based on the soil classes. Therefore, they can be classified as 2-Grade level of microzonation following the guidelines proposed by ISSMGE (1999). In this study, the soil class was derived by the most recent soil class map of Italy published in Forte et al. (2019). This map was developed accounting for the shear-waves measurements distribution for each geological formation identified on a reinterpretation of 1:100.000 geological maps. The results are provided in a 500 × 500 m grid resolution, which is higher than the other available maps such as those in Michelini et al. (2008, Italian ShakeMap).

Generation of GMPE-based ground-shaking fields for the Central Italy Earthquake (24 August 2016) scenario
According to the procedure outlined in Sects. 2.2 and 2.3, N sim = 1,000,000 realizations of the GMPE-based ground shaking fields are generated for the Amatrice Earthquake scenario (Mw= 6.0) providing the peak ground acceleration (PGA) at each building location. It is worth noting that the realizations are generated for the entire portfolio of buildings surveyed (i.e., all the classes together). Depending on the approach followed for considering the stratigraphic amplification, following steps are taken: for approach (a) outlined in Sect. 2.3.2, the ground-shaking field realizations are generated directly by employing the GMPE and considering the specific coefficients embedded in the GMPE for EC8 soil categories for stratigraphic amplification. These realizations are further modified in order to consider the topographic effect; for approach (b) outlined in Sect. 2.3.2, the generated ground-shaking fields on bedrock (i.e., assuming soil type A) are multiplied by both the stratigraphic and topographic amplification/deamplification factors calculated according to Eq. 9 and Table 1, respectively. It should be emphasized that the topographic effects are considered in the same manner for both approaches. Since the uncertainty in the evaluation of these amplification factors is not considered, their application affects only the median M (multiplied linearly by amplification factors) and leaves the covariance matrix Σ unaffected. Finally, in order to obtain the conditional GMPE-based fields, PGA registrations for eighty one accelerometric stations within a 100 km distance of the epicentre (see Fig. 6a and b) are employed in order to update the ground motion fields according to the procedure described in Sect. 2.3.3. Figure 6a shows the location of the 81 stations used for this work. Figure 6b shows the 16th, median, and 84th percentiles of the GMPE (Bindi et al. 2011), considering the stratigraphic (soil B, Landolfi et al. 2011) and topographic (curvature class 3, Silvestri et al. 2016) factors applied to the GMPE, before and after the updating with the available registrations. It should be noted that the results can be sensitive both to the position and the number of accelerometric stations considered. It is worth noting that: 1-theoretically-speaking, the recorded ground shaking information from all the accelerometric stations can be included in the updating; 2-the results are most sensitive to the accelerometric recordings inside the polygon delineating the surveyed buildings; 3-in the case of Amatrice earthquake, only one station is located within the polygon that delineates the surveyed buildings (see Fig. 6a). Therefore, in lieu of enough accelerometric stations within the polygon of the surveyed buildings, a 100 km distance limit from the epicentre was set. Roughly speaking, this marks the distance beyond which the ground shaking intensity is not going to be significant for shallow crustal earthquakes. Figure 6c maps the median PGA values for 25,000 realizations (note that smaller number of realizations are necessary for generating the map in Fig. 6c with respect to those necessary for estimating empirical fragility) rendered by a mesh-grid of 500 m × 500 m resolution. Figure 6d, demonstrates the median of the conditional GMPE-based ground-shaking field for the surveyed buildings in Amatrice (the largest of the seven towns considered).

Empirical fragility assessment for masonry building classes MBC1-MBC4
The empirical fragility curves and the corresponding plus/minus one standard deviation interval for each class of masonry buildings defined in Sect. solid lines (using established colour codes, i.e., green for D2, yellow for D3, orange for D4 and red for D5). It should be noted that the buildings marked as D0 and D1 are grouped together in this study. Therefore, no empirical fragility curves are reported for D1. The fragility curves are calculated based on the conditional ground-shaking fields generated as described in Sect. 2.3, corresponding to a total set of 1,000,000 simulations. Moreover, the expected value plus/minus one standard deviation fragility curves (Eq. 2) are plotted in dashed lines; the binned damage fractiles are plotted as circles. The observed damage fractiles, which are reported only for the sake of comparison with the resulting fragility curves, are obtained by dividing the PGA domain (ranging from smallest value to the maximum PGA value generated in the ground-shaking fields) into equidistant bins-making sure that each bin contains at least one surveyed building. For instance, the 0.50 damage fractile (shown in full circles) for damage grade D i and a given bin of PGA is calculated as the 50th percentile of the ratio of number of buildings in that bin whose damage grade is equal to or exceeds D i to the total number of buildings in the same bin over all the realizations. The damage fractiles for each bin are positioned at the centre of each bin. It can be observed that the fragility curves obtained through the logistic regression procedure (Sect. 2.4) show a very good agreement with the observed damage fractiles. As mentioned before, those ground-shaking fields that fail to produce monotonically increasing fragility curves are discarded as not plausible (see Sect. 2.3.4). Tables S1 and S2 (reported in the electronic supplement to the paper) summarize the equivalent lognormal median and logarithmic standard deviation (dispersion) for the fragility curves, respectively for Copernicus EMS and visual survey damage surveys.

Comparing the results based on ShakeMap (INGV) and on the proposed GMPE-based conditional ground-shaking fields
In this section, the empirical fragility curves obtained using the proposed conditional GMPE-based ground-shaking fields and the fragility curves obtained based on the INGV ShakeMap (http://shake map.rm.ingv.it/shake /index .html, last access 15/12/2018) are compared for the building classes MBC1-4. It is to note that approach (b) as defined in Sect. 2.3.2 is used for consideration of the stratigraphic amplification. Figure 9 shows the fragility curves obtained using the INGV ShakeMap based on both Copernicus and visual survey and plotted with solid and dashed grey lines, respectively. The figure shows the empirical fragility curves (for D2-D5 and based on the standard colour code) obtained by employing the proposed conditional ground-shaking fields based on Copernicus and visual survey (plotted as solid and dashed lines). The consideration of ShakeMap for estimation of ground shaking (as it is common in the literature) shows significant discrepancies between the fragility curves obtained based on Copernicus and visual survey for MBC1 (D2 to D4), MBC2 (D2 to D3), MBC4(D2). This difference can be attributed to the fact that areal ground shaking values according to ShakeMap consisted of a maximum of six distinct levels of ground shaking. In some cases, differences in damage level attribution has led to significant changes in the number of buildings associated to each of the six ground shaking levels. Consequently, this can lead to significant difference in the resulting fragility curves. On the other hand, the fragility curves calculated based on the conditional GMPEbased random fields are based on point-wise estimates of ground shaking. Therefore, they are going to be less sensitive to variations in damage grade assignments (see Figure S4 in the electronic supplementary material for more detailed analysis of MBC1, D2). In fact, the fragility curves obtained by employing the proposed conditional GMPE-based groundshaking fields based on Copernicus and visual survey show better agreement. The empirical fragility curves based on visual survey (as D Damage ) and ShakeMap (as PGA) seem in closer agreement with those obtained through the proposed procedure (MBC1: D2-D5, MC4: D2). The empirical fragility curves obtained based on Copernicus (as damage survey) and the ShakeMap (for ground shaking estimation) seem to overestimate, with respect to the proposed procedure, the structural capacity for all the damage levels and for all the structural classes studied. In other words, the fragility curves obtained by employing the generated conditional ground-shaking field demonstrate higher vulnerability with respect to the fragility curves obtained based on the ShakeMap (paired up with Copernicus) for all four classes (the equivalent lognormal statistics of the ShakeMap-based fragility curves are reported in the electronic supplement to the paper, Table S3).

Stratigraphic amplification: Bindi et al. (2011) GMPE versus Landolfi et al. (2011)
In this section, the empirical fragility curves obtained based on ground shaking fields PGA that consider the effects of stratigraphic amplification as in Bindi et al. It is to note that, for the fragility curves obtained based on Landolfi et al. (2011) stratigraphic factors are the same as those plotted in Fig. 7. It can be observed that the two approaches adopted for considering the stratigraphic amplification lead to a good agreement. No significant improvement is observed, in terms of the relative predictive capacity of the empirical fragility curves, when the Landolfi et al. (2011) model with the ground shaking-dependent non-linear soil column effect is adopted instead of the GMPE-embedded constant coefficients. More specifically, the use of Landolfi et al. (2011) model leads to a slight reduction both in the median and standard deviation with respect to Bindi et al. (2011) model for all the four building classes considered (with the exception of MBC3 where approach (a) fragility curves demonstrate slightly lower dispersion). Therefore, in lieu of site-specific information on stratigraphic amplification, the use of approach (a) based on ITA10-embedded coefficients is recommended. The equivalent lognormal statistics of the approach (a)-based fragility curves are reported in the electronic supplement to the paper, Table S4.   Fig. 11. It is to note that, the Bindi et al. (2011) GMPE furnishes estimates of the geometric mean between maxima of the two horizontal components of the ground shaking IM g.m. . This is while, the national hazard curves Meletti and Montaldo 2007, Project S1) are expressed in terms of the maximum of one arbitrary horizontal component IM arb . Since the empirical fragility curves obtained in this work are based on ground shaking estimates from the GMPE (i.e., IM g.m. ), it is important to ensure they are consistent with the hazard curves in terms of intensity measure (Baker and Cornell 2006). That is, the risk convolution requires that hazard and fragility curves are expressed in terms of the same intensity measure. Based on a set of working assumptions, Baker and Cornell (2006) demonstrate that the expected value of the logarithm of IM g.m. is the same as that of IM arb . However, the logarithmic standard deviation of IM g.m. is lower than that or equal to that of IM arb , depending on the correlation between the maxima of ground shaking for the two horizonal components. A constant correction factor (estimated in Ebrahimian et al. 2019) has been applied to both inter-event and intra-event components of logarithmic standard deviation of the GMPE to make the fragility and hazard curves consistent (for the purpose of risk convolution). Table 5 reports the annual probabilities of exceeding the different damage levels (risk R). R F denotes the risk calculated based on fragility curve F; R F±σ denote the risk calculated based on fragility plus/minus one standard deviation curves, respectively. Except for MBC2 (D5), MBC3 (D3-D5) and MBC4 (D4-D5), the overall agreement observed between the fragility curves is fine (when compared to Rota et al. 2008 fragility curves with flexible floors). The quantified comparison in Table 5 reveals that the annual probability of damage level exceedance estimates rendered by fragility curves of Rota et al. 2008 and the ones proposed in study are quite different (although some agreement in ballpark can be observed). The following factors can be responsible for the lack of a good agreement: differences in building portfolio characteristics, differences in the definitions of the building classes, the methods used for estimating the ground-shaking intensity, and the methods used for fitting the damage data. Furthermore, the results in Table 5 show an overall agreement between visual survey and Copernicus damage gradings paired up the conditional GMPE-based ground shaking fields (the first two rows of the table). Moreover, the results in Table 5 indicate that the fragility curves rendered by pairing up of visual survey and INGV ShakeMap provide un-conservative estimates with respect to both Rota et al. 2008 and conditional GMPE-based ground shaking fields. This provides evidence for the remarkable sensitivity of the empirical fragility curves to the methods used for the estimation of ground shaking level.

Conclusions
Empirical fragility assessment can reveal significant sensitivity to the definition of the building class, building-to-building variability within a given class, observed damage and ground shaking levels for the considered portfolio of buildings for a given class. In such context, accurate estimation of the ground shaking level at the location of buildings and considering the possible correlations between these ground shaking levels is as important as accurate estimation of the damage. The (Bayesian) concept of robust updated fragilities (a.k.a. predictive fragilities) can be applied to consider the uncertainty in the estimation of the ground shaking level at the position of the buildings of the interest. The joint probability distribution of the ground shaking levels-and hence the spatial correlation structure-can be modelled through simulation of plausible realizations of a random ground motion field. This random field is defined based on the (updated) statistics of the GMPE considering the ground shaking levels recorded at the nearby accelerometric stations. The set of plausible ground motion fields is further filtered as those fields which lead to physically meaningful fragility curves for all damage levels. This work uses damage data obtained based on Copernicus-EMS damage grading maps and complemented by visual survey for buildings damaged due to the Amatrice Earthquake of 24th of August 2016 and the aftershocks immediately following it. Using Copernicus EMS and visual survey after the Amatrice Earthquake of 24 August 2016 as data basis, the paper explores several factors that may severely affect empirical fragility assessment.
Effect of stratigraphic amplification It can be observed that explicit consideration of the non-linear soil column amplification by using (deterministic) class-specific coefficients applied to ground shaking at bedrock level (approach b, Sect. 2.3.2) shows very slight improvement with respect to embedded coefficients of ITA10 GMPE (approach a, Sect. 2.3.2). This could be somewhat expected because, technically speaking, both approaches considered in this work can be classified by the same level of micro-zonation (since they are both based on soil class). Therefore, in lieu of site-specific information on stratigraphic amplification, the use of the embedded coefficients of ITA10 is recommended.

Comparison between Copernicus and Visual Survey
The two approaches seem to be in good agreement in terms of the cumulative probability of exceeding a damage state for damage states of D3 and higher. It can be observed that the two approaches for damage estimation lead to comparable results in ballpark in terms of the annual probability of exceeding a specific damage level. Nevertheless, the Copernicus-EMS database does not provide the building classification; therefore, it needs to rely on complementary information regarding the building classification.
Comparison with ShakeMap While Copernicus and visual survey demonstrate overall fine agreement with each other when the proposed procedure is followed for the estimation of the ground shaking levels, the two damage grading datasets lead in some cases to totally different results when paired up with ShakeMap. This can be attributed to the fact that ShakeMap relies on areal estimates of ground shaking and does not consider the uncertainty in the estimation of ground-shaking at the position of buildings of interest. In the context of the survey performed in this work, pairing up Copernicus damage grades with ShakeMap leads to significant underestimation of vulnerability and risk with respect to the conditional GMPE-based ground shaking fields.
Comparison with literature Comparison with the fragility curves reported in Rota et al. 2008 for analogous classes of masonry buildings, shows agreement in the ballparks of the probability of exceeding various damage levels for buildings without tie rods or ring beams (MBC1-2). The results demonstrated are subjected to the following limitations: • In some cases, the adopted survey techniques might not have assigned the correct building class. For example, the ring beams may not be visible from the building façade. • There might be some bias in the ratio of un-damaged buildings to the total number of buildings surveyed. This can lead to a bias in the survey sample size. The potential bias can be attributed to the nature of the survey that tended to identify damaged buildings. • The survey coverage is not full. It is assumed that the surveyed portfolio of buildings provides a sufficient and un-bias coverage of the portfolio of buildings in the seven towns considered. • There was only one accelerometric stations inside the polygon defining the portfolio of buildings considered (see Fig. 6b). • It is assumed that the damage survey is not subjected to error. However, if there is information available for calibrating the accuracy of damage results, the predictive fragilities can be extended to also consider the error in damage survey. • The choice of the GMPE and the corresponding spatial correlation structure between intra-events residuals can affect the results. The results could potentially improve, if site-specific data was available for the stratigraphic amplification. In other words, the fragility curves would be more accurate if based on higher levels of seismic microzonation. • The empirical fragility parameters are estimated without enforcing the ordinal hierarchy in the damage scale.
As a final word, the proposed method provides practical means for considering the uncertainties in the estimated ground shaking level and correlations in the ground motion at the location of the damaged buildings. The statistical updating considering recorded ground motion constitutes a fundamental step in incorporating the local site conditions and the ground motion levels effectively experienced. After being subjected to sufficient validation and testing, the proposed methods can be considered as a potentially valid alternative to the use of ShakeMap for estimation of the ground shaking level, for the purpose of empirical fragility estimation. Moreover, the proposed method can incorporate damage estimation methods using remote sensing techniques for fast fragility estimation.