Physics-Based Ground Motion Simulations for the Prediction of the Seismic Vulnerability of Masonry Building Compounds in Mirandola (Italy)

: The current paper aims at investigating the seismic capacity of a masonry building aggregate in the historical centre of Mirandola based on a reliable ground motion simulation procedure. The examined clustered building is composed of eleven structural units (SUs) mutually interconnected to each other, which are made of brick walls and are characterized by wooden ﬂoors poorly connected to the vertical structures. Non-linear static analyses are performed by adopting the 3Muri software to characterize the seismic capacity of both the entire aggregate and the individual SUs. In this framework, a multi-scenario physics-based approach is considered for the deﬁnition of the seismic input in terms of broadband seismic signals inclusive of source and site effects. Finally, the incidence of the seismic input variability is discussed for the prediction of the global capacity response of the case study building.


Introduction
The repeated occurrence of devastating seismic events has raised the technical-scientific sensitivity towards the issue of the protection of human lives and exposed buildings [1,2]. Generally, existing buildings represent the weakest vulnerable part of urban systems since they are susceptible to damage under seismic phenomena. Throughout the various construction eras, they have preserved techniques and structural details typical of a systematic design procedure conceived to resist, predominantly, gravitational actions only [3,4]. Therefore, existing unreinforced masonry (URM) buildings intrinsically present potential vulnerability factors due to multiple aspects related to the nature of the basic material, geometric complexity, construction heterogeneity and absence of appropriate seismic protection devices, that alter the entire architectural, functional and structural system [5,6].
In general, it was observed that several historical centres are characterized by an ancient building heritage having an inadequate safety level against earthquakes. In the urban centres, the existing construction typologies are frequently grouped in aggregate configurations, and the constructive matrix of structural units (SUs) is the living representation of a well-consolidated design practice that is very vulnerable to seismic actions [7][8][9][10].
Generally, the seismic capacity offered by single SUs is considerably different from that of the entire building aggregate. So, accounting for the highly non-linear response of

The Historical Centre of Mirandola
The historic compartment of Mirandola is one of the oldest urban centres on the Italian peninsula. It developed over the centuries based on a spontaneous urban planning process that gave it the first "Forma Urbis", apparently casual, but responding to accurate settlement details prepared for the creation of roads and public spaces necessary for civil and public life, where important bureaucratic activities were carried out for the entire community. The historic centre of Mirandola is a tangible representation of the spatial configuration of a well-defined urban project starting from the early Middle Ages. The city was built as a fortified city with a well-structured and organized conformation of the historical settlements. The urban conformation of the municipality started from the construction of buildings of public interest (e.g., churches and palaces around the Castle of Pico), conceived with new architectural and functional elements to create a global planimetric structure more meticulous and responsive to the concept of the city. The city was connected to the neighbouring villages while maintaining the conformation of the main roads. In general, the buildings were arranged in a grouped configuration and bordered on the main façade, thus presenting a free front for external view and access, and an opposite front overlooking the internal pertinence area (Figure 1) [30].  . On the left, a bird-eye view of the city, in which the highlighted areas indicate the urban expansion process, is illustrated. On the right, the urban development started from the "ancient village" (XIII Century) to the complete configuration of the citadel (XV Century).
Nowadays, the union of the villages and the organization of new districts have kept the original planimetric configuration of the unchanged octagonal city (Figure 2).
The historic centre is mainly characterized by masonry buildings, erected in a clustered conformation, in which the coexistence of construction heterogeneities, sometimes incongruous, has led to the manifestation of material degradation. The historic centre is mainly characterized by masonry buildings, erected in a clustered conformation, in which the coexistence of construction heterogeneities, sometimes incongruous, has led to the manifestation of material degradation.

Main Features of the Case Study Building Compound
The historic centre of Mirandola is made up of 43 urban areas composed of buildings erected and conceived in the aggregate configuration.
The configuration of the current building stock highlights how the size of the builtup area is correlated to the evolution of the city foreseen in the Recovery Plan approved in July 2001. However, although Mirandola did not have a radial expansion, it is noted that the dimensional variation of the building aggregates is linked to the transition from the first expansion phase of the city with a quadrangular plan to the second one with an octagonal urbanized distribution [30,31]. In particular, the buildings located in the historic centre are mainly characterized by a construction uniformity of vertical structures, made

Main Features of the Case Study Building Compound
The historic centre of Mirandola is made up of 43 urban areas composed of buildings erected and conceived in the aggregate configuration.
The configuration of the current building stock highlights how the size of the built-up area is correlated to the evolution of the city foreseen in the Recovery Plan approved in July 2001. However, although Mirandola did not have a radial expansion, it is noted that the dimensional variation of the building aggregates is linked to the transition from the first expansion phase of the city with a quadrangular plan to the second one with an octagonal urbanized distribution [30,31]. In particular, the buildings located in the historic centre are mainly characterized by a construction uniformity of vertical structures, made of solid brick walls, with the presence of semi-rigid and deformable floors, which represents the most widespread construction typology of the historic centres of the Emilia-Romagna Region [10].
Concerning the case study aggregate, it is an existing complex masonry building with a non-regular shape composed of 11 SUs (Figure 3a). The clustered units forming the aggregate occupy different structural positions. In particular, as reported in Figure 3b, two SUs (called SU1 and SU11) occupy corner positions, nine cells (from SU2 to SU10) are placed in intermediate positions, and only SU7 occupies the heading position.
Buildings 2021, 11, x FOR PEER REVIEW 5 of 25 of solid brick walls, with the presence of semi-rigid and deformable floors, which represents the most widespread construction typology of the historic centres of the Emilia-Romagna Region [10]. Concerning the case study aggregate, it is an existing complex masonry building with a non-regular shape composed of 11 SUs (Figure 3a). The clustered units forming the aggregate occupy different structural positions. In particular, as reported in Figure 3b, two SUs (called SU1 and SU11) occupy corner positions, nine cells (from SU2 to SU10) are placed in intermediate positions, and only SU7 occupies the heading position.    The brick masonry walls have an average thickness ranging from 0.24 m to 0.50 m. The horizontal structures are mainly made up of wooden elements with upper double planking. For these elements, a design vertical load, F d , of 5.5 KNm −2 has been estimated, considering the fundamental load combination at the ultimate limit state (ULS) [30].
Regarding the physical conditions of the building aggregate, it is noted that SUs have common walls simply juxtaposed to each other. In addition, the pushing effects of the roofing systems could trigger the overturning of the main façades. This phenomenon is very common in historical buildings since, in the case of earthquakes, the masonry walls, poorly connected to the structure top level, are subjected to horizontal pushing actions deriving from roof, which usually activate out-of-plane mechanisms.
Concerning the material degradation, along the facades facing the main streets, even if partial detachments of the plaster are evident, they do not compromise the functional statics of the building compound.
The mechanical properties of the structural elements characterizing the case study building are deduced from the indications prescribed in Table C8.5.I of the Italian NTC18 standard [32]. Due to the absence of accurate on-site test procedures, the mechanical characteristics of the masonry are reduced by assuming a confidence factor, FC, equal to 1.35, corresponding to a LC1 (limited) knowledge level.
Specifically, for the mechanical parameters (f m , τ 0 , f v0 ), the minimum values of the intervals are reported in Table C8.5.I of the Italian standard are considered [32], while for the elastic moduli (E and G), the mean values are used. The mechanical parameters of solid brick masonry walls are summarized in Table 1.

Seismological Structural Model
To take into account the entire variability of structural behaviour from elastic to inelastic, up to global collapse [33,34], as well as the record-to-record variability (variability related to the mechanism of the seismic source, path attenuation effects and local site), ground motion records must be properly selected and, if necessary, scaled to a certain seismic intensity level to provide relatively stable estimates of the median values of the damage thresholds for the investigated structure.
For the seismic input definition, with the physics-based approach followed here, the modelling of the propagation of seismic waves for the estimation of the seismic hazard requires a prior characterization of the mechanical properties of the soil layers crossed by the waves generated by a seismogenic source. In this regard, the XeRIS web application, as proposed in [35,36], is used to compute the seismic wavefield responsible for structural damage. The modelling of the ground motion scenarios is based on the modal summation (MS) technique, as described in [37], which considers the vibrational modes of the layered medium excited by a seismic source, providing a rapid computational procedure for an adequate simulation of the ground motion in the case of a far-field event. Therefore, in the case of near-field events, the discrete wave number (DWN) technique is used according to the theory proposed in [38]. Nevertheless, these techniques allow for a very efficient calculation of synthetic seismograms in laterally homogeneous layered models, taking into account anelastic attenuation. Thus, the proposed site-specific geological model has been defined starting from the soil structure provided by the Italian Accelerometric Archive, ITACA, [39] considering a soil class C, then adopting the mechanical layers parameters (thickness, density, P and S waves velocity, and attenuation) discussed in [40], which represent the geological setting of the investigated area.

Fault Scenarios
The calculation of a fault-based ground-shaking scenario allows modelling, at a specific site, the ground motion caused by an earthquake with an extended fault. The computational technique used is based on the theory proposed by [41], where the seismogenic sources are modelled as a set of sub-sources, with a specific scaling source spectra law, also to predict the effects related to the kinematic rupture process (i.e., directivity).
Moreover, some of the kinematic parameters are made by varying the nucleation point, rupture velocity and slip field, and also take into consideration the source stochastic component [42,43].
Four different sources (Mirandola (MIR), Finale Emilia (FE), Veronese (VR) and Ferrara (FE)) are selected from the Web-GIS application DISS (Database of Individual Seismogenic Sources) [44] to compute the seismic input at the selected site. The sources have rectangular geometry of length, L, and width, W, with a variable distance between 1.9 km up to 44 km about the reference site, identified as MR (Lon ( • ): 11.06 and Lat ( • ): 44.88-geographical position of the study aggregate). In particular, Mirandola's source descends solely on the evidence of the recent tectonic activity of the Arc of Ferrara [45]. The 29 May 2012 earthquake in the Emilia-Romagna region of Italy, which was the second main shock of the Emilian sequence, activated this seismogenic source, producing a detectable uplift of the buried Mirandola anticline. The final geometry of the fault adopted is the one validated in the study proposed by [40]. The Finale Emilia seismic source was formed after the occurrence of the 20th May 2012 earthquake, which was the first mainshock of the Emilian sequence giving rise to a detectable uplift of the associated buried anticline. The Veronese source was generated following the 3rd January 1117 Verona earthquake, which was perhaps the strongest event that occurred in the Po Valley. The data collected in both historical and instrumental catalogues have shown that the seismicity of the Po Valley is concentrated mainly along the foothills and the buried thrust fronts of the northern Apennines and the southern Alps [46]. Finally, concerning the existence of the Ferrara's source, it dates back to the earthquake that occurred on 17th November 1570 based on data relating to both the recent tectonic activity of the Ferrara Arch [45], and from the subsoil geological structure [47]. Thus, based on the above considerations, some of the characteristics of the four considered seismic sources are reported in Table 2. Table 2. Some characteristics of the four selected seismogenic sources [40,44]. The acronyms are related to: edi = epicentral distance; depth = focal depth of the earthquake occurred; strike = defines the orientation of the fault counter-clockwise till the North; dip = represents the inclination of the fault plane; rake = the direction the hanging wall moves during rupture measured relative to the fault strike; sre = represents the angle formed between the focal point and the site. Furthermore, a map of the location of the case study building (MR) to source configuration has been presented in Figure 4   Subsequently, starting from the definition of the above-introduced seismic sources, the modelled seismic scenarios consider the moment magnitude, the position of the nucleation point, and the mean rapture velocity, Vs. The simulation proposed is based on the Monte-Carlo technique for the quantification of the distribution of the final slip, considering 100 statistical realizations and assuming the random variability of the nucleation point of coordinates (X, Y) (the maximum allowable variation is ±1/3 of the fault length (X direction) and ±1/3 of the fault width (Y direction)).

ID
The simulation of the source rupture processes has been performed using the PUL-SYN06 algorithm based on the theory presented in [41,48].
For a sake of example, the fault rupture scenarios generated by four realizations (i.e., #1, #50, #90, #100) out of 100 for the Mirandola (MIR) source, are shown in Figure 5. For each simulated rupture scenario, a 2D final slip map function is derived. The sub-source grid is made up of 17 × 9 crosses in which the red dot represents the nucleation point, the dark blue indicates the high slip of the fault, and the white contours (isochrones) represent the time-domain slip rupture propagation of a non-stationary random process. Subsequently, starting from the definition of the above-introduced seismic sources, the modelled seismic scenarios consider the moment magnitude, the position of the nucleation point, and the mean rapture velocity, V s . The simulation proposed is based on the Monte-Carlo technique for the quantification of the distribution of the final slip, considering 100 statistical realizations and assuming the random variability of the nucleation point of coordinates (X, Y) (the maximum allowable variation is ±1/3 of the fault length (X direction) and ±1/3 of the fault width (Y direction)).
The simulation of the source rupture processes has been performed using the PUL-SYN06 algorithm based on the theory presented in [41,48].
For a sake of example, the fault rupture scenarios generated by four realizations (i.e., #1, #50, #90, #100) out of 100 for the Mirandola (MIR) source, are shown in Figure 5. For each simulated rupture scenario, a 2D final slip map function is derived. The sub-source grid is made up of 17 × 9 crosses in which the red dot represents the nucleation point, the dark blue indicates the high slip of the fault, and the white contours (isochrones) represent the time-domain slip rupture propagation of a non-stationary random process. Buildings 2021, 11, x FOR PEER REVIEW 9 of 25

Maximum Credible Seismic Input (MCSI)
For the assessment of the seismic vulnerability of the considered case study, an appropriate calculation procedure is presented to predict the expected hazard for the site of interest.
For earthquake engineering purposes, the maximum credible seismic input (MCSI) represents a reliable estimation of the expected ground shaking level for a specific site, independently of the occurrence rate of earthquakes that could affect the investigated area. In particular, the proposed methodology [37,43] does not use empirical equations, such as GMPE, to derive the intensity measure (e.g., PGA or SA, PGV, PGD engineering parameters), but takes into consideration the seismic scenarios generated by independent seismic sources providing realistic ground motions in the time domain (either as a response spectrum or as a set of seismograms). Nevertheless, the use of source spectra [41] introduces a stochastic element in the proposed methodology; therefore, to define the MCSI spectrum, the uncertainties correlated to the hazard estimation must be considered.
Physics-based synthetic engineering parameters can be calculated through knowledge of the earthquake generation process and the propagation of seismic waves in an inelastic soil profile.
In this regard, the structural and topographical heterogeneities, as well as the influence of the source rupture process (slip distribution and rupturing velocity of the fault), were evaluated using the XeRIS web application, as proposed in [35].
Regarding the physical definition of the MCSI response spectrum, it is described in [29,42] where, for each seismogenic source, n-scenarios (in terms of magnitude, epicentre distance and focal mechanism of the earthquake) have been considered, the obtained nspectral acceleration values (SA) have been compared, and the maximum median has been selected, as depicted in Figure 6.

Maximum Credible Seismic Input (MCSI)
For the assessment of the seismic vulnerability of the considered case study, an appropriate calculation procedure is presented to predict the expected hazard for the site of interest.
For earthquake engineering purposes, the maximum credible seismic input (MCSI) represents a reliable estimation of the expected ground shaking level for a specific site, independently of the occurrence rate of earthquakes that could affect the investigated area. In particular, the proposed methodology [37,43] does not use empirical equations, such as GMPE, to derive the intensity measure (e.g., PGA or SA, PGV, PGD engineering parameters), but takes into consideration the seismic scenarios generated by independent seismic sources providing realistic ground motions in the time domain (either as a response spectrum or as a set of seismograms). Nevertheless, the use of source spectra [41] introduces a stochastic element in the proposed methodology; therefore, to define the MCSI spectrum, the uncertainties correlated to the hazard estimation must be considered.
Physics-based synthetic engineering parameters can be calculated through knowledge of the earthquake generation process and the propagation of seismic waves in an inelastic soil profile.
In this regard, the structural and topographical heterogeneities, as well as the influence of the source rupture process (slip distribution and rupturing velocity of the fault), were evaluated using the XeRIS web application, as proposed in [35].
Regarding the physical definition of the MCSI response spectrum, it is described in [29,42] where, for each seismogenic source, n-scenarios (in terms of magnitude, epicentre distance and focal mechanism of the earthquake) have been considered, the obtained n-spectral acceleration values (SA) have been compared, and the maximum median has been selected, as depicted in Figure 6. uildings 2021, 11, x FOR PEER REVIEW 10 of 25 Figure 6. Simulated maximum credible seismic input (MCSI) for the selected sources. The grey band identifies the MCSI, as defined in [43], which is controlled by three sources: curves 1, 2 and 3 are the median spectra obtained from one hundred realizations of the rupture process for the sources MIR, FE and VR, respectively.

Global Behaviour of the Case Study Aggregate
The case study compound has been numerically analysed using the 3Muri software founded on the Frame Macro-Elements (FME) theory ( Figure 7) [49]. This methodological approach assumes that masonry walls are considered as a set of one-dimensional macroelements (columns, beams, and nodes) mutually interconnected. The resistance criteria of the deformable elements have been assumed based on the requirements of EN 1998-3 [50], which establish maximum threshold values of 0.4% and 0.6% for shear failure and flexural collapse, respectively. The analyses were carried out along the two main global directions of the building, X and Y. The above-introduced simulations were interrupted when the shear strength decay was equal to 20% of the maximum shear capacity, as reported in the NTC18 standard [32].
To this purpose, two different distributions of seismic forces and four combinations (±X, ±Y) have been considered, neglecting the influence of accidental eccentricities in the calculation procedure. The results in terms of capacity curves are presented in Figure 8.  [43], which is controlled by three sources: curves 1, 2 and 3 are the median spectra obtained from one hundred realizations of the rupture process for the sources MIR, FE and VR, respectively.

Global Behaviour of the Case Study Aggregate
The case study compound has been numerically analysed using the 3Muri software founded on the Frame Macro-Elements (FME) theory ( Figure 7) [49]. This methodological approach assumes that masonry walls are considered as a set of one-dimensional macroelements (columns, beams, and nodes) mutually interconnected. The resistance criteria of the deformable elements have been assumed based on the requirements of EN 1998-3 [50], which establish maximum threshold values of 0.4% and 0.6% for shear failure and flexural collapse, respectively. The analyses were carried out along the two main global directions of the building, X and Y.
Buildings 2021, 11, x FOR PEER REVIEW 10 of 25 Figure 6. Simulated maximum credible seismic input (MCSI) for the selected sources. The grey band identifies the MCSI, as defined in [43], which is controlled by three sources: curves 1, 2 and 3 are the median spectra obtained from one hundred realizations of the rupture process for the sources MIR, FE and VR, respectively.

Global Behaviour of the Case Study Aggregate
The case study compound has been numerically analysed using the 3Muri software founded on the Frame Macro-Elements (FME) theory ( Figure 7) [49]. This methodological approach assumes that masonry walls are considered as a set of one-dimensional macroelements (columns, beams, and nodes) mutually interconnected. The resistance criteria of the deformable elements have been assumed based on the requirements of EN 1998-3 [50], which establish maximum threshold values of 0.4% and 0.6% for shear failure and flexural collapse, respectively. The analyses were carried out along the two main global directions of the building, X and Y. The above-introduced simulations were interrupted when the shear strength decay was equal to 20% of the maximum shear capacity, as reported in the NTC18 standard [32].
To this purpose, two different distributions of seismic forces and four combinations (±X, ±Y) have been considered, neglecting the influence of accidental eccentricities in the calculation procedure. The results in terms of capacity curves are presented in Figure 8. The above-introduced simulations were interrupted when the shear strength decay was equal to 20% of the maximum shear capacity, as reported in the NTC18 standard [32].
To this purpose, two different distributions of seismic forces and four combinations (±X, ±Y) have been considered, neglecting the influence of accidental eccentricities in the calculation procedure. The results in terms of capacity curves are presented in Figure 8. From the gotten results, it was observed how the aggregate provides a better structural capacity in the transverse direction Y, rather than in the orthogonal one X. In particular, the ultimate displacements achieved in the positive Y direction (Y + ) are approximately two times higher than the corresponding ones detected in the perpendicular Xdirection referred to the uniform distribution of the forces.
Conversely, in the negative X direction (X − ), an increase of the displacements was observed, compared to the Y direction.
In particular, the static condition (X − 1st Mode) has shown an ultimate displacement, du, equal to −5.18 cm, which is approximately 61% greater than the analogous displacement attained in the direction Y (du = −3.18 cm). Regarding the maximum base shear threshold, it was noted that in the Y direction, the maximum strength is approximately two times greater than the corresponding value in direction X. This discrepancy in terms of global behaviour (e.g., displacements and maximum shear thresholds) in the two analyses directions is due to the irregular planimetric conformation of the entire building conjunctly with a higher presence of masonry resistant area.
Finally, for a global view of the vulnerability threshold and expected damage concerning the examined building, two indices were analysed: the vulnerability index (VI) and the damage index (DI). With this intention, the vulnerability index, VI, is intended as the ratio between the acceleration capacity, PGACi, intended as the peak acceleration that determines the achievement of the ultimate limit state (ULS) condition based on the return periods, Tr, and the corresponding acceleration demand, PGADi, related to the spectral coordinates of the site where the construction is located. In the specific case, the 3 Muri software [49] automatically estimates the maximum spectral acceleration, starting from the definition of the elastic response spectrum (Tr = 475 year) for the site of interest (Mirandola) having assumed a soil class C, as reported in the Italian Accelerometric Archive, ITACA [39]. Consequently, the damage index has been evaluated, taking into account the mathematical formulation proposed in [51]: where µmax and µd are the required and available ductility associated with the structural system, respectively, which can be expressed in terms of displacements; dmax represents the required displacement, evaluated for Ti < TC, based on the area equivalence criterion; and dy and du are the yielding and ultimate displacements associated with the SDoF system, respectively. The numerator of the above-mentioned formulation indicates how the maximum displacement, dmax, deviates from the elastic one, while the denominator provides indications From the gotten results, it was observed how the aggregate provides a better structural capacity in the transverse direction Y, rather than in the orthogonal one X. In particular, the ultimate displacements achieved in the positive Y direction (Y + ) are approximately two times higher than the corresponding ones detected in the perpendicular X-direction referred to the uniform distribution of the forces.
Conversely, in the negative X direction (X − ), an increase of the displacements was observed, compared to the Y direction.
In particular, the static condition (X − 1st Mode) has shown an ultimate displacement, d u , equal to −5.18 cm, which is approximately 61% greater than the analogous displacement attained in the direction Y (d u = −3.18 cm). Regarding the maximum base shear threshold, it was noted that in the Y direction, the maximum strength is approximately two times greater than the corresponding value in direction X. This discrepancy in terms of global behaviour (e.g., displacements and maximum shear thresholds) in the two analyses directions is due to the irregular planimetric conformation of the entire building conjunctly with a higher presence of masonry resistant area.
Finally, for a global view of the vulnerability threshold and expected damage concerning the examined building, two indices were analysed: the vulnerability index (V I ) and the damage index (DI). With this intention, the vulnerability index, V I , is intended as the ratio between the acceleration capacity, PGA Ci , intended as the peak acceleration that determines the achievement of the ultimate limit state (ULS) condition based on the return periods, Tr, and the corresponding acceleration demand, PGA Di , related to the spectral coordinates of the site where the construction is located. In the specific case, the 3 Muri software [49] automatically estimates the maximum spectral acceleration, starting from the definition of the elastic response spectrum (Tr = 475 year) for the site of interest (Mirandola) having assumed a soil class C, as reported in the Italian Accelerometric Archive, ITACA [39]. Consequently, the damage index has been evaluated, taking into account the mathematical formulation proposed in [51]: where µ max and µ d are the required and available ductility associated with the structural system, respectively, which can be expressed in terms of displacements; d max represents the required displacement, evaluated for T i < T C , based on the area equivalence criterion; and d y and d u are the yielding and ultimate displacements associated with the SDoF system, respectively. The numerator of the above-mentioned formulation indicates how the maximum displacement, d max , deviates from the elastic one, while the denominator provides indications about the type of failure mechanisms (e.g., when the ultimate displacement, d u , is comparable to d y , the fragile mechanism is activated, vice versa, ductile ones). Thus, congruently to what has been introduced, the synthesis of the results has been depicted in Figure 9.
Buildings 2021, 11, x FOR PEER REVIEW 12 of 25 about the type of failure mechanisms (e.g., when the ultimate displacement, du, is comparable to dy, the fragile mechanism is activated, vice versa, ductile ones). Thus, congruently to what has been introduced, the synthesis of the results has been depicted in Figure 9.
(a) (b) First, the achieved results presented in Figure 9a have shown that the worst-case scenario is attained for the direction X + (uniform), which corresponds to a C/D ratio equal to 0.67. Similarly, in the direction Y (Y + 1st Mode), the expected vulnerability level is 0.40 (e.g., medium-low level).
Concerning the damage index DI, as reported in Figure 9b, it has been observed that the value of 0.57 is reached for a uniform distribution of forces in the X + direction, since the difference between the maximum displacement demand, dmax, and the corresponding yielding displacement, dy, is influenced by the achievement of a premature plastic displacement threshold.
Subsequently, the vulnerability indexes have been associated with the estimated damage index, DI, and then correlated to the damage thresholds provided by the European Macroseismic Scale, EMS-98 [52,53], as reported in Table 3.  First, the achieved results presented in Figure 9a have shown that the worst-case scenario is attained for the direction X + (uniform), which corresponds to a C/D ratio equal to 0.67. Similarly, in the direction Y (Y + 1st Mode), the expected vulnerability level is 0.40 (e.g., medium-low level).
Concerning the damage index DI, as reported in Figure 9b, it has been observed that the value of 0.57 is reached for a uniform distribution of forces in the X + direction, since the difference between the maximum displacement demand, d max , and the corresponding yielding displacement, d y , is influenced by the achievement of a premature plastic displacement threshold.
Subsequently, the vulnerability indexes have been associated with the estimated damage index, DI, and then correlated to the damage thresholds provided by the European Macroseismic Scale, EMS-98 [52,53], as reported in Table 3. [-] In Figure 10, the validation between the real damages detected after the Emilia-Romagna earthquake and the corresponding ones derived from the mechanical analysis procedure evaluated at the ultimate displacement, du, referring to the ULS (ultimate limit state) combination, has been provided. In Figure 10, the validation between the real damages detected after the Emilia-Romagna earthquake and the corresponding ones derived from the mechanical analysis procedure evaluated at the ultimate displacement, du, referring to the ULS (ultimate limit state) combination, has been provided. More specifically, there is a propensity to damage identified mainly in the spandrel elements, since they have no tensile strength (instead guaranteed by a steel tie-rod or reinforced concrete ring beam) towards the combined effects induced by both bending and shear.
Congruently to what was discussed in the above-introduced comparison, it is possible to appreciate how the simulated structural model provides a global damage level not dissimilar to that really detected in the case study aggregate. More specifically, there is a propensity to damage identified mainly in the spandrel elements, since they have no tensile strength (instead guaranteed by a steel tie-rod or reinforced concrete ring beam) towards the combined effects induced by both bending and shear.
Congruently to what was discussed in the above-introduced comparison, it is possible to appreciate how the simulated structural model provides a global damage level not dissimilar to that really detected in the case study aggregate.

Capacity Response of Building Structural Units
This section analyses the seismic behaviour of the individual structural units in the aggregate configuration using a numerical procedure based on the capacity curve reconstruction process derived from that of the whole aggregate. This methodological approach allows estimation of the propensity to damage of the SUs, considering their interaction in the aggregate configuration and, therefore, accounting for their contribution in terms of both stiffness and mass. Thus, after assessing the seismic response of the whole aggregate along with the main analysis directions, the behaviours of the reference cells included in the aggregate have been evaluated. In particular, the clustered SUs response has been achieved by adopting a step-by-step simulation process by considering as base shear, V b i , the analytical sum of the base reactions of the piers of the reference case study unit and, as displacement, the average value of the displacements of the selected monitoring nodes, D m i , at the top floor. The analytical formulations, adopted for the definition of the MDoF capacity curves, are presented according to the study proposed in [5,54,55]: having assumed |V j | as the absolute value of the base shear force associated with the j-th node at the i-th step, D j as the absolute value of the top displacement associated with the j-th node at the i-th step, and N as the total number of the nodes at the base of the masonry piers. The resulting curves in both analysis directions, X and Y, after a bi-linearization process, have been plotted in Figure 11. both stiffness and mass. Thus, after assessing the seismic response of the whole aggregate along with the main analysis directions, the behaviours of the reference cells included in the aggregate have been evaluated. In particular, the clustered SUs response has been achieved by adopting a step-by-step simulation process by considering as base shear, Vb i , the analytical sum of the base reactions of the piers of the reference case study unit and, as displacement, the average value of the displacements of the selected monitoring nodes, Dm i , at the top floor. The analytical formulations, adopted for the definition of the MDoF capacity curves, are presented according to the study proposed in [5,54,55]: having assumed |Vj| as the absolute value of the base shear force associated with the j-th node at the i-th step, Dj as the absolute value of the top displacement associated with the j-th node at the i-th step, and N as the total number of the nodes at the base of the masonry piers.
The resulting curves in both analysis directions, X and Y, after a bi-linearization process, have been plotted in Figure 11. From the acquired results, it appears that the behaviour of the SUs in the two analysis directions is quite heterogeneous. In particular, the head units have been denoted with black-grey, the corner units in red, and a black scale colour has been adopted for the other intermediate units. The capacity parameters expressed in terms of forces and displacements for each of the analysed SUs are presented in Table 4. From the acquired results, it appears that the behaviour of the SUs in the two analysis directions is quite heterogeneous. In particular, the head units have been denoted with black-grey, the corner units in red, and a black scale colour has been adopted for the other intermediate units. The capacity parameters expressed in terms of forces and displacements for each of the analysed SUs are presented in Table 4.
From the obtained results, it has been observed that, in the X direction, the head unit (SU7) has higher stiffness and strength than the other building units which occupy intermediate and corner positions, respectively. This circumstance has been revealed since the corner and intermediate units are affected by the beneficial effect of mutual structural confinement, while the head unit is constrained only on one side, thus resulting in greater susceptibility to displacements. In terms of displacements, it was noted that the SU7 (head) provides displacements comparable to the SU5 (intermediate), which has a higher mass due to the presence of contiguous SUs.
Moreover, in the Y direction, it has been observed that the SU7 provides a higher contribution in terms of displacement capacity than the other cells due to the interconnected confinement effect provided by the adjacent SUs. Furthermore, by comparing the maximum shear threshold of the SU7 in the two directions, it can be seen that, in the Y direction, the structural cell has an increased shear strength of about 44% compared to the corresponding threshold achieved in the orthogonal direction.

Analysis of the Capacity Parameters
To interpret the different global responses of the examined structural units, the main capacity parameters (stiffness, ductility, shear strength, mass and vibration period) have been plotted and correlated to each other, as reported in Figure 12. From the obtained results, it has been observed that, in the X direction, the head u (SU7) has higher stiffness and strength than the other building units which occupy int mediate and corner positions, respectively. This circumstance has been revealed since the corner and intermediate units are fected by the beneficial effect of mutual structural confinement, while the head unit constrained only on one side, thus resulting in greater susceptibility to displacements. terms of displacements, it was noted that the SU7 (head) provides displacements comp rable to the SU5 (intermediate), which has a higher mass due to the presence of contiguo SUs.
Moreover, in the Y direction, it has been observed that the SU7 provides a high contribution in terms of displacement capacity than the other cells due to the interco nected confinement effect provided by the adjacent SUs. Furthermore, by comparing t maximum shear threshold of the SU7 in the two directions, it can be seen that, in the direction, the structural cell has an increased shear strength of about 44% compared to t corresponding threshold achieved in the orthogonal direction.

Analysis of the Capacity Parameters
To interpret the different global responses of the examined structural units, the ma capacity parameters (stiffness, ductility, shear strength, mass and vibration period) ha been plotted and correlated to each other, as reported in Figure 12. First of all, the vibration period of the single structural units analysed, taking i account their contribution in terms of mass and stiffness, has been appropriately co pared with the simplified formulation derived from the NTC18 standard [32] (see Fig  12a). From the comparison, it has been noted that in the two analysis directions, the p First of all, the vibration period of the single structural units analysed, taking into account their contribution in terms of mass and stiffness, has been appropriately compared with the simplified formulation derived from the NTC18 standard [32] (see Figure 12a). From the comparison, it has been noted that in the two analysis directions, the periods associated with the single SUs are heterogeneous with a higher incidence in the X direction, where they exceed the corresponding period provided by the NTC18 standard [32].
However, it is worth pointing out that the NTC18 standard [32] provides a simplified formulation for estimating the vibration period of buildings considered as isolated structures, without considering the aggregate effect. Thus, it is noted that the SUs 1, 2, 3, 6b, 8, 9 and 10 have a low seismic mass-to-stiffness ratio, presenting a reduced vibration period compared to the other structural cells examined in the same analysis direction. Conversely, in Y direction, the vibration period is, on average, 24% lower compared to the analogue values examined in the orthogonal direction.
This circumstance denotes that, in the transverse direction, the stiffness associated with the single SUs is on average 79% higher than that obtained in the X direction, as depicted in Figure 12b. As far as the displacement ductility is concerned, it has been correlated to the wall resistant area (see Figure 12c) to consider the incidence of the structural resistant area (W R.A ) on the expected ductility (µ). From the results achieved, it was observed that in the X direction, the point cloud is quite homogeneous for 2.0 < µ < 7.0, which corresponds to a wall resistant area W R.A enclosed in the range between 2.60 m 2 and 21.0 m 2 . On the contrary, due to the higher number of wall panels, in the Y direction, the wall resistant area is within 7.40 m 2 and 32.0 m 2 and reduces the corresponding displacement.
Finally, in Figure 12d, the correlation between the maximum shear force and the total volume of SUs is shown. The analysis of the results in both analysis directions has revealed how the linear regression function presents a direct proportionality between the monitored parameters. The regression coefficient (R 2 ) is robust, assuming values greater than 0.5 in both examined cases, taking into account all the building-to-building variability in terms of mass, stiffness, and in-elevation interaction.

Vulnerability Assessment for the Bilinear Elasto-Plastic SUs Models
Following the above-introduced sections, the generation of MCSI spectra was implemented according to the soil structural model and fault scenario selected [56]. In particular, the MCSI has been defined in terms of the response spectrum to accomplish the engineering seismic verification of the case study masonry compound.
With this intent, according to [29,42], the MCSI response spectrum, taking into account the stochastic nature of the algorithm implemented by the XeRIS tool [35], is evaluated at the 95th percentile, based on the number of realizations of the fault rupture process considered. Therefore, by defining the maximum credible seismic input, the maximum contribution in terms of seismic hazard is offered by the Mirandola source (MIR, see Table 4) due to a reduced site-source distance. In these circumstances, the source spectrum derived from the rapture process of the Mirandola source represents the maximum expected MCSI for the generated scenario, enveloping, for periods of interest, the spectra generated by the other sources considered. The derived spectra for both the 50th and 95th percentile (median and maximum), are presented in Figure 13a, with their regularized shape according to the Italian seismic code as illustrated in Figure 13b. Finally, in Figure 13c, the design response spectra for the Mirandola site [32] have been presented, assuming soil class C and Tr = 475 years.
by the other sources considered. The derived spectra for both the 50th and 95th percentile (median and maximum), are presented in Figure 13a, with their regularized shape according to the Italian seismic code as illustrated in Figure 13b. Finally, in Figure 13c, the design response spectra for the Mirandola site [32] have been presented, assuming soil class C and Tr = 475 years. Consequently, the N2 method [57][58][59] has been adopted for the evaluation of the seismic performance of the single SUs in aggregate conditions by adopting the above-introduced MCSI spectrum. The method combines non-linear static analysis and the response spectrum approach to derive the ADRS (acceleration displacement response spectrum) domain, as reported in Figure 14. Moreover, the seismic verification is carried out by estimating the vulnerability index, appropriately evaluated as the ratio between the seismic demand, Di, and the relative capacity, Ci, provided by the studied structures. Consequently, the N2 method [57][58][59] has been adopted for the evaluation of the seismic performance of the single SUs in aggregate conditions by adopting the aboveintroduced MCSI spectrum. The method combines non-linear static analysis and the response spectrum approach to derive the ADRS (acceleration displacement response spectrum) domain, as reported in Figure 14. Moreover, the seismic verification is carried out by estimating the vulnerability index, appropriately evaluated as the ratio between the seismic demand, D i , and the relative capacity, C i , provided by the studied structures.
The results presented have shown how the expected vulnerability is heterogeneous in both analysis directions, depending on the structural position of SUs, as discussed in Sections 3.1 and 3.2. For the sake of example, a detailed application of the N2 method to the head unit is depicted in Figure 15, where also the ADRS obtained using the source spectra is compared with the corresponding spectrum derived from the NTC18 standard [32]. The results presented have shown how the expected vulnerability is heterogeneous in both analysis directions, depending on the structural position of SUs, as discussed in Sections 3.1 and 3.2. For the sake of example, a detailed application of the N2 method to the head unit is depicted in Figure 15, where also the ADRS obtained using the source spectra is compared with the corresponding spectrum derived from the NTC18 standard [32]. From the obtained results, it is noted that, in the transverse direction, there is an increase of both strength (55%) and stiffness (40%) concerning the longitudinal direction response. Subsequently, the results obtained by applying the ADRS spectrum, derived from the NTC18 standard [32], have revealed how the verification in the X direction provides a D/C ratio equal to 0.83 while, considering MCSI at the 50th percentile, the mechanical vulnerability, D/C, is equal to 0.96, with an increase of 14%. In the Y direction, it is From the obtained results, it is noted that, in the transverse direction, there is an increase of both strength (55%) and stiffness (40%) concerning the longitudinal direction response. Subsequently, the results obtained by applying the ADRS spectrum, derived from the NTC18 standard [32], have revealed how the verification in the X direction provides a D/C ratio equal to 0.83 while, considering MCSI at the 50th percentile, the mechanical vulnerability, D/C, is equal to 0.96, with an increase of 14%. In the Y direction, it is noted as D/C| NTC18 = 0.7, while the corresponding value obtained by MCSI is equal to 0.9, with an increase of 22%. Moreover, the ADRS spectrum achieved from the MCSI evaluated at the 95th percentile provides unsafe conditions, since the seismic checks are not satisfied (D/C > 1) in both analysis directions.
In conclusion, the N2 method verifications implemented through the proposed approach to derive seismic spectra have provided results on the safe side, in comparison to those provided by the standard code.

Fragility Assessment
In this section, the mean fragility curves of the different SUs have been derived. The fragility curves define the probability that a given damage state, d Si , is equalled or exceeded by a specific intensity measure. The adopted methodology considers as intensity measure the spectral displacement, S d , according to the following Equation (3) [10,55,60]: where Φ [•] stands for the cumulative normal distribution, while the denominator and β are the median and the standard deviation of the corresponding normal distribution, respectively. The damage thresholds have been assumed according to the theory proposed by [3,60] as a function of the structure yielding and ultimate displacements. Moreover, to reduce the number of uncertainties deriving from the deterministic approach, a single value has been assumed for the standard deviation.
Specifically, the standard deviation has been set as equal to the displacement ductility, µ, multiplied for a fixed correlation coefficient equal to 0.45. Thus, the damage states and the standard deviations have been identified, as reported in Table 5. The expected damage levels are in line with the assumption proposed in the framework of the Risk-UE project [61], which is generally used in many seismic risk studies on European building stocks in prone cities. Therefore, the fragility functions have been derived, taking into consideration the average displacements of the above-mentioned case study SUs, as depicted in Figure 16.

0.45‧ln(µ)
The expected damage levels are in line with the assumption proposed in the framework of the Risk-UE project [61], which is generally used in many seismic risk studies on European building stocks in prone cities. Therefore, the fragility functions have been derived, taking into consideration the average displacements of the above-mentioned case study SUs, as depicted in Figure 16.  The obtained results clearly show how the corner structural unit provides damage probabilities higher than those of the other SUs occupying intermediate and head positions. In the specific case, the corner unit, despite having a number of floors less than the contiguous structural cells, has a reduced lateral constraint guaranteed by the adjacent structural units, so that the torsional phenomena limit the displacement capacity and, consequently, the expected damage is much more pronounced. Regarding the torsional effects, they are strictly correlated by the addition of structural mass rates offered by the contiguous SUs, resulting in a non-uniform distribution of the resistant elements in the plane.
Therefore, the more pronounced eccentricity between the centre of gravity and centre of stiffness triggered the premature failure of the masonry wall panels, confirming the damage level depicted in Section 4.1. In the specific case, the corner unit, despite having a number of floors less than the contiguous structural cells, has a reduced lateral constraint guaranteed by the adjacent structural units, so that the torsional phenomena limit the displacement capacity and, consequently, the expected damage is much more pronounced. Regarding the torsional effects, they are strictly correlated by the addition of structural mass rates offered by the contiguous SUs, resulting in a non-uniform distribution of the resistant elements in the plane. Therefore, the more pronounced eccentricity between the centre of gravity and centre of stiffness triggered the premature failure of the masonry wall panels, confirming the damage level depicted in Section 4.1.

Conclusions
The research work analysed the seismic vulnerability of an unreinforced masonry building aggregate located in the historic centre of Mirandola using a physics-based simulation approach.
The first step aims to characterize the expected hazard for the site of interest, based on the physics-based ground motion simulations. To this end, the web application XeRIS has been adopted. The starting point was the definition of the soil structure profile to model the propagation of the seismic waves for the estimation of the seismic hazard. In this regard, parametric studies were conducted both on the geological structural model and on the configurations of the selected seismic sources.
The proposed structural model contains the specific geological properties of the investigated area. Subsequently, four different sources (Mirandola (MIR), Finale Emilia (FE), Veronese (VR) and, Ferrara (FE)), were selected to derive the MCSI source spectrum by adopting a stochastic seismic scenario that considers all the possible variability related to the moment magnitude, the position of the nucleation point, and the mean slip rapture velocity.
The second step concerned the evaluation of the global seismic behaviour of the examined masonry compound through non-linear static analyses. The different global behaviour (e.g., displacements and maximum shear thresholds) in the two analysis directions is mainly due to the irregular planimetric configuration of the entire building conjunctly to the higher presence of resistant wall areas. Subsequently, concerning the global response, two synthetic parameters, namely the vulnerability index (V I ) and the damage index (DI), were analysed. To this end, the vulnerability index, V I , was mechanically estimated directly by extrapolating the results from the 3Muri software, which once defined the elastic response spectrum for the site of interest. From the results achieved, it was shown that the X + (uniform) direction provides a PGA C /PGA D ratio equal to 0.70 while, in the transverse direction (Y + 1st Mode), the vulnerability stands at 0.40 (e.g., medium-low level). As far as the damage index (DI) is concerned, it has been used for establishing a failure hierarchy based on the analyses performed. In particular, it has been found that the selected parameter varied in the range between 0.1 and 0.6, which mainly corresponds to a ductile (bending) crisis of the structural elements and the corresponding damage level according to the EMS-98 scale.
Consequently, the seismic behaviour of the individual structural units in the aggregate configuration was analysed using a numerical procedure based on the reconstruction process of the capacity curves derived from those of the entire aggregate. This methodological approach allows estimation of the damage propensity of the single SUs, considering their interaction in the aggregate, so taking into account the contribution given by SUs adjacent to the investigated one in terms of both stiffness and mass.
The performed analyses provided important insight into the correlation between shear force and ductility. It has been observed that, by correlating the base shear force and the total volume of the single SUs, there is a direct proportionality among the monitoring parameters. However, the ductility tends to decrease as the total resistant area of walls in Y direction increases up to 35 m 2 .
Consecutively, the seismic verification was implemented using the bilinear elastoplastic curves of SUs. To this purpose, the ADRS domain was adopted for the MCSI spectrum (derived at 50th and 9th percentile). In the first instance, the acquired results showed that the expected vulnerability is heterogeneous in the analysis directions, providing unsatisfactory results using the ADRS at the 95th percentile (worst-case scenario). Furthermore, to propose a suitable example regarding the performance of the SUs towards the seismic action, the investigated head unit was considered as a reference case for conducting the seismic verification in terms of demand and capacity displacements. In this regard, the ADRS obtained using the source spectra were compared with the corresponding spectrum derived according to the Italian design code. From the obtained results, it was noticed that, in the transverse direction, an increase of both strength (55%) and stiffness (40%) concerning the longitudinal direction was found.
The seismic verification in the longitudinal direction, obtained by applying the ADRS domain derived from the Italian NTC18 standard, provided a D/C ratio equal to 0.83, while, considering MCSI at the 50th percentile, the mechanical vulnerability, D/C, was equal to 0.96, with an increase of 14%. Moreover, in the Y direction, it was seen as D/C| NTC18 = 0.7, while the corresponding value obtained by MCSI was equal to 0.9, with an increase of 22%.
This aspect demonstrated how the displacement demands achieved by adopting the simulated scenario provide conforming results concerning the prediction of the code spectrum vulnerability procedure.
Finally, the mean fragility curves were derived for the different SUs investigated. The results showed how the corner structural unit provides higher damage probabilities than the other SUs occupying intermediate and head positions. Firstly, this outcome was due to the reduced lateral constraint guaranteed by the adjacent structural units, and, secondly, by the induced effect of the torsional phenomena strictly correlated to the addition of structural mass rates offered by the contiguous SUs, resulting in a non-uniform distribution of the resistant elements in the plane.
Overall, based on the results presented, this research topic provided complete information concerning the adoption of a specific procedure based on ground motion simulations for the seismic assessment of URM buildings.
For this reason, this study represents an initial first step that can provide clear insights into the causes of collapse observed in masonry structures after destructive earthquakes. In this framework, future developments should be oriented to the simulation of other seismic scenarios characterized by different fault mechanisms, epicentre distances, hypocentre depths, and ground amplification, to have a greater consistency on the effects caused by a possible earthquake on a larger building dataset.

Data Availability Statement:
Data deriving from the current study can be provided to the readers based upon their explicit request.