Typological fragility analysis of masonry buildings in Emilia Romagna region (Italy)

UnReinforced Masonry (URM) buildings constitute the majority of the Italian building stock and represent one of the most vulnerable construction typologies, as emerged during the seismic events occurred in the last decades. The objective of the present work is the proposal of a methodology for the seismic vulnerability and fragility assessment of different typologies of URM buildings. Through the statistical analysis of post-earthquake damage survey forms and of census data, the most prevalent structural typologies were identified in the Emilia Romagna region (Italy), leading to the definition of two reference masonry buildings, built at the beginning of the XX century and representative of the building stock of the study area. These buildings were modeled according to the equivalent frame method and their seismic performance was studied by means of nonlinear static analyses, considering uncertainties on the mechanical properties of masonry, on the intensity of the seismic action and on the displacement demand calculation method. Through the definition of Response Surfaces, analytical fragility curves relative to the investigated structural typologies were determined for the different damage states.


Introduction
The seismic vulnerability of existing UnReinforced Masonry (URM) buildings is a relevant issue for many European countries, because of the age of their building stock, built without a proper seismic approach and often suffering materials deterioration (Ceroni et al. 2012;Parisi and Augenti 2013). At the same time, the safety assessment of historical constructions is complicated by their geometrical complexity, the variability of material properties, the poor knowledge on past events which might have affected the current condition of the constructions and the lack of specific design codes (Lourenço et al. 2011;Indirli et al. 2013;Penna et al. 2014b). Over time, there has been a progressive evolution of the construction techniques, resulting in a great variety of structural typologies, in terms of configuration in plan and in elevation, of mechanical properties of the masonry, of the type of horizontal structural elements, i.e., slabs and roofs, and of their connections to the vertical masonry walls. As a consequence, the study of the seismic behavior of masonry buildings is associated to a large number of uncertainties (Sousa et al. 2016). Therefore, probabilistic methods are recommended to carry out quantitative assessments of structural safety (Parisi and Augenti 2012). In this framework, the seismic structural safety of existing buildings is often expressed in the form of fragility curves (Buratti et al. 2010), which define the probability of exceeding specific damage levels as a function of the intensity of the seismic ground-motion.
Seismic fragility curves can be obtained by means of numerical simulations (Silva et al. 2014;Rota et al. 2014;Simões et al. 2015;Milosevic et al. 2020;Battaglia et al. 2021), statistical analysis of observational damage data (Rota et al. 2008;Ioannou et al. 2021) or engineering judgment (Jaiswal et al. 2011). To calibrate fragility curves through numerical simulations it is necessary to develop nonlinear models that describe the behavior of structures subject to ground-motions, reproducing their failure modes. Given the large number of uncertainties involved in the problem, rigorous probabilistic models are often unpractical and approximate approaches are often preferred (Milani and Benasciutti 2010;Bracchi et al. 2015).
The development of detailed vulnerability models at a territorial scale requires the identification of different building classes or typologies and the definition of the factors mostly affecting their structural response. Indeed, buildings with similar architectural and structural features, located in similar geotechnical conditions, are expected to have similar performances in the event of an earthquake.
The objective of the present paper is the proposal of a procedure for the evaluation of the seismic fragility of different classes of URM buildings, through the definition of two case studies, representative of a relevant part of the building stock of the Emilia Romagna region. Their structural capacity was defined by applying two different displacement demand calculation methods on the capacity curves, obtained from nonlinear static analyses (pushover). The variability of the seismic action was taken into account by adopting response spectra of recorded accelerograms suitable for the displacement demand estimate. The uncertainties related to the structural behavior of the case studies were considered by assuming the mechanical properties of masonry as random variables and by adopting a Response Surface (RS) model (Faravelli 1989;Rajashekhar and Ellingwood 1993;Franchin et al. 2003): a statistical model which allows to define the dependency of the structural capacity on a set of factors, through simplified analytical relationships (Buratti et al. 2010). Fragility curves were then obtained with reference to four damage states, accounting for uncertainties related to the mechanical properties of masonry, to the seismic action and to the displacement demand calculation method.

Identification of the building typologies
In this work, building types were defined by considering information gathered from two databases including information about the residential URM buildings of the Emilia Romagna region: (i) the 2011 Italian census, by the Italian National Institute of Statistics (ISTAT) and (ii) the post-earthquake damage survey, collected by the Italian Department of Civil Protection after the 2012 Emilia earthquake through the AeDES (Accessibility and Damage detection in the post-Seismic Emergency) form (Baggio et al. 2002). Both databases contain information, for each building, about the adopted structural technology (concrete, masonry, etc.…), the construction period, and the number of floors. In addition, the AeDES forms contain further data about the structural components as well as the level and the extent of damage. In more detail, vertical and horizontal structural elements are classified, the quality of materials is visually assessed and the presence of structural features that might affect the seismic behavior (such as steel ties, bond beams or ring beams, irregularities either in plan or in elevation) are reported. The damage states are classified based on a five levels scale: no damage (DS0), slight damage (DS1), medium or heavy damage (DS2-DS3), severe or very heavy damage, even leading to collapse (DS4-DS5).
A statistical analysis of these combined sets of data can lead not only to the identification of the prevalent structural types in the Emilia Romagna region, but also to the definition of the structural characteristics associated with the most severe damage. To obtain a more uniform dataset, the analysis of the census data was limited to the municipalities located in the area stroke by the 2012 Emilia earthquakes, belonging to the provinces of Bologna, Modena and Ferrara. A total of 110,655 URMbuildings was obtained from the census data, while 29,322 buildings, representing a data subset of the census data, were considered by the AeDES forms, which were filled upon request after the seismic events. The graphs presented in Figs. 1 and 2 show the distribution of the buildings according to the construction period and the number of floors for the two data sources. It can be noticed that the biggest amount of buildings were built before 1919 and that they were mainly characterized by the presence of two or three floors above ground.  < 1919 1919-1945 1946-1960 1961-1970 1971-1980 1981-1990 1991-2000 2001-2005 > 2005  Considering the data collected from the post-earthquake surveys (AeDES), it was possible to investigate the damage distribution for the different construction periods: its graphical representation (Fig. 3a) shows that the newer the buildings, the lighter the damage produced. Moreover, it can be observed that the highest percentages of buildings with the most severe damage (either DS4 or DS5) referred to the construction period prior to 1919. The most frequent case, among the presented data, was represented by buildings built before 1919, with either 2 or 3 floors above ground and damage state DS4. For this combination of parameters, the type of vertical and horizontal structural elements was further investigated (Fig. 3b, c). The number of buildings having a regular or irregular masonry bond pattern is reported in Fig. 3b, while the distribution of the different types of horizontal structural elements is shown in Fig. 3c. The majority of the buildings of the considered sub-category were characterized by a regular configuration in plan (69.8%) and by the presence of light-weight thrusting (40.3%) or light-weight non-thrusting (39.3%) roofing elements.  < 1919 1919-1945 1946-1960 1961-1971 1972-1981 1982-1991 1992-2001 1919-1945 1946-1960 1961-1971 1972-1981 1982-1991 1992-2001  Based on the statistical analysis of the databases here presented, two representative building types were considered in this research for the fragility assessment. They are described in Section 3.

Case study buildings
Two existing buildings, having the characteristics identified by means of the statistical analysis presented in Sect. 2, were selected in this research for the determination of analytical fragility curves. These real buildings were built before 1919 and they are representative of the portion of the Emilia-Romagna building stock dating back to the beginning of the XX century. Their plan view is similar and very regular, as shown in Fig. 4, and in elevation they both feature 2 floors and an attic. The first building (BT1) is a brick masonry building with a regular bond pattern, while the second building (BT2) has stone masonry external walls and brick masonry internal walls. Both buildings have deformable timber slabs and roofing elements.
The two buildings were modelled with the software 3Muri ( Fig. 5), based on the equivalent frame modelling approach, to perform nonlinear static (pushover) analyses 1 3 (Lagomarsino et al. 2013;Penna et al. 2014a). All the geometric data needed in the modelling phase and the data needed for the load analysis of the slabs and the roof were directly derived from surveys on the buildings. With the aim of investigating the seismic behavior of different building typologies, three distinct models were considered based on the characteristics of the selected buildings: -BT1_BM: the first building was modeled by considering all the walls as constituted by clay brick masonry and by adopting a shear failure criterion according to the Mann and Müller's theory, which is appropriate for masonries characterized by a regular bond pattern (Mann and Muller 1980); -BT2_SM: the second building was modeled by considering all the walls as constituted by stone masonry and by adopting the Turnšek-Čačovič's shear failure criterion (Turnsek and Cacovic 1971); -BT2_SM-BM: the second building was modeled by considering the external walls constituted by stone masonry and the internal walls constituted by clay brick masonry. The Turnšek-Čačovič's shear failure criterion (Turnsek and Cacovic 1971) was adopted for all the walls for sake of consistency.
The mechanical properties of masonry, in terms of compressive strength (f m ), shear strength (τ 0 or f v0 ), Young's modulus (E) and shear modulus (G), were assigned to the masonry piers and spandrels according to the recommended ranges provided by the Italian Building Code (Ministero delle Infrastrutture e dei Trasporti 2018, 2019) for existing buildings, reported in Table 1. In more detail, concerning the shear strength, the parameter τ 0 controls the diagonal cracking failure mode according to the Turnšek-Čačovič's criterion (1971), while the parameter f v0 controls the stair-stepped diagonal cracking according to the Mann-Müller's criterion (Mann and Muller 1980). Confidence factors, to be used when dealing with existing buildings to reduce the mean mechanical properties of Table 1 according to the knowledge level, and partial safety factors were considered equal to 1.
To account for uncertainties in the structural behavior, nonlinear static analyses were performed by considering several permutations of the mechanical properties of the masonry walls. The permutations were defined by combining the minimum and the maximum values of the parameters f m , τ 0 or f v0 , E and G ( Table 1), according to the Design of Experiments theory (Box and Draper 1987;Khuri and Cornell 1997) with a two-level factorial design (Box and Wilson 1951). A total of 2 n permutations were obtained for each model, being n the number of independent variables, plus an additional permutation in which all the parameters were set to their average value. The Young's and the shear modulus were considered dependent on each other; therefore, their ratios did not change in the permutations. For the models BT1_BM and BT1_SM, constituted by a unique masonry typology, 9 permutations were obtained, while 65 permutations were obtained for the model BT2_SM-BM, given the presence of two distinct masonry typologies.
As an example of the dynamic characteristics of the investigated buildings, Fig. 6 shows the first modes of vibration in the X and Y directions for the permutation with the average values of mechanical parameters, while the fundamental periods are reported in Table 2.

Nonlinear static analysis
For each model, different nonlinear static (pushover) analyses were performed by considering two different distributions of lateral forces applied along the two main horizontal directions (± X and ± Y) and their corresponding accidental eccentricities, i.e., − 5%, 0%, and + 5% of the building width in each direction. According to the Italian Building Code (Ministero delle Infrastrutture e dei Trasporti 2018), the first lateral load distribution is  proportional to the shape of the fundamental mode of vibration in each direction, while the second is uniform along the height. These distributions will be indicated in the following as modal and uniform, respectively. Combining the 2 different lateral load distributions, the 4 loading directions and the 3 values of accidental eccentricity, a total of 24 capacity curves were obtained for each permutation of the material properties. These curves were then converted into capacity curves for equivalent Single-Degree-Of-Freedom (SDOF) systems whose base shear (F b *) and displacement (d c *) were computed from the capacity curves from the Multi-Degree-Of-Freedom (MDOF) model (F b , d c ), as: where Γ is the modal participation factor for the force distribution considered (Fajfar 1999(Fajfar , 2000. The curves for the equivalent SDOF systems were then approximated by means of elastic-plastic curves. The bilinearization technique was based on energy equivalence principles, while the ultimate displacement was defined, according to the indications of the Italian Building Code for the collapse limit state (Ministero delle Infrastrutture e dei Trasporti 2018, 2019), as the minimum between the displacement associated with 20% reduction of the base shear after the maximum of the capacity curve and the displacement corresponding to the collapse of the masonry piers belonging to any wall considered relevant for the safety of the construction on any floor of the building. In Fig. 7 and Fig. 8, a selection of the capacity curves (for the MDOF model) obtained from the models BT1_BM and BT1_SM is reported, considering the permutations in terms of shear strength, for which the other mechanical parameters were equal to their maximum value. Directions + X and + Y are considered separately. As expected, by increasing the shear strength value, the shear capacity increased. For the model BT2_SM-BM, three sets of capacity curves are presented in Fig. 9, considering the permutations in terms of shear strength for both the stone masonry and the brick masonry. In the selected permutations, all the other mechanical parameters were set at their maximum value. Variations of both shear strength values (i.e., for brick and stone masonry) led to an increased capacity, in terms of maximum base shear and maximum displacement. Results of the pushover analyses showed, in general, that the masonry piers were failing mainly for shear (see Fig. 10 in Fig. 7 Capacity curves for the model BT1_BM with permutations in terms of shear strength f v0 : a + X direction, b + Y direction which representative failure modes are reported) and this was coherent with the effect that a variation in the shear strength had on the capacity curves.

Damage state definition
As evidenced in a recent work (Cattari and Angiolilli 2022), the debate about the definition of damage levels for masonry buildings is still open. Several methods exist, some of which based on a multiscale approach, some others based on the definition of specific thresholds for a unique parameter, e.g. roof drift, interstory drift ratio, etc. (Lagomarsino and Cattari 2015). Previous works, starting from the proposal by Calvi (1999), defined the damage states in terms of specific displacement values on the pushover curves (Cattari et al. 2004(Cattari et al. , 2014Lagomarsino and Giovinazzi 2006). Given the regularity and simplicity of the investigated buildings, this method was here adopted, instead of the multiscale approach, which is more suitable when dealing with large and complex constructions (Lagomarsino and Cattari 2015).
According to the adopted damage scale, an expected structural performance is assigned to each damage state, as described in Table 3. The reported damage states were associated

Seismic displacement demand
Given a capacity curve and a seismic action in the form of an elastic response spectrum, it is possible to estimate the maximum displacement that a structure will attain during a seismic ground-motion, also referred to as performance point. According to both the Italian Building code (Ministero delle Infrastrutture e dei Trasporti 2018) and the Eurocode 8 (CEN 2004), two approaches can be adopted to evaluate the seismic displacement demand. The first accepted methodology is known as N2 Method, and it is based on empirical R-µ-T relationships (Fajfar and Fischinger 1988;Vidic et al. 1994;Fajfar 1999Fajfar , 2000. The (2)  Table 3 Definitions of damage limit states (Calvi 1999) Damage state The building collapses and cannot be repaired either from a practical or economical point of view second procedure is the Capacity Spectrum Method, defined in ATC-40 (Applied Technology Council 1996) and included in FEMA-247 (FEMA 1999). It is based on the definition of an equivalent viscous damping ratio related to the dissipative capacity of the structure under consideration. In the present research, the uncertainty related to the choice of the displacement demand calculation method was considered. Therefore, the two different approaches were employed in the seismic demand computation.

N2 method
The first method adopted for computing the displacement demand was the N2 Method as modified by Guerrini et al. (2017). These authors derived R-µ-T relationships for shortperiod URM buildings by considering two MDOF systems that were representative of two limit-configurations: flexural dominated structures with slender piers and shear dominated structures characterised by squat piers. From both models, equivalent SDOF systems were calibrated so that the hysteretic demand would coincide between the associated MDOF and SDOF systems. The dissipative hysteretic capacity per cycle was quantified by means of the Jacobsen's equivalent damping ratio ξ hyst. (Jacobsen 1930). The SDOF response was then evaluated through several hysteretic models, featuring different slopes of the initial acceleration vs displacement curve, different elastic periods T 0 , masses m and viscous damping coefficient values (ξ vd ). Capacity curves for each model were transformed into equivalent elastic-plastic relationships through equivalence of the underlying areas of the base shear vs displacement diagram, starting from the origin of the axes up to the collapse point. The collapse displacement was identified as described in Section 4.1. The secant stiffness was evaluated in correspondence with the 70% of the maximum base shear, while the corresponding elastic period T * was computed as: where m and k are the mass and the stiffness of the SDOF system, respectively.
(3) The force reduction factor R was instead obtained from the elasto-plastic curve by following the relationship: where d e = a e ·(T/2π) 2 is the elastic displacement demand, a e = S e (T * ) is the pseudo-spectral acceleration associated with the elastic period T * and a damping ratio ξ = 5%, and F e = m·a e is the corresponding elastic force. Then, d y = a y ·(T/2π) 2 is the displacement demand at yielding, a y is the pseudo-spectral acceleration at yielding and F y = m·a y is the associated yielding force. Guerrini et al. studied the response of a variety of SDOF systems by means of nonlinear dynamic analyses, considering many combinations of elastic periods, pseudospectral acceleration or force-reduction factor, and they obtained the formulation for the displacement demand. In particular, for force reduction factor R ≥ 1, Guerrini et al. proposed the following equation for evaluating the inelastic displacement d max : where T hyst , a hyst , b and c are parameters calibrated by Guerrini et al. (Guerrini et al. 2017) by analysing the dynamic response of the first SDOF systems defined.
The Capacity Spectrum Method (CSM), first proposed by Freeman (Freeman 1994(Freeman , 2004, estimates the displacement demand through elastic spectra scaled by viscous damping coefficients associated to the hysteretic dissipation. The comparison between demand and capacity, in this method, was implemented in the Acceleration-Displacement Response Spectra (ADRS) plane. To express the capacity curve in terms of acceleration vs displacement, the relationships presented in Eq. (6) were used: where d c * and F b * are the displacement and the base shear of the SDOF capacity curve, M is the total mass of the building, α 1 is the effective modal mass ratio for the force distribution considered, while φ 1 is the amplitude of the force distribution at the control point. After defining an initial tentative value of the performance displacement d pi , the associated dissipated energy was obtained and the equivalent viscous damping ξ eq was computed as well according to the expression: where k = 0.33 was taken from ATC-40 (Applied Technology Council 1996) with reference to poor existing buildings and short shaking duration (Structural Behaviour Type C). This parameter was used to calculate the reduced elastic spectrum. Then, the intersection between the scaled elastic spectrum and the capacity curve represented the performance point, d pi+1 . If this point was close to the initial guess d pi , then d pi+1 was taken as the , otherwise the iterative process continued until convergence of the solution and d pi+1 was considered as a new tentative displacement.

Definition of seismic actions
In the present work, seismic actions were represented through the response spectra of recorded accelerograms, which were selected based on spectral compatibility with elastic target response spectra with the shape adopted by the Italian Building Code. To find the target spectra for the accelerogram selection, the capacity curves obtained for the models with a unique masonry typology (BT1_BM and BT2_SM) were considered, selecting the permutation in which all the mechanical parameters were set to their average value ( Table 1). By means of the N2 method, they were used to determine the return period of the code spectra that led to displacement demands equal, on average, to the displacement capacities associated to each of the four damage states under consideration (see Section 4.2). The elastic spectra in Fig. 12 represent the target for the accelerogram selection procedure, with values of the return period reported in the legend. It was checked that the target spectra obtained for the model BT2_SM could be used as target spectra also for the model BT2_SM-BM. Clearly the target spectra are characterized by return periods much higher than the ones which are usually considered in design applications, because, on one hand, the damage states considered are more severe that those conventionally adopted by building codes and, on the other hand, mean values of material strengths are used. The procedure adopted allowed to identify response spectra that are consistent with the uniform hazard spectra while limiting the need for scaling. Spectral compatibility was searched for a range of periods including the elastic and secant periods of vibration of the studied structures, specifically from 0.1 to 0.7 s. Recorded accelerogram response spectra were selected from the PEER Ground motion database (Pacific Earthquake Engineering Research Center), precisely referring to the NGAWest2 project. To account for variability in the seismic action, 40 accelerograms were selected for each damage state and for each model. The following criteria were introduced in the selection: (i) restrictions in terms of soil category were considered, selecting a range of V S,30 values that allowed to include only recordings on soil class B and C (i.e., 180 m/s ≤ V S,30 ≤ 800 m/s), as defined by the Italian Building Code; (ii) Fig. 12 Target elastic spectra for four damage limit states: a BT1_BM, b BT2_SM and BT2_SM-BM impulsive records were excluded; (iii) a maximum of 4 records per earthquake was used. After this preliminary selection, the remaining accelerograms were sorted according to the value of their average root-mean squared deviation (RMSD) from the target spectrum and the 40 accelerograms with the lowest RSMD value were considered. The accelerograms selection was carried out by using the PSARotD50, as defined by Booree et al. (2006), namely the 50th percentile of the values obtained from the rotation of the acceleration horizontal components from 0° to 180°. As an example, in Fig. 13 the accelerograms selection for the model BT1_BM is reported for each damage state: the grey and red curves represent the spectra from the selected accelerograms and the average ones, respectively, while the black curves represent the target spectra ( Fig. 12a). The complete list of the selected accelerograms for all the models and damage states is reported in the Appendix: for each registration, the earthquake name and year are reported together with the recording station name, the moment magnitude, the Joyner-Booree source-to-site distance and the value of V S,30 .
The pseudo-acceleration spectra for two orthogonal horizontal directions were used in the analysis. Recordings, however, came from different sources and events and thus the horizontal components did not always correspond to the N-S and E-W directions, which was the most common notation. To avoid confusion, in the following paragraphs, the components will be referred to as H 1 and H 2 as to identify two generic perpendicular horizontal components. Due to the unknown orientation of the typological buildings with respect to the two acceleration components of the seismic Fig. 13 Accelerogram selection related to the BT1_BM model for the damage states: a D 1 , b D 2 , c D 3 , d D 4 1 3 actions, both components were adopted for the definition of the structural capacity, as described in Section 6.1.

Definition of the structural capacity
To carry out the fragility assessment of the typological structural models, it was first necessary to evaluate the ground-motion intensity ( PGA (j) i ) associated to the attainment of the j-th damage state for the i-th accelerogram. These PGA values, referred to as structural capacity, were computed by means of an iterative procedure. In particular, by finding, for the response spectrum of the i-th accelerogram, the scaling factor SF (j) i that led to a displacement demand equal to the displacement capacity D j , associated to the j-th damage state. Then the structural capacity in terms of PGA was computed as: where PGA i is the PGA of the i-th accelerogram under consideration. This procedure was repeated for each permutation of the values of the mechanical properties, for each pushover analysis (i.e., for the different force distributions, directions and eccentricities) and for each of the two horizontal acceleration components of the ground-motion records and for both the displacement demand calculation methods (Guerrini et al. (2017) and CSM). A limit was imposed on the values of the scaling factors; namely, all the accelerograms associated to a scaling factor SF (j) i lower than 0.5 or greater than 6 were neglected in the successive steps of the present work.

Response surface method
Response Surface (RS) models were used to define empirical relationships between structural capacity in terms of PGA and the random variables associated to the mechanical properties of the case study buildings. RSs were fitted to the observed data, i.e., the values of the structural capacity obtained from the different analyses carried out, aggregating the PGA (j) i values for both the components of all accelerograms selected for the damage state D j under consideration. Additionally, given the significant variability of the response variable, PGA (j) i , homogeneous k-th data groups were considered to define the RSs, i.e. separating the results according to the load distribution (either uniform or modal), loading direction (± X or ± Y) and displacement demand calculation method (Guerrini et al. (Guerrini et al. 2017) or CSM) considered in the analyses. A total of sixteen RSs for each damage state of each model was obtained. The RS functional forms adopted for the k-th data group of the models BT1_BM, BT2_SM and BT2_SM-BM were, respectively: , c 1 (j,k) , …, c 7 (j,k) coefficients, associated, for each model, to the k-th data group and to the j-th damage state, and the standard deviation of the error term were estimated by means of linear regressions. Sections of the Response Surfaces for the model BT2_SM are presented in Fig. 14. Among the sixteen RSs calibrated for this model, these sections were obtained by considering the data group with a uniform load distribution, a loading direction + X and the Guerrini et al. Method. In more detail, Fig. 14a, b show the effect of a variation in terms of shear strength on the response, i.e., PGA for the damage states D 1 and D 4 , respectively. These sections were evaluated by fixing the compressive strength of the stone masonry to its maximum value and by varying the Young's modulus, considering the minimum and the maximum values. It can be observed that a variation in terms of shear strength influenced the response, i.e., lower shear strength values corresponded to lower PGA values, and that this effect was more significant for the damage state D 4 , as expected. Sections of the Response Surfaces for the model BT2_SM-BM, accounting for the presence of two masonry typologies, are also reported in Fig. 15 for the combination with a uniform distribution, a loading direction + Y and the Guerrini et al. Method. The effects on the structural response given by variations in terms of shear strength for both the materials, i.e., stone masonry (Fig. 15a) and clay brick masonry (Fig. 15b), are here presented for the damage state D 4 only, and it was observed that they did influence the PGA values. In more detail, the effect on the structural response seemed more evident for the clay brick masonry, but it should be observed that a wider range for the shear strength was considered with respect to stone masonry, according to Table 1.
The RS models provided analytical estimates of the values of the structural capacity in terms of PGA, as a function of the random variables associated to the material properties, representing also the uncertainty related to the structural response to different groundmotions. Therefore, the RS models can be efficiently adopted to compute the exceedance probability of a damage state D j as a function of the seismic demand, defined in terms of PGA.

Derivation of fragility curves
The random variables used to represent uncertainties on material properties were assumed to follow lognormal distributions, with the parameters reported in Table 4. These were derived from the CNR DT212/2013 Guidelines (CNR-DT212/2013 2014), while the values referring to fv 0 were obtained in a previous research (Ferretti et al. 2016(Ferretti et al. , 2019, devoted to the investigation of the mechanical properties of existing masonry buildings, similar to the ones here considered and located in the Emilia Romagna region. Once the unknown coefficients in Eqs. (9-11) were determined, it was straightforward, for the k-th data group of each model, to compute the mean and the variance of the PGA capacity for the D j damage state. This was done by considering Eq. (12), Eq. (13) or Eq. (14) for the models BT1_BM, BT2_SM and BT2_SM-BM, respectively: Assuming that the natural logarithm of the PGA capacity is normally distributed with mean ln(PGA (j,k) i ) and standard deviation ln(PGA (j,k) i ) , the probability of exceedance for the k-th data group of each model and for each damage state D j was computed as:

Fragility curves
Depending on the loading direction, the force distribution, and the displacement demand calculation method, sixteen sets of fragility curves were obtained for the considered models. Few examples of the fragility curves obtained using the Guerrini et al. method (14) ln(PGA  (Guerrini et al. 2017) are presented in Figs. 16 and 17 for the model BT1_BM and BT2_ SM, respectively. The influence of the force distribution and of the loading direction on the fragility curves for brick and stone masonry and for the damage states D 2 and D 4 was evidenced. The most significant differences were observed for the damage state D 4 . For this representation, positive and negative X and Y directions were considered together.
In the present paper, the approach proposed by Shinozuka et al. (2000) to combine fragility curves was adopted. In particular, being M the total number of k-th data groups, it was considered: Given that the k-th data groups were equally populated, P h = 1/M was the weight assigned to each group. As suggested in the cited work (Shinozuka et al. 2000), the combined fragility curves may be assumed to be lognormal with respect to the mean ln(PGA (j) ) and the variance 2 ln(PGA (j) ) can be estimated on the basis of the variances of the single fragility curves, see Eqs. (12)(13)(14)16).
The described approach was first applied considering separately the fragility curves obtained from the two displacement demand computation methods (M = 8). In this way, combined fragility curves were obtained, two for each model, independent from the loading direction and the shape of the load distribution. As an example, the fragility curves for the model BT1_BM are reported in Fig. 18. It can be observed that the two sets of curves   Fig. 19 for the model BT1_BM at the damage limit state D 4 , for which, in most cases, the lowest scale factors were obtained through the Capacity Spectrum Method. Similar observations can be coherently applied to all other damage limit states and to all the structural models. Therefore, it can be highlighted that the PGA capacity computed from the Capacity Spectrum Method, according to Eq. (8), was systematically lower compared to the same value obtained by implementing the Guerrini et al. method. The displacement demand calculation method, as observed, represent one of the different sources of uncertainties involved in the problem. Then, with the objective of obtaining fragility curves independent from the displacement demand calculation method, as well as from the loading direction and the shape of the load distribution, the Shinozuka's approach (Shinozuka et al. 2000) was eventually applied combining the fragility curves obtained for all the k-th data groups of each model, for each damage state (M = 16). The parameters (μ ln , σ ln ) of the combined fragility curves are reported in Table 5, while the curves are presented in Fig. 20 for the three typological buildings. Figure 21 shows an attempt of comparing the fragility curves with others obtained in the literature, both analytically and empirically, for similar building typologies. To this aim, the curves obtained for the model BT1_BM (i.e., brick masonry typology) are considered; the other models are not considered because their features do not find a direct match in the  literature. In the comparison, the analytical fragility curves were obtained for unreinforced clay brick masonry buildings, characterized by a regular plan and elevation configuration and by the presence of rigid concrete floors (Cattari et al. 2014;Lagomarsino and Cattari 2014). The definition of damage states (D 1 to D 4 ) was compatible with the one adopted in the present work with the addition of further drift limits. Differences between the fragility curves can be observed especially for the damage state D 2 and can be attributed to differences in the structural configuration of the considered buildings (e.g., deformable vs rigid slab) and to possible discrepancies in the identification of the different damage states, in particular due to drift limits. The empirical curves were derived from a recent work (Ioannou et al. 2021), focused on RC and masonry buildings of the Emilia-Romagna Region. In this case, the comparison was only made for damage states D 1 and D 3 , since limited observational data were obtained for D 4 and D 5 in the cited work. It can be noticed that analytical models feature an higher fragility than empirical models, this is due to the criteria adopted for the definition of the damage states and to the modeling assumptions. Actual buildings normally have a higher dissipative capacity than numerical models because these latter introduce various simplifications of the mechanical behavior of materials and structural elements and do not consider the contribution of non-structural elements.

Conclusions
This paper presented a methodology to analytically assess the fragility of typological unreinforced masonry buildings. To this aim, two case study buildings, representative of part of the building stock of the Emilia Romagna region, were identified: their construction details and geometry were selected from a statistical analysis of both the 2011 census data and the post-earthquake surveys database from the 2012 Emilia earthquake. Based on the selected case studies, three tri-dimensional models were considered to investigate the behaviour of different building typologies, i.e., with stone masonry, brick masonry or both. The capacity of these structures was computed by means of pushover analyses and damage states were identified according to literature indications. Moreover, several permutations of the mechanical parameters were implemented to account for the variability of the masonry mechanical properties, according to the Design of Experiments theory. Local failure modes were not analysed.
To consider the variability of the seismic action, ground motion recordings were selected from the PEER Ground motion database. The displacement demand was then calculated by means of two different displacement demand calculation methods, derived both from building codes and available literature: the N2 Method, as modified by Guerrini et al. (2017), and the Capacity Spectrum Method. Response Surfaces were fitted to the observed data to estimate the structural capacity, in terms of PGA, as a function of the mechanical properties of masonry, assumed as random variables. It was observed that the parameter mostly affecting the structural response was the masonry shear strength, coherently with the failure modes obtained in the nonlinear static analyses.
Fragility curves were then obtained for homogenous groups of data, aggregating the results, in terms of PGA capacity, according to the displacement demand calculation method, the load distribution and the loading direction. These sets of curves were then combined, obtaining overall lognormal fragility curves for each building and for each damage limit state. A critical comparison with analytical and empirical fragility curves was also presented for the brick masonry typology.
In conclusion, the proposed methodology allowed to thoroughly consider the uncertainties related to the mechanical properties of masonry as well as to the seismic action and the displacement demand calculation method. Based on the use of Response Surface to perform the fragility assessment, it can be successfully applied for the determination of typological fragility curves for URM buildings, as demonstrated by its application to the case studies presented in this research.