Nonlinear Solid–Fluid Coupled Seismic Response Analysis of Layered Liquefiable Deposit

: A seismic response analysis of layered, liquefiable sites plays an important role in the seismic design of both aboveground and underground structures. This study presents a detailed dynamic site response analysis procedure with advanced nonlinear soil constitutive models for non ‐ liquefiable and liquefiable soils in the OpenSees computational platform. The stress ratio con ‐ trolled, bounding surface plasticity constitutive model, PM4Sand, is used to simulate the nonline ‐ ar response of the liquefiable soil layers subjected to two seismic ground motions with different characteristics. The nonlinear hysteretic behavior of the non ‐ liquefiable soil under earthquake exci ‐ tations is captured by the Pressure Independent Multi Yield kinematic plasticity model with a von Mises multi ‐ yield surface. The soil elements are modelled utilizing the solid–fluid fully coupled plane ‐ strain u ‐ p elements. The seismic response of the layered liquefiable site in terms of the de ‐ velopment of excess pore water pressure, acceleration, ground surface settlement, and stress– strain and effective stress path time histories under two representative earthquake excitations are investigated in this study. The numerical results indicate that both the characteristics of ground motions and the site profile have a significant influence on the dynamic response of the layered liquefiable site. Under the same intensity of ground motion, the loose sand layer with a 35% rela ‐ tive density is more prone to liquefaction and contractive deformation, which causes irreversible residual deformation and vertical settlement. The saturated soil layer can effectively filter the high ‐ frequency components and amplify the low ‐ frequency components of ground motions. Moreover, the liquified site produces a 40% post ‐ earthquake consolidation settlement after the excess pore pressure dissipation. test simulations using OpenSees. The results showed that PM4Sand can effectively capture the stress–strain response, strength degradation during cyclic load ‐ ing, and the number of cycles to liquefaction for sand with a relative density of 35% and 75%.


Introduction
Earthquake-induced ground failure, resulting from liquefaction of loose sand, silt, and gravelly soil deposits, has caused tremendous damage to engineering structures and infrastructure. Ground failures due to the lateral spreading of soil liquefaction at sites on horizontally stratified or mildly sloping ground impose significant threats to both aboveground and underground engineering structures. Dramatic damage caused by liquefaction has been reported in past historical earthquakes [1], such as Niigata, Japan 1964, Japan 2011, Christchurch, New Zealand 2011, and Palu, Indonesia 2018. During the 2011 Great East Japan Earthquake, the settlement in the Tokyo Bay area was generally in the range of 30 to 50 cm and settlement at some severely liquefied sites was as large as 70 to 100 cm. This macroscopic phenomenon means that the soil strata generate substantial consolidation settlement due to the generation and dissipation of excess pore water pressure (EPWP), which may lead to significant damage to engineering structures. Therefore, it is critical to study the dynamic characteristics of the sites that are susceptible to seismic liquefaction [2].
Numerical simulation is an important tool to help estimate the site response and assess the seismic safety of engineering sites. Many numerical studies on the seismic response of saturated soil sites have been carried out in recent years. Liyanathirana et al. [3] proposed a coupled effective stress method by combining the Hardin-Drnevich model with the pore pressure growth model proposed by Seed et al. [4]. Then, the proposed method was validated against the seismic acceleration records of the Port Island downhole array during the 1995 Hanshin earthquake in Japan. Ziotopoulou et al. (2012) [5] numerically simulated the seismic liquefaction process of the two sites utilizing three different constitutive models after calibrating the parameters of the nonlinear constitutive models against field test data. The results indicated that dilatancy is an important factor in capturing appropriate ground motions in liquefiable sites. Ramirez et al. (2018) [6] systematically compared the response characteristics of the liquefiable site using two different constitutive models of sand (i.e., Pressure Dependent Multi Yield02 model and SANISAND model), in terms of settlement, acceleration, and excess pore water pressure, under seismic loading using two numerical computational platforms, OpenSees and FLAC. By comparing the results of centrifuge-free field tests, it was discovered that the two constitutive models can predict most seismic responses for the liquefiable site, but there is still a need for improvement in simulating soil settlement. Tropeano et al. (2019) [7] developed a simplified pore pressure model based on the stress-based approach in a 1D computer code to carry out effective stress analyses of the liquefiable site. The code provides reliable predictions by comparing the test results of ideal soil profiles, a centrifuge test, and seismic events recorded on a well-instrumented test site. Based on the physical modeling and numerical simulation, Adampira and Derakhshandi (2020) [8] evaluated the seismic behavior of liquefiable sublayers and their effect on the seismic ground response. The study showed that the thicker and shallower liquefiable sublayers can effectively reduce the energy of seismic waves and earthquake-induced forces and the thinner liquefiable sublayers at greater depths could increase seismic intensity and reduce settlement. Chen et al. (2021) [9] proposed a loosely coupled nonlinear effective stress method for site seismic response analysis. In this method, the nonlinear hysteresis model of soil was combined with an EPWP generation model characterized by cyclic shear-volume strain coupling, which established the coupling relationship between the degradation of cyclic stiffness and the generation of EPWP associated with earthquake events. This method was then used to simulate the seismic response of non-liquefied and liquefied sites with a downhole seismic array in Japan. The simulated results are in good agreement with the recorded spectral acceleration curves of the non-liquefied and liquefied sites. Wang et al. (2016) [10] established a three-dimensional finite element dynamic analysis model of piles in a liquefiable site based on the OpenSees software platform and conducted numerical simulations for the centrifuge shaking table tests of pile in liquefiable horizontal foundations. Special attention was mainly dedicated towards the ground acceleration, EPWP, and the bending moment of the pile shaft during the shaking. However, the settlement and pore water pressure dissipation process of the soil strata after the earthquake were not mentioned in the study.
In view of the above, the dynamic response of saturated soils under earthquake excitations is a complex mechanical process that includes the nonlinear hysteretic behavior of the soil and the coupled volumetric-distortional deformation. The former is composed of a decrease in the shear modulus associated with an increase in energy dissipation due to the damping. The latter is related to plastic irreversible strain under high confinement. The effects of liquefaction on geotechnical structures are increasingly being evaluated in practice using nonlinear time history analyses, but a 'loosely coupled' approach that predicts pore pressure by adopting relationships used in combination with total stress constitutive models (e.g., pore pressure prediction based on accumulated strains or stresses) is often used in the previous studies [3,9,11]. This approach fails to ef-fectively consider the influence of dilatancy and the post-liquefaction shear deformation of the soil strata. In this study, by concentrating on the effects of liquefaction, a 'fully coupled' approach is adopted to evaluate the nonlinear seismic response of liquefiable multi-layered soils under earthquake loading and post-earthquake site response. The advanced soil constitutive model is used to characterize the complicated nonlinear dynamic behavior of sand. Based on the numerical analyses, the dynamic response characteristics associated with the soil liquefaction are described in detail in this study.

Liquefiable Site Model
The Open System for Earthquake Engineering Simulation (OpenSees) [12] is an open-source software framework that can perform the pseudo-static and dynamic analyses with the finite element method and add additional constitutive models and algorithms based on user needs. Moreover, OpenSees has good capability for a nonlinear analysis of soil and structure, and can perform a comprehensive analysis from the material stress-strain level to the section level to the element level, and finally, the structure level. This computational platform, OpenSees, has been employed to simulate the nonlinear behaviors of soil-structure interaction systems under static and dynamic loading by many researchers [13][14][15][16]. Figure 1 shows the schematic diagram of a generic site consisting of three soil layers, following existing studies [17,18] without the initial static shear stress. The soil domain was discretized using four-node quadrilateral elements with u-p formulations [19] such that water pressure and soil deformation could be computed simultaneously. The overall depth and width of the soil deposit are 21 m and 50 m, respectively. The soil deposit mainly consists of sand with relative density, Dr = 35% and 75% and crust clayed soil. The water table is assumed to be coincident with the ground surface level. Therefore, the ground surface is supposed to be a free drainage surface, while the lateral and base boundaries are impermeable. Horizontal kinematic constraints are introduced to the nodes on two side boundaries to ensure the same horizontal movement of the two nodes at the same burial depth and effectively simulate the shear deformations of the soil layers under upward propagation of in-plane waves. The analyses consist of two steps: a gravity application step and a dynamic excitation step. Gravity is firstly applied to the model with the nodes fixed at the base to generate the distributions of the initial effective confinement and pore pressure in the soil model for the site before the application of the seismic excitation. The permeabilities for all of the elements in the gravity analysis were initially set to 1.0 m/s to facilitate the development of hydrostatic conditions in the model. Subsequently, the seismic input motion is applied to the base of the soil model as a force-time history at the node, which shares equal degrees-offreedom with the Lysmer-Kuhlemeyer [20] dashpot. This dashpot depicts the finite rigidity of the underlying bedrock by using a bedrock shear wave velocity of 760 m/s for the soil profiles and a mass density of 2200 kg/m 3 . Meanwhile, the nodes at the base of the soil model are only constrained in the vertical direction and are left free to move in the horizontal direction.
In the dynamic analysis, the constant average acceleration Newmark integrator method with the parameters γ = 0.6 and β = 0.3025 is used in order to resolve the integration in time. The energy increment test with a tolerance of 1.0 × 10 −6 was utilized to control convergence checking. For each time step, the solution is obtained using the Krylov-Newton algorithm with Krylov subspace acceleration [21]. Rayleigh damping with an initial damping ratio of 3% is applied for the purpose of stabilizing the numerical results in dynamic numerical analyses. For the soil deposit, the first-and second-order natural frequencies of the site are 0.46 Hz and 1.42 Hz, which were obtained through model analysis. The damping coefficients for the mass matrix and stiffness matrix are calculated by selecting the damping ratio of 3% and the ground natural frequencies are equal to 0.124 and 0.0008 by numerical model analysis, respectively.

Constitutive Model and Input Parameters
The selection of a constitutive model to describe the behavior of both fully or partially liquefiable soils is critical in the numerical analysis. Many researchers have proposed various types of constitutive models to investigate the cyclic behavior of soil subjected to dynamic loading [22,23].
In this study, the critical-state, bounding surface dynamic elastoplastic constitutive model, PM4Sand model, which was developed by Boulanger and Ziotopoulou [22], is used to simulate the nonlinear response of the saturated sand layers, wherein the model could capture the dilatancy, nonflow liquefaction, permanent shear strain accumulation, and continuously update the pressure dependent shear modulus. These features play an important role for modeling the undrained response of sands subjected to cyclic loading. The PM4Sand model parameters in Table 1 demonstrate the loose sand and dense sand of two relative densities, which have been extensively calibrated to emphasize the realistic modeling of liquefaction behavior [24]. By calibrating three primary input parameters, the shear modulus coefficient, G O , apparent relative density, Dr, contraction rate parameter, and hpo, reasonable approximations of the desired dynamic behavior of soil can be obtained from the numerical analyses, including pore pressure generation and dissipation, limiting strains, and cyclic mobility, in which G O is a parameter controlling the small-strain shear modulus defined through the following equation: where Gmax is the small-strain shear modulus; pA is the atmospheric pressure; and p is the current mean effective stress.
Default values are assigned to the rest of the secondary parameters. Meanwhile, the nonlinear hysteretic behavior of the clay soil under cyclic loading was simulated using the Pressure Independent Multi-Yield kinematic plasticity model [25] with a von Mises multi-yield surface, in which the plasticity appears only in the deviatoric stress-strain response. The volumetric stress-strain response in the linear-elastic response is independent of the deviatoric stress, and it is insensitive to the confining pressure variation. The critical mechanical parameters of the crust clay are listed in Table 2.  As shown in Figure 2, the cyclic response of sand consolidated at the initial vertical stress of 100 kPa is evaluated based on the stress-controlled undrained cyclic direct simple shear (DSS) test simulations using OpenSees. The results showed that PM4Sand can effectively capture the stress-strain response, strength degradation during cyclic loading, and the number of cycles to liquefaction for sand with a relative density of 35% and 75%.

Input Ground Motions
In this study, two acceleration time histories were used as input motions in the horizontal direction of the elastic bedrock for the numerical analyses. The two motions were recorded during the 1940 Imperial Valley earthquake at the El-Centro station (termed as El-Centro motion in this study) and the 1995 Kobe earthquake at the KJMA station (termed as Kobe motion in this study), respectively. The El-Centro motion, with an amplitude of 0.348 g and the main shaking duration being about 26 s, has the obvious characteristics of the middle-to far-field ground motions. The frequency content of the El-Centro motion concentrates in a range of 0~10 Hz. The Kobe motion, with an amplitude of 0.83 g and the main shaking duration being about 10 s, has typical characteristics of near-fault ground motion with a velocity pulse and the concentration range of the excitation frequency is 0~5 Hz. Figure 3 presents the acceleration time histories, Arias intensity time histories, Fourier spectra, and the energy flux time histories of the two input motions after being scaled to a peak bedrock acceleration of 0.30 g. It is worth noting that the energy flux for the scaled El-Centro motion is 2218 J/m 2 /s, which is more than twice that of the scaled Kobe motion (991 J/m 2 /s). The more detailed information of the two selected input motions is introduced, as shown in Table 3.

Numerical Analysis and Results
It can be seen in Figure 1 that the computational domain is symmetrical in x = 25 m. For brevity, in this study, the results of the seismic dynamic analysis in x = 25 m are taken as the representatives to exhibit the dynamic responses of the layered liquefiable site.

Nonlinear Stress Response
Under earthquake excitations, the loose sand exhibits nonlinear behavior. The liquefaction in loose sand is closely related to the build-up of EPWP and the reduction of effective stresses. Figures 4 and 5 present a nonlinear shear stress-strain relation and effective stress path along the soil depth under the El-Centro motion. As shown in Figure  4, the shear stress and strain of the crust clay layer exhibit a linear relationship, indicating that the soil shear stiffness remains the same during the earthquake excitation. For the loose sand layers at the burial depths of 3 m, 10 m, and 15 m, the hysteresis loops of the shear stress and strain gradually flatten, and the shear stiffness gradually decreases under cyclic loading until the saturated loose sand layer reaches the liquefaction initiation state. The cyclic shear strain in the saturated sand layer increases rapidly after the initial liquefaction, wherein the maximum single amplitude shear strain γSA reaches 22.3% in the upper loose sand layer. The maximum single amplitude shear strain γSA in the loose sand layer reaches 2.5% at the interface of the loose sand and dense sand layers, which indicates that the numerical model used in this paper can effectively simulate the cyclic shear deformation characteristics of the saturated sand layer after liquefaction initiation. For the dense sand layer, it is more difficult to liquefy and soil shear stiffness remains almost unchanged, which can be attributed to the surrounding higher effective confining pressure and the soil densification in the dense sand layer. Figures 6 and 7 illustrate the shear stress-strain relationships and effective stress path under the Kobe motion, which are, in general, similar to those of the El-Centro case in Figures 4 and 5.
The responses of the soil layers at different burial depths under the Kobe motion are smaller than those subjected to the El-Centro case. The maximum single amplitude shear strain γSA reaches 16.5% at the interface of the crust clay and loose sand layers, and 3.7% at the interface of the loose sand and dense sand layers. The difference may be attributed to the fact that the El-Centro motion is a characteristic far-field ground motion. The main shaking of the El-Centro motion lasts longer with similar peak acceleration amplitudes, thereby resulting in a higher degree of nonlinear shear deformation at the site. However, the Kobe motion has the feature of a near-field ground motion with fewer repeated cycles, which leads to the smaller shear deformation of the site.

Progressive Liquefaction of Site
Compared with non-liquefied soils, the accumulation of EPWP during an earthquake is a typical seismic response of liquefiable soils. EPWP is a critical indicator for the initiation of liquefaction and the stress state in saturated soils. Figure 8 displays the EPWP along the soil depth at different moments under two input motions. The initial effective vertical stress, σ v0 = γ h, is also plotted as a reference, where γ′ is the submerged unit weight of the loose sand and h is the burial depth. When EPWP at a specific depth reaches the initial effective vertical stress line, the soil loses its shear strength and behaves like a liquid. It can be seen from Figure 8a that due to the liquefiable loose sand layer with an overlying crust clay layer with extremely low permeability, the crust clay layer prevents the pore water from dissipating upward and causes EPWP to gradually accumulate at the interface between the non-liquefied crust clay layer and the liquefied loose sand layer. Therefore, the upper layer of saturated loose sandy soil firstly reaches the liquefaction initiation state within 3 s, and then the soil liquefaction gradually develops downward. The whole loose sand layer reaches the initiation of liquefaction at about 20 s. As shown in Figure 8b, it is observed that most of EPWP are developed during 6~20 s of the Kobe motion, which corresponds to most of its Arias intensity accumulation. The entire loose sand layer liquefies in about 10 s, after which the pore pressure begins to dissipate and the degree of soil liquefaction is significantly less than in the case of the Kobe Motion. It is also worth noting that the excess pore water pressure of the dense sand layer is much lower than the initial effective stress under both input motions. Figure 9 illustrates the shear strain distribution along the soil depth at different moments under two input motions. Large shear strain mainly concentrates in the loose sand layer. Figure 9a demonstrates the maximum shear strain accumulation along the interface between the crust clay and loose sand layers after reaching an initial liquefaction at approximately 2 s until it reaches the maximum shear strain value (γSA = 22.3%) under the El-Centro motion. In comparison with the aforementioned position, the shear strain (γSA = 16.5%) of loose sand at the interface between the loose and dense sand layers is smaller. The difference may be attributed to the energy flux of two input motions. The distribution pattern of the shear strain of the Kobe case is similar to the El-Centro case, as shown in Figure 9b.   Figure 10 presents the acceleration time histories of the free-field at different depths, and it can be seen from the figure that the acceleration response patterns at the depths of 15 m and 18 m are basically consistent with the input acceleration time histories. However, the acceleration amplitudes of the loose sand layer at 3 m and 9 m depths are significantly reduced compared with the input acceleration, thereby showing that the liquefied soil layer dissipates seismic energy due to soil softening and increased damping from the transmission of seismic waves. Besides, the liquefied soil layer prevents the upward propagation of seismic waves, leading to the reduced seismic response of the crust clay layer at the ground surface. Also, as shown in Figures 11 and 12, due to the filtering effect of the liquefied loose sand layer on the high frequency components of ground motion, the acceleration time curves of the soil in the loose sandy soil layer are smoother and the high frequency acceleration spikes are not remarkable in Figure 10. By comparing the input acceleration Fourier spectra with the acceleration Fourier spectra of the loose sand layer and the crust clay layer in Figures 11 and 12, it can be found that the predominant acceleration frequency decreases from 2.5 Hz~4.0 Hz to 0.5 Hz~1.5 Hz due to the liquefaction of the loose sand layer. The obtained numerical results are in good agreement with the findings from the experiment [25].

Settlement Deformation of the Ground Surface
Many studies [26][27][28] have shown that the post-earthquake consolidation settlement of the site with the dissipation of the EPWP is close to the liquefaction subsidence of the site during the earthquake, which cannot be ignored. The PM4Sand constitutive model can better reflect the process of consolidation settlement by reducing the value for the elastic shear and bulk modulus after the earthquake in OpenSees, where the relationship between the dissipation of EPWP and the post-earthquake recovery of soil strength can be well simulated. Figure 13 compares the results of the layered liquefiable site obtained using Open-Sees after the main shaking under the two ground motions. For the El-Centro case, the settlement of the ground surface is 11.46 mm, of which the settlement of 5.3 mm is generated during the post-earthquake pore pressure dissipation period, thereby accounting for 31% of the total settlement. In the Kobe case, a total settlement of 14.60 mm is generated at the ground surface, and 9.77 mm occurs during the post-earthquake pore pressure dissipation period, thereby accounting for 40% of the total settlement.

Further Comparison and Discussion
To further evaluate the efficacy of the fully coupled nonlinear effective stress analysis approach utilized in this study, the numerical results were compared with those obtained using the one-dimensional site response analysis program DEEPSOIL V7.0 [29], which adopts a loosely coupled nonlinear dynamic effective stress analysis approach for layered saturated soil. For comparison purposes, the MKZ (Matasovic-Konder-Zelasko) nonlinear constitutive model loosely coupled with the pore water pressure model proposed by Vucetic and Dobry [30] in DEEPSOIL was selected for the effective stress analysis of the nonlinear seismic response of liquefiable sites in this manuscript. The material parameters of different soil layers are listed in Table 4. Moreover, the pore pressure ratio time histories of the liquefiable site under the El-Centro motion and Kobe motion were obtained from OpenSees and DEEPSOIL, as shown in Figure 14. It can be seen that liquefaction occurs in the soil at −3 m and −15 m, and the moment to reach liquefaction is very close for both analysis methods. However, the DEEPSOIL only approximately simulates the overall EPWP accumulation process without remarkable oscillation phenomenon in the pore pressure ratio time histories. The oscillation phenomenon of EPWP has been frequently observed in liquefiable sites during the previous numerical simulations and dynamic experiments [31][32][33]. Therefore, the fully coupled nonlinear effective stress analysis approach in OpenSees can more realistically capture the development of EPWP in the liquefiable soil during earthquakes.
The acceleration time histories and spectral accelerations at different burial depths are depicted in Figures 15 and 16. It can be seen that although the peak accelerations at the same burial depth of the liquefiable site are close to each other using two numerical approaches, the results from OpenSees using the fully coupled nonlinear effective stress analysis approach consists of more high frequency contents in the acceleration time histories. The loosely coupled nonlinear effective stress analysis approach in DEEPSOIL fails to properly capture the high-frequency component of ground motions during the seismic propagation in the liquefiable site. Moreover, significant differences between the spectral accelerations and the predominant periods of the soil responses at three different burial depths using two numerical approaches. The spectral acceleration values calculated by OpenSees in short periods are significantly larger than those calculated by DEEPSOIL. Overall, the fully coupled dynamic effective stress analysis method employed in this study can better simulate the development of EPWP in the fully-saturated soil during the earthquake and wave propagation of high-frequency components of seismic motions at a liquefiable site.

Conclusions
In this study, the fully-coupled nonlinear dynamic analyses of a layered liquefiable site were performed by adopting the OpenSees finite element program. The dynamic behavior of the liquefiable soils was modeled by an extensively validated soil constitutive model based on the critical state soil mechanics and bounding surface theory for sandy soil. The dynamic behavior of the non-liquefiable soils was modeled by a von Mises multi-surface model, which can characterize the nonlinear hysteretic behavior of the crust clay soil under cyclic loading. The excess pore water pressure, acceleration, shear stress-strain, and effective stress path in a liquefiable site were investigated in detail during and after the earthquakes. The following conclusions can be drawn from this study.
(1) Under strong earthquake shaking, the loose sand can be liquefied quickly, accompanying the EPWP build-up. Numerical results show that liquefied sand loses the shear strength and behaves like a fluid, and shear waves can't be transmitted in liquefied soil. This phenomenon proves that the PM4Sand constitutive model is capable of describing the post-liquefaction characteristics of soil. (2) The acceleration amplitude of the soil significantly decreases after soil liquefaction, which indicates that the liquefied soil hinders the propagation of seismic waves from the base of the numerical model to the ground surface. The peak accelerations of the ground surface under the El-Centro motion and the Kobe motion attenuate to 64% and 42% of the input ground motion, respectively. (3) The low permeability of crust clay prevents the upward dissipation of excess pore water pressure in the liquefiable sand, thus gradually accumulating excess pore water pressure at the interface of the crust clay and loose sand layers. The large postliquefaction flow deformation occurs at the upper layer of the loose sand layer. (4) The total settlement increases by approximately 40% when considering the postearthquake pore pressure dissipation and reconsolidation. (5) Compared with the loosely coupled effective stress analysis method, the fully coupled effective stress analysis method can better simulate the propagation of highfrequency ground shaking in layered saturated soils and the nonlinear seismic response characteristics of liquefiable soil layers.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: All data generated or analyzed during this study are included in this published article.

Conflicts of Interest:
The authors declare no conflict of interest.