Seismic Risk Assessment Using Stochastic Nonlinear Models

: The basic input when seismic risk is estimated in urban environments is the expected physical damage level of buildings. The vulnerability index and capacity spectrum-based methods are the tools that have been used most to estimate the probability of occurrence of this important variable. Although both methods provide adequate estimates, they involve simpliﬁcations that are no longer necessary, given the current capacity of computers. In this study, an advanced method is developed that avoids many of these simpliﬁcations. The method starts from current state-of-the-art approaches, but it incorporates non-linear dynamic analysis and a probabilistic focus. Thus, the method considers not only the nonlinear dynamic response of the structures, modeled as multi degree of freedom systems (MDoF), but also uncertainties related to the loads, the geometry of the buildings, the mechanical properties of the materials and the seismic action. Once the method has been developed, the buildings are subjected to earthquake records that are selected and scaled according to the seismic hazard of the site and considering the probabilistic nature of the seismic actions. The practical applications of the method are illustrated with a case study: framed reinforced concrete buildings that are typical of an important district, the Eixample , in Barcelona (Spain). The building typology and the district were chosen because the seismic risk in Barcelona has been thoroughly studied, so detailed information about buildings’ features, seismic hazard and expected risk is available. Hence, the current results can be compared with those obtained using simpler, less sophisticated methods. The main aspects of the method are presented and discussed ﬁrst. Then, the case study is described and the results obtained with the capacity spectrum method are compared with the results using the approach presented here. The results at hand show reasonably good agreement with previous seismic damage and risk scenarios in Barcelona, but the new method provides richer, more detailed, more reliable information. This is particularly useful for seismic risk reduction, prevention and management, to move towards more resilient, sustainable cities.


Introduction
Several ongoing research projects are focused on the mitigation of seismic risk. One of these is the Kairos project [1], which is aimed to maintain and increase the resilience and sustainability of communities against earthquakes. One research area in the Kairos project is the development of advanced numerical tools to assess the seismic risk of structures. The current capacity of computers, combined with advanced statistical approaches, means that new probabilistic frameworks can be developed for estimating the stochastic nonlinear response of civil structures submitted to seismic actions. Vargas-Alzate et al. (2019) [2] presented the development through implementation of the structures. An algorithm based on the Monte Carlo and Latin hypercube sampling techniques has been developed to perform the probabilistic simulations. Ruaumoko software is used for the structural analyses [15].

Description of the Building Typology
Because of the low seismic hazard in Barcelona, buildings constructed before the last decade of the twentieth century were not subjected to any seismic design code. Hence, in most of the buildings in the city, seismic loads were not considered. For instance, most of the RC structures built in Barcelona have waffle slabs that are directly connected to the columns, i.e., the strong-column-weak-beam concept is not satisfied. This design practice avoids progressive collapse of RC structures due to a cascade effect created by column failure in the lower levels. In addition, as the columns are designed solely to withstand gravity loads, their cross sections tend to be smaller than those in seismic areas. This makes the fundamental periods of these buildings longer than those of analogous structures in areas where the seismic hazard is relevant. Figure 1 shows a sketch of an RC waffle slab typical of buildings in Barcelona.
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 21 been developed to perform the probabilistic simulations. Ruaumoko software is used for the structural analyses [15].

Description of the Building Typology
Because of the low seismic hazard in Barcelona, buildings constructed before the last decade of the twentieth century were not subjected to any seismic design code. Hence, in most of the buildings in the city, seismic loads were not considered. For instance, most of the RC structures built in Barcelona have waffle slabs that are directly connected to the columns, i.e., the strong-column-weakbeam concept is not satisfied. This design practice avoids progressive collapse of RC structures due to a cascade effect created by column failure in the lower levels. In addition, as the columns are designed solely to withstand gravity loads, their cross sections tend to be smaller than those in seismic areas. This makes the fundamental periods of these buildings longer than those of analogous structures in areas where the seismic hazard is relevant. Figure 1 shows a sketch of an RC waffle slab typical of buildings in Barcelona. These features make the buildings' vulnerability high, even in a low-to-moderate earthquake. This means that, despite the low-to-moderate seismic hazard, seismic risk in Barcelona is expected to be significant, mainly due to a highly vulnerable environment and the high value of exposed property and goods. Moreover, most of the buildings in this city are made from unreinforced masonry, which is an even more vulnerable structural typology than RC buildings.

Building Distribution
RC buildings in the Eixample district are used as a testbed. Figure 2 shows the distribution of the number of stories of these buildings, the spatial distribution of RC buildings and the seismic zonation of soils. Data for buildings were provided by the City Council's Municipal Institute of Information Technology (IMI). In previous seismic risk studies in Barcelona [5][6][7][8], buildings were grouped into three classes: low-, mid-and high-rise, called RCL, RCM and RCH, with 1 to 3, 4 to 6 and over 6 stories, respectively. Notably, most of the buildings in this district range from 6 to 11 stories, because buildings with 5 or fewer stories in Barcelona are mainly made from unreinforced masonry. The 3 buildings with over 20 stories were excluded from the simulation, because they are special cases and probably have RC wall cores. These features make the buildings' vulnerability high, even in a low-to-moderate earthquake. This means that, despite the low-to-moderate seismic hazard, seismic risk in Barcelona is expected to be significant, mainly due to a highly vulnerable environment and the high value of exposed property and goods. Moreover, most of the buildings in this city are made from unreinforced masonry, which is an even more vulnerable structural typology than RC buildings.

Building Distribution
RC buildings in the Eixample district are used as a testbed. Figure 2 shows the distribution of the number of stories of these buildings, the spatial distribution of RC buildings and the seismic zonation of soils. Data for buildings were provided by the City Council's Municipal Institute of Information Technology (IMI). In previous seismic risk studies in Barcelona [5][6][7][8], buildings were grouped into three classes: low-, mid-and high-rise, called RCL, RCM and RCH, with 1 to 3, 4 to 6 and over 6 stories, respectively. Notably, most of the buildings in this district range from 6 to 11 stories, because buildings with 5 or fewer stories in Barcelona are mainly made from unreinforced masonry. The 3 buildings with over 20 stories were excluded from the simulation, because they are special cases and probably have RC wall cores.

Modeling Considerations
Typical RC frame members of buildings designed in low seismic hazard areas have structural deficiencies against horizontal dynamic loads. For example, they have poor longitudinal and transversal reinforcement and/or inadequate sizes of cross-sections. These shortcomings may be responsible for significant strength and stiffness degradation when non-linear dynamic behavior occurs. To allow for strength degradation, the elastic limits in interaction diagrams must be reduced as a function of the ductility and the number of cycles of inelastic behavior. Figure 3 shows a function for reducing the yield strength due to the ductility reached by a structural element. In this figure, µ1, µ2 and µ3 are related to the ductilities at which the strength degradation begins (µ1) and stops (µ2), and when the strength achieves 0.01 of its initial value (µ3). R1 is the residual strength, which is a fraction of the initial yield strength whose value has been considered equal to 0.4. Because of the low performance of the structural elements, µ1, µ2 and µ3 are in the order of 1.1, 2.5 and 5, respectively. Higher values would be expected in a seismic protected urban environment, due to earthquake resistant design and better seismic performance. In addition, a reduction factor of strength per cycle of inelastic behavior equal to 0.9 has been considered. To allow for stiffness degradation, the degrading bilinear hysteresis was used [16]. The α degrading factor is equal to 0.4. The yielding surfaces are defined by the bending moment-axial load interaction diagram for columns and bending moment-curvature for beams. Note that in the case of beams, an equivalent flat beam element, whose depth is equal to that of the waffle slab, has been considered. A frame element, where its nonlinear behavior is concentrated at both ends, has been used for modeling the column and the equivalent beam elements [17].

Modeling Considerations
Typical RC frame members of buildings designed in low seismic hazard areas have structural deficiencies against horizontal dynamic loads. For example, they have poor longitudinal and transversal reinforcement and/or inadequate sizes of cross-sections. These shortcomings may be responsible for significant strength and stiffness degradation when non-linear dynamic behavior occurs. To allow for strength degradation, the elastic limits in interaction diagrams must be reduced as a function of the ductility and the number of cycles of inelastic behavior. Figure 3 shows a function for reducing the yield strength due to the ductility reached by a structural element. In this figure, µ 1 , µ 2 and µ 3 are related to the ductilities at which the strength degradation begins (µ 1 ) and stops (µ 2 ), and when the strength achieves 0.01 of its initial value (µ 3 ). R 1 is the residual strength, which is a fraction of the initial yield strength whose value has been considered equal to 0.4. Because of the low performance of the structural elements, µ 1 , µ 2 and µ 3 are in the order of 1.1, 2.5 and 5, respectively. Higher values would be expected in a seismic protected urban environment, due to earthquake resistant design and better seismic performance. In addition, a reduction factor of strength per cycle of inelastic behavior equal to 0.9 has been considered. To allow for stiffness degradation, the degrading bilinear hysteresis was used [16]. The α degrading factor is equal to 0.4. The yielding surfaces are defined by the bending moment-axial load interaction diagram for columns and bending moment-curvature for beams. Note that in the case of beams, an equivalent flat beam element, whose depth is equal to that of the waffle slab, has been considered. A frame element, where its nonlinear behavior is concentrated at both ends, has been used for modeling the column and the equivalent beam elements [17].

Preliminary Statistical Considerations
When a probabilistic approach is used to characterize exposure, two types of random variables must be distinguished: those that refer to the building stock (for instance, at neighborhood, district or city levels) and those that refer to intrinsic properties of a single building. Consideration of these types of random variables depends on whether the risk assessments are oriented to earthquake scenarios in urban environments (structural classes related to the building typology matrix [BTM]) or they concern individual buildings. Modeling the random variation in building-to-building structural characteristics within a structural typology is standard practice. This modeling must reflect the epistemic uncertainty and how many structural types are grouped together into a structural typology. Usually, structures sharing similar features are grouped into the same structural typology.
One of the most used properties to classify a structure within a typology is the number of stories. For instance, in the Risk-UE project [10,11], RC buildings were classified as low-rise (1-3 stories), midrise (4-6 stories) and high-rise (>6 stories). Of course, many other features are used to identify a structural typology. This classification scheme simplifies the characterization of the exposure when the seismic damage is estimated at urban level. However, using the approach presented in Vargas-Alzate et al. (2019), [2], the response of buildings can be analyzed without grouping them into height categories as low-, mid-and high-rise buildings, since many building models can be generated, varying the numbers of stories as necessary. Thus, the structural typologies in this study are classified as RCT where T is the number of stories.
According to the distribution presented in Figure 2, structures whose number of stories range from 1 to 14 should be analyzed. Two approaches can be taken to consider uncertainties: (i) generate the exact number of building models according to the exact/deterministic distribution presented in Figure 2 (ii) generate the same fixed number, k, of building samples in each building class defined by the number of stories, i.e., k samples of 1 story, k samples of 2 stories, and so on. Following the first approach, the expected distribution of damage could be analyzed directly, in the case of a single earthquake. This first approach would also permit an analysis of how much the distribution of damage would vary with the frequency content of the seismic actions.
In both approaches, the uncertainties of the parameters can be considered when building samples are generated. However, in the first case, when the number of buildings with a specific number of stories is small, for instance, buildings with 4 stories and others in Figure 2, the uncertainty cannot be adequately represented because of the low number of samples. Instead, the second

Preliminary Statistical Considerations
When a probabilistic approach is used to characterize exposure, two types of random variables must be distinguished: those that refer to the building stock (for instance, at neighborhood, district or city levels) and those that refer to intrinsic properties of a single building. Consideration of these types of random variables depends on whether the risk assessments are oriented to earthquake scenarios in urban environments (structural classes related to the building typology matrix [BTM]) or they concern individual buildings. Modeling the random variation in building-to-building structural characteristics within a structural typology is standard practice. This modeling must reflect the epistemic uncertainty and how many structural types are grouped together into a structural typology. Usually, structures sharing similar features are grouped into the same structural typology.
One of the most used properties to classify a structure within a typology is the number of stories. For instance, in the Risk-UE project [10,11], RC buildings were classified as low-rise (1-3 stories), mid-rise (4-6 stories) and high-rise (>6 stories). Of course, many other features are used to identify a structural typology. This classification scheme simplifies the characterization of the exposure when the seismic damage is estimated at urban level. However, using the approach presented in Vargas-Alzate et al. (2019) [2], the response of buildings can be analyzed without grouping them into height categories as low-, mid-and high-rise buildings, since many building models can be generated, varying the numbers of stories as necessary. Thus, the structural typologies in this study are classified as RC T where T is the number of stories.
According to the distribution presented in Figure 2, structures whose number of stories range from 1 to 14 should be analyzed. Two approaches can be taken to consider uncertainties: (i) generate the exact number of building models according to the exact/deterministic distribution presented in Figure 2 (ii) generate the same fixed number, k, of building samples in each building class defined by the number of stories, i.e., k samples of 1 story, k samples of 2 stories, and so on. Following the first approach, the expected distribution of damage could be analyzed directly, in the case of a single earthquake. This first approach would also permit an analysis of how much the distribution of damage would vary with the frequency content of the seismic actions.
In both approaches, the uncertainties of the parameters can be considered when building samples are generated. However, in the first case, when the number of buildings with a specific number of stories is small, for instance, buildings with 4 stories and others in Figure 2, the uncertainty cannot be adequately represented because of the low number of samples. Instead, the second approach allows uncertainties to be considered in a more adequate manner. In this study, the second approach was preferred.

Building to Building Variation
Since the distribution of the number of stories is known, other properties related to the number of spans, N sp , the story height, H st , and the span length, S l , are considered random variables. Thus, N sp follows a uniform, discrete distribution in the interval [3,6]. H st is distributed uniformly in the interval [2.8, 3.2] m. S l is distributed uniformly in the interval [4,6] m. These ranges are typical in the RC buildings in Barcelona and they were estimated after analyzing several blueprints and visiting some buildings in the district. To define the cross sections of structural elements of the first story, several functions are used that are based on design criteria. Again, values observed in RC buildings of the Eixample district were used to validate the cross sections' dimensions in the structural models. Based on these models, the width, W c , and depth, D c , of the columns of the first story will depend on the number of stories, N st , and the span length, S l . The following equation was used to calculate W c or D c : where c i are coefficients that may be adjusted depending on the data distribution of the analyzed area. For this study, c 1 = 0.45, c 2 = 0.05, c 3 = 0.01 and c 4 = −0.6. Φ 1,0 is the standard normal distribution. Note that the columns are not necessarily square, that is, one random sample is generated for the width and one for the depth of the columns according to Equation (1). Figure 4 shows W c for the first story columns according to the conditions described above. For upper stories, the size of columns decreases systematically by 5 cm every three stories. The values generated are rounded to the nearest multiple of 5 cm. Note that Equation (1) may produce very small or even negative values for buildings whose N st is lower than 5, which makes no sense (blue dots). This issue has been solved by adding the condition that W c or D c cannot be lower than 0.25 m (red dots).
Sustainability 2020, 12, x FOR PEER REVIEW 6 of 21 approach allows uncertainties to be considered in a more adequate manner. In this study, the second approach was preferred.

Building to Building Variation
Since the distribution of the number of stories is known, other properties related to the number of spans, Nsp, the story height, Hst, and the span length, Sl, are considered random variables. Thus, Nsp follows a uniform, discrete distribution in the interval [3,6]. Hst is distributed uniformly in the interval [2.8, 3.2] m. Sl is distributed uniformly in the interval [4,6] m. These ranges are typical in the RC buildings in Barcelona and they were estimated after analyzing several blueprints and visiting some buildings in the district. To define the cross sections of structural elements of the first story, several functions are used that are based on design criteria. Again, values observed in RC buildings of the Eixample district were used to validate the cross sections' dimensions in the structural models. Based on these models, the width, Wc, and depth, Dc, of the columns of the first story will depend on the number of stories, Nst, and the span length, Sl. The following equation was used to calculate Wc or Dc: where ci are coefficients that may be adjusted depending on the data distribution of the analyzed area. For this study, 1 = 0.45, 2 = 0.05, 3 = 0.01 and 4 = −0.6. Φ1,0 is the standard normal distribution. Note that the columns are not necessarily square, that is, one random sample is generated for the width and one for the depth of the columns according to Equation (1). Figure 4 shows Wc for the first story columns according to the conditions described above. For upper stories, the size of columns decreases systematically by 5 cm every three stories. The values generated are rounded to the nearest multiple of 5 cm. Note that Equation (1) may produce very small or even negative values for buildings whose is lower than 5, which makes no sense (blue dots). This issue has been solved by adding the condition that Wc or Dc cannot be lower than 0.25 m (red dots). In Barcelona, it is very common that the first story of RC buildings to be higher than the upper stories, so most of the buildings have soft first stories. This is common practice because first stories in the Eixample district are normally intended for commercial use. To consider this effect, the height of the first story of the simulated models has been multiplied by a random value that is obtained from the Gaussian distribution shown in Figure 5. To assign the steel percentage of the columns, ρc, a continuous Gaussian distribution is assumed whose mean and standard deviation values are 1% and 0.1%, respectively. In Barcelona, it is very common that the first story of RC buildings to be higher than the upper stories, so most of the buildings have soft first stories. This is common practice because first stories in the Eixample district are normally intended for commercial use. To consider this effect, the height of the first story of the simulated models has been multiplied by a random value that is obtained from the Gaussian distribution shown in Figure 5. To assign the steel percentage of the columns, ρ c , a continuous Gaussian distribution is assumed whose mean and standard deviation values are 1% and 0.1%, respectively. The waffle slabs have been modeled by means of flat beam elements whose width, Wb, has been fixed to 0.7 m and whose depth, Db, depends on the span length, Sl. The following stepwise function has been used to calculate Db: In this equation, m stands for meter units. Figure 6 depicts the depth of the beams as a function of the span length. No random term was considered for the beams. To assign the steel percentage of the beams, ρb, a continuous Gaussian distribution is assumed whose mean and standard deviation values are 0.6% and 0.05%, respectively. In fact, either the coefficients or the type of equation considered for assigning the cross sections to columns and beams can be modified, so that they better represent the urban environment studied. The values presented herein are representative for the structural conditions of RC buildings in Barcelona.  The waffle slabs have been modeled by means of flat beam elements whose width, W b , has been fixed to 0.7 m and whose depth, D b , depends on the span length, S l . The following stepwise function has been used to calculate D b : ( In this equation, m stands for meter units. Figure 6 depicts the depth of the beams as a function of the span length. No random term was considered for the beams. To assign the steel percentage of the beams, ρ b , a continuous Gaussian distribution is assumed whose mean and standard deviation values are 0.6% and 0.05%, respectively. In fact, either the coefficients or the type of equation considered for assigning the cross sections to columns and beams can be modified, so that they better represent the urban environment studied. The values presented herein are representative for the structural conditions of RC buildings in Barcelona. The waffle slabs have been modeled by means of flat beam elements whose width, Wb, has been fixed to 0.7 m and whose depth, Db, depends on the span length, Sl. The following stepwise function has been used to calculate Db: In this equation, m stands for meter units. Figure 6 depicts the depth of the beams as a function of the span length. No random term was considered for the beams. To assign the steel percentage of the beams, ρb, a continuous Gaussian distribution is assumed whose mean and standard deviation values are 0.6% and 0.05%, respectively. In fact, either the coefficients or the type of equation considered for assigning the cross sections to columns and beams can be modified, so that they better represent the urban environment studied. The values presented herein are representative for the structural conditions of RC buildings in Barcelona.   Several random variables should be considered when the seismic behavior of a group of buildings is modeled. In this study, the live loads, LL, the superimposed loads, SL, the compressive strength of the concrete, fc, the tensile strength of the steel, fy, the elastic modulus of the concrete, Ec, and the elastic modulus of the steel, Es, were considered random variables. A continuous Gaussian distribution is assumed for these variables. The mean values, µ, and standard deviations, σ, are shown in Table 1. The correlation between pairs of random variables has also been considered. The Monte Carlo method allows modeling of complex systems with a large number of random parameters [18]. This method has been used to generate random samples of buildings according to the considerations described above. The algorithm implemented in Vargas-Alzate et al. (2019) [2], has been adapted to this end. Figure 7 shows one hundred 7-story building models generated with this algorithm. This figure shows the variation in geometrical properties of the models within the considered intervals. Several random variables should be considered when the seismic behavior of a group of buildings is modeled. In this study, the live loads, LL, the superimposed loads, SL, the compressive strength of the concrete, fc, the tensile strength of the steel, fy, the elastic modulus of the concrete, Ec, and the elastic modulus of the steel, Es, were considered random variables. A continuous Gaussian distribution is assumed for these variables. The mean values, µ, and standard deviations, σ, are shown in Table 1. The correlation between pairs of random variables has also been considered. The Monte Carlo method allows modeling of complex systems with a large number of random parameters [18]. This method has been used to generate random samples of buildings according to the considerations described above. The algorithm implemented in Vargas-Alzate et al. (2019), [2], has been adapted to this end. Figure 7 shows one hundred 7-story building models generated with this algorithm. This figure shows the variation in geometrical properties of the models within the considered intervals.  Figure 8 shows the first three periods, T1, T2 and T3, of the simulated models as a function of building height. The fundamental period would be expected to increase with building height. However, for building models whose number of stories is equal or lower than 4, with a height lower than 15 m, the increment in fundamental period follows a different pattern to buildings with over 4 stories. This can be explained considering that the size of the column elements does not change as long as the number of stories is lower than 4 (see Figure 4).  Figure 8 shows the first three periods, T 1 , T 2 and T 3 , of the simulated models as a function of building height. The fundamental period would be expected to increase with building height. However, for building models whose number of stories is equal or lower than 4, with a height lower than 15 m, the increment in fundamental period follows a different pattern to buildings with over 4 stories. This can be explained considering that the size of the column elements does not change as long as the number of stories is lower than 4 (see Figure 4).

Uniform Hazard Spectrum
Seismic hazard is probably the main source of uncertainty in assessments of expected seismic damage in urban environments. PSHA allows dealing with this uncertainty to produce an explicit description of the distribution of future earthquakes that may occur at a specific site [4]. One of the most important outcomes of a PSHA is the uniform hazard spectrum (UHS). The 'uniformness' of this spectrum means that the rate of it being exceeded is the same for every ordinate. The reciprocal of this rate of exceedance is the return period. Thus, if the rate of exceedance is known for a wide range of periods of vibration, then the UHS can be calculated. This spectrum is considered an envelope of the intensity measure (IM) under consideration. One of the most common variables when calculating the UHS is spectral acceleration. Salgado et al. (2015) [14] recently performed a PSHA for Barcelona. These authors obtained the UHS at bedrock level for several return periods (Figure 9a).
(a) (b) Figure 9. Uniform hazard spectra for return periods for the Eixample district, at bedrock level [14] (a) and with soil effects (b).

Uniform Hazard Spectrum
Seismic hazard is probably the main source of uncertainty in assessments of expected seismic damage in urban environments. PSHA allows dealing with this uncertainty to produce an explicit description of the distribution of future earthquakes that may occur at a specific site [4]. One of the most important outcomes of a PSHA is the uniform hazard spectrum (UHS). The 'uniformness' of this spectrum means that the rate of it being exceeded is the same for every ordinate. The reciprocal of this rate of exceedance is the return period. Thus, if the rate of exceedance is known for a wide range of periods of vibration, then the UHS can be calculated. This spectrum is considered an envelope of the intensity measure (IM) under consideration. One of the most common variables when calculating the UHS is spectral acceleration. Salgado et al. (2015) [14] recently performed a PSHA for Barcelona. These authors obtained the UHS at bedrock level for several return periods (Figure 9a).

Uniform Hazard Spectrum
Seismic hazard is probably the main source of uncertainty in assessments of expected seismic damage in urban environments. PSHA allows dealing with this uncertainty to produce an explicit description of the distribution of future earthquakes that may occur at a specific site [4]. One of the most important outcomes of a PSHA is the uniform hazard spectrum (UHS). The 'uniformness' of this spectrum means that the rate of it being exceeded is the same for every ordinate. The reciprocal of this rate of exceedance is the return period. Thus, if the rate of exceedance is known for a wide range of periods of vibration, then the UHS can be calculated. This spectrum is considered an envelope of the intensity measure (IM) under consideration. One of the most common variables when calculating the UHS is spectral acceleration. Salgado et al. (2015) [14] recently performed a PSHA for Barcelona. These authors obtained the UHS at bedrock level for several return periods (Figure 9a).
(a) (b) Figure 9. Uniform hazard spectra for return periods for the Eixample district, at bedrock level [14] (a) and with soil effects (b). Figure 9. Uniform hazard spectra for return periods for the Eixample district, at bedrock level [14] (a) and with soil effects (b).

Soil Considerations
The UHS presented in Figure 9a is calculated at bedrock level. However, the Eixample district lies on soil type 2, according to the description given in Cid et al. (2001), [19] (see Figure 2). For this soil type, the Vs 30 is 324 m/s. Thus, the UHS presented in Figure 9a has been multiplied by the soil amplification factor function shown in Figure 10. The UHS with soil effects (Figure 9b) have been used in this study, as they represent the seismic hazard of the Eixample district. The UHS presented in Figure 9a is calculated at bedrock level. However, the Eixample district lies on soil type 2, according to the description given in Cid et al. (2001), [19] (see Figure 2). For this soil type, the Vs30 is 324 m/s. Thus, the UHS presented in Figure 9a has been multiplied by the soil amplification factor function shown in Figure 10. The UHS with soil effects (Figure 9b) have been used in this study, as they represent the seismic hazard of the Eixample district.

Selecting Strong Ground Motions
Once the expected seismic hazard is known, the next point is to obtain a group of acceleration records compatible with the UHS. This task can be done either for the entire range of periods or for selected intervals, depending on the structures. In this study, a group of earthquakes has been selected using a modified version of the procedure presented in Vargas-Alzate et al. (2013) [20]. In that study, the authors proposed selecting from a strong ground motion database a group of acceleration records that have a mean of their spectra similar to a target spectrum. To achieve this, they processed the database of Ambraseys et al. (2004) [21]. Then, an algorithm was developed to sort and select accelerograms. This algorithm was based on minimizing the mean square error to obtain a representative group of accelerograms for several target spectra.
Vargas-Alzate (2013) [22] showed that this procedure allows the selection of earthquakes that are compatible with whichever design spectra is proposed in Eurocode 8, CEN (2004) [23]. Notably, in the procedure proposed by Vargas-Alzate et al. (2013) [20], the PGAs of the database records and the PGAs of the target spectra were previously normalized to '1'. This normalization helped to match any spectral shape. However, this normalization could also lead to records from earthquakes of low magnitude being scaled up to represent strong motion records from earthquakes of high magnitude and vice versa. Thus, normalization of PGAs has been omitted herein, and the selected accelerograms naturally fit the remaining conditions presented in Vargas-Alzate et al. (2013) [20]. Based on this procedure, one hundred earthquake records have been selected for each return period. In Figure 11a, the spectra of 100 records, their mean spectrum and the UHS for a return period, Tr, of 475 years are presented. In this figure, the good agreement between the mean and the UHS can be seen.

Scaling
The mean spectrum of the selected earthquakes fits the UHS well. However, in order to fit the UHS in periods close to the fundamental period, T1, the earthquake records are scaled by equating

Selecting Strong Ground Motions
Once the expected seismic hazard is known, the next point is to obtain a group of acceleration records compatible with the UHS. This task can be done either for the entire range of periods or for selected intervals, depending on the structures. In this study, a group of earthquakes has been selected using a modified version of the procedure presented in Vargas-Alzate et al. (2013) [20]. In that study, the authors proposed selecting from a strong ground motion database a group of acceleration records that have a mean of their spectra similar to a target spectrum. To achieve this, they processed the database of Ambraseys et al. (2004) [21]. Then, an algorithm was developed to sort and select accelerograms. This algorithm was based on minimizing the mean square error to obtain a representative group of accelerograms for several target spectra.
Vargas-Alzate (2013) [22] showed that this procedure allows the selection of earthquakes that are compatible with whichever design spectra is proposed in Eurocode 8, CEN (2004) [23]. Notably, in the procedure proposed by Vargas-Alzate et al. (2013) [20], the PGAs of the database records and the PGAs of the target spectra were previously normalized to '1'. This normalization helped to match any spectral shape. However, this normalization could also lead to records from earthquakes of low magnitude being scaled up to represent strong motion records from earthquakes of high magnitude and vice versa. Thus, normalization of PGAs has been omitted herein, and the selected accelerograms naturally fit the remaining conditions presented in Vargas-Alzate et al. (2013) [20]. Based on this procedure, one hundred earthquake records have been selected for each return period. In Figure 11a, the spectra of 100 records, their mean spectrum and the UHS for a return period, Tr, of 475 years are presented. In this figure, the good agreement between the mean and the UHS can be seen.

Scaling
The mean spectrum of the selected earthquakes fits the UHS well. However, in order to fit the UHS in periods close to the fundamental period, T 1 , the earthquake records are scaled by equating the average of the UHS with the average of the spectrum of each record, called AvSa, in an interval around T 1 . The low and high limits of this interval are 0.1 T 1 and 1.8 T 1 . Figure 11b shows an example of this scaling procedure for a period of vibration equal to 1 s using the spectra of the accelerograms selected for Tr = 475 years. This scaling procedure is often used to estimate fragility functions for a single structure, since AvSa tends to be highly correlated with the MIDR, which is an engineering demand parameter (EDP) that is also highly correlated with seismic damage (Eads et al. (2015) [24] and Kohrangi et al. (2017) [25]). the average of the UHS with the average of the spectrum of each record, called AvSa, in an interval around T1. The low and high limits of this interval are 0.1 T1 and 1.8 T1. Figure 11b shows an example of this scaling procedure for a period of vibration equal to 1 s using the spectra of the accelerograms selected for Tr = 475 years. This scaling procedure is often used to estimate fragility functions for a single structure, since AvSa tends to be highly correlated with the MIDR, which is an engineering demand parameter (EDP) that is also highly correlated with seismic damage (Eads et al. (2015) [24] and Kohrangi et al. (2017) [25]).

NLDA
NLDA is a kind of structural analysis that can be used to simulate the nonlinear dynamic response of a structure subjected to strong ground motions. NLDA can be used to calculate the envelope of several structural response variables such as the displacement at the roof, the MIDR and the global damage index, amongst others. NLDA is the most rigorous existing numerical tool to estimate the dynamic response of a structure. However, in the authors' view, the use of this procedure is not as widespread in the engineering community as it should be. This is mainly because the building modeling and characterization of the seismic hazard by means of acceleration time histories are not easy tasks. Moreover, NLDA is time consuming when it is used to perform probabilistic estimations of seismic damage at urban level. In any case, the extended use of probabilistic NLDA, i.e., considering uncertainties, should be a priority in the earthquake engineering community. In this study, probabilistic NLDA has been used to analyze the seismic response of the structures in detail.

Park and Ang Damage Index
The global damage index of Park and Ang, DIPA, can be calculated from the output of an NLDA [26]. When the building goes into nonlinear behavior, DIPA becomes positive and can be obtained as the weighted sum of the damage indices at element level, DIE: where the weights, , are the ratios of the hysteretic energy dissipated by each element, E, to the total hysteretic energy dissipated by the structure [27] and DIE is defined by the following equation:

NLDA
NLDA is a kind of structural analysis that can be used to simulate the nonlinear dynamic response of a structure subjected to strong ground motions. NLDA can be used to calculate the envelope of several structural response variables such as the displacement at the roof, the MIDR and the global damage index, amongst others. NLDA is the most rigorous existing numerical tool to estimate the dynamic response of a structure. However, in the authors' view, the use of this procedure is not as widespread in the engineering community as it should be. This is mainly because the building modeling and characterization of the seismic hazard by means of acceleration time histories are not easy tasks. Moreover, NLDA is time consuming when it is used to perform probabilistic estimations of seismic damage at urban level. In any case, the extended use of probabilistic NLDA, i.e., considering uncertainties, should be a priority in the earthquake engineering community. In this study, probabilistic NLDA has been used to analyze the seismic response of the structures in detail.

Park and Ang Damage Index
The global damage index of Park and Ang, DI PA , can be calculated from the output of an NLDA [26]. When the building goes into nonlinear behavior, DI PA becomes positive and can be obtained as the weighted sum of the damage indices at element level, DI E : where the weights, λ i , are the ratios of the hysteretic energy dissipated by each element, E, to the total hysteretic energy dissipated by the structure [27] and DI E is defined by the following equation: where µ m and µ u are the maximum and ultimate ductilities, respectively, β is a non-negative parameter that considers the contribution of cyclic loading on structural damage, E h is the dissipated hysteretic energy, F y is the yield load and δ y is the yield displacement.

Damage States
Damage states, DS K , represent the level of affectation of a structure after a catastrophic event. DS K can vary from null, which indicates that the structure was not damaged at all (i.e., K = 0), to collapse, which indicates that the structure has lost its capacity to withstand gravity loads. In the middle are other DS K associated with different damage grades. Obviously, damage also affects the functionality of the structure so that these intermediate DS K are strongly related to the resilience level of an urban environment, since they are strongly correlated with the downtime of a structure after an earthquake.
In earthquake engineering, there are many types and definitions of DS K , which are generally related to EDPs, like MIDR, displacement at the roof, maximum story acceleration, amongst many others. One of the most common EDP when an NLDA is used is the MIDR. However, MIDR only indicates the maximum damage observed at the most affected story. This means that the damage level observed in other stories is not considered within this EDP. In contrast, DI PA provides an estimation of the overall damage which is more strongly correlated with the general health of the structure. Thus, the DI PA has been used as an EDP to define the DS K of the structures analyzed after an earthquake. The limits presented in Table 2 are used to define the DS K thresholds for a structure. These limits are classical values of the damage index of Park and Ang [27]. In addition, Table 2 presents the MIDR that produces similar damage probabilities to those shown later on in Section 4.

Simulation Scheme
After an explanation of how the exposure and the hazard have been characterized and how the damage has been defined, the next step is to perform probabilistic simulations. Recall that five scenarios and fourteen building types have been defined. For each typology, one hundred building samples have been generated considering uncertainties. Moreover, for each scenario, one hundred accelerograms matching the corresponding UHS were selected and scaled. Then, to perform the probabilistic analysis for each building type and for each earthquake scenario, the following steps are followed. STEP 1: for each building sample, one of the one hundred available accelerograms is randomly selected. Then, the NLDA is performed to obtain the DI PA and the corresponding damage state, according to the definitions in Table 2. As a result of this first step, one hundred damage states were obtained corresponding to the building type and earthquake scenario. STEP 2: the probabilities of each damage state K and typology T, P DS K,T , are obtained, simply by counting the buildings with damage state DS K . Then, DSm T , which can be interpreted as the mean damage state for each typology, is defined by the following equation: where n is the number of non-null damage states, which is 5 in this case. Thus, DSm T may be used as a general indicator of the damage expected within a structural typology. In order to obtain an overall indicator of the expected damage in the Eixample district, the following equation can be used: In this equation, GDSm is the global indicator of the expected damage and N T represents the number of buildings of typology T. Recall that this procedure is iterated for the five earthquake scenarios taken from Salgado et al. (2015) [14]. Thus, for each scenario, 1400 NLDAs have been performed (one hundred building samples per fourteen structural typologies). This procedure can easily be performed with a modified version of the algorithm presented in Vargas-Alzate et al. 2019 [2]. Tables S1-S5 found in the Supplementary present the probabilities of occurrence of each DSm T , as well as the mean damage state for each typology, DSm T , for the return periods presented in Figure 9 and for the structural typologies RC 1 to RC 14 . Each Supplementary Table also contains a graphical representation of the DSm T , spatially distributed building by building.
GDSm can be used to denote how strongly an earthquake affects a specific place. Hence, this variable can be compared with the EMS-98 intensity values when some damage exists. For instance, for the return period of Tr = 31 years, a few buildings are slightly damaged. Therefore, this is equivalent to intensity VI in EMS-98 (slightly damaging). In general, equivalence can be calculated between GDSm and EMS-98 by assuming a linear variation of GDSm from 0 to 5 in the interval VI to XII of EMS-98. This equivalence is only valid for GDSm above 0. When GDSm is equal to zero, this variable cannot describe the EMS-98 intensities below VI since they are not related to structural damage. Thus, for the return periods 225, 475, 975 and 2475 years, the EMS-98 intensity values are VI, VII, VIII and X, respectively. The intensity values were obtained after rounding those calculated by linear transformation to the nearest integer. This is because the EMS-98 scale only admits discrete values.

Level 2 Method, LM2
Lantada et al. (2010) [6] and Irizarry et al. (2011) [8] estimated the seismic risk of Barcelona using LM1 and LM2 respectively, which were developed in the Risk-UE project. A comparison between these approaches can be found in Lantada et al. (2009) [7]. In LM2, RC buildings of Barcelona were classified into three types: low-, mid-and high-rise, called RCL, RCM and RCH, respectively. RCL included buildings with 1 to 3 stories, RCM 4 to 6 and RCH more than 6 stories. Only one specific building was modeled for each height range, assuming that it was representative of all the buildings in the category. Thus, buildings were characterized by bilinear capacity spectra. Figure 12a shows the three bilinear capacity spectra, which represent the mean seismic behavior of these structural typologies.
According to LM2, fragility curves can be obtained from the bilinear capacity spectra. For a specific damage state DS K , the fragility function provides the probability that this damage state is attained or exceeded. A detailed description of the calculation of fragility functions in LM2 can be found in Milutinoviç and Trendafiloski (2003) [28]. Similar to the damage state definitions presented in Table 2, LM2 proposes 4 non-null damage states identified as slight (DS 1 ), moderate (DS 2 ), severe (DS 3 ) and complete (DS 4 ). Notice that the complete damage state in LM2 groups DS 4 (Complete) and DS 5 (Collapse) from Table 2. Figure 12b,c, and d depict the fragility functions for RCL, RCM and RCH, respectively, for the four non-null damage states. Recall that the fragility function for the null (DS 0 ) damage state is equal to one. Sustainability 2020, 12, x FOR PEER REVIEW 14 of 21 The probability of occurrence of each DSK can be obtained from the fragility curves, given a spectral displacement value [28]. For a specific earthquake scenario, there are several methods to obtain the expected spectral displacement, given a capacity spectrum and a response spectrum function. The expected spectral displacement is generally known as the performance point, pp. The most simplified assumption, and one of the most used in practice, is to estimate the pp based on the equal displacement approximation (EDA) rule. EDA is a well-known empirical rule for the assessment of the non-linear behavior of structures exposed to strong ground motions. This procedure states that the predicted inelastic displacement response of oscillators is often very similar to the predicted elastic displacement response of oscillators with the same period (ATC-40) [29]. In this study, the EDA rule is used to estimate the pp for the capacity spectra of Figure 12a and the UHS of Figure 9b. Once the pp is known, the probability of occurrence of DSK, represented in the fragility curves, can be obtained. Tables S1-S5 found in the Supplementary also show the probabilities of occurrence of each DSK for the return periods presented in Figure 9 and for the structural typologies RCL, RCM and RCH. Note that subscript T, which stands for building type in Equation (4), has been omitted when referring to the probabilities of DSK in LM2.

Comparison of the Results
In this section, the results obtained with LM2 and LM3 are compared and discussed. However, in LM2, the structural typology classification brings together buildings with different numbers of stories. Therefore, to compare both methods, the following is necessary. (i) The probabilities of the DSK obtained with LM3 must be grouped in terms of number of stories. Hence, RC1, RC2 and RC3 will be aggregated so that the results are equivalent to those obtained for the RCL class in LM2. This The probability of occurrence of each DS K can be obtained from the fragility curves, given a spectral displacement value [28]. For a specific earthquake scenario, there are several methods to obtain the expected spectral displacement, given a capacity spectrum and a response spectrum function. The expected spectral displacement is generally known as the performance point, pp. The most simplified assumption, and one of the most used in practice, is to estimate the pp based on the equal displacement approximation (EDA) rule. EDA is a well-known empirical rule for the assessment of the non-linear behavior of structures exposed to strong ground motions. This procedure states that the predicted inelastic displacement response of oscillators is often very similar to the predicted elastic displacement response of oscillators with the same period (ATC-40) [29]. In this study, the EDA rule is used to estimate the pp for the capacity spectra of Figure 12a and the UHS of Figure 9b. Once the pp is known, the probability of occurrence of DS K , represented in the fragility curves, can be obtained. Tables S1-S5 found in the Supplementary also show the probabilities of occurrence of each DS K for the return periods presented in Figure 9 and for the structural typologies RCL, RCM and RCH. Note that subscript T, which stands for building type in Equation (4), has been omitted when referring to the probabilities of DS K in LM2.

Comparison of the Results
In this section, the results obtained with LM2 and LM3 are compared and discussed. However, in LM2, the structural typology classification brings together buildings with different numbers of stories. Therefore, to compare both methods, the following is necessary. (i) The probabilities of the DS K obtained with LM3 must be grouped in terms of number of stories. Hence, RC 1 , RC 2 and RC 3 will be aggregated so that the results are equivalent to those obtained for the RCL class in LM2. This grouping Sustainability 2020, 12, 1308 15 of 21 should be performed for RCM and RCH as well. (ii) Complete (DS 4 ) and collapse (DS 5 ) damage states in LM3 should be aggregated to make the results compatible with those corresponding to complete (DS 4 ) damage state of LM2. In the following, this will be referred to as LM3*, whenever this grouping is made. Figure 13 shows comparisons of the probabilities of occurrence of DS K , named P DS K , for the UHS presented in Figure 9b by using LM2 ((a), (c) and (e), left side) and LM3* ((b), (d) and (f), right side). Again, subscript T has been omitted in reference to the probabilities of the DS K in LM3*. grouping should be performed for RCM and RCH as well. (ii) Complete (DS4) and collapse (DS5) damage states in LM3 should be aggregated to make the results compatible with those corresponding to complete (DS4) damage state of LM2. In the following, this will be referred to as LM3*, whenever this grouping is made. Figure 13 shows comparisons of the probabilities of occurrence of DSK, named , for the UHS presented in Figure   The results depicted in Figure 13 can also be represented and spatially distributed building by building using a GIS. For instance, Figure 14 shows the probability of occurrence of the damage state moderate, P DS 2 for the return period of 475 years by considering the LM3. Figure 13a,b (first line) shows that the probability of the complete damage state in low-rise structures, P DS 4 , is always higher in LM3*. In these structures, the probability of collapse increases because of the low redundancy level. That is, the fewer the number of nodes in the structure, the fewer the number of failure modes. Therefore, the beginning of plasticization can quickly trigger the collapse of the structure. Note that P DS 4 strongly affects the scenarios of victims and the economic costs, which in turn are crucial to the development of strategies for earthquake risk assessments, prevention and management. The results depicted in Figure 13 can also be represented and spatially distributed building by building using a GIS. For instance, Figure 14 shows the probability of occurrence of the damage state moderate, 2 for the return period of 475 years by considering the LM3. Figure 13a,b (first line) shows that the probability of the complete damage state in low-rise structures, 4 , is always higher in LM3*. In these structures, the probability of collapse increases because of the low redundancy level. That is, the fewer the number of nodes in the structure, the fewer the number of failure modes. Therefore, the beginning of plasticization can quickly trigger the collapse of the structure. Note that 4 strongly affects the scenarios of victims and the economic costs, which in turn are crucial to the development of strategies for earthquake risk assessments, prevention and management. From the probabilities presented in Figure 13, the mean damage grade, DSm, can be obtained according to Equation (4). Figure 15 depicts the evolution of DSm as a function of the return period by using both LM2 and LM3* for RCL, RCM and RCH building classes. The above indicates that, on average, the LM2 simplified method provides realistic results. However, it does not reveal, for instance, the differences between the expected damage in a six-story and a twelve-story building. LM3 has better accuracy and resolution and should therefore be preferred. From the probabilities presented in Figure 13, the mean damage grade, DSm, can be obtained according to Equation (4). Figure 15 depicts the evolution of DSm as a function of the return period by using both LM2 and LM3* for RCL, RCM and RCH building classes. The results depicted in Figure 13 can also be represented and spatially distributed building by building using a GIS. For instance, Figure 14 shows the probability of occurrence of the damage state moderate, 2 for the return period of 475 years by considering the LM3. Figure 13a,b (first line) shows that the probability of the complete damage state in low-rise structures, 4 , is always higher in LM3*. In these structures, the probability of collapse increases because of the low redundancy level. That is, the fewer the number of nodes in the structure, the fewer the number of failure modes. Therefore, the beginning of plasticization can quickly trigger the collapse of the structure. Note that 4 strongly affects the scenarios of victims and the economic costs, which in turn are crucial to the development of strategies for earthquake risk assessments, prevention and management. From the probabilities presented in Figure 13, the mean damage grade, DSm, can be obtained according to Equation (4). Figure 15 depicts the evolution of DSm as a function of the return period by using both LM2 and LM3* for RCL, RCM and RCH building classes. The above indicates that, on average, the LM2 simplified method provides realistic results. However, it does not reveal, for instance, the differences between the expected damage in a six-story and a twelve-story building. LM3 has better accuracy and resolution and should therefore be preferred. The above indicates that, on average, the LM2 simplified method provides realistic results. However, it does not reveal, for instance, the differences between the expected damage in a six-story and a twelve-story building. LM3 has better accuracy and resolution and should therefore be preferred.

Discussion and Conclusions
A principal issue in seismic risk assessments concerns the capacity of civil structures to withstand the strong ground motions produced by earthquakes. When these ground motions produce catastrophic events, the resilience and sustainability of affected communities are stretched and compromised [1]. Note that a significant issue for global economic stability concerns the resilience of civil structures, since increasing globalization is causing a redistribution of seismic risk. This means that areas without seismic hazards or with low-seismic risk will be affected if they are economically connected to areas where the seismic risk is high. This is another reason why the estimation and mitigation of seismic risk must be a worldwide effort.
Advanced knowledge of the expected seismic behavior of civil structures will contribute to correct quantification, and indeed to the mitigation of seismic risk. In this respect, the current and increasing capacity of modern computers means that a tremendous amount of numerical analyses can be performed that only a few years ago were costly and time-consuming. This enhanced capability of computer technology allows consideration of the non-linear time history response of sophisticated structural models and the uncertainties related to a number of input variables. A couple of decades ago, this kind of structural analysis was unaffordable, even for a single model of a building. Along with these technological advances, conceptual bases, probabilistic methods and programming tools have been developed to tackle the nonlinear stochastic response of a structure. These advanced tools will contribute to the planning of optimal strategies for reducing seismic risk. In this study, a new method for estimating the probabilistic seismic damage of structures at urban scale has been presented (LM3). In this method, the seismic hazard is characterized by means of actual earthquake records and the exposure by using MDoF systems. The results are compared with those obtained using LM2, which is a simplified method for assessing expected seismic damage of structures.
The main differences between LM2 and LM3 are in the characterization of seismic actions (hazard) and buildings (exposure). In LM2, seismic actions are characterized by response spectra and buildings are defined by means of capacity spectra. In LM3, comprehensive acceleration time histories (earthquake records) and realistic structural models are used. In LM2, the expected seismic response of a structure is obtained by means of the EDA rule, which crosses the capacity and response spectra. In LM3 the expected seismic response is calculated by means of the NLDA. The results obtained with both methods, LM2 and LM3, were compared and presented in Section 4. Several conclusions can be drawn from these comparisons. For instance, for low return periods (31 and 225 years), the expected mean damage states, DSm, are similar with both methods in RCM and RCH structures. In RCL structures, DSm tends to be higher when the LM2 method is used. This proves that for RCL structures, LM2 is more conservative than LM3*. In mid-rise structures, DSm tends to be lower when the LM3* method is used. Conversely, in high-rise structures, DSm tends to be slighter when LM3* is used. This fact can be attributed to the participation of higher modes as well as to the nonlinear response of the structure, which are considered in a more exhaustive manner in LM3. Another factor that may have a significant impact on this increment is that LM3 contemplates p-delta effects. In any case, the differences between LM2 and LM3* for RCM and RCH are lower than those observed in the analysis of RCL structures.
In LM2, related typologies are essentially defined through capacity spectra. Therefore, it would be interesting to compare the capacity spectra obtained with both methods. Figure 16 presents a comparison between the capacity spectra obtained with models generated in this study and those proposed in LM2 for Barcelona. A total of 1400 capacity spectra are obtained, 100 for each structural typology, as in the LM3 calculations. Figure 16 also shows the mean capacity spectra, after grouping the generated models as in LM2. It can be seen that these mean curves are similar to those originally proposed for Barcelona, except for RCH structures. In the case of RCL structures (see Figure 16a), the fundamental period associated with the mean capacity spectrum is lower than the one observed in LM2 (blue line). This may be because in LM3 modeling, for the upmost story, gravity loads are one third of those considered for intermediate stories. In low structures, this reduction significantly decreases the fundamental period of the models. As the EDA rule is applied to estimate the performance points, which in turn will be used to define the probabilities of each damage state, the elastic part of the capacity spectra is fundamental to perform a good estimation of the seismic response. Moreover, due to the low yielding point assigned to the capacity spectrum in LM2, the mean damage state for low performance points is higher than the one obtained with LM3. These aspects could explain the differences observed in Figure 13 when the LM2 and LM3* results are compared for RCL buildings. In the case of RCM structures, the differences between the capacity spectra obtained with both methods are insignificant (at least in the elastic part). This is reflected in Figure 13 when similar estimations are provided by both approaches (LM2 and LM3*). In the case of RCH structures, the mean capacity spectrum differs considerably from that proposed in LM2. Interestingly, both approaches (LM2 and LM3*) provide similar estimations of the DSm as can be observed in Figure 13. In any case, Figure 16c presents the mean capacity spectra (green line) that provide a better approximation to the LM2 capacity spectra for RCH structures. A comparison between all the capacity spectra using the LM3 approach can be seen in Figure 16d. This figure reflects the high scattering inherent to the exposure. performance points, which in turn will be used to define the probabilities of each damage state, the elastic part of the capacity spectra is fundamental to perform a good estimation of the seismic response. Moreover, due to the low yielding point assigned to the capacity spectrum in LM2, the mean damage state for low performance points is higher than the one obtained with LM3. These aspects could explain the differences observed in Figure 13 when the LM2 and LM3* results are compared for RCL buildings. In the case of RCM structures, the differences between the capacity spectra obtained with both methods are insignificant (at least in the elastic part). This is reflected in Figure 13 when similar estimations are provided by both approaches (LM2 and LM3*). In the case of RCH structures, the mean capacity spectrum differs considerably from that proposed in LM2. Interestingly, both approaches (LM2 and LM3*) provide similar estimations of the DSm as can be observed in Figure 13. In any case, Figure 16c presents the mean capacity spectra (green line) that provide a better approximation to the LM2 capacity spectra for RCH structures. A comparison between all the capacity spectra using the LM3 approach can be seen in Figure 16d. This figure reflects the high scattering inherent to the exposure.  In LM2, structures with over six stories are grouped as RCH, which may be too general. In this respect, LM3 increases the resolution of the results in estimations of seismic damage. The cost of repairing a 14-story structure can be much higher than the cost of repairing a 7-story one. For this reason, it is important to develop and characterize new typologies that represent the exposure of an urban environment in a more realistic, exhaustive manner. These new typologies must be assessed considering the current state of the art regarding the modeling of structures and the characterization of seismic action. All these aspects have been considered in the development of LM3.
LM3 also improves the accuracy of the results as it considers the contribution of higher modes in the non-linear dynamic response of structures without simplifying the behavior of MDoF through In LM2, structures with over six stories are grouped as RCH, which may be too general. In this respect, LM3 increases the resolution of the results in estimations of seismic damage. The cost of repairing a 14-story structure can be much higher than the cost of repairing a 7-story one. For this reason, it is important to develop and characterize new typologies that represent the exposure of an urban environment in a more realistic, exhaustive manner. These new typologies must be assessed considering the current state of the art regarding the modeling of structures and the characterization of seismic action. All these aspects have been considered in the development of LM3.
LM3 also improves the accuracy of the results as it considers the contribution of higher modes in the non-linear dynamic response of structures without simplifying the behavior of MDoF through single degree of freedom systems. The method presented herein, LM3, can be extended to other structural typologies. Thus, structural typologies related to steel and masonry can be redefined and recalculated considering the aspects presented in this study. This effort should be made to obtain a better estimation of the expected seismic damage and risk at an urban level.
Regarding the reliability of the probabilities of the DS K obtained with LM3 and presented in this study, several simulations have been repeated. The average error between results of different simulations is approximately 2%. AvSa has been used to scale the earthquake records to the UHS. It would be of interest to review whether other intensity measures reduce this error. Pinzon et al. (2019) [30] recently analyzed the correlation between several IMs and the MIDR. They found that there are IMs which tend to be strongly correlated with the MIDR than, for instance, the more traditional IM spectral acceleration at the first-mode period of the structure, Sa(T 1 ). Moreover, the results are sensitive to the definition of damage limits presented in Table 2. It will be of interest to analyze the damage observed in real structures after an earthquake, to calibrate and improve these limits.
Currently, the simulation scheme allows for obtaining the probabilities of DS K by considering several structural typologies and earthquake records, the latter conditioned to a single UHS. However, parallelization tools could be used in such a way that several return periods can be assessed simultaneously.
In the current simulation scheme, participation of non-structural elements was not included, since 2D models are being simulated. Hence, the infills that affect the dynamic behavior of the structure are mainly those located at the facades. In 2D modeling approaches, one should look for a frame that characterizes the structure under study in the best possible way. Thus, considering the infills in this singular frame could overestimate the contribution to the global stiffness of these elements. Moreover, as this contribution to global stiffness disappears at low deformation values of the structure and generally makes the structure less vulnerable, the calculation was performed without considering the contribution of these elements. However, further research will be conducted to perform the calculation presented herein using 3D models. This approach will take into consideration not only the participation of non-structural elements but also torsion, which is a key issue in assessments of a structure's seismic behavior.
Finally, the results presented here will be open access and can be downloaded from a repository administrated by the Polytechnic University of Catalonia (UPC), the host institution of the KaIROS project. Link to the repository: http://kairoseq.upc.edu.