Probabilistic Seismic Assessment of Existing Masonry Buildings

The evaluation of seismic performance of existing masonry buildings is a critical issue in assessing the seismic vulnerability of the built environment. With this aim, non-linear static analysis is commonly used, but results are influenced significantly by the collapse criteria adopted, as well as by the assumptions about material properties and drift capacity of masonry walls. A methodology for the probabilistic assessment of the seismic risk index is proposed by means of an original non-linear pushover type algorithm developed by the authors. The main sources of uncertainties related to masonry parameters and their influence on seismic risk indices are identified by means of sensitivity analysis. Response surfaces for the seismic risk indices are thus defined through general polynomial chaos expansion in order to quantify the uncertainties in the resulting seismic risk index. Finally, a seismic performance classification is presented to help stakeholders to manage risks and define priorities for seismic retrofit. The methodology together with the outcomes is illustrated for a set of existing masonry buildings that are part of the school system in the Municipality of Florence.


Introduction
Seismic vulnerability of the built environment has been confirmed and emphasized by the heavy consequences of recent earthquakes. For this reason, the evaluation of seismic performance of existing buildings is increasingly relevant, especially in historical towns characterized by constructions mostly erected without following any specific seismic provisions.
In order to set-up intervention strategies as well as to identify priorities for seismic retrofit and strengthening, also in view of future planning and optimal allocation of available resources in seismic regions, city planners and owners of buildings need a suitable classification of construction in terms of seismic vulnerability. In this context, masonry buildings are topical since they may cover not only relevant and strategic buildings, such as schools, public buildings, hospital facilities and even large industrial buildings, but also include a large part of the architectural heritage to be preserved.
The paper focuses on the development of a probabilistic methodology for the seismic risk assessment and management of existing masonry buildings.
In the framework of an important research project regarding the seismic assessment of more than 80 school masonry buildings in Florence (I), funded by the Municipality of Florence, the authors developed an appropriate methodology in line with the most updated structural codes [1] to classify the seismic performance of masonry buildings, aiming to define priorities in seismic upgrading as well as to identify suitable retrofitting techniques.
Assessment of the seismic performance of existing masonry buildings is the subject of several up to date studies concerning both the modelling strategies and the evaluation of masonry properties.
Indeed, understanding and modelling the mechanical behavior of masonry elements is a longstanding challenge in civil engineering [2,3].
Depending on the field of application and the level of complexity, several modeling strategies can be envisaged, like macro-modeling, which is the most common method in engineering practice, simplified micro-modeling, and detailed micro-modeling, where masonry units and mortar are modelled explicitly.
Whatever the approach adopted, in any case results of the assessment are very sensitive to the masonry's mechanical parameters [4]. Since mechanical parameters of similar masonry typologies can vary in a very wide range; it is not surprising that contrasting conclusions can be derived about the assessment of a given construction simply modifying the assumed values of their relevant mechanical properties. Moreover, despite several "standardized" masonry typologies that can be identified in existing buildings, depending on the material, on the block shape and on its texture, the evaluation of relevant mechanical parameters associated with them is often an open question.
Experimental in situ tests are also limited and characterization of masonry properties is often obtained by visual inspection methods [5] combined with engineering judgment.
A first attempt to consider the uncertainty of masonry properties in the analysis of existing masonry buildings is discussed in [6] where a Bayesian updating procedure is proposed for the calibration of material properties to be used in the structural analysis.
In this study, a method for the analysis of masonry structures is implemented in a probabilistic framework in order to take into account the uncertainties related to the modelling of the shear behavior of masonry walls, aiming to provide a probabilistic assessment of the seismic performance.
Starting from preliminary historical inquiries and in situ tests and surveys, the proposed procedure relies on an ad hoc developed original, reliable and robust, method for the non-linear static seismic analysis of masonry building [7], as summarized in Section 2. Section 3 describes the 11 masonry buildings considered in the present study, while in Section 4 the main sources of uncertainties related to the shear behavior of masonry walls are discussed with special reference to a sensitivity analysis concerning the propagation of uncertainties in the definition of seismic risk index. A probabilistic seismic risk assessment is then proposed leading to the seismic performance classification presented in Section 5.

Assessment of the Seismic Performance of Existing Masonry Buildings
The structural assessment of existing buildings and infrastructures requires the preliminary acquisition of all the relevant information regarding the construction and the action that it withstands during its life [8]. Inter alia, to define the actual performance of the construction, particularly significant are data concerning the original configuration and the static scheme of the structure, the sequence of structural modifications, alterations and interventions that occurred during its life, the actual material properties, diagnosis and prognosis of existing crack patterns, failures, settlements extent of degradation and damage.
Focusing on old masonry structures in seismic zones, it must be stressed that they were often built disregarding any seismic design rules and without particularly technical competent workmanship; in Italy, for example, this situation was very common at the end of the 1970s. For this reason, masonry buildings are often highly vulnerable, thereby explaining why they frequently collapse during strong earthquakes.
Furthermore, to estimate the actual performance of existing masonry buildings, specialized methods, based on appropriate mechanical models, should be adopted for their static and seismic analysis.
One of the first methods for the seismic resistance verification of unreinforced masonry building named POR and considering floors infinitely rigid in the horizontal plane, was introduced by Tomažević after the Friuli earthquake in 1976 [9,10]. This method was adopted by the Italian regulations published in 1978 [11] and 1981 [12] for the repair of masonry buildings after earthquake events.
For a consistent assessment of the seismic performance, robust analysis methods but also a proper characterization of the behavior of masonry walls under horizontal loads are needed. Commonly, beside local analyses, non-linear seismic static analyses are performed, supplemented, if necessary, by dynamic linear analyses, mainly devoted to check the order of magnitude of the results. In non-linear static analysis, masonry buildings are mostly modelled by means the so-called equivalent frame model [13], but the resulting structural scheme is usually very complicated and the analysis needs to be "driven" step by step in order to obtain reliable results.

The E-PUSH Program
As already said, an original method for non-linear static seismic analysis of mono or multi-story buildings masonry building has been developed by the authors and is presented in [7]. The method, which allows taking into account the ductility of the whole structure, recovers some basic assumptions of the classical POR method, but without its limitations and inaccuracies. In addition, the so called E-PUSH software program requires very simple structural models, so that, by contrast with commercial pushover programs, the input results in being extremely robust and nearly independent of the user. As usual, floors are considered, in turn, "flexible" or "rigid", setting their in-plane stiffness to zero or infinity, respectively, according its deformations in the horizontal plane.
Basic assumptions of the method are: 1. only lateral stiffness of masonry walls is taken into account, disregarding transverse stiffness; 2.
horizontal force-lateral displacement diagram of the wall is elastic-plastic, with the plastic plateau limited by the elastic drift δ e and the ultimate drift δ u . 3.
verification is performed in terms of seismic capacity and demand using the acceleration-displacement response spectra (ADRS) taking into account the story displacements.
The step by step procedure considers increasing values of the given horizontal forces, which are distributed in the plan layout to the resistant shear walls according their stiffness. In short, at each step of the procedure, the inter-story displacement δ of each shear wall is evaluated and three alternative conditions are checked: i.
δ < δ e , the wall is still in the elastic range; ii. δ e ≤ δ ≤ δ u , the wall is in the plastic range, it achieved its shear force resistance and its equivalent stiffness is reduced; iii.
δ > δ u , the wall sustains only vertical loads and its shear resistance and its lateral stiffness are set to zero.
The lateral forces are increased until the base shear resistance reduces to about 80% of the relative maximum base shear resistance or when there is the collapse of all the walls pertaining to a same floor, defining in this way the capacity curve of the whole structure.
The verification procedure, which follows the classical steps of a pushover analysis, is based on the so-called N2 method developed by the University of Ljubljana, described in [14]. The N2 method combines the non-linear static analysis of the multi-degree of freedom model with the response spectrum analysis of an equivalent single degree of freedom system and in the proposed procedure has been implemented according to the following steps: 1.
the non-linear capacity curve of the structure is transformed in an equivalent bi-linear elastic-plastic curve as described in [14]; the curve is characterized by the maximum force F y,eq , calculated averaging the maximum base shear in the non-linear capacity curve and the base shear corresponding to the attainment of the elastic limit; by the yield displacement δ y,eq , evaluated as the ratio between F y,eq and the effective stiffness of the structure K E ; and by the ultimate displacement δ u,eq which is the maximum displacement in the capacity curve; 2. the bi-linear, force-displacement capacity curve, F − δ, of the multi-degree of freedom (MDOF) system is converted in an acceleration-displacement, S a − S d , capacity diagram, for an equivalent single degree of freedom (SDOF) system, according to the following formulae where Γ is the mass participation factor given by and m * is the mass of the equivalent SDOF system.
In Equation (3) m j are the story masses and φ j are the normalized displacements in the considered direction. The elastic period T * of the idealized bi-linear system is thus determined as: where S d,y is the yield displacement and S a,y is the corresponding acceleration.
3. in order to obtain the demand diagram, the elastic design spectrum defined in the Italian Building code [15] from the standard pseudo acceleration-natural period, S ae − T is converted into the pseudo acceleration-displacement format S ae − S de through: 4. the capacity spectrum and demand spectrum curves are plotted in the same graph, to define displacement demand. If the capacity curve intersects the demand curve, the displacement demand is assumed equal to the intersection point d t . Otherwise, the displacement demand is determined starting from the intersection of the radial line corresponding to the elastic period T * of the structure with the elastic design spectrum defining the acceleration demand S ae (T * ). A reduction factor R µ is defined as the ratio between the acceleration demand S ae (T * ) and the yield acceleration S a,y.
Said T C the characteristic period of the ground motion, if T * ≥ T C the ductility demand µ, defined as the ratio between the displacement demand and the yield displacement, µ = d t /S d,y , is equal to R µ , so that if instead T * < T C , the ductility demand µ and the displacement demand d t are calculated according to the formulas given in [10]: Buildings 2019, 9, 237 5 of 20 5. the seismic performance of the structure is evaluated comparing the displacement demand d t with the ultimate displacement defined by the capacity curve d c .
The procedure is iterated for different reference period of the seismic demand until it results d t = d c , finally allowing the evaluation of the seismic risk index I R of the structure, where PGA C is the peak ground acceleration resisted by the structure, i.e., when d t = d c , and PGA D the design peak ground acceleration provided by the Italian Building Code [15]. An example of the output of the proposed procedure in terms of ADRS for the two relevant directions of the seismic force is reported in Figure 1, referring to one of the school masonry buildings in Florence, described in the following.
where is the peak ground acceleration resisted by the structure, i.e., when = , and the design peak ground acceleration provided by the Italian Building Code [15].
An example of the output of the proposed procedure in terms of ADRS for the two relevant directions of the seismic force is reported in Figure 1, referring to one of the school masonry buildings in Florence, described in the following.
The curves in Figure 1, derived from the E-PUSH output for the given case study, include: • the bi-linear capacity curve of the structure (red solid line); • the design spectrum (blue solid curve) for a return period of 712 years ( ), referring to the ultimate limit state for life safety ( ) of occupants in a given location (Florence Municipality in this case); • the design inelastic spectrum (blue dashed curve); • the displacement demand (scarlet dashed line) to be compared with the ultimate displacement (green dashed line); • the elastic spectrum (green solid curve) for the return period TRC consistent with the capacity of the structure ( = ) and the corresponding inelastic spectrum (green dashed curve). The inelastic spectra obtained with the N2 method are based on the commonly adopted "equal displacement rule" for structures in the medium and long period range, * ≥ , and on the procedure proposed by Vidic et al. [16] for the short period structures. The equal displacement rule has been used quite successfully for many years [14] and many statistical studies [16][17][18] have confirmed its applicability in case of structures on firm site, characterized by an average shear wave velocity greater than 180 m/s, with fundamental period in the medium and long period range. For structures characterized by lower fundamental period, the equal displacement rule does not work and < . Formulae for were calibrated by Vidic et al. in [16] based on a statistical study, and then adopted in the N2 method, assuming the transition period equal to , thereby disregarding its dependency on the ductility . Anyhow, in case the ductility demand is < 4 this simplified assumption leads to conservative results. The curves in Figure 1, derived from the E-PUSH output for the given case study, include: • the bi-linear capacity curve of the structure (red solid line); • the design spectrum SLV e (blue solid curve) for a return period of 712 years (T RD ), referring to the ultimate limit state for life safety (SLV) of occupants in a given location (Florence Municipality in this case); • the design inelastic spectrum SLV a (blue dashed curve); • the displacement demand d t (scarlet dashed line) to be compared with the ultimate displacement d c (green dashed line); • the elastic spectrum (green solid curve) for the return period T RC consistent with the capacity of the structure (d t = d c ) and the corresponding inelastic spectrum (green dashed curve).
The inelastic spectra obtained with the N2 method are based on the commonly adopted "equal displacement rule" for structures in the medium and long period range, T * ≥ T C , and on the procedure proposed by Vidic et al. [16] for the short period structures. The equal displacement rule has been used quite successfully for many years [14] and many statistical studies [16][17][18] have confirmed its applicability in case of structures on firm site, characterized by an average shear wave velocity greater than 180 m/s, with fundamental period in the medium and long period range. For structures characterized by lower fundamental period, the equal displacement rule does not work and R µ < µ. Formulae for R µ were calibrated by Vidic et al. in [16] based on a statistical study, and then adopted in the N2 method, assuming the transition period equal to T C , thereby disregarding its dependency on the ductility µ. Anyhow, in case the ductility demand is µ < 4 this simplified assumption leads to conservative results.
The influence of the mathematical model adopted for damping and the level of damping ratio in the definition of the R µ factor has been investigated by Vidic et al. in [16]. It is shown that R µ increases with decreasing damping coefficient; for example, when ξ decreases from 5% to 2%, R µ increases of about 10% considering the mass-proportional model and of about 20% for instantaneous-stiffness-proportional damping. Consequently, when ξ < 5%, adoption of R µ factor corresponding to 5% damping is safe-sided.
In structural codes the response spectra are usually based on a conventional value ξ = 5% of the viscous damping ratio; damping effects, if any, are taken into account suitably scaling the response spectrum. In the already cited Italian Building Code [15], in agreement with corresponding clauses of Eurocode 8, EN1998-1 [19], the damping correction factor η as a function of ξ is given by: Under earthquake excitations, damping increases as the level of damage increases and values of ξ up to 10% to 20% can be reached as discussed in [20].
To analyze the effects of varying viscous damping coefficients on the seismic risk index, a parametric investigation has been carried out for the given case study. In the sensitivity analysis, the response spectra, evaluated according the Italian Building Code, have been modified by varying ξ in the range 5% to 20% by 2.5% steps. The results of the sensitivity study are diagrammatically illustrated in Figures 2 and 3, for the two relevant directions, x and y, of the seismic excitation. More precisely, Figure 2, analogous to Figure 1, synthetizes in the ADRS plane results for different values of the viscous damping coefficient ξ (ξ = 5% ; ξ = 10% ; ξ = 15% ; ξ = 20%), while Figure 3 shows the variation of seismic risk indexes as a function of the viscous damping coefficient ξ. In the case of ξ = 5%, the outputs in Figures 2 and 3 obviously correspond to those already illustrated in Figure 1. Figure 3 also indicates the correction factors η obtained introducing in Equation (10) the effective viscous damping coefficients ξ e calculated according the formulation proposed in [20]: where δ y is the displacement of the structure at the first crack and δ u its ultimate displacement. Evidently, according to Equation (11), ξ e values can be different in the two directions, like in the present case study, where ξ e,x ≈ 17.1% and ξ e,y ≈ 14.1%, so that also correction factors are different in the two directions. In the present case, Equation (10) provides η ≈ 0.67 in the x direction and η = 0.72 in the y direction.
As expected, the seismic demands are reduced and a significant increase in seismic risk indexes is obtained considering the reduced damping correction factors η: in fact, the seismic risk index in the x direction becomes bigger than one, increasing by about 46%; while in the y direction it increases by about 41%.
It must be considered that such high values of the damping coefficients correspond to heavy damages of the masonry walls: most likely, they can be achieved by modern brick masonry walls characterized by high ductility, but for historical or ancient masonry walls, the assumption of such high values of damping should be adequately justified.     An analytical alternative formulations of damping coefficient is provided, for example, in [21] where, on the basis of some experimental results, ξ e is expressed as a function of the ductility µ by: Assuming a ductility factor µ equal to 1.5, as recommended for unreinforced masonry in the guidelines for the application of the POR method [12], in the case study, Equation (12) leads to an equivalent damping of ξ e ≈ 8.4%. This damping value, leading to a correction factor η ≈ 0.86, appears to be more consistent with the actual behavior of historical masonry and leads to results in terms of seismic performance more in line with the usual conservative assumption ξ e = 5% (see dashed green line in Figure 3). In this case, the increment of the seismic risk index is about 14.8% in the x direction and about 15.3% in the y direction.

Validation of the E-PUSH Program
The proposed pushover analysis method has been validated, focusing on the ability to reproduce the seismic response of masonry building prototypes tested in the laboratory and for which experimental load displacement curves were available. An example of these structures is the full scale two-story unreinforced masonry (URM) building prototype tested at the Georgia Tech laboratory [22].
The prototype was a two-story unreinforced brick masonry structure (see Figure 4) with timber floor and roof diagrams; the in plan dimensions of the building were 7.32 m by 7.32 m and the story heights were 3.6 m and 3.54 m for the first and the second story, respectively.  During the experimental campaign, two series of in-plane wall tests were carried out (one parallel to Walls 1 and 2; and one parallel to Walls A and B) controlling the horizontal displacements, and adopting a displacement profile based on the first vibration mode. A complete description of the structure and the loading tests can be found in [22].
The masonry structure has then been modeled and analyzed with the E-PUSH program in order to obtain the capacity curves to be compared with the experimental ones presented in [22]. In Figure 5, the 3D layout of the shear resistant walls modelled in E-PUSH is illustrated.
The experimental force-displacement ( − ) curves obtained for walls A, B, 1, and 2 have been compared with those evaluated analyzing the two-story URM structure with the E-PUSH program.
The results are illustrated in Figure 6 for walls A and B, and in Figure 7 for walls 1 and 2. The satisfactory agreement substantially confirms the reliability of the proposed approach. The building was composed of four URM masonry walls connected at the corners, labeled Walls A, B, 1, and 2 in [22], designed with different thickness and opening ratios to represent typical masonry walls.
During the experimental campaign, two series of in-plane wall tests were carried out (one parallel to Walls 1 and 2; and one parallel to Walls A and B) controlling the horizontal displacements, and Buildings 2019, 9, 237 9 of 20 adopting a displacement profile based on the first vibration mode. A complete description of the structure and the loading tests can be found in [22].
The masonry structure has then been modeled and analyzed with the E-PUSH program in order to obtain the capacity curves to be compared with the experimental ones presented in [22]. In Figure 5, the 3D layout of the shear resistant walls modelled in E-PUSH is illustrated.   Results fit very satisfactorily both in terms of ultimate resistance, as well as in terms of ultimate displacement, in all considered case studies. Detailed comparisons of the outcomes for two relevant The experimental force-displacement (H − δ) curves obtained for walls A, B, 1, and 2 have been compared with those evaluated analyzing the two-story URM structure with the E-PUSH program.
The results are illustrated in Figure 6 for walls A and B, and in Figure 7 for walls 1 and 2. The satisfactory agreement substantially confirms the reliability of the proposed approach.   The E-PUSH program has been also tested on real case studies, such those presented in the next section, and for which no experimental curves are available, comparing the outcomes of the algorithm in terms of capacity curves with those obtained by means of a widely used commercial pushover analysis software, Aedes PCM [23].
Results fit very satisfactorily both in terms of ultimate resistance, as well as in terms of ultimate displacement, in all considered case studies. Detailed comparisons of the outcomes for two relevant case studies are reported in [4] and [7], referring to a two-story and a four-story masonry building, respectively.  Results fit very satisfactorily both in terms of ultimate resistance, as well as in terms of ultimate displacement, in all considered case studies. Detailed comparisons of the outcomes for two relevant case studies are reported in [4] and [7], referring to a two-story and a four-story masonry building, respectively.

Challenges in the Assessment of the Seismic Performance
As mentioned before, a proper evaluation of the seismic performance requires a suitable definition of masonry mechanical parameters, duly considering that masonry is a non-homogeneous material characterized by an inelastic and anisotropic behavior. The identification of the appropriate force-displacement response of in-plane loaded walls is then particularly crucial as it significantly influences the results [4]. Moreover, depending on the material as well as on the texture, several masonry types can be detected in existing buildings, characterized by mechanical parameters varying in a very wide range. Anyhow, stiffness identification [24,25] and evaluation of drift capacity of in-plane loaded walls [26,27] are still critical issues [28].
Flat jacks, diagonal compression and shear compression tests, in situ or in laboratory conditions, are the subject of a wide literature concerning the shear resistance of masonry walls. But their outcomes are often contradictory; in fact, the values of mechanical parameters depend on the test procedure and may vary significantly for the same class of masonry or even for the same wall [24,25]. In this respect, guidelines for application of the Italian Building Code [29] provide some guidance, about the classification of masonry in terms of quality and typology.
The probabilistic seismic assessment procedure proposed here has been applied to 11 case studies suitably chosen among the masonry buildings under investigation. For the sake of the study, masonry has been mainly classified based on visual inspection and double flat-jack compressive test results, duly taking into account epistemic and aleatoric uncertainties affecting the force displacement response of in-plane loaded walls.

Case Studies
The aforementioned set of buildings, numbered from B1 to B11, are representative of a rather long historical period, since they were built between 1900 and 1980. The set includes: 3 single-story buildings (B3, B6 and B8), 4 two-story building (B1, B2, B7 and B10), 3 three-story buildings (B4, B5 and B11) and one four-story building (B10).
The buildings, mainly built with irregular stone and brick masonry characterized by different quality and workmanship skill, have been selected in order to cover wide ranges in terms of structural typologies, dimensions (volume varies between 3700 m 3 and 45,000 m 3 ) and inter-story heights (height varies between 3 m and 5 m).
The 3D arrangements of the seismic resistant walls of the investigated buildings is synthesized in Figure 8.

Identification of Uncertain Parameters
As discussed in the previous paragraph, the shear behavior of masonry walls is characterized by significant uncertainty regarding the stiffness identification, the displacement capacity and the ultimate strength.
In non-linear static analysis, the hysteretic behavior of masonry walls subjected to constant vertical load and cyclic horizontal loads is idealized by means of a bi-linear resistance envelope.
As anticipated, the bi-linear envelope is characterized by an initial elastic slope, defined by the lateral stiffness k and by a plastic plateau, bounded by the elastic inter-story drift and by the ultimate inter-story drift , corresponding to the shear resistance of the wall (see Figure 9). Assuming both ends are clamped, the lateral stiffness k of a masonry wall can be calculated, taking into account shear deformations and bending effects [30], as:

Identification of Uncertain Parameters
As discussed in the previous paragraph, the shear behavior of masonry walls is characterized by significant uncertainty regarding the stiffness identification, the displacement capacity and the ultimate strength.
In non-linear static analysis, the hysteretic behavior of masonry walls subjected to constant vertical load and cyclic horizontal loads is idealized by means of a bi-linear resistance envelope.
As anticipated, the bi-linear envelope is characterized by an initial elastic slope, defined by the lateral stiffness k and by a plastic plateau, bounded by the elastic inter-story drift δ e and by the ultimate inter-story drift δ u , corresponding to the shear resistance of the wall H Rd (see Figure 9). where h is the inter-story height of the wall, l its length, A the area of its cross-section, and E and G are the modulus of elasticity and the shear modulus of masonry respectively. If diagonal tension shear failure governs, the shear resistance of the wall can be derived from where is the compressive stress induced in the wall by the seismic load combination, is the shear strength of masonry and b is the shear resistance factor. The shear resistance factor depends on the aspect ratio of the wall, ℎ⁄ . When ℎ⁄ ≥ 1.5, it can be assumed b = 1.5 [30].
Concerning the ultimate displacement, unreinforced masonry is characterized through the so called "ductility". This property is far from being a ductility in a conventional sense [31], but rather a synthetic parameter accounting for the relative slip along crack surfaces that parts of the wall elements can sustain without significant shear stress loss. Wall ductility depends on normal stress [31,32], on boundary conditions [13], as well as, via the aspect ratio ℎ⁄ , of the wall geometry.
In the Eurocode 8-Part 3 [33] as well as in the guidelines for application of the Italian Building Code [29], the ultimate drift is defined as a percentage of the inter-story height of the wall, depending on the failure mode. For shear failure, the ultimate drift is set to Excluding geometric variables, uncertain parameters characterizing the definition of the force-displacement curve are: • the shear modulus G and the elastic modulus E of masonry; • the shear strength ; • the ultimate displacement .
For the purposes of the present study, suitable pdf has been derived for each relevant parameter, considering that the outcomes of seismic assessment are strongly dependent in non-trivial way, beside the adopted failure criterion, also on the probabilistic distributions of the material properties [34].
Following an approach similar to that proposed in [6], a normal distribution has been associated to the shear strength . The mean value of , ̅ , has been estimated combining in situ test results with the indications reported in [29], taking into account the masonry class and its qualitative features, while its standard deviation has been derived in such a way that 90% of the experimental data lies within the range of values given in [29] for the identified masonry class. In case of stone masonry, it has been obtained a coefficient of variation, COV = ̅ = 0.14 ⁄ . Regarding the shear modulus, due to belittle consideration influence of anisotropy and cracking of masonry, it seems that values provided in almost all national codes [29] as well as in Eurocode 6 [35] overestimate the actual values as reported in [24,36]. Based on a critical analysis of a Assuming both ends are clamped, the lateral stiffness k of a masonry wall can be calculated, taking into account shear deformations and bending effects [30], as: where h is the inter-story height of the wall, l its length, A the area of its cross-section, and E and G are the modulus of elasticity and the shear modulus of masonry respectively. If diagonal tension shear failure governs, the shear resistance of the wall H Rd can be derived from where σ 0 is the compressive stress induced in the wall by the seismic load combination, τ k is the shear strength of masonry and b is the shear resistance factor. The shear resistance factor depends on the aspect ratio of the wall, h/l. When h/l ≥ 1.5, it can be assumed b = 1.5 [30]. Concerning the ultimate displacement, unreinforced masonry is characterized through the so called "ductility". This property is far from being a ductility in a conventional sense [31], but rather a synthetic parameter accounting for the relative slip along crack surfaces that parts of the wall elements can sustain without significant shear stress loss. Wall ductility depends on normal stress σ 0 [31,32], on boundary conditions [13], as well as, via the aspect ratio h/l, of the wall geometry.
In the Eurocode 8-Part 3 [33] as well as in the guidelines for application of the Italian Building Code [29], the ultimate drift is defined as a percentage of the inter-story height of the wall, depending on the failure mode. For shear failure, the ultimate drift is set to Excluding geometric variables, uncertain parameters characterizing the definition of the force-displacement curve are: • the shear modulus G and the elastic modulus E of masonry; • the shear strength τ k ; • the ultimate displacement δ u .
For the purposes of the present study, suitable pdf has been derived for each relevant parameter, considering that the outcomes of seismic assessment are strongly dependent in non-trivial way, beside the adopted failure criterion, also on the probabilistic distributions of the material properties [34].
Following an approach similar to that proposed in [6], a normal distribution has been associated to the shear strength τ k . The mean value of τ k , τ k , has been estimated combining in situ test results with the indications reported in [29], taking into account the masonry class and its qualitative features, while its standard deviation σ τk has been derived in such a way that 90% of the experimental data lies within the range of values given in [29] for the identified masonry class. In case of stone masonry, it has been obtained a coefficient of variation, COV = σ τk /τ k = 0.14.
Regarding the shear modulus, due to belittle consideration influence of anisotropy and cracking of masonry, it seems that values provided in almost all national codes [29] as well as in Eurocode 6 [35] overestimate the actual values as reported in [24,36]. Based on a critical analysis of a large database of experimental test results [25], empirical relationships can be found for shear modulus and elastic modulus: where α 1 = G/τ k , is a normally distributed uncertain parameter with mean value 1500 and coefficient of variation 0.3, while α 2 = G/E, is again normally distributed with mean 0.15 and coefficient of variation 0.2. For the ultimate displacement δ u , a database of test results is collected and analyzed in [26] and [27], mirroring the high variability of this parameter. A median value δ u = 0.47% h is reported in [26] in case of shear failure, which is in accordance with the value of δ u = 0.4% h given in Eurocode 8-Part 3 [33] and the Italian Building Code. Since in [26] a minimum value of 0.14% h is given for δ u , δ u can be assumed normally distributed with mean value equal to 0.4% h and coefficient of variation 0.2, so that the minimum value 0.14% h corresponds to the 0.1 fractile.
The statistical models adopted to model uncertain parameters are summarized in Table 1, in terms of mean values (X) and the coefficient of variation (COV).

Sensitivity Analysis
Sensitivity analysis has been performed on the aforementioned set of case studies in order to quantify the relevance of each parameter in the evaluation of the seismic risk index.
Response surface methods based on general polynomial chaos expansions [37] have been used to represent uncertain parameters and their propagation through the models.

Response Surface via generalized Polynomial Chaos Expansion
Let Q : Ω → R nq the vector of input random parameters [ τ k ; G/τ k ; G/E; δ u ], described in the previous paragraph, characterizing the mechanical model, said Ω the set of possible events. The epistemic uncertainty coming from our lack of knowledge is modelled by this randomness, at which prior probability distributions π(q) have been assigned in Table 1, based on professional expertise. For mathematical convenience we also introduce maps F : Ω → R n . These maps associate the set of random variables to a set of mutually independent random variables Z i , Z i : Ω → R n , characterized by standard normal distributions π(z i ).
Then, the non-linear static analysis algorithm E-PUSH can be viewed as the forward model G, a deterministic solver that takes as input a given set of parameters q and provides a unique response vector u u = G(q), where u is a vector that gathers the response quantities, in this case the seismic risk index I R.
Performing uncertainty quantification analysis consists in running the forward model G many times to quantify the uncertainty of the output and may require an high computational cost. In order to speed up the computations, a proxy model for the predicted measurable u replacing the G map can be of great use. Taking advantage of functional approximations of the random variables by means of the generalized polynomial chaos (gPC) expansion, the construction of response surfaces is significantly facilitated and a computationally cheap proxy model is obtained [38]. For an extensive review of this topic, please refer to [39].
Without going here into much mathematical detail, assuming that U = G(Q) has a finite variance, the U N approximation of the response reads: In Equation (19) Φ i (Z) = φ i 1 . . . φ i n are the multivariate gPC basis functions, some orthogonal polynomials of total degree less than or equal to N, which is the degree at which the expansion is truncated; Z is a Gaussian vector whose corresponding orthogonal polynomials are the Hermite polynomials, which are orthogonal with respect to the Gaussian measure;û i are the coefficients of the PCE, i = (i 1 , . . . , i n ) ∈ N n 0 is a multiindex (with |i| = i 1 + . . . + i n ). The computation of the coefficients can be carried out with different methods such as the interpolation/regression, and the pseudo-spectral projection. Further details about this topic are given in [39,40]. When the gPC expansion of a given forward model is available, one has in fact an analytical representation of u in terms of z, with the advantage that for any realization of the random vector Q the response u can be easily evaluated by first mapping q to z and then evaluating the gPC expansion without much computational expense.
The response surface is then obtained by running few times the forward model in correspondence of the integration points determined via Gaussian quadrature rule for the random parameters. In this specific case, a degree six polynomial expansion has been used (2401 integration points) and the accuracy of the response given by the surrogate model has been assessed checking that it fits the exact solution.
As an example, the response surfaces obtained for seismic risk index of the first case study (Building 1) are reported in the four diagrams in Figure 10: • in terms of G/τ k and δ u , assuming mean values for τ k and G/E; • in terms of G/τ k and τ k , assuming mean values for δ u and G/E; • in terms of τ k and δ u , assuming mean values for G/τ k and G/E; • in terms of G/τ k and G/E assuming mean values for τ k and δ u .

Evaluation of Sobol Indices
As demonstrated in [41], it must be remarked that, besides providing a complete representation of the random response of a model, the gPC expansion also allows an analytical derivation of Sobol indices [42,43], so avoiding Monte Carlo simulations, which are much more computationally demanding.
Sobol indices, which measure the quota of the total variance referable to the uncertainties of each input parameters, are recognized to be accurate descriptors of the sensitivity of the model since they do not assume any kind of linearity or monotonicity of the model [43].
Once the gPC representation of the model is available, the Sobol indices S q can be analytically computed as where D PC is total variance of the response.
points) and the accuracy of the response given by the surrogate model has been assessed checking that it fits the exact solution.
As an example, the response surfaces obtained for seismic risk index of the first case study (Building 1) are reported in the four diagrams in Figure 10: • in terms of / and , assuming mean values for and / ; • in terms of / and , assuming mean values for and / ; • in terms of and , assuming mean values for / and / ; • in terms of / and / assuming mean values for and .  Referring to the 11 case studies, the results of the sensitivity analysis, carried out considering both relevant direction of application of the seismic action, are summarized in Figure 11. These results demonstrate that the most relevant model parameters influencing the total variance are ultimate displacement δ u (≈54% of the total variance) and the G/τ k ratio (≈36% of the total variance) and, to a lesser extent, the shear strength τ k (≈9% of the total variance), while the influence of G/E ratio is negligible, as it represents only the 1% of the total variance (see Figure 12). As demonstrated in [41], it must be remarked that, besides providing a complete representation of the random response of a model, the gPC expansion also allows an analytical derivation of Sobol indices [42,43], so avoiding Monte Carlo simulations, which are much more computationally demanding.

Evaluation of Sobol Indices
Sobol indices, which measure the quota of the total variance referable to the uncertainties of each input parameters, are recognized to be accurate descriptors of the sensitivity of the model since they do not assume any kind of linearity or monotonicity of the model [43].
Once the gPC representation of the model is available, the Sobol indices can be analytically computed as where is total variance of the response. Referring to the 11 case studies, the results of the sensitivity analysis, carried out considering both relevant direction of application of the seismic action, are summarized in Figure 11. These results demonstrate that the most relevant model parameters influencing the total variance are ultimate displacement (≈54% of the total variance) and the / ratio (≈36% of the total variance) and, to a lesser extent, the shear strength (≈9% of the total variance), while the influence of / ratio is negligible, as it represents only the 1% of the total variance (see Figure 12).
It must be underlined that in the considered case the sum of contributions to the total variance associated with and / is not particularly sensitive on the scattering of . For example, when the coefficient of variation of the shear strength increases from 0.14 to 0.2, the contribution of raises to 19% of the total variance, while the contribution of / decreases to 30%., but the total variance associated to both contributions remain almost the same (≈49% vs. ≈45%).     It must be underlined that in the considered case the sum of contributions to the total variance associated with τ k and G/τ k is not particularly sensitive on the scattering of τ k . For example, when the coefficient of variation of the shear strength increases from 0.14 to 0.2, the contribution of τ k raises to 19% of the total variance, while the contribution of G/τ k decreases to 30%., but the total variance associated to both contributions remain almost the same (≈49% vs. ≈45%).

Uncertainty Quantification
From the response surfaces previously defined, statistics of the seismic risk index can be finally easily computed quantifying the uncertainty in the results and providing an important information for the performance classification of the investigated structures.
In Table 2 and in the box-plot in Figure 13, the results in terms of median values, 5-th and 95-th percentiles are reported for the 11 case studies considering the most relevant direction of application of the seismic action.

Seismic Performance Classification
As already underlined, a robust and reliable performance classification of buildings is of paramount importance in decision-making strategies about seismic retrofit.
In order to introduce a suitable performance indicator, a performance classification, based on the probabilistic risk assessment performed in the previous section, is proposed here.

Seismic Performance Classification
As already underlined, a robust and reliable performance classification of buildings is of paramount importance in decision-making strategies about seismic retrofit.
In order to introduce a suitable performance indicator, a performance classification, based on the probabilistic risk assessment performed in the previous section, is proposed here.
Seismic risk categories are identified considering the median value (I R,50 ) and the 5th percentile (I R,05 ) of the seismic risk index. In particular, five categories (A-E) with increasing magnitude of vulnerability are defined according the interval reported in Table 3. Table 3. Seismic performance classification-definition of categories. Allowing us to quantify and to communicate risks, the seismic performance classification defined in Table 3 could represent an important tool for stakeholders to undertake decision about seismic retrofit strategies, helping them to prioritize the interventions for the most vulnerable structures.
Focusing the attention on the buildings investigated, the classification proposed in Table 3 leads to the classification: single-story buildings B3, B6 and B8 are in class A, buildings B9 and B11 are in class C, and buildings B1, B2, B4, B5, B7 and B10 result in being in class D (see Figure 14).

Conclusions
In this paper, a methodology for the probabilistic assessment of the seismic risk index of existing masonry buildings is proposed.
Starting from a simplified non-linear pushover-type algorithm developed by the authors, the procedure allows us to quantify the propagation of uncertainties related to masonry parameters in the resulting seismic risk index.
The procedure is based on the definition of response surfaces for the seismic risk indices IR through gPC-expansion; in this way, quantification of uncertainties in the resulting seismic risk

Conclusions
In this paper, a methodology for the probabilistic assessment of the seismic risk index of existing masonry buildings is proposed.
Starting from a simplified non-linear pushover-type algorithm developed by the authors, the procedure allows us to quantify the propagation of uncertainties related to masonry parameters in the resulting seismic risk index.
The procedure is based on the definition of response surfaces for the seismic risk indices I R through gPC-expansion; in this way, quantification of uncertainties in the resulting seismic risk index is easily obtained, reducing the computation effort, and Sobol sensitivity indices can be computed analytically from the gPCE coefficients.
The main sources of uncertainties, related to masonry material parameters and in-plane displacement capacity, have been investigated and, by means of sensitivity analysis, their impact in the evaluation of seismic risk indices has been quantified. The results confirm that the most relevant parameters are related to definition of the ultimate displacement δ u and the shear modulus G, highlighting the necessity for further investigations on the shear behavior of masonry walls.
Finally, a seismic performance classification has been presented, defining five seismic risk categories on the basis of the median value (I R,50 ) and the 5th percentile (I R,05 ) of the seismic risk index.
The outcomes have been illustrated for a set of existing masonry buildings that are part of the school system in the Municipality of Florence, demonstrating how the presented procedure is very promising and provides a robust and reliable tool to support decisions about seismic retrofit and the priority of interventions for existing masonry buildings.
Funding: This research received no external funding.