Seismic behaviour of piles in non-liquefiable and liquefiable soil

This paper investigates the nonlinear soil–pile–structure interaction employing three-dimensional nonlinear finite element models verified with the results of large-scale shaking table tests of model pile groups-superstructure systems. The responses of piles in both liquefiable and non-liquefiable soil sites to ground motion with varying intensities were evaluated considering both kinematic and inertial interaction. The calculated piles and soil responses agreed well with the responses measured during the shaking events. The numerical models correctly predicted the different pile deformation modes that were exhibited in the experiments. The finite element analysis was then employed to perform a parametric study to evaluate the kinematic and inertial effects on the piles' response, considering different ground motion Intensity and piles characteristics. It was found that the bending moment of piles in the liquefiable site increases significantly, compared to the non-liquefiable site, due to the loss of lateral support of the liquified soil, and the maximum bending moment occurs at the interface between the loose and dense sand layers. The inertial interaction contributes the most to the bending moments at the pile top and the interface between the top clay and liquefied loose sand layers. For piles with a larger diameter, the bending moment due to kinematic interaction increases significantly, and the bending moment distribution corresponds to short (rigid) pile behaviour. In addition, the piles at the saturated site displace laterally as a rigid body during strong ground motions because the pile base loses the lateral support due to the soil liquefaction. Finally, the kinematic interaction effect becomes more significant for piles with higher elastic modulus.


Introduction
Analysis of soil-pile-structure interaction (SPSI) is essential in designing significant structures, such as high-rise buildings, bridges, and infrastructure supported by deep foundations installed in loose sand. In regions characterized by significant seismic hazard, the pile foundation could experience severe damage during large earthquakes, which could lead to catastrophic failure of the supporting structure (Cubrinovski et al. 2016;Cubrinovski et al. 2014;Dash et al. 2009;Finn and Fujita 2002;Palermo et al. 2011;Wotherspoon et al. 2012). Sites with saturated loose sand layers are susceptible to liquefaction due to the build-up of pore water pressure (PWP). The potential of liquefaction and consequent damage have attracted significant research efforts to examine the complex behaviour of SPSI in liquefiable soil (Holzer et al. 1989;Kramer 1996). Moreover, if the soil layers are sloped, the ground could move laterally, causing lateral spreading, which adds supplementary lateral forces on piles (Abdoun and Dobry 2002;Ebeido et al. 2019;Ishihara and Cubrinovski 1998;Su et al. 2016). Failures of pile foundations and the supporting structures during recent earthquakes (e.g., 1989 Loma Prieta earthquake, 1995 Kobe earthquake, and 1999 Chi-Chi earthquake) underscored the importance of understanding and evaluating seismic SPSI effects (Boulanger et al. 1995;Chu et al. 2006;Su et al. 2020;Sugimura et al. 1995).
The difficulty of analyzing the SPSI in liquefiable soil is related to the involvement of PWP, the inertial force from the interaction between the superstructure mass with the pilessoil system, and the kinematic interaction between the pile and soil. The superstructure will generate an inertial force when it oscillates, which is transmitted eventually as lateral force and bending moment on the pile cap. The inertial force subjects the piles to additional lateral and axial loads, increasing the bending moment on the pile shaft. Excess pore water pressure (EPWP) is generated during the shaking, which reduces the axial pile capacity, and can lead to an excessive post-earthquake settlement. The loss of soil resistance along a section of the pile shaft due to liquefaction can cause buckling and axial instability of the pile and deterioration in its bending stiffness, which may initiate plastic yielding (Bhattacharya and Bolton 2004b;Boulanger et al. 1999;Veletsos and Meek 1974). Pile buckling is a common failure mechanism in liquefiable soil owing to the loss of support from the surrounding soil, eventually causing a plastic hinge in a pile. The buckling failure mechanism can result from installation errors or weak support of soft soil Bhattacharya and Bolton 2004a, b;Bhattacharya et al. 2005;Bhattacharya and Madabhushi 2008).
Several researchers scrutinized the response of the SPSI using physical and numerical modelling techniques. Model tests for evaluating the SPSI, such as centrifuge and 1 g shaking table tests, provide valuable insights into different aspects of seismic SPSI and generate a propitious database that can be utilized for verifying advanced numerical models to investigate the complex phenomena further. In scaled modelling experiments, typical measurements include the responses of soil, piles, and superstructures (Abdoun et al. 2003;Kagawa et al. 2004). Yasuda et al. (2000) conducted large-scale shake table tests to investigate the behaviour of piles in liquefiable soil, considering the potential of ground flow. Tokimatsu et al. (2005) studied the effect of inertial and kinematic forces on pile stresses based on the results of an extensive shaking table testing program for piles installed in both dry and saturated sand.  performed large-scale shake table tests on two flexible piles and two stiff piles subjected to ground displacement. The study showed that the key parameters that controlled the pile response were the lateral ground displacement, ultimate pressure from the surface layer, and stiffness reduction in the liquefied soil. The behaviour of a low-cap pile group in liquefied dense sand was examined using shake table tests . Ebeido et al. (2019) performed large-scale shake table experiments to study the response of singles pile and pile groups to liquefaction-induced lateral spreading.
This experimental database provides valuable information for calibrating advanced numerical models for high-performance computational simulation tools. Several approaches model the SPSI, such as beam on the nonlinear Winkler foundation (BNWF), elastic beam theory, and numerical modelling. The BNWF model is the simplest and is most commonly used, which estimates the pile response to lateral ground deformation due to liquefaction by simulating the piles as beam elements and the soil as longitudinal and transversal springs. These simplified models can capture the pile bending moments and reduce the lateral soil resistance due to liquefaction. They incorporated empirical p-y curves to describe the relationship between load and deformation (Ashour and Ardalan 2012;Cubrinovski and Ishihara 2004;Liyanapathirana and Poulos 2005;O'Rourke et al. 1994;Su et al. 2016), where the spring stiffness can allow non-proportional relations between the load and the displacement. Also, Valsamis et al. (2012) conducted a study to explore the SPSI using the elastic beam theory, where the liquefied soil stiffness was ignored to account for liquefaction. The simplicity of this method allowed examining the SPSI using the elastic equations of the deflection and the bending moment, and it produced a general overview of the mechanism. However, these methods approximated the reality of the SPSI response, and a sophisticated simulation was essential.
Recent research efforts adopted complex 3D finite element models to investigate the nonlinear soil behaviour by simulating the whole SPSI system to properly account for the interaction between the different soil layers with different stiffnesses, piles and superstructure. However, this analysis requires complex mathematical and computational modelling efforts, which require extensive specialized computational tools (Kampitsis et al. 2015). Meanwhile, recent developments in numerical simulation, including sophisticated and realistic constitutive models for soil behaviour, provide powerful tools for analyzing the seismic SPSI. Consequently, several researchers investigated the pile behaviour in liquefied soil employing numerical models (Chae et al. 2004;Elgamal and Lu 2009;Maheshwari et al. 2004;Yang and Jeremić 2003;Hussein and El Naggar 2021a). Also, numerical studies overcome the disadvantages of experimental simulation such as scaling and hence can be used to validate experimental observations and study different aspects of SPSI further considering a wide range of soil and pile parameters (Cheng and Jeremić 2009;Hussein and El Naggar 2021b). Cubrinovski et al. (2008) studied the effect of non-liquefied crust layer on the bending moment of piles by performing three-dimension effective stress analysis. Assimaki and Shafieezadeh (2013) performed nonlinear finite element models to investigate the SPSI in liquefied soil-induced lateral spreading and validated the numerical model with physical tests. He et al. (2017) calibrated a numerical model employing results of large-scale shaking table experiments to investigate the effect of soil permeability on triggering liquefaction and the liquefaction-induced lateral force. Su et al. (2020) simulated the response of the SPSI system behind a quay wall in liquefiable soil-induced lateral spreading and calibrated the model using the results obtained from large shake table tests.
The studies reviewed above explored the response of piles in either dry soil or saturated and liquefiable soil individually. None of these studies investigated the SPSI response considering the exact configuration of piles in similar site conditions but with liquefiable versus non-liquefiable soil to allow for direct comparison of the obtained responses and delineating the performance aspects and failure mechanisms for each case. It is imperative to separate the kinematic and inertial effects considering liquefaction conditions. Xu et al. (2020) performed shaking table tests to investigate the seismic response of the SPSI in dry and saturated sands, employing a similar test configuration (soil profile and pile group) to enable direct comparison of the nonlinear response for non-liquefiable and liquefiable soil cases. They applied different shakings including, 0.05 g sine wave beat motion to induce linear soil-pile responses and the Wenchuan earthquake scaled motion of 0.3 g to induce nonlinear soil-pile responses. Due to the ever-existing limitations of experimental programs, these experiments were limited to one soil profile and a limited range of earthquake shaking.
To further extend the observations from this comparative study, the present paper develops rigorous 3D nonlinear finite element models to investigate the dynamic response of the SPSI considered in the shaking table testing program to evaluate the effects of different soil and pile parameters (Xu et al. 2020). First, the paper used the experimental results to calibrate and validate the numerical models and explore the different failure mechanisms. Then, the validated models were employed to implement a parametric study to explore the kinematic and inertial effects on the SPSI. Xu et al. (2020)   a potential disturbance of the soil profile during installation, and all the soil parameters were therefore assumed as average values. The reported soil disturbance in the experiments was a conspicuous illustration of the low values of the measured soil parameters, especially in the test with water. It is worth noting that a sensitivity analysis was conducted to examine the geotechnical soil parameters variability while calibrating the material model to capture a reasonable SPSI response.

Description of shake table experiments
A 2 × 2 pile group was installed in the soil bed, and the pile heads were connected by a 0.25 m thick square cap (0.8 m × 0.8 m), and the piles were spaced at five times the pile diameter. The concrete compressive strength of the pile cap was 27.6 MPa, and the concrete compressive strength of the piles was 17.1 MPa. The pile cap was reinforced with bars 5 mm in diameter spaced at 80 mm across the cap. The piles had diameter, d = 0.1 m and length, L = 1.65 m, with the pile toe located 10 cm above the base. Each pile was reinforced with six longitudinal steel bars with a diameter of 4 mm. The top 1.40 m of the pile was reinforced by 1.2 mm diameter spiral stirrups with a 10 mm pitch, while the rest of the pile had stirrups with a 20 mm pitch. A model superstructure was bolted on top of the pile cap, which comprised two rigid steel blocks weighing 410 kg each, connected by a flexible steel column. Both the liquefiable and non-liquefiable soil experiments had the same configuration for soil layers, pile group, and superstructure, but the soil was saturated in the liquefiable case.

Finite element model
The 3D soil-pile-structure interaction numerical model was established using the Open-Sees Framework (Mazzoni et al. 2006), an open-source computational platform that aids in developing applications to simulate geotechnical structures for static and seismic analysis . The post-processing of these analyses was produced using GID (Shinozuka et al. 2000). . Thus, the element size was adequate as the maximum frequency was well above the dominant frequency of the Wolong motion, and the soil discretization would guarantee proper propagation of the ground motion through the soil elements. A finer mesh could cause redundant elements within the wavelength and possibly numerical problems (Zerwer et al. 2002). The annular zone around the pile was divided into six elements (each 0.05 m) to simulate the pile diameter and the surrounding 10 cm soil interface element, as shown in Fig. 3. Figure 4. shows a 3D view for the soil section discretization. The model had 3933 nodes and 3168 elements.

Defining elements and materials
The soil elements were modelled utilizing Brick u-p Element with 8-node hexahedral linear isoperimetric elements . The element has four degrees of freedom: three for the solid displacement in the three perpendicular directions (u), while the fourth is for the fluid pressure (p). This element simulates the dynamic response of the  solid-fluid fully coupled material, and its ability to capture the dynamic response of solids was verified for soil with permeability less than 10 -3 in the complete frequency range (Zienkiewicz et al. 1980(Zienkiewicz et al. , 1999. The constitutive material model PressureDependMulti-Yield02 (PDMY02) was used to simulate the nonlinear response of the sand layers. The material response is elastic-plastic, and it is sensitive to pressure change. The behaviour of the clay layer is simulated employing the PressureIndependMultiYield elastic-plastic material model, in which the plasticity appears only in the deviatoric stress-strain response; however, 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 Yang and Elgamal 2002;Yang et al. 2008).
The interpretation of the soil model geotechnical parameters (e.g., low stain shear modulus, low strain bulk modulus, soil density, and friction angle) was based on the documented measured parameters from the shaking table test program and correlations from the literature. The sand minimum and maximum voids ratios were estimated using the correlations with its particles' mean diameter and uniformity coefficient (Sarkar et al. 2019). The voids ratio was then used to calculate the low strain shear modulus at each depth through Eq. (1) (Das and Ramana 2011). A sensitivity analysis was implemented to optimize the geotechnical soil parameter selection (i.e. friction angle and the small strain shear modulus), which captures the appropriate SPSI behaviour; Li and Motamed (2015) conducted a similar sensitivity calculation.
where e is the voids ratio, pr is the effective confining pressure in kPa, A = 3230, B = 2.97, and C = 0.5.
The soil parameters for the constitutive models, controlling the liquefaction-induced cyclic shear strain, were selected based on the recommended values of the model ) and correlations with the parameters measured from the shaking table test. The sand with different relative densities (i.e., Dr = 60%, and 90%) was calibrated to approximate the potential cyclic resistance under simple shear loading constraints. The material model constant parameters were modified and adjected to capture the liquefaction triggering response for the sand, and the calibration was based on the field results of standard penetration test (SPT), and cone penetration test (CPT) summarised by Idriss and Youd (1997) and Robertson (1985), respectively. The corrected SPT blow number was correlated with the sand relative density using Eq. (2) (Skempton 1986), and the normalized cone resistance was calculated using Eq. (3) (Idriss and Boulanger 2007).
where N1, 60 is the corrected SPT blows, Dr is the relative density, q c1n is the CPT normalized cone resistance. Figure 5 illustrates the results of the FEA calibration, where the liquefaction was initiated after 15 uniform cycles. The sand with Dr = 60% was calibrated to reach a 4% amplitude shear strain after liquefaction initiation (EPWPR = 1.0), and the cyclic stress ratio (CSR) was 0.215. Finally, the dense sand (Dr = 90%) was calibrated to reach 1% shear strain after 15 uniform cycles, and the CSR was 0.47. The calibrated soil parameters were summarized in Tables 1 and 2. Table 1 provides the parameters of the constitutive models for the different soil layers, while Table 2 presents the parameters used to account for liquefaction in the constitutive model. A thin-layer element was introduced to simulate the soil-pile-structure interaction, with a thickness that ranges from 0.01 to 0.1 m depending on the relative stiffness between the two adjacent elements (Desai et al. 1984). For efficient seismic analysis, the element is modified by incorporating the Rayleigh damping to simulate energy dissipation (Miao et al. 2016). A thin element 0.1 m thick (to account for the high pile-soil relative stiffness) with Rayleigh damping was used surrounding the pile (as shown in Fig. 3) to simulate the pile-soil interface. The shear and normal stiffnesses and shear parameters of the interface element were reduced by a factor of 0.70 as recommended by Ghalibafian (2006). The soil elements within the annular zone around the pile were discretized to 0.05 m, corresponding to the size of the pile elements and the interface thickness (Lu et al. 2011).
The pile and the superstructure were simulated using 3D nonlinear displacement-based beam-column elements, while the cap and rigid links elements were modelled as elastic beam-column elements (Mazzoni et al. 2006). A fibre section was implemented to capture the elastoplastic behaviour of the pile cross-section, including cover, concrete core, longitudinal reinforcement, and spiral stirrups, as shown in Fig. 6a, b. The uniaxialMaterial Concrete02 constitutive material model was used to describe the concrete part (the coreconfined and outer-cover-unconfined) (Kent and Park 1971;Scott et al. 1982), while the Giuffrè-Menegotto-Pinto material model (uniaxialMaterial Steel02) was used to simulate the behaviour of the steel part (Filippou et al. 1983;Giuffrè 1970).
Tables 3 and 4 illustrate the parameters used to describe the concrete and steel models, respectively. Thes superstructure cross-section shown in Fig. 6c was also described by a fibre section, and its material behaviour was simulated utilizing the constitutive Calibration of PDMY02 against field observation (Idriss and Youd 1997) and (Robertson 1985)-( } vo = 100 kPa): a relation between shear stress and number of cycles (Dr = 60%); b relation between shear strain and number of cycles (Dr = 60%); c relation between deviatoric stress confinement pressure (Dr = 60%); d relation between shear stress and shear strain (Dr = 60%); e relation between CSR and SPT blow count N 1,60 ; f Relation between CSR and CPT Normalized cone resistance model axialMaterial Hardening. Rigid beam-column link elements were used to fill the space of the pile within the soil, with elastic stiffness ten thousand times the elastic stiffness of the pile (E Irigid = 10 4 E Ipile ) (Su et al. 2017); the rigid links were normal to the pile cross-section ). The pile cap was modelled with high flexure stiffness, and the total mass of the cap was assigned to the cap centre (i.e., depth of 0.125 m). The rigid link and the pile cap were simulated as elastic material. The pile head had the same degrees of freedom as the pile cap, i.e., the same displacement (Li and Motamed 2017). Table 5 illustrates the properties of the parameters used to describe the pile cap and rigid link elements. The two superstructure masses (410 kg each) were  assigned at 2.0 m and 3.0 m for the bottom and top masses, respectively, above the ground level and top of the pile cap (Xu et al. 2020). Table 6 illustrates the properties of the parameters used to describe the constitutive superstructure model.  1 3

Input motion
The shake table experiments comprised harmonic sine waves with an amplitude of 0.05 g and scaled Wolong ground motion record (Li et al. 2008), as shown in Figs. 7 and 8. The Wolong peak acceleration was scaled to 0.3 g, and the 220-s earthquake motion was scaled to 70 s. The peak acceleration of actual motion was reduced because the ground motion dominant frequency was close to the natural frequency of the liquefied soil test, which posed a risk to the safety and integrity of the experiment setup should higher intensity were to be used. The dominant frequency of the motion was 4.69 Hz, which matched the motion used in the tests. A harmonic wave with a peak amplitude of 0.05 g was also applied to study the linear response. It is worth noting that during the test, the shake table actuator could not produce the same small amplitude for both liquefied and non-liquefied tests, but in the numerical study, the same motion was used for the two different sites.

Boundary conditions and analysis stages
To simulate the physical model experiment, staged analysis was adopted to capture the initial conditions before the seismic loading was applied (Wang et al. 2016). Different boundary conditions were applied to the model at different stages, as explained below.
1. In the first stage, only the soil block was implemented using the u-p brick elements. The soil base nodes were fixed against the movement in the gravity direction only, and the outer base points were fixed against the translation along x and y-directions. This stage represented an elastic stage, and the soil behaviour was assumed to be linear elastic. 2. In the second stage, the state of soil was updated to account for plasticity, and the initial free-field stresses were obtained from a plastic gravity analysis. The base boundary condition was fixed in all directions for all nodes. The 4 th degree of freedom of the soil elements defines the conditions of the pore pressure. Therefore, boundary conditions were specified to constrain the nodes above the water lever (dry nodes) to account for the built-up of the PWP. 3. The third stage comprised installing the structural elements (piles, rigid links, pile cap, and superstructure). The three degrees of freedom of the pile and rigid link nodes and the soil interface elements were tied together. Structure element masses were defined along with the mass of the superstructure; then plastic gravity analysis was performed to calculate the consolidation and settlement due to the pile group installation. This analysis would not capture the effect of the pile installation, such as stress relaxation for drilled piles or soil densification for driving piles, but these are minimal in the current analysis (Wang et al. 2016). The response during gravity steps was evaluated using the Newmark integration method with the time integration parameters of · = 1.50 and · = 1.0. 4. In the seismic analysis stage, a shear beam boundary was applied to the node at both sides of the model at each depth to ensure that the movements of the nodes in the horizontal direction were equal at the same height (Su et al. 2020). The permeability coefficients of the different soil layers were updated to account for the undrained behaviour during the seismic motion. The base motion was implemented to the fixed base of the model using the UniformExcitation command (Mazzoni et al. 2006) with a time step of 0.01 s, and the seismic analysis was performed.
For the dynamic steps, the solution was obtained using the modified Newton-Raphson approach (KrylovNewton) (Carlson and Miller 1998;Mazzoni et al. 2006). The initial tangent stiffness of the system was used for all steps to achieve an energy increment test less than 10 −4 . Finally, the Rayleigh damping was defined by a damping ratio of 2% for the dry site and 3% for the saturated site, with an estimated first mode frequency of 1.99, adopted to simulate the energy dissipation and enhance the stability of the numerical analysis. A low stiffness-proportional coefficient of 0.0006 and 0.0008 was used for the dry and saturated sites, respectively.

Excess pore water pressure build-up
The excess pore water pressure ratio (EPWPR) is expressed as the ratio between the measured variation in pore water pressure, (EPWP), to the initial vertical effective stress } vo ( EPWPR = EPWP } vo ) at a specific depth. When EPWPR reaches one, the soil loses its shear strength and behaves like a liquid. During the shaking table tests, sensors were placed at different depths through the potentially liquefiable soil to monitor the generation and dissipation of the water pressure. Figure 9a, c compares the calculated and measured EPWPR for both the 0.05 g sine beat wave motion and the 0.3 g Wolong ground motion. Figure 9b, d shows a rise in the EPWPR at the free field points and the points near the piles, which can be used to evaluate the soil-pile interaction effect on pore water pressure propagation. Figure 9a shows that the calculated maximum EPWPR was less than 0.4 at all depths, which was slightly higher than the measured EPWPR during the weak 0.05 g sine wave motion. The variation herein could be related to the inability of the actuators to produce the sine wave at exactly 0.05 g. The calculated and measured EPWPR values indicate that the weak shaking did not trigger liquefaction in the loose sand. Figure 9c shows that calculated and measured EPWPR for the Wolong motion reached 1.0 (i.e., liquefaction occurred) at all points up to a depth of 0.9 m. At depth 1.5 m, the calculated EPWPR was 0.67, and the measured value during the experiment was 0.5. The maximum EPWPR of 0.6 at the base of the dense sand layer is attributed to the densification of the sand and the higher overburden pressure. Both experimental and numerical results depicted the triggering of liquefaction in the period of 8 s to 13 s of Wolong motion. However, the experimental results for EPWPR were slightly less than the calculated values at the lower depths after 23 s of the input motion. These observations indicate that the numerical model correctly captured the build-up of PWP and liquefaction. Figure 10 shows a 3D visualization of the EPWP distribution for the Wolong motion at 15.21 s. Figure 10 also shows the soil deformation (scaled up 70 times) after reaching the (a) (b) (c) (d) Fig. 9 Calculated excess pore pressure ratio at saturated test: a free field points for 0.05 g sine beat motion; b near pile points for 0.05 g sine beat motion; c free field points for 0.3 g Wolong ground motion; d near pile points for 0.3 g Wolong ground motion liquefaction state. The calculated EPWP at the free-field point and near the pile are compared in Fig. 9. For the points near the pile in Fig. 9b, d, the behaviour was oscillatory during the PWP build-up, which corresponded to the first 5 s with 0.05 g acceleration and the peak accelerations throughout the Wolong motion between seconds 6 and 25. The oscillatory response is attributed to the vibration of piles, which developed then dissipated the PWP at the pile-soil interface as a gap opened and closed, causing significant increases and decreases in EPWPR, especially during Wolong motion, and it was also observed in the experiments, especially for the top points. However, it was not observed at larger depths because the gap was limited to the surficial layer. Similar behaviour was also reported by Dash and Bhattacharya (2015) and Wang et al. (2016).

Free-field soil acceleration amplification and deformation
Shape acceleration arrays (SAA) were used to measure the deformations and accelerations of the piles and soil at specified points (Xu et al. 2020). The measured accelerations revealed amplification of the input motion as the seismic waves propagated towards the surface. Figure 11 compares the measured and calculated profiles of the acceleration amplification in terms of the acceleration amplification factor (AAF) for the dry and saturated soil tests during the 0.05 g sine beat wave and the 0.3 g Wolong events. For the 0.05 g sine beat wave event, the measured AAF at the surface was 2.2 and 7 for the saturated and dry sites, respectively. Meanwhile, the calculated AAF was 2.3 and 6.4 at the exact locations, as shown in Fig. 11a. The lower AAF value in the saturated site is due to soil softening because of EPWPR and increased damping, which dissipated more energy than the dry site.

3
The experimental and numerical results for the Wolong ground motion presented in Fig. 11b show AAF values of 1.0 for the points in the dense sand; however, the topsoil layers exhibited different responses. The motion amplified through the top two layers for the dry site, and the measured and calculated AAF values were 1.32 and 1.41. On the other hand, the saturated site exhibited de-amplification of the acceleration with AAF of 0.64 and 0.73 for the experimental and numerical results. The de-amplification of the acceleration is attributed to the strength degradation of the liquefied soil. Although the AAF was less than one throughout the liquefied sand and the top crust clay, the AAF increased through the crust clay.
The calculated lateral soil displacements are compared with the measured values in Fig. 12. The measured and calculated lateral displacements increased from bottom to top. For the 0.05 g sine beat wave, the measured lateral displacements increased from 2 to 8 mm for the saturated site and up to 6 mm for the dry site, while the maximum calculated displacement was 5.0 mm for both sites. For the Wolong ground motion, the calculated maximum lateral displacement was 10 mm and 12.3 mm for the saturated and the dry sites, which is less than what was measured in the experiments. The calculated and measured lateral displacements were higher in the dry sand site than the saturated sand site, attributed to the increased energy dissipation of the liquified soil. It is worth noting that the increase of the small strain shear modulus or the friction angle of the liquefiable layer caused a remarkable reduction in the lateral displacement, and this observation was concurrent with the finding of Li and Motamed (2015); while a marginal change was observed with the variation in the dense sand layer. Overall, the large deformations and acceleration amplification in the dry sand tests and the liquefaction in the saturated sand tests were predicted fairly accurately by the developed finite element model. Figure 13 shows the lateral displacement time histories of the dry and saturated sites during Wolong ground motion. The lateral displacement of the saturated site had a longer duration, as shown in Fig. 13b, indicating a more extended vibration period of the soil profile. Figure 14 exhibits the 3D visualization of the lateral deformation (scaled up 70 times) for both the dry and saturated sites at the second 11.71 of the Wolong events. For the dry site, the soil block lateral deformation was distributed along with the soil profile almost uniformly, while the lateral deformation increased at a higher rate (almost concentrated) within the liquified soil layer compared to the dense sand at the bottom and the clay crust at the top.
Several factors could explain the difference between the experiment results and the numerical results (e.g., test setup and execution, test measurements, soil parameter interpretation, numerical model assumptions, mesh size, damping ratio, scaled ground motion, time step, and mass and proportional stiffness coefficient). The reports of a high soil disturbance during installing pile groups and the actuators' inability to produce the sine wave at exactly 0.05 g impose uncertainties in the measured parameters (Xu et al. 2020). These variabilities could provide some explanation for the differences between the experimental and numerical results. However, the experimental and numerical AAF and the lateral

Pile acceleration amplification
The SAA was used to measure the pile's accelerations and deformations during the experiments. Correspondingly, the pile accelerations and deformations were calculated at the exact locations. Figure 29 compares the numerical and experimental pile AAF values during the 0.05 g sine wave and Wolong earthquake ground motions for the dry and the saturated sites. The trends and the AAF values obtained from the numerical model agree well with the experimental results. For the 0.05 g sine wave, the calculated AAF increased from 1.0 at the base to 4.3 and 2.1 at the pile top for the dry and saturated sites, while the maximum experimental values were 3.3 and 2.52 for the dry and saturated sites. Like the soil, the difference between the calculated and measured values could be mainly attributed to the inability of the actuators to produce the sine wave at exactly 0.05 g (the actual peak acceleration was less than 0.05 g). For the Wolong motion, the calculated and experimental AAF in the dry site increased from bottom to top with a maximum value of 1.99 and 1.83 at the pile top. Meanwhile, the AAF was one in the dense soil and decreased through the liquified soil to a minimum calculated and measured values of 0.64 and 0.7 at 1.0 m below surface.
Comparing the soil and pile AAF values during the 0.05 g sine wave (Figs. 11a, 15a), it is observed that the pile AAF was less than that for the soil. On the other hand, the Wolong motion caused more prominent oscillation in a pile compared to the soil in the dry site, while for the saturated site, the soil acceleration was de-amplified more than the pile's acceleration as noted from Figs. 11b, 15b.

Pile lateral displacement and vertical settlement
The lateral pile deflection was significantly affected by the ground motion intensity and the site condition. For the weak motion, the pile's maximum lateral displacement was less than 1.5 mm for both dry and saturated sites. The pile lateral displacement increased during the large Wolong motion as shown in Fig. 16, i.e., 7.37 mm and 5.47 mm for the dry and saturated sites, respectively. The pile group in the saturated site, however, exhibited a more significant rotational response. Also, the pile base lateral displacement was 0.0 mm during the dry tests and was − 0.65 mm during the saturated test, and that indicates the pile group rotated in an almost rigid body mode, as opposed to flexural deflection at the top during the dry test. Figure 17 displays the calculated vertical settlement time history of piles one and two in the dry and the saturated sites for both ground motions. The piles' static settlement was 0.35 mm and 0.5 mm at the dry and saturated sites. During the shaking events, the pile group settled more: a gradual increase in settlement of piles was observed during the 0.05 g sine wave ground motion to a maximum of 1.8 mm and 1.6 mm for the dry and saturated sites, respectively. For the Wolong ground motion, the settlement increased to 3.3 mm for the dry site and 6.1 mm for the saturated site. It is noted that a significant difference in the settlement was observed after the 11.7 s when liquefaction occurred. Again, a longer duration of response cycles of the pile-soil system is observed after the soil is liquefied. It is also noted from Fig. 17d that the piles experienced uneven settlement, which corresponded to the rigid body rotational response of the group in the saturated site. This rigid body rotation of the pile-cap structure system was also observed in the experiments. These results clearly indicate that pile groups exhibit different performance characteristics and failure mechanisms in non-liquefiable and liquefiable sites subjected to the same excitation. It is also noted that the settlement of the pile group during strong shaking events is attributed to the densification of loose sand and the settlement associated with liquefaction.

Pile bending moment
Strain gauges were attached to the pile reinforcement to measure the pile's strains, which were then used to calculate the bending moment (Xu et al. 2020). Correspondingly, the time history of the pile's bending moments and the results are compared with the measured values. The comparison demonstrates that the calculated bending moments were slightly larger than the experimental results. For the 0.05 sine beat wave, the measured maximum bending moments were 40 N m at the saturated site and 12 N m at the dry site, while the calculated values were 52 N m at the saturated site and 18 N m at the dry site. The difference may be attributed to the higher acceleration used in the numerical model than the actuator input motion during the experiments. The interface factor which controls the reduction of the soil shear strength at the interface elements could also affect the calculated bending moment. For the Wolong motion, the measured maximum bending moments were 380 N m at the saturated site and 28 N m at the dry site, while the calculated values were 420 N m at the saturated site and 31 N m at the dry site (i.e., the difference is about 5%). The maximum bending moment along the shaft varied for piles in the dry site, indicating an inflection point, and residual bending moment existed at the end o shaking, as shown in Fig. 18a. The bending moment of the pile in the saturated site was much higher than that in the dry site due to the degradation of the pile's lateral support, but there was no residual moment after the shaking, as shown in Fig. 18b.
Generally, the results indicated that the numerical model captured the pile bending moment characteristics correctly and depicted the different deformation modes for both saturated test-0.05 g motion; c dry test-Wolong earthquake; d saturated test-Wolong earthquake the saturated and dry sites. Nevertheless, the experimental results showed that piles did not have any macroscopic failure because the bending moments were less than the yield moment; there were no signs of plastic hinges, which the numerical model also predicted.

Inertial and kinematic effects on piles response
Piles seismic performance is dominated by the coupling between the inertial forces and the kinematic effects from the ground movement interacting with the pile shaft. This topic was investigated through physical tests, cross-correlation analysis techniques and pseudo-static analysis methods (Abdoun and Dobry 2002;Tokimatsu et al. 2005;Wang et al. 2016). However, complete 3D finite element analysis can clarify the ambiguity of the complex conjugation between the kinematic and inertial forces, especially when the liquefaction phenomenon is involved. Plastic hinges are created through the pile shaft because of the excessive internal bending moment induced by the kinematic and the inertial force coupling. The location of these hinges depends on the intensity of ground motion (kinematic effect), masses on the pile top (inertial effect), and pile configuration. A parametric study was conducted to explore the variation of the kinematic and inertial effects on the pile response.
In the following sections, piles along the initial movement direction are denoted as Leading Piles, while piles in the opposite direction are referred to as Trailing Piles. The results of the bending moments shown in the following graphs are the maximum values through the time histories, and they are normalized by the yield bending moment for the pile cross-section with d = 0.1 m. The calculated bending moment is compared to the results from the shaking table configuration. Table 7 summarized the parameters used in the parametric study for the saturated and dry cases.

Effect of kinematic interaction-excluding superstructure mass
The superstructure was removed to explore the kinematic effect of the ground motion on the pile response, and the response was calculated for the Wolong motion scaled to 0.3 g, 0.4 g,0.5 g and 0.6 g. The analyses were performed for both the saturated and dry sites. Figure 19a, b shows the pile normalized bending moments (PNBM) at the saturated site, excluding the superstructure mass. The PNBM increased for all piles as the PA of input motion increased; however, the bending moments remained below the section yield (a) (b) Fig. 18 Pile bending moment time history for Wenchuan earthquake: a dry test; b saturated test moment capacity (maximum PNBM = 0.89 at the peak 0.6 g acceleration ground motion). The significant increase in bending moment occurred at the intersection of the loose layer with other layers. Abdoun et al. (2003) modelled steel-driven piles using centrifuge tests and reported similar observation, which was attributed to the shear discontinuity effect.
The results revealed that PNBMs were reduced through the liquified layer to almost half (i.e., at a depth of 0.6 m). The pile diameter affects pile stiffness and, correspondingly, the kinematic pile-soil interaction. Thus, analyses were conducted to study the kinematic effect (without superstructure) of piles with d = 0.2 m and d = 0.3 m at the saturated site, as shown in Fig. 20. The piles were spaced at 4d, and the pile longitudinal reinforcement was maintained constant at 0.96% of the concrete section. For the 0.3 g input motion used in the tests, the maximum bending moment of the pile with d = 0.1 m was 420 N m, while the maximum bending moment for the pile with d = 0.2 m was 1273 N m; and for the pile with d = 0.3 m  was 3719 N m. The maximum bending moment increases from the base to a maximum value at the top (i.e., no inflection point), which indicates short (rigid) pile behaviour due to the low slenderness ratio, (L/d) = 8.25 and 5.5 for the piles with d = 0.2 m and 0.3 m, respectively. The maximum bending moment occurred at the pile top and increased exponentially as the pile diameter increased due to the higher section stiffness. For example, the maximum bending moment during the 0.6 g ground motion was 2552 N m for d = 0.2 m and was 8372 N m for d = 0.3 m. Despite the increase of the bending moments, no plastic hinges occurred because of the increase in pile bending capacity (PNBM was 54% and 50.7% for d = 0.2 m and 0.3 m, respectively). Figure 21 shows the lateral displacement of piles with d = 0.1 m, d = 0.2 m, and d = 0.3 m. As expected, the pile lateral displacement increased as the peak acceleration of ground motion increased. For a pile with d = 0.1 m, the pile rotated at a point just above its tip, and the lateral displacement increased to a maximum at the pile head for the ground motions with peak accelerations of 0.3 g, 0.4 g, and 0.5 g. However, the pile displaced laterally in almost rigid body movement for peak acceleration of 0.6 g because the dense sand layer liquefied and lost its shear strength. It was observed that the pile diameter also influenced the lateral displacement. For piles with d = 0.2 m and d = 0.3 m in Fig. 21b, c, due to the increased pile rigidity, it displaced laterally almost as a rigid body with some rotation, resulting in maximum displacement at the pile head. The rigid body movement was particularly pronounced for the strong shaking as the dense soil layer liquified.
Several models were implemented to explore the kinematic effect on the SPSI in the dry site. Like the saturated site, the superstructure mass was removed, and the ground motion was scaled to 1.33, 1.67, and 2 times the 0.3 g Wolong ground motion. Figure 22a demonstrates that the kinematic effect on the bending moment of piles (without superstructure) in the dry site was relatively small. However, the bending moment increased almost proportionally through the piles as the ground motion intensity increase; and the maximum bending moments occurred at the pile top and near the  Figure 22b shows that the pile lateral displacement also increased from zero at the pile base to the maximum value at the pile top, and as expected, was increased as the ground motion acceleration increased.

Effect of kinematic interaction-including superstructure mass
To explore the kinematic effect of the SPSI during the existence of the internal force at the saturated site, the FEMs were updated to account for the superstructure mass used during the shaking table test series. Figure 23 presents the bending moment for piles with d = 0.1 m, where several sections of the pile approached the yield state. The bending moment of pile cross-section at the interface of clay and sand layers exceeded the yield moment (i.e., plastic hinge developed) for both the leading and trailing piles, while the cross-section at the intersection between the loose and dense sand layers approached yield state with a maximum PNBM of 0.8 (cracks initiated). This increase is attributed to the inertial forces and the weak resistance of the clay layer. It is also noted that the bending moment increased consistently through the pile shafts because of the inertial effect. These results demonstrate that when the piles lost the lateral support after reaching the liquefaction state, the bending moment increased significantly, especially along the liquefied part. Two significant parameters define the pile stiffness: the cross-section dimensions and the material elastic modulus. Thus, the effect of pile material elastic modulus on kinematic interaction was investigated in this section by considering the response of steel pipe pile, with the same diameter as the test pile (i.e., 10 cm external diameter) and the wall thickness of 8 mm. The pile elastic modulus was assumed to be 2.1E8 MPa, and the steel unit weight was 78 kN/m 3 . The pipe pile section modulus is 49 cm 3 , and the bending moment capacity is 18E3 N m, assuming a yield strength of 355 MPa.
A series of analyses were performed by varying the peak ground acceleration, and the results of the maximum bending moments through the time history are shown in Fig. 24. The maximum bending moment during the 0.3 g Wolong earthquake was 1098 N m at the leading pile and 1497 N m at the trailing pile, while the maximum bending moment recorded during the shaking table test (as well as the calculated) for the 0.1 m reinforced concrete pile was 420 N m, which represents an increase of 356% due to the pile rigidity.  Figure 24a shows that the bending moment increases with higher ground motion acceleration, and the bending moment for the 0.6 g case is more than twice that for the 0.3 g case. Figure 24b shows that the trailing piles bending moments increased with higher ground motion, but this increase was located throughout the loose sand, and it became prominent at the intersection of loose and dense sand layers due to the shear discontinuity. The minor reduction of the bending moment in the clay soil is attributed to the deterioration of the soil passive resistance with higher ground motion acceleration. The kinematic interaction for the steel pile was evaluated at the dry site. Figure 25 illustrates that the increase of the pile stiffness caused a higher bending moment; the maximum bending moment during the 0.3 g Wolong motion increased by almost 50% over the bending moment for the concrete pile (from 31 to 45.4 N m). It is also found that the bending moment increased by 300% from the 0.3 g motion to the 0.6 g motion. Figure 25 shows that the trailing pile experienced a higher bending moment than the leading pile, which is attributed to the high passive resistance of the soil against the lateral displacement of the piles. Figure 26a, b presents the steel pile lateral displacement at the saturated and the dry sites. As expected, the pile lateral displacement decreased at both the dry and the saturated sites, which was attributed to the increase of the pile stiffness.

Inertial effects
Several failure mechanisms could be attributed to the inertial force, including the development of plastic hinges due to excessive bending moments, buckling failure and the loss of foundation capacity. Figure 27 shows the response of piles in the saturated site considering different superstructure mass defined in terms of mass factor; MF = considered mass/test mass. Figure 27a, b shows that the leading and trailing piles experience larger PNBM as MF increases. The calculated maximum bending moments for MF ≥ 2 exceeded the yield moment and caused plastic hinging at the clay-sand interface. The pile used 88% of its bending capacity at the MF = 1, while it used 136% of its bending capacity at MF = 5 (i.e., a 53% difference). The excessive increase of the bending moment at the pile top is due to the significant increase of the acceleration at the pile top, and the superstructure mass increases the fixation of the pile top. The inertial effect was investigated for the piles at the dry site subjected to the Wolong motion with 0.3 g peak acceleration and varying superstructure mass. Figure 28 shows that the PNBMs were small, and the bending moment increased at the intersection of the layers. The rise of the superstructure mass increased the PNBM from 0.06 at MF = 1 to a The minor effect on the bending moment is attributed to the higher soil lateral subgrade reaction at the dry site, which degrades at the saturated site with higher EPWPR; until the soil liquefies (i.e., EPWPR = 1) and the pile becomes laterally unsupported. The effect of the inertial interaction on the pile's lateral displacements is presented for the saturated and dry sites in Fig. 29. Unlike the kinematic effect, the difference between the pile lateral displacement at the saturated and dry sites was small.
To investigate the effect of the pile diameter on the inertial effect at the saturated site, several models have been conducted using pile with d = 0.20 m. Figure 30 shows the maximum bending moments during the time history at the saturated site supporting different superstructure mass during the 0.3 g Wolong motion. Like the kinematic effect, the pile bending moment increased with a larger pile cross-section and a higher superstructure mass. Nevertheless, there were no signs of plastic hinges since the maximum bending moment was 3021 N m at MF = 5, which is less than the section bending capacity of 4750 N m. The distribution of the pile bending moment displayed in Fig. 30, increased from zero at the pile base to a maximum value at the pile top, indicates a rigid pile response. On the other hand, the bending moment displayed in Fig. 30 shows somewhat flexible pile behaviour, with maximum bending moments at the intersections between the layers. The pile with d = 0.2 m had a maximum PNBM of 41% and 64% for MF = 1 and 5, respectively.

Piles vertical settlement response to inertial forces
During the static analysis, the increase of the superstructure mass caused higher settlement, and the settlement values at time zero in Fig. 31a show that the settlement increased from 0.53 mm at MF = 1 to 2.28 mm at MF = 5 (430% increase). On the other hand, the variation of pile settlement with an increase in superstructure mass changing during the dynamic motion was not significant as the static and kinematic effects mainly governed it. The dense sand layer did not liquefy during the 0.3 g ground acceleration (maximum EPWPR was 0.6 from the FEA results), so it preserved its strong resistance (stiffness) to support the higher masses. The pile group rotation in Fig. 31b also indicates that mass variation has a minor effect. The lateral and vertical deformation response for the pile with d = 0.2 m was similar to the pile with d = 0.1 m. Therefore, the serviceability conditions are significantly controlled by the kinematic force and marginally by the inertial forces.

Summary
Based on the validated FEMs, the performance of the SPSI is found to be controlled by the soil parameters (e.g. the relative density, small strain shear modulus, friction angle, and density), as well as the representation of the interaction between the pile and the adjacent soil. Brandenberg et al. (2013) emphasized the importance of adopting a 3D continuum with proper size elements close to the pile to simulate the interface element, but it could be computationally expensive. Nevertheless, this paper implemented a full 3D FEA with a thin layer interface element between the pile and the soil as proposed by Desai et al. (1984). The obtained numerical results are aligned with findings from several studies. For example, Dash and Bhattacharya (2015) illustrated that the EPWPR decreases near the pile due to the water dissipation at the soil-pile interface, which was observed in the current study. Also, the large period of the soil lateral displacement at the liquefiable site case (shown in Fig. 13)  the calculated pile lateral displacement indicated longer vibration response as shown in Fig. 16. The calculated pile responses revealed different failure modes for piles at the dry and saturated tests. The piles exhibited permanent deformation at the intersection between the clay and the loose sand at the dry site, while the pile group was exposed to a rigid rotation at the saturated site due to the loss of the lateral soil support. Also, at the saturated site, the maximum moment occurred during the highest acceleration cycles, which triggered liquefaction, then the pile bending moment reduced to its initial value, similar to observations reported by Ebeido et al. (2019). The pile maximum bending moments were observed at the intersection between the different layers owing to the shear discontinuity, which was consistent with Abdoun and Dobry (2002); Abdoun et al. (2003) centrifuge test results. All these observations confirm Martin and Chen (2005) suggestion that the relative stiffness between the pile and soil predicts the failure mode.
The capability of the FEA to capture the physical model behaviour was encouraging to explore the kinematic and inertial effects separately. The obtained results illustrated that piles with larger flexural stiffness (EI) experienced larger bending moments. Furthermore, short (rigid) pile behaviour was observed for piles with larger diameters due to the low slenderness ratio, and the pile group experienced rigid body movement during the strong shaking when the dense soil layer liquified. However, this paper revealed a marginal effect of the pile diameter on the lateral displacement, especially on the lateral displacement. Therefore, the findings from this paper confirm the complexity of the SPSI response and underscore the importance of the further investigation.

Conclusion
Several 3D finite element analyses, using the OpenSees platform, were performed to address the dynamic response of the soil-pile group-superstructure interaction through non-liquefiable and liquefiable soil states. The numerical modelling was based on four shaking table experiments at the China Academy of Building Research aiming to distinguish between the behaviour of the soil and piles in different site conditions during different intensity motions. Overall, a reasonable dynamic response agreement between the experiment results and the numerical analysis was achieved for the pile and the soil.
The shaking table test results validated the numerical models. The validated models were then employed to study the kinematic and inertial effects on the soil-pile interaction. A comprehensive parametric study was conducted by varying the pile diameter, pile material, intensity of ground motion, and the superstructure masses. The main findings from the study are summarized below.
(1) As expected, the pile bending moments increase as the peak acceleration of ground motions. The increase at the saturated site is significant relative to the dry side due to the loss of the lateral support of the liquified soil. The maximum bending moment occurs at the interface of the loose and dense sand layers.
(2) The inertial interaction due to the superstructure mass contributes the most to the bending moments at the pile top and the interface of the top clay and loose sand layers. These bending moments increase as the supported mass increases.
(3) Piles with larger diameters experience a significant increase in bending moment due to kinematic interaction, and the bending moment distribution corresponds to short (rigid) pile behaviour due to the low slenderness ratio. (4) As expected, the pile lateral displacement increases as the peak acceleration increases for both the dry and saturated sites. However, the piles displace laterally as a rigid body at the saturated site during strong ground motions as the pile base loses the lateral support due to the soil liquefaction, especially for more enormous diameter piles. (5) The increase of the pile material elastic modulus causes a dramatic increase in the bending moment due to the kinematic interaction, but the lateral displacement decreases due to increased pile stiffness. (6) The variation of dynamic pile settlement with the superstructure mass is marginal given that the dense soil layer at the base does not liquefy.