1. Introduction
Historical earthquakes (Niigata, Japan 1964, Dagupan City, Philippines 1990, Chi-Chi, Taiwan 1999, Japan 2011, Kocaeli, Turkey 1999 and Christchurch, New Zealand, 2011) demonstrated the importance of estimating liquefaction-induced damages, mainly in terms of permanent lateral spread and settlements. In this regard, liquefaction potential assessment has been based on empirical data (e.g., [
1,
2,
3,
4]) and many relationships between liquefaction resistance and soil parameters have been proposed [
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15]. Among them, Halder, A. [
7] applied SPT results to a second-moment statistical analysis as a development of the empirical procedure introduced by Robertson, P.K [
12]. In the Italian background, Crespellani, T. [
16,
17,
18] proposed several approaches for the assessment of liquefaction potential. In this regard, Di Ludovico, M. [
19] developed empirical fragility curves to analyse liquefaction-induced effects during the 2012 Emilia Romagna earthquake. Moreover, Kafali, C. [
20] estimated the probability of liquefaction as a function of earthquake load and SPT resistance by logistic regression analyses, while other researchers proposed an advanced first-order second-moment (AFOSM) technique to calculate the probability of failure, [
9]. Juang, C.H. [
8] performed an extensive series of sensitivity analyses using the first-order reliability method associated with mapping functions to characterise the uncertainties in the performance function.
However, as pointed out by Juang, C.H. [
8], such methodologies aimed to assess ground and foundation deformation but rarely to determine liquefaction effects on existing structures and the subsequent liquefaction-induced structural damages. In this regard, the determination of liquefaction potential has become a fundamental issue in the assessment of the induced damage to civil engineering structures by applying several methodologies. A class of approaches was based on the use of artificial neural networks (ANN) models [
21,
22,
23,
24,
25,
26] that may estimate the relationships between the soil and earthquake characteristics with the liquefaction potential, requiring no prior knowledge of the form of the relationships. However, the main disadvantage of such approaches consists of the deterministic definition of soil parameters that cannot consider many uncertainties. Recently, the development of numerical simulations of nonlinear dynamic soil–structure-interaction (SSI) analyses provided probabilistic-based insights into the liquefaction-induced effects (e.g., [
27,
28,
29,
30,
31,
32,
33,
34,
35,
36,
37,
38]).
Furthermore, damage prediction is essential for the assessment of economic losses, and fragility curves are interesting representations of its cumulative distribution, providing helpful information on the probability that specific levels of damage have risen or exceeded. Fragility curves have been applied to estimate structural damages due to ground shaking (e.g., [
39,
40]) as a way of expressing the probability of exceeding specific limit states. In this regard, probabilistic-based approaches based on analytical fragility curves are widely developed in the literature for structures and infrastructures [
41,
42,
43,
44,
45,
46]. Recently, Liang, H. [
47] proposed one of the most recent applications of fragility curves to study the seismic vulnerability of a high concrete arch dam based on seismic instability and taking the friction coefficients and cohesions as random parameters. In another contribution, Caverzan, A. [
41] applied analytical fragility curves to study the influence of dynamic SSI effects on the assessment of the seismic vulnerability of buildings. In particular, Shinozuka, M. [
48] developed a statistical procedure for reliability analyses by applying two-parameter lognormal distribution functions to develop empirical and analytical fragility curves for bridges. The Monte Carlo approach was also applied by Kafali, C. and Hwang, H.M. [
20,
49].
Fragility curves may include many sources of uncertainties, such as the ones connected with the soil characterisations and properties. For example, Hwang, H.M. [
49] considered uncertainties, such as viscous damping ratio, strength and stiffness of structural materials, Shinozuka, M. [
48] included variability in ground motions, and Popescu, R. [
50] developed fragility curves that account for the variability of soil properties during liquefaction. In this regard, the study of liquefaction effects on structural fragility curves was the object of several contributions, such as Aygün, B. [
51] who performed a coupled 3D bridge-foundation system with 2D heterogeneous soil strata finite element (FE) models to monitor multiple bridge failure mechanisms in the presence of 1D soil site amplifications and liquefaction effects. Recently, Fotopoulou, S. [
52] presented a numerical framework for the vulnerability assessment of low-code Masonry Reinforced Conrcete (MRF RC) buildings subjected to liquefaction-induced differential displacements.
In the previous contributions, fragility curves were used to include liquefaction to demonstrate its role in the assessment of structural vulnerability. The novelty of the presented study instead consists of proposing the application of fragility curves to estimate liquefaction potential and applying it to two case studies. Analytical fragility curves are here developed to estimate the probability of liquefaction in terms of lateral spread and settlement with 3D numerical models, first for free field (FF) conditions, and then considering SSI effects.
2. Liquefaction-Potential Curves
The fragility curve has been proposed as a well-established approach in earthquake engineering to consider uncertainties in demand and capacity. They represent graphical relationships between the conditional probability of exceeding predefined damage states (DS) for given levels of EDP (engineering demand parameters) and are used to estimate the amount of damage for a particular level of shaking defined with intensity measures (IM). For example, fragility curves, based on HAZUS-MH MR5: Technical Manual [
43], consist of considering four states (slight, moderate, extensive, complete), and they display the probability that the response of a system exceeds a given threshold limit given a certain excitation level, as applied in Ranjbar, P.R.’s work [
46]. Lateral spread displacements and settlement at the surface were chosen as EDPs to assess the seismic effects of liquefaction on FF and SSI systems and thus developing analytical fragility curves to define liquefaction-induced damage. Peak Ground Acceleration (PGA) is selected as the IM, scaled to multiple levels from zero to 1.0 g with steps of 0.1 g, following the well-credited approach of incremental dynamic analysis (IDA) proposed by Hwang, H.M. [
49]. Thus, a number of 10 × 7 IDA nonlinear analyses were performed (with Opensees PL) for each case study (total 140). The input motions (
Figure 1 and
Table 1) were selected from the Pacific Earthquake Engineering Research (PEER) Center Next generation Attenuation (NGA) database (
http://peer.berkeley.edu/nga/), ([
53,
54] on the base of Eurocode prescriptions) to significantly affect the dynamic characteristics of the system and was applied at the base of the models, along the x-axis (longitudinal direction). In particular, IDA requires a series of nonlinear, dynamic time history analyses for the selected ensemble of ground motions with increasing intensity levels to cover the entire range of the response, from elastic behaviour through liquefaction condition. Therefore, the limit state was chosen as the liquefaction condition (ru = 1). Pore pressure ratio (ru), defined previously in [
5], considers the pore pressure as a fraction of the vertical stress (including the weight of ponded water, if the material is submerged) [
32]. This study assumes that all uncertainties in the fragility curves can be represented by lognormal distributions, and, thus, only two parameters were needed to plot the curves, the logarithmic mean (µ) and standard deviation (β) of the lognormal seismic intensity measure (PGA). In particular, there are several reasons that justify the use of lognormal distribution: (1) its simplicity in approximating an uncertainty quantity that must take on a positive value, using only an estimate of the central value and uncertainty, (2) the fact that it has been widely used for several decades in earthquake engineering and (3) it is the distribution that assumes the least knowledge (the only variables are the mean and logarithmic standard deviation). Therefore, IDA analysis results were represented to build linear regression and define the mean and log-standard deviation values. Moreover, the probability of liquefaction (ru = 1) for a specific intensity (PGA) was then calculated as
where
ϕ is the standard normal cumulative distribution function.
3. Finite Elements Models
Finite element models were performed by OpenSeesPL [
23,
35], from the Pacific Earthquake Engineering Research (PEER) Center. 3D finite models were developed to capture many outputs (e.g., the mechanism of excess pore pressure, accelerations, displacements, settlements, tilts, and inter-story drifts) that are otherwise impossible to be assessed with 1D numerical simulations. The study here proposed was divided into two steps (
Table 2), while 3D numerical models are shown in
Figure 2 and
Figure 3.
Step 1. A 3D numerical model was developed to build liquefaction potential curves in FF conditions. In particular, the soil performed in [
33] was considered. The 3D soil models consisted of a 100 × 100 × 20 m mesh, performed with 8000 20-node BrickUP elements and 9163 nodes to simulate the dynamic response of solid-fluid fully coupled material [
54,
55]. Finite elements were defined on the wavelengths of the seismic signal and the maximum frequency above which the spectral content of the input was considered negligible (15 Hz). The dimension of the elements was increased from the structure to the lateral boundaries that were modelled to behave in pure shear and located far away from the structures, as previously applied in [
32,
33]. For each BrickUP element, 20 nodes described the solid translational degrees of freedom, while the eight nodes on the corners represented the fluid pressure 4 degrees-of-freedom (DOF). For each node, DOFs 1, 2 and 3, represent solid displacement (µ), and DOF 4 describes fluid pressure (p), which was recorded using OpenSees Node Recorder [
55,
56] at the corresponding integration points. The soil model applied the two-phase material µ -p formulation (where µ is the displacement of the soil skeleton and p is pore pressure) and simulated with PressureDependMultiYield02 (PDMY02) model [
55,
56], based on the multi-yield-surface plasticity framework and calibrated with centrifuge tests and consolidated undrained cyclic triaxial tests (VELACS No. 40–58), more details in [
32,
33].
Table 3 shows the adopted parameters, such as the low-strain shear modulus and friction angle, as well as parameters to control the contraction (c1), dilatancy (d1 and d2) and the level of liquefaction-induced yield strain (l1, l2 and l3), [
55,
56]. The backbone curve is shown in
Figure 4. The penalty method was adopted for modelling the boundary conditions (tolerance: of 10
−4), chosen as a compromise to be large enough to ensure strong constrain conditions but not too large to avoid problems associated with conditioning of the system of equations. Base boundaries (20 m depth) were considered rigid in all directions, and the input motions were applied at the base as acceleration time histories [
32,
33]. Vertical direction (described by the 3rd DOF) was constrained in correspondence with all the boundaries (at the base and lateral), while longitudinal and transversal directions were left unconstrained at lateral boundaries to allow shear deformations. The definition of the mesh dimensions followed the approach already adopted [
32,
33,
53,
54] and, to verify proper simulation of FF conditions at lateral boundaries, accelerations at the top of the meshes were compared with the FF ones, that were found to be identical, confirming the effective performance of the mesh.
Step 2 consisted of considering a complete system (soil + foundation + structure),
Figure 3. The benchmark structure was calibrated to represent residential buildings and already applied in [
33]. In this regard, a two-storey concrete structure was selected (floor height: 3.40 m, total heights: 6.80 m), with a 4 × 2 column scheme (4 in the longitudinal direction (8 m spaced) and 2 in the transversal direction (10 m spaced)). Some simplifications consisted of considering a shear-type behaviour (with plan and vertical regularity) and concentrating the seismic masses at each floor level.
Table 4 shows the dynamic characteristics of the structure, assumed to be linear elastic and modelled with elastic-beam column elements (
Table 5). A concrete shallow foundation (28.4 m × 34.4 m × 0.50 m) was considered to represent recurrent typologies for residential buildings. The foundation was performed as rigid by linking all the nodes at the base of the columns together and to boundaries, by applying the equaldof [
55,
56]. The foundation was designed by calculating the eccentricity (the ratio between the overturning bending moment at the foundation level and the vertical forces) in the most detrimental condition of minimum vertical loads (gravity and seismic loads) and maximum bending moments. The foundation was modelled with a Pressure Independent MultiYield [
55,
56] as an equivalent concrete nonlinear hysteretic material with a Von Mises multi-surface kinematic plasticity model to simulate the monotonic or cyclic response of materials whose shear behaviour was insensitive to the confinement change,
Table 6. The first 0.5 m of soil around the foundation was constructed of an infill layer defined by PressureDependMultiYield, and
Table 7 shows the adopted parameters. The number of yield surface was equal to 20.
Figure 5 shows the Backbone curve.
It is worth noting that due to the large amount of nonlinear dynamic analyses, the following simplifications were introduced to improve the computation speed: (1) The foundations were modelled as rigid slabs. (2) Deformations inside the foundation, weak soil, and intermediate nodes along the height of the superstructure were calculated during the analyses but not saved and memorised as results. (3) Concrete material was modelled as equivalent linear elastic. (4) The structure was modelled as a shear-type building by introducing rigid diaphragms. (5) The soil was modelled as an equivalent one-layered with constant shear wave velocity and (6) the base of the mesh was considered fixed (with no absorbing boundaries). These assumptions allowed the maintenance of the computation of response for the performed earthquake records in less than 72 h by using the parallel computation procedure.
4. Liquefaction-Potential Curves for Free-Field
In this section, step 1 (FF conditions) is considered. 3D nonlinear dynamic analyses (total 70) were performed, and the maximum longitudinal displacements (named lateral spread, LS) and vertical displacements (named settlements, ST) were calculated for each.
Figure 6 shows the results for input motion 2 (100%). For each analysis, the LS and ST at the surface in correspondence with the centre of the mesh were considered at with liquefaction condition (ru = 1) and defined on scaled PGA values (
Figure 7 and
Figure 8) to evaluate the variation of demand and capacity with ground motion characteristics. In particular, standard deviation (β) defines the dispersion about the mean value, and it significantly affects the slope of the curves. In particular, higher values of β increase the probabilities at low IMs and, at the same time, they decrease the probability at high IMs. The mean value (µ) defines the central tendency of the curves affecting the translation of them inside the range of values. Thus, bigger values of the mean correspond to a higher probability of reaching or exceeding a particular damage state. In particular,
Table 8 shows the comparison between the values of β and µ that affect the difference between the two liquefaction-potential (LP) curves. LP probabilities at low PGA were bigger in correspondence with ST-FF (β = 0.624) than with LS-FF (β = 0.595). The mean values (µ) also were significantly different.
Figure 9 shows the LP (liquefaction potential) curves for both the considered parameters (LS and ST), showing that both the curves represented a severe increase in the LP at a low value of PGA (low intensities). Comparing the two curves, it is worth noting that when lateral spread was considered, there was a significant difference of LP in the range of PGA = 0.20–0.30 g. For example, given an input PGA = 0.26 g, LP = 0.575 and 0.284 (around 50%), respectively, for ST-FF and LS-FF. If other levels of PGA were considered, this difference decreased significantly (at least to zero). In particular, the difference decreased gradually as the severity of the motion increased until around 0.80–0.90 g, where the curves were almost equivalent. For PGA values less than 0.20 g, referring to settlement was still more conservative than considering lateral spread. The reason for the difference is due to the fact that settlements are more sensitive to inertial forces connected with the vertical stresses (especially the weight) and start at a low value of PGA, while lateral spread needs a higher level of intensity to occur since they are connected with shear mechanisms.
5. LP Curves for SSI
Step 2 is shown in this section, and the role of the structure is discussed by calculating LS and ST for the 70 3D nonlinear dynamic analyses.
Figure 10 shows the results for input motion 2 (100%) to give an example of the soil behaviour. LS and ST were calculated at the surface in correspondence with the centre of the foundation that was demonstrated (see [
47]) being mostly affected by the presence of the structure (mainly in terms of vertical stresses). LS and ST were considered in correspondence with liquefaction condition (ru = 1) and plotted in
Figure 11 and
Figure 12.
Table 9 shows the comparison between the values of β and µ that affects the difference between the two LP curves. In particular, LP probabilities at low PGA were significantly bigger in correspondence with ST-SSI (β = 0.653) than with LS-SSI (β = 0.624). The mean values (µ) were significantly different.
Figure 13 displays LP curves for both the considered parameters (LS and ST) showing that ST-SSI curve represented a severe increase in the LP at a low value of PGA (for earthquakes with low intensity). Comparing these two curves, it is worth noting that when lateral spread was considered, there was a significant difference of LP in the range of PGA = 0.10–0.25 g. For example, given an input PGA = 0.175 g, LP = 0.666 and 0.195 (around 30%), respectively, for ST-SSI and LS-SSI. If other levels of PGA were considered, this difference decreased significantly (at least to zero). In particular, the difference decreased gradually as the severity of the motion increased until around 0.70–0.90, where the curves were almost equivalent. These curves show the role of the structure in aggravating the vertical stresses (due to the structural weight). In particular, the low-rise structure was selected to consider the effect of a rigid and heavy superstructure on the inertial interaction between the structure and the soil that is herein demonstrated to be significantly detrimental in increasing LP values. These results confirm that the damage is mainly governed by the rigid rotations and settlements, and liquefaction works as a natural isolation system against the transmission of inertial seismic actions on the superstructure, as shown in [
19].
The comparison between FF and SSI cases are shown in
Figure 14 and
Figure 15 that plot LP curves for both the considered parameters (LS and ST). It is worth noting the detrimental effects of the presence of the structure in both cases. The role of SSI was particularly significant when the settlement was considered as a damage parameter with a severe increase in the LP values (for a range of PGA that was around 0.10–0.20 g). For example, given an input PGA = 0.148 g, LP = 0.621 and 0.337 (around 54%), respectively, for ST-SSI and ST-FF. These curves show that considering the settlement as a damage parameter is more conservative in assessing the role of the SSI.