Seismic vulnerability assessment and fragility analysis of pre-code masonry buildings in Portugal

Despite the fact that in recent years Portugal has not seen the occurrence of high-magnitude earthquakes, it remains threatened by these events due to its geographic location. Since the 1960s, reinforced concrete has been the most used material for new constructions; however, the historic urban centers are dominated by old unreinforced masonry (URM) buildings, which techniques and construction materials have evolved since the Great Lisbon earthquake that occurred in 1755 (Mw = 8.5). Given the presence of these buildings in areas of significant seismicity, extensive research is needed to assess the seismic risk and define mitigation policies. This kind of studies is often supported by empirical methods and based on expert judgment due to the high variability of the building stock and lack of information. The main purpose of this work is: (i) to provide analytical fragility curves, supported by nonlinear static analysis, for the entire population of old masonry buildings, built before the introduction of the first design code for building safety against earthquakes (RSSCS) in 1958; (ii) define vulnerability curves to be used by the technical community for seismic assessment of pre code URM buildings. The characterization of the building stock geometry and material properties is based on information previously collected, which was essential to define representative archetypes and typologies.


Introduction
Over the years, masonry structures have shown evidence of good behavior under vertical static loads. However, its characteristics, such as the high specific mass, low tensile and shear strength, make the use of this heterogeneous material unsuitable in earthquake prone areas, e.g., Andradiva-Greece (2008), L'Aquila-Italy (2009), Emilia-Romagna-Italy (2012), Umbria-Italy (2016), Abruzzo-Italy (2017). Although Portugal has not been the target of high magnitude earthquakes in recent years, it remains susceptible, due to this geographical location, as it occurred in the past (Oliveira 1986): the 1755 Lisbon earthquake (M w = 8.5), 1909 Benavente earthquake (M w = 6.3), the 1969 Algarve earthquake (M w = 7.8), Azores 1980 (M w = 7.2) and 1998 (M w = 5.8). These events caused significant damage in the affected regions, and particularly on masonry constructions (Correia et al. 2015).
The Portuguese building stock in historic urban centers is predominantly constituted by old unreinforced masonry (URM) residential buildings (INE 2012). Their characteristics are the result of different periods of construction and construction practice due to the available materials, existing techniques, and society needs. In general, four typologies, described in Sect. 3, can be identified in the urban centers: "Pre-Pombalino" -before 1755 Lisbon earthquake (Domingos 2010); "Pombalino" -1755 to 1870 (Meireles and Bento 2012;Lopes et al. 2014); "Gaioleiro" -1870 to 1930 (Candeias 2008;Mendes 2012;Simões 2018) and "Placa" -1930to 1960(Lamego 2014Ferrito et al. 2016;Milosevic 2019;Bernardo et al. 2022). It is also worth pointing out that no impact of earthquake has been considered in their design as the first Portuguese seismic design regulation appeared only in 1958 (RSCCS 1958).
In the last decades, the performance of buildings under seismic action has received special attention due to the interest in the conservation of heritage and protection of human life. The seismic risk at a national scale was evaluated by (Sousa 2006) and , and by other authors at urban scale, e.g.: Coimbra (Vicente et al. 2011), Faro (Vicente et al. 2014), Seixal (Ferreira et al. 2013). Most of these studies employed statistical data and expert opinion combined with empirical methods to derive fragility and vulnerability functions to characterize the building stock. related to vulnerability and seismic risk assessment, the first part of this work aims to provide analytical fragility curves for the population of old (pre-code) URM buildings in Portugal with rigid and flexible floor diaphragms. The second part of the work consists in deriving mechanistic vulnerability functions for seismic assessment of these typologies, in compliance with the specifications provided by the NP EN 1998-3:2017 -Portuguese version of Eurocode 8-part 3 (EN1998-3:2017), hereinafter EC8-3, -thus can be employed by the technical community to evaluate the seismic performance of a pre-code masonry building in Portugal. Therefore, the behavior of the buildings analyzed is only limited to the inplane mechanisms according to the current version of EC8-3. Although the analyses carried out in the present work have been derived to account only for the in-plane mechanisms to be compliant with the current version of EC8-3, which assumes these as being prevented from occurring, several works regarding out-of-plane fragility functions for masonry buildings can be found in literature, emphasizing the importance of such mechanisms, namely in buildings with flexible floor diaphragms (Costa 2012;Ceran and Erberik 2013;Ferreira et al. 2015;Simões et al. 2020;Sumerente et al. 2020;Giordano et al. 2021;Jaimes et al. 2021).
In the framework of the present study, the development of the fragility/vulnerability curves, involves the following steps: (i) generate a synthetic database of 18.000 masonry buildings up to 5 stories high, including different archetypes based on statistical information previously collected and different material properties to cover the variability found in the literature; (ii) estimate the in-plane seismic behavior of the entire database through displacement-based nonlinear static methods; (iii) derivation of fragility functions for the capacity expressed by the maximum interstorey drift. The fragility curves proposed are only related to the deformation capacity of the buildings in order to be applied in seismic risk studies or safety assessment; (iv) computation of the seismic demand and derivation of vulnerability curves expressed by the interstorey drift (EDP) for different levels of seismicity in compliance with EC8-3.

Seismic fragility curves: background
Traditionally, fragility curves are probabilistic relationships introduced to reflect the stochastic nature of seismic input, which is increased by the structural nonlinear behavior response, expressing the probability of exceeding a given limit state related to the behavior for a given seismic intensity measure of the input ground motions.
Fragility curves can be classified in four categories (Porter et al. 2007;Frankie et al. 2013;Asteris et al. 2014;Lagomarsino and Cattari 2014;Pitilakis et al. 2014;Kappos and Papanikolaou 2016): (i) empirical, which are based on post-earthquake observations and very useful for calibration/validation of fragility curves computed from other methods (Colombi et al. 2008;Azizi-Bondarabadi et al. 2016;Del Gaudio et al. 2017;Rosti et al. 2021); (ii) expert elicitation/judgement, where expert opinion is used to predict a relationship between the damage and the seismic intensity level when data is limited or not available; (iii) analytical, obtained from numerical approaches using simplified mechanical models or based on non-linear analyses (D'Ayala 2005; Rota et al. 2010;Lagomarsino and Cattari 2014;Giordano et al. 2021); (iv) hybrid, combining the different aforementioned methods (Kappos et al. 2006;Maio et al. 2015;Sandoli et al. 2021).
Currently, analytical fragility curves have been more used given the increased computational capacities nowadays, allowing to derive parametric studies (e.g., simulating strengthening and retrofitting solutions) and controlling the associated uncertainties. In general, the probabilistic description of analytical fragility functions is obtained by computing several Incremental Dynamic Analyses (Vamvatsikos and Cornell 2002) for representative sets of ground motions normalized to the same seismic intensity measure. For risk analysis purposes, the dispersion of the fragility distributions is increased to cope with other uncertainties affecting the structural performance. The randomness in the seismic input is tackled by the dispersion of attenuation laws in seismic hazard analysis.
A literature review reveals several works regarding the development of analytical fragility curves, in particular, for masonry buildings, which is not a straightforward task given the wide variability found in this kind of constructions, related to the different types of masonry, techniques used, structural systems adopted and state of maintenance. For example, Lagomarsino and Cattari 2014 proposed the construction of analytical fragility curves using simplified mechanical models (Displacement-based Vulnerability method) based on geometrical and mechanical parameters, loading scenario and correction factors for the evaluation of the acceleration at yielding and fundamental period. Alternatively, analytical fragility curves can be directly derived using nonlinear static or dynamic analysis on prototype buildings (Rota et al. 2010), which are more suitable to account for specific constructive details, such as addressed by Milosevic et al. 2020 andSimões et al. 2020 for a set of prototype Portuguese masonry buildings. Recently, Lagomarsino et al. 2021 also proposed a heuristic vulnerability model to derive fragility curves for masonry buildings, based on the expertise that is implicit in the European Macroseismic Scale, with fuzzy assumptions on the binomial damage distribution, and calibrated by damage observation in Italy.
The "Pre Pombalino" buildings, constructed before the 1755 earthquake, are recognized by their irregular geometry, reduced dimensions, narrow facades, high density of walls and few openings to the exterior. They usually are four stories high and are constituted by poor-quality masonry walls supporting the timber floors. The "Pombalino" buildings emerged after the Lisbon earthquake and are particularly known by the improvements in the anti-seismic conception in that period. They usually have up to five stories high and regular geometry. This typology was standard in the building construction practice for more than one century. On the other hand, the "Gaioleiro" buildings represent a downgrade when compared to the previous typology, with the adoption of more simplified construction techniques and the use of low-quality materials which was promoted by the rapid expansion of the urban centers and the housing demand. This typology is significantly more vulnerable, from a seismic point of view, compared with the previous one. Finally, "Placa" buildings emerged before the enforcement of the first seismic-code in 1958 and introduced the use of lightly reinforced concrete slabs at the floors level. The high mass of the RC slabs and the low strength capacity of the load bearing walls to horizontal forces results in an unsatisfactory structural seismic performance.

Geometry characterization and definition of archetypes
The geometry characterization comprises the information gathered through detailed drawings from the original projects and collected from municipal archives, for a population of 100 old (pre-code) masonry buildings. This data represents the geometry for the most typical masonry buildings built before the decade of 1960 and described in the previous section. The geometric parameters obtained, such as plan dimensions, height of the stories, openings ratio, interior walls density, walls thickness and type/thickness of floors, were statistically characterized and described in Bernardo et al. 2021.
Based on this information, 9 archetypes -A1, A2, A3, B1, B2, B3, C1, C2 and C3 -with different configurations and up to 5 stories high, were generated for the subsequent analyses. Figure 2 presents the layout of these archetypes. Table 1 summarizes the statistical information for the geometry and the parameters adopted (underlined) to represent the archetypes.
The archetype B2 ( Fig. 2) represents the mean size configuration (12.6 × 12.1 m). The plan dimensions for the remaining were derived from the mean archetype, considering a dispersion equals to one standard deviation: L x = 12.6 ± 5.0 m and L y = 12.1 ± 4.1 m. The total area ranges approximately between 60.0m 2 to 285.0m 2 . The layout for the arrangement of the partitions/interior walls follows, in a reasonable manner, the typical size of the compartments for theses typologies (3 × 3 up to 4 × 5 m), representing a mean value for the interior walls' density equal to 0.054.
Regarding the walls thickness, considering the enormous variability in the type of material, arrangement and absence of information in the documentation gathered, were considered the mean thickness for the facades and side walls. For the interior/partition walls the most common value (mode) were adopted, which are representative for more than 60% of the buildings collected (Bernardo et al. 2021). The intrinsic variability in the wall's characteristics (type of masonry, morphology and arrangement) was computed in the material mechanical properties uncertainty, carried out in Sect. 4.2. For the remaining variables, the mean values adopted are mentioned in Table 1.

Material properties selection
The definition of material mechanical properties was based on literature review. Table 2 presents a background on the types of masonry that can be identified in the Portuguese building stock and the suggested properties by the latest version of EC8-3 (see Candeias et al. 2020).
Taking into account the wide range of masonry mechanical properties, the uncertainty was propagated through Monte Carlo simulations (Fryer and Rubinstein 1983). For that purpose, two groups of buildings were considered with different material properties: Type I -typologies with good quality masonry (e.g., regular and squared masonry, brick masonry with cement lime mortar) and Type II -typologies with poor quality masonry (e.g., ruble stone masonry, brick masonry with lime mortar). Given the differences of the interior/ partition walls (e.g., tabique, 1 frontal walls, 2 perforated brick masonry) when compared with the exterior walls (e.g., solid masonry bricks, stone masonry or concrete blocks), two sub categories were defined: Type I-1 and Type II-1, to represent the properties of exterior walls, and Type I-2 and Type II-2 for the interior/partition walls. A set of 100 samples for each typology -Type I and Type II -were generated to describe the material variability, attaining an error of around 5% (95% confidence interval) for the material population   Table 3 resumes the distributions adopted for the random variables, which are plotted in Fig. 3.

Modelling strategy and general assumptions
Considering the previous geometric statistical information, tridimensional multi degree of freedom models (MDOF) were developed to simulate the nonlinear response of the buildings. For this purpose, an equivalent frame modeling strategy available in the research version (2.1.104) of TreMuri software (Lagomarsino et al. 2013) was used. Some of the features of the model include the accurate representation of the principal in-plane failure mechanism, including the stiffness and strength degradation, such as bending rocking, diagonal shear and sliding. The software was originally developed for frame-type analysis of the entire URM buildings whereby the response is governed by the in-plane behavior of the walls. On the other hand, the current version of EC8-3 also does not include the out-of-plane mechanisms or assumes these as being prevented from occurring. Hence, the behavior of the buildings analyzed is only restricted to in-plane mechanisms.

Macroelement model validation
The macroelement model is defined by a set of mechanical parameters at macroscopic scale that should be representative of an average of the masonry panel properties: Young's modulus-E, shear modulus-G, density-ρ, compressive strength-f c , cohesion-τ 0, friction coefficient-μ, and by two phenomenological parameters related to the shape of nonlinear shear constitutive modelc t that expresses the shear deformability in the inelastic range, wherein the amplitude in the inelastic displacement is proportional to the product G Ct ; and β s that controls the slope of the softening branch in the post-peak region (Lagomarsino et al. 2013).
The validation of the macroelement model was performed using the results of an inplane quasistatic experimental test carried out on a full-scale masonry panel (see Fig. 4a), with aspect ratio of 1.325:1 (2.65 × 2.00 × 0.25 m-height, length and thickness), made of solid clay bricks, extracted from an old masonry building. The panel was tested under double-bending boundary conditions with an axial vertical constant force of 200 kN during the test. The horizontal displacement, with increasing amplitude cycles, was applied at the base and restrained at the top.
The macroscopic mechanical parameters for the macroelement compatible with the wall tested were taken from laboratory destructive tests on small specimens (Marques 2020), including both clay bricks and mortar, in order to obtain the individual material properties, and on which were applied a strain-based first-order homogenization process to provide the effective properties for the masonry, as detailed in (Bernardo et al. 2020). The nonlinear shear deformation parameters (G Ct and βs) were directly obtained through the calibration against the experimental test. Table 4 lists the macroelement mechanical properties and the shear parameters calibrated.
The numerical hysteresis curve was computed by means of sequential monotonic pushover analyses, imposing the same loading history and boundary conditions equivalent to those mentioned above in the description of the experimental test. The criteria used for the validation of the nonlinear macroelement model consisted in approximating the initial stiffness and maximum strength values to those obtained in the experimental test. The comparison of the numerical hysteresis curve with the experimental response is shown in Fig. 4, where a shear failure mechanism is clearly observed Fig. 4b). Although a reasonable prediction of the total hysteretic energy dissipation is achieved, it is important to point out the differences obtained between the numerical and experimental models, namely the higher shear strength in the initial loops ( Fig. 4c) of the numerical model. These results can be explained by the constitutive model for the shear-dominated response implemented in the nonlinear macrolement, which is only able to reproduce the softening behavior, while in the experimental results an initial post elastic stiffening is observed. Furthermore, the friction observed between the masonry wall and the setup device, namely up to 100 mm of cumulative displacement (Fig. 4d), may also influence the post elastic stiffening in the experimental hysteretic response. However, note that the variability considered in the randomness of materials properties (see Table 3) for the pushover analyses subsequently performed, include the values obtained in the experimental test validation.

Floor diaphragms modelling
Floor diaphragms were modelled as a two-dimensional orthotropic membrane element, defined by four nodes with two displacement degrees of freedom each, and characterized by the equivalent mechanical properties: equivalent thicknesst eq , modulus of elasticity of the  diaphragm in the principal direction -E 1 -and perpendicular direction -E 2 , shear modulus -G, that influence the horizontal force distribution between walls, and Poisson ratioʋ. Two types of floor diaphragms, rigid and flexible, were considered. Rigid diaphragms were modelled by RC slabs and assuming a load distribution in the walls proportional to the influence area. In this case, was assumed a good connection between the walls to ensure equal planar displacements at floor level, simulated through rigid links beams. For flexible diaphragms was adopted a typical timber floor, constituted by timber sheathing and timber joists perpendicular to the facades. In this case, the load was distributed by the main timber joists and the connections between walls were modelled through equivalent elastic link beams at the floor level to simulate medium to weak connections, according to (Simöes et al. 2018).
Tables 5 and 6 and summarize the mechanical properties for the membrane elements (rigid and flexible) and for the connection between walls, respectively.

Representative buildings models
Considering the assumptions related to geometry layout (Fig. 2) and discussed in Sect. 4.1, were modelled 45 archetypes of buildings (A1, A2, A3, B1, B2, B3, C1, C2 and C3), up to 5 stories high, as shown in Fig. 5. Attending to materials variability (Table 2) and different type of floor diaphragms (rigid and flexible), a population of 18,000 buildings was generated, based on the modelling assumptions previously described.
With regard to gravity loads applied, are followed the prescriptions of Eurocode 8 (CEN 2004), combining the nominal values of permanent loads G with the quasi-permanent live loading Ѱ E Q. The permanent loads are defined by the self-weight of the masonry (Table 3), while the other gravity loads are presented in Table 7. The live loads depend on the building category, which is assumed for domestic and residential purpose (category A). 6 Numerical analysis and seismic behavior

Methodology for seismic assessment
In order to develop the fragility curves and vulnerability functions compatible with the seismic verification purpose by the current version of EC8-3, the methodology to evaluate the seismic behavior of the buildings follow the recommendations of the standard for the in-plane global safety verification using nonlinear methods. The general methodology of EC8-3 uses a performance-and displacement-based approach to assess the safety level of a given structure. Regarding the required performance levels for Portugal, implicitly related to the seismic hazard, three limit states (LS) are defined to assess the structural performance of an existing building, depending on its importance class: Damage Limitation (DL), Significant Damage (SD) and Near Collapse (NC). The return periods (RP) prescribed by EC8-3 for these LS are defined in accordance with levels of protection, having values of 73, 308 and 975 years for the DL, SD and NC limit states, respectively, corresponding to probabilities of exceedance of 50%, 15% and 5% in 50 years. For residential buildings (importance class II) the safety verification is only mandatory for the SD limit state.

Definition of seismic action for Portugal
The seismic action at ground surface was modelled according to EC8-1 through the elastic response spectrum S e (T) for a given return period, where T is the vibration period of the linear single degree of freedom system (SDOF). Considering the seismicity of the Portuguese regions, two main seismic sources can be distinguished by magnitude, epicenter, event duration and frequency content: (i) interplates scenario, with offshore epicenters, high magnitude, long duration and lower frequency content; (ii) intraplate scenario, mainly occurred inland with moderate magnitudes, short duration and higher frequency content. In order to avoid the inconsistencies found in both scenarios, two main shapes for the elastic spectrum are defined in the National Annex of EC8-1, resulting in the two distinct seismic zonings (type 1 and type 2) presented in Fig. 6. The parameters to describe the respective horizontal elastic spectrum for are summarized in Table 8.

Nonlinear static analysis
The buildings capacity was predicted from nonlinear numerical static analysis (pushover), through monotonic horizontal forces, which requires particular attention on the choice of load pattern. In general, design codes propose to assume two load patterns (e.g., uniform and triangular) to simulate the distribution of inertial forces during the seismic loading on the deformed shape of the buildings. However, the deformed shape depends on the damage on the building and may change during the loading scenario. To overcome this limitation, (Galasco et al. 2006) propose an adaptive pushover algorithm for masonry buildings, wherein the load pattern is proportional to the displacement shape in the previous step. This approach revealed to be more suitable for masonry buildings with rigid floors, comparing with time-history analysis. However, for flexible floors, given the local mechanisms and the mass participation of each single wall in the vibration mode, that approach does not provide significant improvements. In that case, the pseudo-triangular load pattern is better suited to assure that all mass is mobilized (Lagomarsino and Cattari 2015). For the present study, adaptive pushover analyses with inverse triangular first ratio pattern were adopted for buildings with rigid floors and an inverted pseudo-triangular for flexible floors. The control node was selected at the top level and the shear was measured on the base up to  reaching 20% decay of the maximum shear strength (NC limit state), as recommend by the EC8-3 for a global safety verification (Fig. 7). Figures 8 and 9 present the capacity curves normalized for spectral acceleration S a and spectral displacement S d for the archetypes A1, A3, B2, C1 and C3, with 3 to 5 stories high. The red and blue dots correspond to the maximum shear strength and the ultimate displacement for the NC limit state, respectively. The results are presented for the seismic action in the longitudinal direction (parallel to the facade walls), which revealed to be more critical for the geometry layout defined by exterior lateral walls without openings. To support the discussion, Fig. 7 presents de median values obtained in terms of maximum S a and maximum interstorey drift θ for all archetypes and number of floors.
The values of S a observed are lower for typology with weaker properties (type II), as expected, and flexible diaphragms (FD), which do not allow to exploit the entire capacity of the walls as in rigid diaphragms (RD). For RD, the S a decreases with the number of stories and is inversely proportional to the mass; in FD small decreases are also observed but is not so evident. In fact, buildings with RD have much more capacity in load redistribution and tend to be governed by global mechanisms, contrary to what happens on FD. A closer inspection of the results allows to identify uniform values of Sa, independently of the archetype, which reflects a particularity of these buildings related to the similar ratio between the density of walls and area in plan. Moreover, for archetype A1, B1, and C1, with longer facades (L x > L y ) and lower mass, slightly higher values are reached, evidencing the importance of facade walls compared to interior walls.
Concerning interstorey drift θ, the influence of archetype layout is more notorious in archetypes A3, B3, and C3, in particular for FD, with greater number of interior walls. This fact leads to an increase of the θ provided by the reduction of the participation of each wall in the global behavior and by the greater redistribution of the horizontal forces at floor level. It is also worth pointing out that the values of θ is higher as the number of floors increases, similar to a vertical cantilever, which does not occur in FD, which in turn reach approximately constant higher values independently of the building height.

Seismic performance-based assessment
The seismic performance of the buildings was evaluated following the N2 Method iterative procedure recommended in the Appendix B of EC8-1. The response of the structure  is determined from the intersection of the capacity curve with the seismic demand spectrum in acceleration-displacement response spectrum (ADRS) format. The method employs the dynamic properties of the structure to convert the multi degree of freedom (MDOF) capacity curves into a bilinear equivalent SDOF system, assuming an elasticperfectly plastic force-displacement relationship and incorporates the inelastic response spectrum based on structure's ductility.
In this study, the seismic performance was evaluated for each building, adopting the response spectrum for the seismic action offshore (seismic action type 1) and inland (seismic action type 2) defined in previous Sect. 6.2 for Portugal, an equivalent viscous damping equal to 5%, soil amplification factor corresponding to a soil type A, B and C, and a wide range of return periods (RP) up to 5000 years. Taking the 475-years RP as the reference period RP ref , the acceleration on the ground a g can be scaled to other RP as suggested in the comments to the Portuguese version of EC8-1: a g = a gr (RP ref ∕RP) −1∕k , where a gr is the acceleration for the RP ref and k take the values of 1.5 and 2.5 for seismic action in mainland offshore and onshore, respectively, and 3.6 for Azores. These coefficients assume the mean values obtained in the counties and a first-order power law approximation for the seismic hazard.
As an example, Fig. 10 shows the performance of an archetype B2 building with four stories high for all seismic zones, 308-years RP and soil type B. Figures 11 and 12, as an example, show the box-whiskers plots of the performance points achieved, expressed as a function of interstorey drift θ, for buildings with four stories high, all typologies, RP = 308 years, seismic zones 1.3 (ag = 1.13 m/s 2 ) and 2.3 (ag = 1.43 m/s 2 ) and soil type B. Boxes limits were defined by the quartiles with one standard deviation. Horizontal red line in each box represents the median θ and the square dots the mean θ. Vertical lines (whiskers) extending from each box represents the 5 th and 95 th percentiles, and the dots outside the lower and upper ranges the outliers.
Analyzing Fig. 11 one can readily see that the demand in seismic zone 1.3 is roughly around 0.15-0.25% for RD and 0.40-0.55% for FD. For seismic zone 2.3 (Fig. 12) is approximately 0.10-0.15% and 0.25-0.30% for RD and FD, respectively. In a general way, the differences in the interstorey drift are more significant between typologies than between archetypes. Another point is the higher dispersion obtained in buildings with weak properties (Type II), with more emphasis in RD.

Fig. 10
Seismic performance (red point) for a building with four stories, archetype B2 and type I-rigid typology, considering 308-years return period, soil type B and all seismic zones defined for Portugal 7 Vulnerability and fragility analysis

Limit states
In the present work, the limit state (LS) thresholds follow the global scale criterion as a percentage of maximum shear (F max ) defined in the capacity curves, which are in line with the LS proposed by the EC8-3. Therefore, the LS are expressed in terms of spectral displacement (S d ) as a function of the base shear measured, according to criteria indicated in Table 9. The LS1 was defined at the yielding point of the idealized capacity curve and LS2 and LS3 correspond to the peak of maximum shear and post-peak range, respectively.

Fragility functions
The fragility curves proposed in this section were derived from empirical cumulative distribution functions (CDF) for the data analyzed, directly obtained from the nonlinear response of the buildings, considering the limit states indicated in Table 9. Therefore, the fragility curves presented are independent of the seismic action, or spectrum format, and represent the capacity exceedance probability conditioned on a value of demand for the three limit states adopted.
For this purpose, the archetypes were grouped and analyzed by number of stories, typology (Type I and Type II) and type of floor diaphragm rigid (RD) or flexible (FD).

Fig. 11
Seismic response in terms of interstorey drift for buildings with four stories high, 308-years return period (SD limit state), soil type B and seismic zone 1.3 The best cumulative analytical function was fitted to the data based on Kolmogorov-Smirnov tests, and follow, in a reasonable manner, a LogNormal (LN) and Weibull (W) distribution for typologies Type I and Type II, respectively. Figure 13 presents the proposed fragility curves, expressed as a function of interstorey drift C , by number of floors, typology and type of floor diaphragm. Table 10 summarizes the median values of C and dispersion C for the analytical fragility functions proposed for the buildings' capacity. Note that, the values of dispersion achieved in this section include the randomness in material properties defined for each typology and the variability in the geometry layout of the archetypes, allowing to account both variables in the capacity of the buildings for seismic risk studies or seismic assessment (Fig. 14).
Bar chart of Fig. 15 presents a comparison of the moments proposed for the different typologies, type of floor diaphragm and number of stories, to support the discussion of results. For the mean values C (%) computed, the main differences are observed between typologies and type of floor diaphragm: Type I/RD 0.09 to 0.12 (DL), 0.37 to 0.47 (SD), 0.45 to 0.70 (NC); Type II/RD 0.07 to 0.09 (DL), 0.22 to 0.31 (SD), 0.38 to 0.70 (NC); Type I/FD 0.21 to 0.40 (DL), 0.65 to 0.91 (SD), 0.70 to 1.09 (NC); Type II/ FD 0.23 to 0.33 (DL), 0.51 to 0.75 (SD), 0.70 to 0.91 (NC). As can be noticed, high values of C are attained for structures with FD and mostly for type I typology. In contrast, the type II-rigid presents minor drifts, followed by the type I-rigid. In general, buildings with good quality masonry (Type I) present higher values of C , except for one storey height with FD. Fig. 12 Seismic response in terms of interstorey drift for buildings with four stories high, 308-years return period (SD limit state), soil type B and seismic zone 2.3  Thus, in general, the values of dispersion attained are higher for typology Type II. For the same typology, minor differences between the number of floors are noticed, except for buildings with one story. For the SD limit state, the dispersion has a slight range Fig. 13 Analytical fragility proposed curves for the buildings' capacity expressed by C

Table 10
Moments of the analytical fragility curves proposed for the buildings' capacity a, b Are, respectively, the scale and shape parameter for the Weibull distribution No. of stories between the number of floors, comparing to the other limit states, excluding the Type II-flexible typology with relatively similar dispersion. Finally, for the DL seems to be some trend for lower dispersion in low-rise buildings, in contrast to the higher dispersion computed in taller buildings for the NC limit state. Finally, to exemplify the damage level attained for the LS previously discussed, the structural damage mapping for a given building in the database is visualized in Fig. 14 for the main facade. The propagation of the damage depicted in Fig. 14 corresponds to the most common pattern observed in the numerical models for the adopted LS values. The level of shear damage suffered by each structural element is distinguished by colors, based on the definition of the shear damage parameter s implemented in the macroelement model (Penna et al. 2014): the parameter equals zero when the element is undamaged; equals 1 in correspondence to the attainment of the maximum shear strength, and exceeds 1 in the post-peak (softening) phase. The results obtained for the LS proposed in EC8-3 can be qualitatively compared to the EMS-98 damage grades, allowing to establish a correspondence between the DL, SD and NC with, respectively, the D2 (slight), D3 (moderate), and D4 (heavy) proposed by the EMS-98 damage grading scale.

Seismic vulnerability assessment
In this section are proposed mechanistic seismic vulnerability functions as a relationship between the interstorey drift (taken as the EDP) given a rate per year of seismic input, instead of a traditional intensity measure, and in compliance with the seismic action defined by EC8.
In order to derive these vulnerability functions, the synthetic database of buildings was subjected to different levels of seismicity up to 5000-years return period, RP = {1, 10, 20, 50, 73, 95, 225, 308, 475, 975, 1100RP = {1, 10, 20, 50, 73, 95, 225, 308, 475, 975, , 2475, 5000}, allowing to define the demand for a given seismic intensity level. The analyses were performed for buildings foundations on ground types A, B and C and all seismic zones for Portugal. The results obtained were grouped by different number of stories, building typology/material type and type of floor diaphragms, as performed in previous sections. Figure 16 show the results through the histograms and the PDF fitted by a LogNormal distribution for the 308-years RP, different typologies, four stories high, seismic zone 1.3 and ground type B. Based on these results computed for the entire range of RP considered, analytical curves were derived using a first-order power law function best fitted to the data using a nonlinear least square method (Levenberg-Marquardt algorithm) (Moré 1978). The dispersion in the seismic demand D was computed over the entire range of RP by the standard deviation of the logarithmic error between the analytical function fitted and the empirical data.
Taking as an example the four stories high buildings, Fig. 17 shows the computation of the vulnerability curves for this particular case in seismic zone 1.3 (region of Lisbon and offshore scenario) and ground type B. The vulnerability functions are expressed as a function of the interstorey drift demand D for the 16, 50 and 84% quantile up to 5000-years return period. The variability in the seismic demand is also shown for the 73-, 308-and 975-years RP, corresponding to the limit states defined in Sect. 7.1, respectively, DL, SD, and NC. The grey dot plot histogram corresponds to the seismic performance of the buildings' database when subjected to a certain level of seismicity corresponding to a discrete value of RP.
Tables 11 and 12 present the regression parameters a and b and the dispersion D for the vulnerability curves proposed for buildings up to five stories high, seismic zones 1.3 and 2.3 (see Fig. 6 b), soil types A, B and C. The seismic demand in terms of interstorey drift D (%) is expressed by the rate per year ≅ 1∕RP (input parameter), instead of return period (RP), and can be computed as D = a • b . Figures 18 and 19 show the graphical computation of the analytical vulnerability curves corresponding to the 50%, 84% and 95% quantiles for the case of four-stories high buildings, located in the region of Lisbon (seismic zone 1.3 and 2.3) and ground type B.  Tables 11 and 12, allowing to relate the seismic intensity level (rate per year ) to the values of the Interstorey drift demand D . Furthermore, seismic demand was also computed in terms of spectral acceleration demand S a (Figs. 18 and 19) and related with D . The correct interpretation of Figs. 18 and 19 should follow the following sequential reading order: (i) definition of the rate per year to be evaluated (input variable on the y-blue axis); (ii) definition of the quantile value for the seismic demand; (iii) reading the projection of in terms of interstorey drift demand D (x-axis) or the corresponding value of spectral acceleration demand S a (y-red axis), according to the quantile defined in (ii). Note that, in the case of S a , the values become constant due to the assumption of a bilinear capacity curve for the determination of seismic performance, which implicitly leads to the same value of spectral acceleration (performance point) in the plastic range, independently of the seismic intensity level.
Finally, to demonstrate the usefulness of these vulnerability curves, it is assumed, for instance, a residential building (ordinary building) with four stories high, poor-quality masonry and flexible diaphragm. Adopting the seismic assessment for the SD limit state, corresponding to a rate per year of approximately 3.25 × 10 -3 (308-years RP -input variable), the seismic demand in terms of interstorey drift D is equal to 0.65% and 0.32% for the 95% quantile, seismic zones 1.3 and 2.3, respectively (see Figs. 18 and 19 for Type II-flexible buildings). Regarding the building capacity, expressed as a function of C , the   . 19 Seismic demand for buildings with four stories high, soil B, and seismic zones 2.3 median value of 0.51% was assumed, according to Table 10 (type II-flexible and SD limit  state). Thus, given these assumptions, this particular case verifies the seismic safety for the zone 2.3. Alternatively, conventional pushover analysis could be performed to estimate the building capacity (instead of using Table 10), or even another quantile could be assumed, depending on the knowledge level. Naturally, in the case of using the S a demand curves proposed in Fig. 18 and Fig. 19, the building's capacity should be also evaluated in terms of spectral acceleration capacity.
Recently, relevant studies in line with the present work have been carried out by Milosevic et al. 2020 andSimões et al. 2020 to characterize typical prototypes of pre-code masonry buildings in the region of Lisbon, providing fragility curves for those typologies. The values obtained in these studies for the dispersion in the capacity C ranges, approximately, between 0.04 to 0.20 (Milosevic et al. 2020) and 0.05 to 0.30 (Simões et al. 2020). Regarding the dispersion in the demand D , the values obtained in the previous studies vary, approximately, from 0.18 to 0.26 (Milosevic et al. 2020) and from 0.27 to 0.51 (Simões et al. 2020). It is also important to point out that these values were computed for a prototype with three stories high (Milosevic et al. 2020) and four different types of geometry and five stories high (Simões et al. 2020), considering in both studies a set of scaled response spectra compatible with the seismic action defined in the EC8 for a 475year return period.
Comparing the results achieved in the aforementioned studies, the values of C in the present work vary from 0.21 to 0.58 (three stories high building) and from 0.20 to 0.50 (five stories high building). The discrepancy found between the results obtained in this study and the ones compiled in the previous ones may result from the differences in the aleatory uncertainties considered for the mechanical properties -the initial coefficient of variation assumed to describe the uncertainty in the material properties is around half of the one adopted in the present study. Furthermore, the results presented in this study also account for the variability in the geometry of the buildings, which is not addressed in the previous studies, in particular the one carried out by (Milosevic et al. 2020). Concerning the results for the seismic demand, the present study proposes values for the D ranging from 0.12 to 0.43 (three stories high building) and from 0.11 to 0.36 (five stories high building). However, these values cannot be compared with the previous ones. Note that, in this study, the seismic demand corresponds to the variability in the seismic hazard computed for a wide range of return periods, while in the previous studies the D was only computed for the 475-year return period considering different seismic inputs. Comparing the results obtained in the present work with others reported in the literature, a value of C = 0.30 is proposed by HAZUS for pre-code buildings and a total dispersion varying from 0.37 to 0.80 is achieved by Douglas et al. 2015

Final comments and conclusions
Seismic risk studies are of extreme importance in regions with moderate to high seismicity, such as Portugal, to estimate losses and establish policies for risk mitigation. This kind of studies requires knowledge about the building stock, which is often characterized by empirical methods and expert opinion when performed at large scale. The main purpose of the present paper was to derive analytical fragility and vulnerability curves in compliance with the seismic action defined in EC8, which can be used to conduct more detailed seismic risk studies or employed for seismic assessment of pre-code masonry buildings in Portugal. The proposed fragility curves are not linked to a ground motion intensity, or spectrum format, and can be applied in a more general context to characterize the capacity of the building stock considering the randomness in the material properties and the variability in the geometry. The vulnerability functions derived in this work are presented as relationships between the interstorey drift (EDP) given a rate per year of seismic input instead of a traditional intensity measure.
The development of the proposed fragility and vulnerability curves considered a synthetic database of 18.000 masonry buildings, based on statistical information previously collected about the geometry, allowing to define nine archetypes (A1, A2, A3, B1, B2, B3, C1, C2, C3), which were further combined with a wide range of material properties (Type I -good quality and Type II -poor quality) and type of floor diaphragm (rigid and flexible).
The performance-based assessment supported by nonlinear static analyses and carried out on the synthetic database: (i) the main differences in terms of capacity (maximum spectral acceleration -S a ) are noticed for different typologies rather than specifically between archetypes due to the similar ratio between the density of walls and area in plan, as observed in the data collection (Bernardo et al. 2021); (ii) the facade length along the direction of the seismic action, compared to the increase in area of interior walls, leads to slightly higher values of S a in some typologies (A1, B1 and C1); (iii) the efficiency of load redistribution in buildings with rigid diaphragms (RD) allows to explore more capacity from the walls and higher values of S a can be reached in comparison to flexible diaphragm (FD); (iv) the interstorey drift attained is higher in structures with FD; (v) slightly higher values of interstorey drift were obtained in archetypes with higher number of interior walls (A3, B3 and C3) parallel to the seismic action direction.
Based on the previous analyses, analytical fragility curves for the buildings' capacity were derived for both typologies and type of floor diaphragm. Among the results gathered, the following stand out: (i) buildings with good quality materials and flexile diaphragm (Type I FD) reach higher drift values (up to 0.41%DL, 0.91% SD and 1.09% NC). In contrast, smaller drifts (up to 0.09%DL, 0.31% SD and 0.70% NC) are attained for structures with poor quality masonry and rigid diaphragm (Type II RD); (ii) the dispersion achieved is higher in Type II FD buildings (up to 0.62 DL, 0.60 SD and 0.50 NC) and smaller in Type I RD buildings (up to 0.33 DL, 0.26 SD and 0.32 NC); (iii) in general, drift values and respective dispersion seems to be higher with the increase in number of floors.
Regarding vulnerability functions proposed for the seismic assessment, they are compatible with the methodology prescribed in EC8 and account the variability in the nonlinear seismic response of the buildings generated. This paper presents an example of application for the region of Lisbon. These curves can be used by the technical community to evaluate the seismic performance of pre-code masonry buildings in Portugal, up to five stories high and foundations on a soil type A, B or C.
Although the fragility and vulnerability curves are not presented in conventional format, the former can be used in a more general context to characterize the masonry precode building stock in Portugal, including the geometry and material variability, or can be converted to the conventional format using the vulnerability curves proposed. The latter can also be useful for scenario loss estimation, by converting the EDP (interstorey drift) in losses, or to define safety coefficients in seismic assessment or strengthening design. This study is also fundamental for code calibration purposes as it shows the relationship between the seismic hazard characteristics of Portugal with masonry buildings structural performance, taking the representative existing pre-code masonry building stock.