Seismic Fragility for a Masonry-Inﬁlled RC (MIRC) Building Subjected to Liquefaction

: Historical earthquakes have documented that lateral spread and settlements are the most signiﬁcant damages induced by soil liquefaction. Therefore, assessing its effects on structural performance has become a fundamental issue in seismic engineering. In this regard, the paper proposes to develop analytical fragility curves of a Masonry-Inﬁlled RC (MIRC) structure subjected to liquefaction-induced damages. In order to reproduce the nonlinear cyclic behavior (dilation tendency and the increase in cyclic shear strength) due to liquefaction, nonlinear hysteretic materials and advanced plasticity models were applied. The ﬁndings herein obtained in terms of seismic fragility of the MIRC building subjected to liquefaction may be implemented as guidelines or code provisions.


Background
Liquefaction-induced damages were shown to be the principal cause of disruption of functionality, economic losses, direct and indirect costs for communities during many earthquakes (Niigata, Japan, 1964; Dagupan City, Philippines, 1990; Chi-Chi, Taiwan, 1999; Japan, 2011; Kocaeli, Turkey, 1999; and Christchurch, New Zealand, 2011). During these earthquakes, considerable evidence of liquefaction was produced in terms of ground effects (such as flow slides, lateral spreading, ground oscillations, sand boils), as well as severe soil structure interaction (SSI) damages (such as sinking or tilting of heavy structures, failing of retaining systems, slumping of slopes and settlements of buildings). Therefore, the assessment of liquefaction effects has become the fundamental object of pre-and post-earthquake analyses. The study of such effects has been developed with many methodologies.
First of all, the effects of liquefaction were studied on empirical data in terms of settlements [1,2] and lateral spread [3,4]. Other contributions (e.g., [5,6]) proposed comprehensive methods based on relationships between selected soil parameters and liquefaction potential. However, as pointed out by [7], these approaches explored liquefaction susceptibility under free-filed conditions, without considering the damage to structures. Bullock et al. [8,9] applied semi-empirical methods to predict liquefaction consequences.
Secondly, some contributions proposed numerical simulations to assess the effects of liquefaction. For example, Refs. [10,11] investigated the building response to liquefaction. Numerical simulations were applied to study the effects of liquefaction to shallow-founded buildings by [12,13] and by [14] that concentrated on the effects of several earthquakes (Kocaeli 1999, Maule 2010, Christchurch 2011 and Tohoku 2011. Ref. [15] applied the Finite Different Method to study the bearing capacity degradation of shallow foundations on layered liquefiable soils. Other contributions performed numerical simulations to study countermeasures to reduce liquefaction-induced effects. For example, Ref. [16] considered soil compaction and/or increase of permeability below and around the applied surface load. Other methods were also investigated, such as stone columns that were the object of [17] and a preloading technique that was explored in [18]. Recently, Ref. [19] performed tridimensional (3D) numerical simulations to assess the effects of liquefaction

Case Study
In order to study the complex mechanism of liquefaction on a MIRC structural configuration, OpenSees PL [37,38] was applied to reproduce the system (soil + foundation + structure), as shown in Figures 1 and 2

Soil Model
The 3D mesh (118.4 × 124.4 m) was described with 31860 nodes and 35868 BrickUP, consisting of 8-nodes isoparametric elements that allow to consider both the displacements (longitudinal, transversal and vertical: degree of freedom: 1, 2 and 3) and the pore pressure (degree of freedom: 4) to simulate the dynamic response of a solid-fluid fully coupled material (Biot theory, [38]). The dimensions of the elements inside the mesh followed the previous contributions [19,20] and the performance of the lateral boundaries

Soil Model
The 3D mesh (118.4 × 124.4 m) was described with 31860 nodes and 35868 BrickUP, consisting of 8-nodes isoparametric elements that allow to consider both the displacements (longitudinal, transversal and vertical: degree of freedom: 1, 2 and 3) and the pore pressure (degree of freedom: 4) to simulate the dynamic response of a solid-fluid fully coupled material (Biot theory, [38]). The dimensions of the elements inside the mesh followed the previous contributions [19,20] and the performance of the lateral boundaries

Soil Model
The 3D mesh (118.4 × 124.4 m) was described with 31860 nodes and 35868 BrickUP, consisting of 8-nodes isoparametric elements that allow to consider both the displacements (longitudinal, transversal and vertical: degree of freedom: 1, 2 and 3) and the pore pressure (degree of freedom: 4) to simulate the dynamic response of a solid-fluid fully coupled material (Biot theory [38]). The dimensions of the elements inside the mesh followed the previous contributions [19,20] and the performance of the lateral boundaries was verified Appl. Sci. 2021, 11, 6117 4 of 16 by comparing the accelerations at the top of the mesh with those obtained under Free Field (FF) conditions to reproduce wave mechanisms. Mesh discretization was derived from and considering 120 m/s as the lowest shear wave velocity, and 10 Hz as maximum frequency, following suggestions from [39]. The backbone curves of the infill and foundation soil are shown in Figures 3 and 4, while the material parameters in Tables 1 and 2, by considering the values adopted in [40,41]. Particular attention to boundary conditions was considered: (1) absorbing boundaries were applied at the base to dissipate the radiating waves; (2) in order to model the elastic half-space below the mesh, the base nodes were set free to move along longitudinal and transversal directions, while vertical direction was fixed; (3) lateral nodes were constrained to simulate pure shear by applying period boundaries and to ensure free field conditions. At the lateral nodes, penalty method (tolerance: 10 −4 ) was adopted to avoid problems with equations system conditions [20].

Structural Model
The benchmark MIRC structure was calibrated to be representative of a reside 3-storey concrete building with masonry walls and previously modeled [20] as Buil 1 (B1) that is built up with three floors (3.4 m each, total height: 10.2 m), 4 column longitudinal direction (8 m spaced) and 2 ones in transversal direction (10 m spaced) columns and the beams were modeled with fiber sections, by applying Concrete02 m rial for the core and the cove (Figure 5a,b, Table 3). The unloading slope (at maxim strength) and the initial slope are related with a ratio of 0.1. The bars were represente

Structural Model
The benchmark MIRC structure was calibrated to be representative of a residential 3-storey concrete building with masonry walls and previously modeled [20] as Building 1 (B1) that is built up with three floors (3.4 m each, total height: 10.2 m), 4 columns in longitudinal direction (8 m spaced) and 2 ones in transversal direction (10 m spaced). The columns and the beams were modeled with fiber sections, by applying Concrete02 material for the core and the cove (Figure 5a,b, Table 3). The unloading slope (at maximum strength) and the initial slope are related with a ratio of 0.1. The bars were represented by Steel02 material (see Table 4) and 0.01 was chosen as the ratio between post-yield tangent and initial elastic tangent ( Figure 6). Following [38], the transition parameters were: R0 = 15, CR1 and CR2 = 0.925 and 0.15, respectively, as suggested by [38]. Figure 7 shows the diagram moment-curvature of the section. Table 5 shows the characteristics of the elements in terms of cross sections and reinforcement details. Diagonal elasticBeamColumn elements [38] were used to model the masonry walls (both longitudinally and transver-

Structural Model
The benchmark MIRC structure was calibrated to be representative of a residential 3storey concrete building with masonry walls and previously modeled [20] as Building 1 (B1) that is built up with three floors (3.4 m each, total height: 10.2 m), 4 columns in longitudinal direction (8 m spaced) and 2 ones in transversal direction (10 m spaced). The columns and the beams were modeled with fiber sections, by applying Concrete02 material for the core and the cove (Figure 5a,b and Table 3). The unloading slope (at maximum strength) and the initial slope are related with a ratio of 0.1. The bars were represented by Steel02 material (see Table 4) and 0.01 was chosen as the ratio between post-yield tangent and initial elastic tangent ( Figure 6). Following [38], the transition parameters were: R0 = 15, CR1 and CR2 = 0.925 and 0.15, respectively, as suggested by [38]. Figure 7 shows the diagram moment-curvature of the section. Table 5 shows the characteristics of the elements in terms of cross sections and reinforcement details. Diagonal elasticBeamColumn elements [38] were used to model the masonry walls (both longitudinally and transversally) and with selected properties taken from Table C8A.2.1 [42] and described in Table 6. In this regard, this may be considered a limited assumption since the interaction between frame and infill walls depends on the model adopted to simulate the infill wall behavior. In particular, this is due to the complexity of the problem and the computational burden required to complete the analyses, generally simplified models are employed to account for the nonlinear behavior of the superstructure, as shown in [43].

Foundation Model
The foundation slab (28.4 m × 34.4 m, 0.5 m) was chosen in order to represent the most vulnerable typology to earthquakes, since its strength depends only on contact pressures, without frictional resistance (as it occurs for piles foundations). The rigid concrete slab was modeled by applying equaled of [38] to connect the nodes at the base of the columns with those of the soil domain. In order to simulate the interface between the columns and the slab, horizontal rigid links were defined, following [41]. The foundation slab was modelled elastically with an equivalent material that simulates the concrete and implementing the Pressure-Independent Multi-Yield model (PIMY), Table 7 (more details in [38]). The design of the foundation consisted in assessing the eccentricity for the most severe condition: minimum vertical loads (gravity and seismic loads) and maximum bending moments.

Methodology
The seismic vulnerability of structures has been defined with several approaches and, in the last decades, the probabilistic-based approaches have been particularly developed. Among these, developing of analytical fragility curves has been one of the most applied in order to consider several sources of uncertainties in the demand and capacity definitions. Many contributions proposed fragility curves for RC buildings [44][45][46], masonry buildings [47,48] and for different typologies of structures [49] and infrastructures [50].
Fragility curves graphically represent the relationship between the probability of exceedance of a determined level of damage and the intensity of different seismic scenarios. In particular, Ref. [29] proposed to consider four limit states (slight, moderate, extensive, complete). It is worth noting that the development of fragility curves helps to account for several sources of variabilities, such as the dispersion due to the record-to-record and the intra-model variability and that relative to the definition of the damage states. In addition, most recent developments in probabilistic approaches of the seismic assessment have been proposed, such as [51,52].
In the present paper, PGA is selected as Intensity Measure (IM), scaled to multiple levels from zero to 1.0 g, with steps of 0.1 g, by applying the approach called Incremental Dynamic Analysis (IDA), by [53]. The seismic vulnerability was assessed herein by considering (1) the potential of liquefaction at foundation level in terms of settlement and (2) the structural performance in terms of top-floor drift.
Ten input motions (Figure 8) from the PEER NGA database (http://peer.berkeley. edu/nga/, accessed on 25 June 2021) were selected by following the Eurocode 8, Part 1 prescriptions [54] in order to be compatible with the code-specified spectrum for the lifesafety limit state (corresponding to an earthquake with a return period of 475 years, lat.: 42.333 N, 14.246 E., S = 1.50714, Tb = 0.257 s, Tc = 0.770 s, Td = 2.538, Cc = 2.02847 [54]). Figure 8 shows the selected elastic response spectra (5% damping) and the mean spectral ordinates. Note that the mean spectrum (calculated among the selected ones) is included between the −10% and +30% of the code-based spectral shape, as prescribed in [54]. Table 8 shows the characteristics of the selected input motions in terms of Magnitude (Mw), closest distance, peak-to-peak velocity (PPV), pulse period, peak ground acceleration (PGA), peak ground velocity (PGV) and peak ground displacement (PGD). The inputs are divided into 5 Near-Field (NF) and 5 Far-Field (FF) motions, in order to explore the importance of source distance in earthquake scenarios. In addition, the set was selected to have significant dispersion of intensities to significantly affect the dynamic characteristics of the system.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 16 have significant dispersion of intensities to significantly affect the dynamic characteristics of the system. Therefore, IDA nonlinear analyses were performed and since they were particularly challenging and time-consuming, it was necessary to apply four subsequent sub-steps [19,20]. In step 1, the soil alone was loaded and linear properties (weight, shear and bulk modulus) were considered. The presence of the water (0.00 m depth) was introduced in step 2. In step 3, the structure was added and structural loads were applied. The soil properties were changed from elastic to plastic (in 25 load steps). Finally, in the last steps, the dynamic analyses were performed and Newton Line Search algorithm was used to increase the speed of the solution [38]. In order to save computational time, the longitudinal components of the input motions were considered only and applied at the base of the model.
The study assumes that the probabilistic quantities were defined in terms of lognormal standard deviation (β), the median (µ) of the lognormal seismic intensity measure (PGA). Therefore, linear regressions were built up from the performed IDA analyses and the probability of exceedance for a specific intensity (PGA) was calculated as: where: P is the probability of the structural damage (D) exceeding the i-th damage state (C); ϕ is the standard normal cumulative distribution function.    Therefore, IDA nonlinear analyses were performed and since they were particularly challenging and time-consuming, it was necessary to apply four subsequent substeps [19,20]. In step 1, the soil alone was loaded and linear properties (weight, shear and bulk modulus) were considered. The presence of the water (0.00 m depth) was introduced in step 2. In step 3, the structure was added and structural loads were applied. The soil properties were changed from elastic to plastic (in 25 load steps). Finally, in the last steps, the dynamic analyses were performed and Newton Line Search algorithm was used to increase the speed of the solution [38]. In order to save computational time, the longitudinal components of the input motions were considered only and applied at the base of the model.
The study assumes that the probabilistic quantities were defined in terms of lognormal standard deviation (β), the median (µ) of the lognormal seismic intensity measure (PGA). Therefore, linear regressions were built up from the performed IDA analyses and the probability of exceedance for a specific intensity (PGA) was calculated as: where: P is the probability of the structural damage (D) exceeding the ith damage state (C); Φ is the standard normal cumulative distribution function.

Results
In this section, the developed analytical fragility curves are presented and discussed, by considering the potentiality of liquefaction (with the settlement at the foundation level) and the structural fragility by considering the inter-story drifts, respectively, for the soil and for the structure. In particular, the structural damage was considered conservatively defined with the limit states proposed for High and Moderate Code Seismic concrete moment frame (see [29]): Slight, Moderate, Extensive and Complete Damage (corresponding with drift: 0.33%, 0.58%, 1.56% and 4%), named here LS1, LS2, LS3 and LS4, respectively. More details are available in [43].

Soil Results
Liquefaction condition was defined when the pore pressure ration (ru) reaches the unit value (ru = 1), ru being the ratio between the pore pressure increase and the initial vertical effective overburden stress. Activation of liquefaction was assessed at the center of the foundation since the presence of the structure mostly affects this point, as demonstrated in [20]. Liquefaction potential assessment is shown here by considering the profile of the maximum settlement for the most severe input motion (RRS), Figure 9. The linear regression is plotted and R 2 coefficient was chosen to indicate the quality of correlation, and calculated as: where: n is the number of results; µ x and µ y are the mean for the x and y values; σ x and σ y are the standard deviation for the y values.   In this case, R 2 = 0.7702, meaning a good correlation between the variables (PGA and settlements). The slope of the fragility curve depends on the lognormal standard deviation (β): high values increase the probabilities at low IMs, decreasing the probability at high IMs. The central tendency of the fragility curve instead is defined by the median value (µ): bigger values mean higher probability of exceeding a specific damage state. This study considers the serviceability level SLS1 (1/25), as a requirement for all structure of importance level 2 or above, following the New Zealand code (NZS 1170(NZS .0:2002. The values of β and µ (0.9254 and 0.98 g, respectively) are calculated among the results of the analyses ( Figure 10). Therefore, the developed fragility curve is shown in Figure 11, showing that the curve represents a severe increase for low intensities. For example, given an input of PGA = 0.175 g (considerably low value of intensity), the probability of exceedance is 0.666. These values confirmed what is shown in [20]: rigid and heavy superstructure are significantly detrimental for liquefaction, due to the aggravating vertical stresses.    Figure 12 shows the relationships between the base shear force and the longitudinal displacement (lateral spread) at the central column of the model for input motion RRS, considered as a reference. It is worth noticing the large displacements due to the cyclic mobility induced by liquefaction in the soil. It is important to consider that such high nonlinear mechanism was reproduced by the 3D numerical model that is shown to be able to reproduce the tendency for dilation and thus the large shear excursions during the strong shake. Then, the top-floor drifts were calculated as the ratio between the longitudinal displacements and the total height of the building. Figure 13 shows the results of the IDA analyses in terms of PGA and drift and the corresponding value of R 2 (equal to 0.8178 and  Figure 12 shows the relationships between the base shear force and the longitudinal displacement (lateral spread) at the central column of the model for input motion RRS, considered as a reference. It is worth noticing the large displacements due to the cyclic mobility induced by liquefaction in the soil. It is important to consider that such high non-linear mechanism was reproduced by the 3D numerical model that is shown to be able to reproduce the tendency for dilation and thus the large shear excursions during the strong shake. Then, the top-floor drifts were calculated as the ratio between the longitudinal displacements and the total height of the building. Figure 13 shows the results of the IDA analyses in terms of PGA and drift and the corresponding value of R 2 (equal to 0.8178 and calculated from formula (2)) that shows the quality of the linearization among the obtained data. Horizontal lines represent the limit states LS1, LS2, LS3 and LS4 (longitudinal drift: 0.33%, 0.58%, 1.56% and 4%, respectively), whose exceedance is probabilistically connected with a certain amount of structural damage. Table 9 shows the comparison between the values of β and µ for the four considered limit states: β varies relatively, meaning that the level of uncertainties (the slope of the curves) is similar for every limit state. The analytical fragility curves are shown in Figure 14 for LS4, LS3, LS2 and LS1. In particular, for a relatively low value of PGA equal to 0.20 g, the probabilities are: 0.0384, 0.182, 0.692 and 0.942, respectively, for LS4, LS3, LS2 and LS1, showing that PE values vary significantly. The difference between LS1 and LS2 curves decreases for PGA bigger than 0.40 g, while LS3 and LS4 curves are different at all the intensities. The resulted fragility curves may be compared with those derived by [28], where 2-and 3-storey buildings were considered. In particular, the herein resulted values of the median values are relatively close to the 2-storey results for LS1 and LS2. For LS3 and LS4, the calculated values are more similar to those resulted for 4-storey floors. However, if we compare these results with those presented in [28] qualitatively, it is possible to see that the presented curves shift toward more fragile behaviour because of the presence of soft soils that increase the vulnerability of the system (soil + structure) to liquefaction-induced effects, due to soil nonlinearity, as shown in [27]. Overall, it is worth noting that the proposed numerical models represent the complex phenomenon of liquefaction by reproducing dynamic nonlinear effects. This makes the developed analytical fragility curves conservative, especially at low intensities where the studied building was found to be particularly vulnerable. The findings are limited to the specific case study but the 3D numerical model may be applied to assess other cases with different characteristics of the soil and of the structure. This will be the object of future and more extended studies.

Structural Results
Sci. 2021, 11, x FOR PEER REVIEW 13 of represent the complex phenomenon of liquefaction by reproducing dynamic nonline effects. This makes the developed analytical fragility curves conservative, especially low intensities where the studied building was found to be particularly vulnerable. T findings are limited to the specific case study but the 3D numerical model may be appli to assess other cases with different characteristics of the soil and of the structure. This w be the object of future and more extended studies.

Conclusions
The paper developed fragility curves for a MIRC structure subjected to liquefaction, investigating the role of soil deformability both for the soil (in terms of settlements at the foundation level) and for the structure (in terms of drifts at the top level). Fragility curves show (1) the role of SSI on the seismic vulnerability of the system (soil + structure) and that (2) structural drifts that affect the dynamic behaviour are mainly due to deformations in the soil. The proposed FEM model was demonstrated to reproduce the 3D mutual in-

Conclusions
The paper developed fragility curves for a MIRC structure subjected to liquefaction, investigating the role of soil deformability both for the soil (in terms of settlements at the foundation level) and for the structure (in terms of drifts at the top level). Fragility curves show (1) the role of SSI on the seismic vulnerability of the system (soil + structure) and that (2) structural drifts that affect the dynamic behaviour are mainly due to deformations in the soil. The proposed FEM model was demonstrated to reproduce the 3D mutual interaction between the soil, the foundation and structure. The presented outcomes are limited to the considered conditions (in terms of soil profile and direction of input motion), but may potentially be useful for design and pre-and post-earthquake assessments. In particular, the model of infill RC structures can be refined by considering more sophisticated material models and computational modeling strategies. This will be the object of future studies.