Response of soil–pile–superstructure–quay wall system to lateral displacement under horizontal and vertical earthquake excitations

The excess pore water pressure generated under the combined vertical and horizontal earthquake excitations is generally greater than that under the unidirectional excitation, which may further affect the dynamic response of soil-pile-superstructure-quay wall (SPSQ) system. This study conducted a coupled effective-stress analysis for the response of SPSQ system to liquefaction-induced lateral displacement under horizontal and vertical earthquake excitations. To accurately describe the lateral displacement behind the quay wall, the liquefaction behavior of soils was described by a modified generalized plasticity model with parameters properly determined from laboratory tests. The numerical model of the SPSQ system was validated against the previously reported centrifuge model tests. The results show that the earthquake excitation with abundant low-frequency component causes greater residual lateral displacements of pile cap (ux,RC) and quay wall (ux,RQ) than that with abundant high-frequency component. The ux,RC and ux,RQ tend to increase as the peak vertical acceleration increases for both vertical and inclined piles. Increasing the fixity of the quay wall significantly reduces the influence of vertical ground motion on the ux,RC and ux,RQ. The polarity reversal of the vertical ground motion also has a significant effect on the ux,RC and ux,RQ and should be considered in the analyses. In addition, the lateral displacement of liquefied ground causes an amount of rotation of the cap of the inclined pile, probably leading to the potential shear failure at the head of piles. Special consideration should be given to the earthquake-induced shear force at deeper location because the shear force could be significantly greater than that at the head of both inclined and vertical piles in certain cases.


Introduction
Liquefaction-induced lateral displacement of ground was frequently observed at the backfill of the quay wall in the recent earthquakes (Wang et al. 2021;Zhou et al. 2020;Zhuang et al. 2018;Bray and Dashti 2014;Liu and Tesfamariam 2012;Bhattacharya and Madabhushi 2008). The phenomenon referred to as "liquefaction-induced lateral displacement" during earthquakes is commonly defined as the seismic lateral ground deformation that are associated with liquefaction and could cause collapse or failure of slopes and infrastructures (Towhata 2008;Zhang et al. 2004;Ledezma and Bray 2010;Zhuang et al. 2020). Ishihara et al. (1997) reported that the liquefaction-induced lateral displacement of ground caused the instability of the quay walls during the 1995 Kobe earthquake. The same phenomenon was observed during the 1999 Chi-Chi earthquake (Sumer et al. 2002) and the 1999 Kocaeli earthquake (Chen et al. 2000). Moreover, one significant insight gained from 1995 Kobe earthquake is that the lateral displacement is a primary reason for the failure of piles (Ishihara et al. 1997). Thus, the liquefaction-induced lateral displacement of ground poses a serious threat to the safety of pile-supported structures behind the quay wall. This further necessitates adequate study on the various factors affecting the response of soilpile-superstructure-quay wall (SPSQ) system to the lateral displacement during earthquakes Bhattacharya 2014, 2016;Wang et al. 2019), aiming at improving the seismic design of pile-supported structures near the quay wall. Plenty of physical model tests have been carried out on the SPSQ system subjected to liquefaction-induced lateral displacement. These model tests investigated various influencing factors, including the pile-to-wall distance, thickness of overlaying non-liquefiable layer, shaking direction, the mass of the pile cap, input motion, wall type (i.e., either sheetpile or gravity-type quay wall), and the pile rake angle (Sato 1994;Horikoshi et al. 1998;Sato et al. 2001;Tanimoto et al. 2003;Tazoh et al. 2005;Tang et al. 2014;Liu et al. 2017;Su et al. 2016;Motamed et al. 2009Motamed et al. , 2013. It is overall confirmed that the kinematic loading resulted from the lateral displacement is the major cause of large lateral displacement of piles in the SPSQ system (Bouckovalas and Chaloulos;2014). In particular, Horikoshi et al. (1998) found that the input motion significantly affects the lateral displacement of piles, while the effect of frequency content of input motion is not clarified. Tazoh et al. (2008Tazoh et al. ( , 2010 reported that the inclined piles have smaller bending moments but greater axial forces than vertical piles and are more effective than vertical piles in restricting lateral displacement of superstructures subjected to the lateral displacement after the failure of quay wall. However, Gerolymos et al (2008) and Capatti et al (2020) found both the beneficial and detrimental effects of inclined piles subjected to liquefaction-induced lateral displacement, but these were not detailed or mechanistically justified. In addition to the physical model tests on SPSQ systems, a number of experimental model tests have also been performed on the influence of liquefaction-induced lateral displacement of ground on the soil-pile system in the past two decades (Abdoun et al. 2002;Boulanger et al. 2003;Brandenberg et al. 2005;Madabhushi 2008, 2012;Haeri et al. 2012;Ebeido et al. 2018Ebeido et al. , 2019Ueda et al. 2019;Zhang et al. 2020). These studies also contribute to the comprehensive understanding of the seismic behavior and the engineering design of piles subjected to liquefaction-induced lateral displacement.
In parallel to the experimental study, numerical methods have been used to investigate other influential factors on the SPSQ system with lateral displacement, such as internal friction angle, shear modulus, hydraulic conductivity of soils, the elasticity and material nonlinearity of pile, the combined effect of kinematic load and inertial load, and the effectiveness of different mitigation strategies (Cubrinovski et al. 2008;Tasiopoulou et al. 2013;Tang et al. 2015;Li and Motamed 2017;Su et al. 2018;Li et al. 2019;Chatterjee et al. 2019;Zhao et al. 2022). Particularly, Rajeswari and Sarkar (2021) reported that the vertical ground motion has no additional significant influence on the behavior of inclined pile groups embedded in the gently sloping ground which has been encountered with lateral displacement due to horizontal ground motion. However, it has been confirmed that the vertical ground motion could cause significant fluctuation of excess pore water pressure (EPWP) (Yang et al. 2002) and further aggravate either the post-earthquake settlement of liquefiable soils (Tsaparli et al. 2016) or the pile settlement embedded in saturated sand deposits (Xu et al. 2021a, b). In addition, the vertical ground motion was also reported to have a significant influence on the seismic response of bridge-soil-foundation system (Wang et al. 2013;Dehghanpoor et al. 2019Dehghanpoor et al. , 2021. Thus, it is necessary to clarify how the vertical input motion with the different frequency contents affects the dynamic response of the SPSQ system in lateral displacement of liquefied ground. However, limited work has been carried out on this issue. In this study, a coupled dynamic effective-stress analysis model was used to simulate the centrifuge tests for the SPSQ system affected by liquefaction-induced lateral displacement of ground, i.e., soil flow in original papers of Tazoh et al. (2008Tazoh et al. ( , 2010. The generalized plasticity model (Zienkiewicz et al. 1999) was used to describe the behavior of liquefiable soils. The simulation results were compared with the test results, including the lateral displacements of pile cap and quay wall, EPWP and acceleration at the typical soil layer, and bending moment and axial force of pile. In addition, the effects of frequency content of input motion, vertical excitation component, and pile rake angle on the seismic response of SPSQ system were discussed in detail.
2 Centrifuge modeling of the SPSQ system Tazoh et al. (2005Tazoh et al. ( , 2008Tazoh et al. ( , 2010 conducted a series of centrifuge model tests to study the seismic behavior of SPSQ system. The Case 12 reported by Tazoh et al. (2010), with the scale factor N of 30, was selected as a benchmark model for the numerical verification in this study. As shown in Fig. 1, the soil container in Case 12 had two parallel test models: Model A and Model B, in which the pile groups with supported superstructures were installed behind the quay wall. The difference between the two parallel models was the length of the quay wall: the tips of the quay-wall are attached and unattached to the base of the soil container for Model A and Model B, respectively. The embedded depth of the quay wall is a primary concern in the engineering practice. Normally, the degree of fixity of the quay wall depends on the relative rigidity between the quay wall and soil (Broms 1964). For unfavorable geological areas, the embedded depth should be sufficiently large to prevent against the deformation of quay wall. In such condition, the displacement is minimal at the tip of the quay wall (Kamal et al. 2020). To experimentally simulate large embedded depth of quay wall, the tip of quay wall was fixed in Model A for simplicity. The tips of the group piles were fixed in both of the two models.
There were four soil layers in the container. The first layer was above the ground water level and the second layer was liquefiable. Both of these two layers consisted of No. 8 Silica sand (D r = 50%). It was reported that No. 8 Silica sand has e max = 1.333 and e min = 0.703, where e max and e min are the maximum and minimum void ratios, respectively (Cai et al. 2002). The third and the fourth layers were non-liquefiable and composed of Toyoura sand 1 3 (D r = 90%) and No. 3 Silica sand (D r = 90%), respectively. The thicknesses of four layers were 2.1 m, 3.6 m, 2.4 m, and 0.9 m in the prototype, respectively.
In the centrifuge test, the pile group consisted of four steel pipe piles arranged in a 2 × 2 square pattern with a pile spacing of 0.8 m in the prototype scale. The pile had a length of 8.1 m, outside diameter of 0.3 m, and wall thickness of 6 mm. The cap of pile had a dimension of 1.5 m × 1.5 m × 0.9 m. The superstructure had a dimension of 1.5 m × 1.5 m × 1.2 m and was placed on four columns with 1.8 m in length. The superstructure was modeled by a lumped mass in the simulation. The details of the model setup can be found in the literature (Tazoh et al. 2005(Tazoh et al. , 2008(Tazoh et al. , 2010. In the centrifuge model test of Tazoh et al. (2010), the horizontal input motion was applied at the bottom of the model, with a peak value of 0.26 g, as shown in

Numerical model of the SPSQ system
Both 2D plane-strain model and 3D model have been demonstrated to be able to successfully simulate the behavior of the pile group subjected to the liquefaction-induced lateral displacement of ground (Armstrong et al. 2013;Su et al. 2018). In this study, 2D planestrain model was developed for the SPSQ system because of its efficiency and better convergence in calculating the large liquefaction-induced deformations of soil and pile (Li and Motamed 2017;. Figure 1 shows the finite element meshes for Model A and Model B with the locations of typical sensors, including accelerometers (i.e. A-AG2, B-AG2), displacement sensors for pile cap (i.e. A-DC, B-DC) and quay walls (i.e., A-DQ, B-DQ), pore pressure transducers (i.e., A-PP2, B-PP2), and strain gauges (i.e., A-BS1, A-BS4, B-BS1, B-BS4). The soil, the superstructure, and the pile cap were defined by eight-node quadrilateral element. The cap and superstructure were modeled as elastic material with Young's modulus of 206 GPa and Poisson's ratio of 0.3. The pile, connecting column, and quay wall were modeled as elastic Euler-Bernoulli beam with the material parameters given in Table 1. The interaction of the pile cap with the soil was modelled by merging the nodes at the interface between them. All the beam element nodes were merged with the nodes of quadrilateral elements of soils, which means that besides the displacements, the excess pore water pressure also varied continuously at the position of the pile (Xu et al. 2013b(Xu et al. , 2021a. For the simulation of soil-structure (i.e., quay wall, pile, and pile cap) interface, thin-layer isoparametric element can be employed. This method has been demonstrated to be effective for simulating seismic behaviors of quay walls (Alyami et al. 2009), underground structures (Liu and Song 2005), and piles (Xu et al. 2021a, b) in liquefiable soils. When the liquefaction is initiated, the soil reaction surrounding the pile is definitely decreased (Xu et al. 2021a, b) and the total horizontal stress acting on the quay wall is slightly increased (Alyami et al. 2009). Though the Euler-Bernoulli beam has the degrees of freedom that are incompatible with the 2D plane-strain element of soil, the kinematics of pile and soil at the interface can still be reasonably simulated (Potts et al. 2001), which is due to the fact that the pile rotation has much less significant effect on soils compared with the pile deflection (Potts et al. 2001).
The following equations for estimating the equivalent stiffness of the pile and the connecting column in 2D plane-strain model were recommended (Xu et al. 2013b): (1) A ps = nA 0 B where A and I represent the area and inertia moment of the cross section of piles, respectively. The subscripts "ps" and "0" indicate the values in the 2D plane-strain model and in prototype, respectively. B and n are the pile cap dimension and the total pile number in a row in the out-of-plane direction, respectively. Table 1 gives A ps and I ps for the pile and connecting column in this study. The estimated value of I ps was 1.8 × 10 -5 m 4 for quaywall; however, this value was modified to obtain comparable lateral displacements of pile cap and quay wall with respect to the experimental measurement (see Table 1). The effect of the modification of I ps of quay wall is discussed in Sect. 3.5. Based on the numerical results, the internal force (F) of individual group pile were given by (Xu et al. 2021a, b): where F is a representative of each of bending moment (M), shear force (S), and axial load (N) of individual group pile; F ps represents the internal force (i.e., M ps , S ps , and N ps ) of group pile in the 2D plane-strain model.
The precision of the dynamic finite element method can be guaranteed when the mesh size Δ is not larger than a critical value of c s ∕10f cut , where c s is the wave velocity and f cut is the cutoff frequency (Bao et al. 2020). The average velocity of the shear wave was estimated to be 275 m/s and f cut is taken as 30 Hz to consider the effect of frequency content in the next section, resulting Δ ≤ 0.92 m. The shear wave velocity was estimated from the shear modulus that can be determined using the relative density of sand (Xu et al. 2013a). In this study, the mesh size ranged from 0.25 to 0.63 m and thus met the requirement of accuracy.

Soil constitutive modeling
The generalized plasticity model, i.e., Pastor-Zienkiewicz III (PZ III) model (Pastor et al. 1990), with some modifications by Cai et al. (2002) was adopted for the soils. The modified PZ III model, which has been demonstrated to be effective in calculating the liquefaction-induced lateral displacement of ground (Alyami et al. 2009), was used to characterize the liquefaction behavior of sandy soils. The stress increment and strain increment vectors (dσ and dε) are related through an elastoplastic matrix D ep : where D e is the elastic matrix, H L/U is the plastic modulus, n is the loading direction vectors, and n gL/U is the plastic flow direction vector. The subscripts L and U indicate loading and unloading, respectively.
The shear and bulk modulus (G es and K ev ) are given by: (2) where G es0 , K ev0 , m s , and m v are model parameters, p' is the mean effective stress, and P a is atmospheric pressure. n gL/U and n are expressed as where d g is the dilatancy and d f is a loading scalar. d g and d f are related to the stress ratio (η): where α g , α f , M gc , M fc , and C are model parameters, and Θ is Lode angle. The loading plastic modulus H L is expressed as where H 0 , β 0 , β 1 , and γ are model parameters, η f is defined as a stress ratio which acts as the limit between admissible internal states and impossible external states, and ξ is the accumulated plastic deviatoric strain. The unloading plastic modulus H U is given by where H U0 and γ U are model parameters. Table 2 presents the model parameters for the four soil layers in the container. The model parameters for No. 8 Silica sand (D r = 50%) were calibrated from monotonic drained and cyclic undrained triaxial tests reported by the authors (Cai et al. 2002;2007). The calculation process was completed using UWLC element test program to obtain stress-strain response curves from rate-type and incremental constitutive equations during generalized loading (Xu et al. 2013b(Xu et al. , 2019, where UWLC is a fully coupled dynamic effective stress finite element (FE) analysis software program and is described in the next section. Specifically, the parameters (G es0 , K ev0 , C, m s , m v , M gc , α g , M fc , α f , H 0 , β 0 , β 1 , γ) were calibrated from the monotonic drained triaxial test data, as shown in Fig. 3, where q is the deviatoric stress, ε v is the volumetric strain, and ε a is the axial strain. The unloading modulus parameters (H U0 , γ U ) were calibrated from the experimental CSR-N L curves of No. 8 Silica sand (D r = 50%) in Fig. 4, where CSR is the cyclic stress ratio and defined as CSR = d 2 � c , d is the amplitude of sinusoidal cyclic axial loading, ′ c is the effective confining pressure, N L is the number of loading cycles to cause liquefaction. The liquefaction criterion normally uses either the double amplitude (DA) axial strain or the ratio of EPWP to ′ c for sandy The CSR-N L curve for No. 8 Silica sand (D r = 50%; data from Cai et al. (2002)) and Toyoura sand (D r = 90%; data from Yamamoto et al. (2009)) soils. For No. 8 Silica sand, the occurrence of 5% double amplitude (DA) axial strain was adopted as the liquefaction criterion in the experiment (Cai et al. 2002). For Toyoura sand (D r = 90%), the model parameters were calibrated using the data of monotonic and cyclic undrained triaxial tests reported by Yamamoto et al (2009) (see Figs. 4 and 5), where the liquefaction criterion is defined as the condition when the ratio of EPWP to ′ c reaches 1.0 in the experiment.
The model parameters of No. 3 Silica sand were taken as those used for Toyoura sand because the experimental data of No. 3 Silica sand were not available, assuming that both No. 3 Silica sand and Toyoura sand have the same liquefaction strength. This assumption was reasonable because both No. 3 Silica sand and Toyoura sand had the same relative density which dominates the liquefaction strength (Ishihara 1977).

u-p formulations
The fully coupled dynamic effective stress finite element (FE) analysis software program UWLC (Forum 8 Co. Ltd. 2005;Xu et al. 2013bXu et al. , 2019Xu et al. , 2020Xu et al. , 2021Cai 2020) was used for numerical analyses. UWLC stands for the initial of the family name of four developers, i.e., Ugai, Wakai, Li, and Cai, of the program in Gunma university. The program consists of three parts: (1) a program used to determine the parameters of constitutive laws through laboratory tests, (2) a program to calculate the initial stress before the earthquake occurs, and (3) a fully coupled dynamic effective stress analysis program to calculate the seismic responses. In the dynamic analysis, the solid and fluid phases of soils were governed by the u-p formulations of Biot's dynamic consolidation theory (Biot 1956;Zienkiewicz et al. 1999): where M is the mass matrix of soils, u, ̇ , and ̈ are the displacement, velocity, and acceleration vectors of soil skeletons, respectively, C is Rayleigh damping matrix, K is the stiffness matrix of soils, Q is the coupled matrix, p is the pore water pressure vector, ̇ is the rate of pore water pressure vector, f u is the external load vector, H is the seepage matrix, S is the compression matrix, and f p is the external load vector for pore water. For the earthquake analysis, seismic input motions were transferred to the external load vector of each node of elements. Accordingly, the inertial excitation of each calculation domain was considered in numerical analyses.
M, K, Q, S, and H are expressed as (Zienkiewicz et al. 1999) where N u and N p are the shape function matrices for solid and fluid phase, respectively, n is porosity, ρ s and ρ f are the density of solid phase and fluid phase, respectively. B is the strain matrix, D is the matrix evaluated from constitutive law, m is defined as [1,1,1,0,0,0] T , k is the permeability matrix, 1∕Q * = n∕K f + ( − n)∕K s , K f is the bulk modulus of fluid phase, K s is the average bulk modulus of the soil skeleton, = 1 − K f ∕K s . For soils, α = 1 is simply used because solid grain is incompressible relative to pore water. For the soils below the ground water level, the quadrilateral elements had degrees of freedom of both displacement and pore pressure. Otherwise, the quadrilateral elements were set to have only displacement degree of freedom.
The Rayleigh damping matrix is given by The coefficients α and β are coefficients related to damping ratio ξ and angular frequency ω: These two coefficients are calculated by selecting a damping ratio ξ and two frequencies (f 1 and f 2 ) or periods (T 1 and T 2 ) outside which damping is larger than the damping ratio (Clough and Penzien 1993). According to engineering and computation experience (Wakai and Ugai 2004;Xu et al. 2013b), α and β calculated to be 0.172 and 0.00174, respectively, using T 1 = 0.2 s and T 2 = 2.0 s with constant damping ratio of 0.03. This is because most dominant period of site response varies between these two periods (Wakai and Ugai 2004). According to the computation experience, the variation of viscous damping is small between these two periods. Moreover, the viscous damping could have significant effect on the response of the system when the soils are in the elastic state; however, it turns to have much smaller effect when the soils reach the plastic state. This is because the material damping could be much larger than the viscous damping in such condition, e.g., the material damping might reach ~ 20% when the strain of soils reaches ~ 0.1% (Xu et al. 2016;Zhao et al. 2018). To this end, the Newmark-beta method (Newmark 1959) was used for solving the dynamic coupled hydro-mechanical equations. The Newton-Raphson procedure was used to implicitly solve the nonlinear system at each time step which is taken as 0.001 s in this study (Zienkiewicz et al. 1999).

Boundary conditions
The numerical analysis was performed in two stages in UWLC: the initial stress analysis and the subsequent earthquake dynamic analysis. The initial stress analysis can provide geostatic conditions for the dynamic analysis. In the initial stress analysis, the nodes had zero horizontal motion but were allowed to move vertically at two lateral sides of the model. The nodes at the base of the model had zero displacement in any direction.
The issue of reflecting waves at the boundary is important and challenging, especially in the shaking table model test. Therefore, laminar shear box is commonly used in the shaking table tests (Tazoh et al. 2005(Tazoh et al. , 2008(Tazoh et al. , 2010. In general, only the horizontal seismic excitation is input on the base to simulate the vertical propagation of shear waves, and the shear soil box, accordingly, is able to allow the soils move horizontally freely at two lateral sides. Therefore, the reflection effect of shear wave at the lateral boundaries is weak in this case. At the following stage of earthquake dynamic analysis, the shear beam boundary condition was applied to two lateral sides in the horizontal direction (Su et al. 2020). Thus, the nodes in the proposed numerical model in this study were allowed to move horizontally at two lateral sides of the model, which was also consistent with the experimental setup. Also, the case with a rigid boundary model container was also analyzed in the verification process, imaginably considering that two lateral sides of the model were horizontally constrained. The results showed that the effect of boundary conditions is insignificant in these two models for two reasons: (1) the wave reflection is very weak in soft soils, especially in liquefied soils; (2) the boundary is sufficiently far from structures and has insignificant effect on the lateral displacement of structures, as pointed out in the literatures (Alyami et al. 2009;Liu et al. 2020). Moreover, both the horizontal and vertical input motions can be applied at the bottom of the model. In the simulation of centrifuge model tests, the input motion was only applied in the horizontal direction, as shown in Fig. 2. Figure 6 shows the influence of the modification of I ps of quay wall on the time history of quay wall and pile cap. The I ps modification generally decreased the lateral displacements of pile cap and quay wall for Model A, resulting an agreement between the simulation and the measurement (see Fig. 7a). For model B, the I ps modification had insignificant effect on the lateral displacement of pile cap, while it had an abundant increase in the lateral displacement of quay wall near the end of loading (see Fig. 6b). In general, the calculated lateral displacements were generally smaller than the measurements for Model B (see Fig. 7b). This deviation is due to the fact that severe cracks were developed at the surface behind the quay wall of Model B (see Fig. 7c), which was not satisfactorily captured by the finite element analysis. However, the numerical simulations show that instability of the quay wall of Model B near the end of loading, resulting in the residual lateral displacement of the quay wall close to the measurement (see Fig. 7b). The reason for the I ps modification was that the quay wall is also placed as a boundary in the 2D plane-strain model and its bending behavior is significantly influenced by the value of I ps . Moreover, the nonlinearity of pile and quay wall was not considered in this study for two reasons: (1) the yield strength of stainless steel that used for piles and quay wall were not reported (Tazoh et al. 2005(Tazoh et al. , 2008(Tazoh et al. , 2010; (2) Grade 304 stainless steel is the most common stainless steel (Smith Jr et al. 2022) and has a bending yield strength of ~ 581 MPa (Ren et al. 2019) which is larger than the numerical calculations. This might be another reason for modifying I ps . Moreover, To verify the accuracy of the model in predicting the internal forces of pile, the measured bending moments and axial forces of pile were first calculated from the bending (ε M ) and axial (ε N ) strains measured in Tazoh et al. (2010), respectively:

Verification of the model
where E is Young's modulus of pile and R is the radius of pile.
Then, the predicted M and N were obtained from Eq. (12) using the results from numerical simulation. Figure 10a shows that the predicted bending moments were generally greater than the experimental results for Model A. This was due to the slight overestimation of the lateral displacements of pile cap in the simulation (see Fig. 7a). However, a satisfactory agreement between the predicted and measured bending moments was achieved for Model B (see Fig. 10b). Figure 11 indicates that the predicted transient axial forces were generally smaller than the measured values. Nevertheless, the predicted residual axial force was in a good agreement with the experimental result. To this end, accurate simulation for the seismic response of each structure (i.e., pile and quay wall) and soil is a very challenging task, especially for two models using a unique set of input parameters. Specifically for this study, even after the best efforts, a few simulated results are still different to some extent from those measured in the model tests, which is intentionally and honestly reported for future researches. The above comparisons show the overall capability of the proposed 2D plane-strain model for predicting the seismic behavior of the SPSQ system subjected to the lateral displacement of liquefied ground.

Case analyses and discussions
To investigate the influence of vertical input motion on the dynamic response of the SPSQ system, two sets of input motions (i.e., LSST 7 and LSST 12) were selected from the Hualien large-scale seismic test (LSST) data reported by Zeghal et al. (1995). Particularly, LSST 7 has abundant low-frequency components up to 6 Hz for both horizontal and vertical motions, while LSST 12 has abundant high-frequency components up to 15 Hz and 30 Hz for horizontal and vertical motions, respectively (Xu et al. 2021a, b). LSST12 has a peak horizontal acceleration (PHA) of 0.22 g and a ratio (PVA/PHA) of 0.37, where PVA is the peak vertical acceleration. The PHA and PVA/PHA of LSST7 are equal to 0.088 g and 0.38, respectively. The SeismoSignal (SeismoSoft 2011) was employed to make baseline corrections for the input motion. Figure 12 shows the corrected horizontal and vertical input motions. For liquefaction studies, the PHA is a primary variable to calculate the equivalent cyclic stress ratio (Seed and Idriss 1971). Thus, the horizontal component of LSST 7 was multiplied by a scale factor so that it has the same PHA with that of LSST12. Accordingly, the impact of the frequency content of input motion can be highlighted by comparing the results under the excitation of both LSST 7 and LSST 12. To isolate the impact of the vertical ground motion, the PVA/PHA was increased up to 0.67. The influence of polarity of the vertical component of ground motion was also investigated. Moreover, the cases with inclined piles were schematically plotted in Fig. 13 and were used to study the effect of the pile rake angle. In this study, θ = 10° was adopted because it is commonly used in the engineering practice (Escoffier 2012;González et al. 2020).
Both Model A and Model B were used as the benchmark model for the parametric studies. Table 3 gives a total of 21 cases in the parametric studies. The case ID represents the name of the model, the type of input motion, PHA, and PVA/PHA in succession; for example, B_P10_L7_0.22_0.67 represents the case with Model B and θ = 10° under the excitation of LSST7 that has PHA = 0.22 g and PVA/PHA = 0.67. In addition, the polarity of vertical ground motion was reversed in B_P00_L7_0.22_( − 0.67) and B_P00_L12_0.22_ ( − 0.67) to investigate the influence of polarity of the vertical ground motion, where the minus sign represents the reversal.  Figure 14 illustrates the influence of frequency content of input motion on the time histories of lateral displacements of pile cap (u x,C ) and quay wall (u x,Q ) for Model B. The earthquake excitation with abundant low-frequency component (i.e. LSST7 in Fig. 12) caused greater residual lateral displacements of pile cap (u x,RC ) and quay wall (u x,RQ ) than that with abundant high-frequency component (i.e. LSST12 in Fig. 12). This is because the former input motion caused larger area of liquefaction than the latter input motion, as demonstrated in Fig. 15a and b. Figure 15c further shows the comparison of time histories of   Fig. 1) between the two cases. A higher EPWP was induced by the earthquake excitation with abundant low-frequency component in the soils between the pile and quay wall. Similar results that greater EPWP obtained from the motion with abundant low-frequency component was also reported by Panaghi et al. (2019) and Xu et al. (2021a, b). This phenomenon can be explained based on the dynamic coupled hydro-mechanical equations (see Eqs. 22 and 23): (1) According to the Newmark-beta integration scheme, a relatively larger horizontal residual displacement is usually obtained from a low-frequency earthquake excitation with the same PHA, resulting in the larger plastic strain of soil;

Influence of frequency content of input motion
(2) The high-frequency earthquake excitation generally induces larger velocity in the computational domain and thus results in larger Rayleigh damping force applied to the system. Figure 16a shows the influence of frequency content of input motion on the residual lateral displacement for both Model A and Model B. Note that the only difference between Model A and Model B was the fixity of the quay wall (see Fig. 1). In general, increasing the fixity of the quay wall significantly decreased both u x,RQ and u x,RC . Moreover, the relative difference between either the u x,RQ or u x,RC under low-frequency and high-frequency earthquake excitations was much smaller for Model A (i.e., A_P00_L7_0.15_0.0 and A_P00_ L12_0.22_0.0) than for Model B (i.e., B_P00_L7_0.15_0.0 and B_P00_L12_0.22_0.0). This result implies that increasing the fixity of the quay wall also significantly reduced the influence of frequency content of input motions. Figure 16b also shows the comparison of the earthquake-induced maximum internal forces of upstream pile in the aforementioned four cases, where M ps,max , S ps,max , and N ps,max are the maximum M ps , S ps , and N ps of the pile during earthquakes, respectively. As expected, the low-frequency earthquake excitation induced much greater M ps,max , S ps,max , and N ps,max than those under the high-frequency earthquake excitation.

Influence of vertical ground motion
The influence of vertical ground motion on the response of SPSQ system was illustrated in Figs. 17,18,19. It is shown that both the u x,RC and u x,RQ increased as PVA/PHA increased in the case of Model B. This is likely because the coupling effect of vertical ground motion and horizontal ground motion led to an additional increase in the EPWP in the backfill of quay wall and in the soils in front of quay wall (see Fig. 18). More significant increase in the EPWP was induced by the vertical motion in front of quay wall (see Fig. 18b), compared with that in the backfill of quay wall (see Fig. 18a). As PVA/PHA increased up to 0.67, the u x,RC and u x,RQ increased respectively by 5% and 17% under the low-frequency earthquake excitation, and they increased respectively by 9% and 37% under the highfrequency earthquake excitation. This indicates that the influence of the vertical ground motion was also dependent on the frequency content of input motion. Specifically, the effect of the vertical input dominated by high-frequency components was more significant on u x,RC and u x,RQ . Thus, the adverse effect of vertical ground motion on the piled structures behind the quay wall should be considered in engineering practice. Figure 17b shows that M ps,max and S ps,max tended to increase as PVA increases because larger PVA induced greater u x,RC (see Fig. 17a). It is also observed that N ps,max increased with the increasing PVA under the low-frequency earthquake excitation but decreased as PVA increases under the high-frequency earthquake excitation. The reason for this result is explained as follows. As shown in Fig. 11, the horizontal excitation can also cause a significant vertical fluctuation on the axial force of piles. The maximum axial force is expected to occur, if the pile cap vibration induced by the vertical ground motion follows the same direction with that caused by the horizontal ground motion. Otherwise, the vertical ground In general, u x, RC , u x,RQ , and N ps,max show an increasing trend with the increasing PVA in the case of Model A (see Fig. 19). However, the increase in the u x,RC and u x,RQ was much smaller for the cases for Model A than for those for Model B, as PVA/PHA increased up to 0.67 (see Figs. 17 and 19). Specifically, u x,RC and u x,RQ approximately increased by 3% and 2% under the low-frequency earthquake excitation, respectively, and they approximately increased by 6% and 8% under the high-frequency earthquake excitation, respectively. Note that the tip of quay wall was fixed in Model A (see Fig. 1). Thus, it can be concluded that increasing the fixity of the quay wall can significantly reduce the influence of vertical ground motion on u x,RC and u x,RQ . Figures 20 and 21 show the effect of polarity of the vertical ground motion on u x,C and u x,Q under low-frequency and high-frequency earthquake excitations, respectively. The polarity reversal of vertical ground motion may either increase or decrease u x,RQ , as shown in Figs. 20a and 21a, respectively. However, the u x,RQ in the case of polarity reversal remained much larger than that in the case when no vertical ground motion was considered (see Fig. 21a). The polarity reversal generally decreased u x,RC in this study (see Figs. 20b and 21b). This is because the vertical component acts combined with the horizontal component and thus affects the inertial force of the superstructure and pile. Thus, the polarity reversal of the vertical ground motion should be considered in the analyses.

Influence of the pile rake angle
In order to investigate the influence of the pile rake angle on the SPSQ system under combined vertical and horizontal ground motions, six additional cases are analyzed in this study (see Fig. 22). Similar to the cases for vertical pile, u x,RC , u x,RQ , and N ps,max also generally increased as PVA increased. As the pile rake angle (θ) increased from 0° to 10°, u x,RC and u x,RQ were decreased by approximately 79% and 48%, respectively. However, the reduction of M ps,max was less significant than that of the displacement (i.e., u x,RC and u x,RQ ), as the decrease percentage was approximately 43% and 24% under the low-frequency and highfrequency earthquake excitations, respectively. It is indicated that the pile rake angle had more significant effect on the lateral deflection than on the bending moment. Figures 23  and 24 show the envelope of the internal forces of piles, where M ps,peak , S ps,peak , and N ps,peak represents the absolute peak value of the time history of M ps , S ps , and N ps at each node of the pile segment, respectively. As shown in Figs. 23a and 24a, the depth of zero moment increased with the increasing of pile rake angle. The positive M ps,max was observed near the interface between the unsaturated layer and the liquefied layer. Moreover, the difference between the positive M ps,max and negative M ps,max was not significant for the inclined pile and was thus much smaller than that for the vertical pile, indicating an advantage of using inclined pile subjected to the lateral displacement of ground. Note that the sign conventions for the beam element with positive nodal displacements in finite element method were referred to Xu et al. (2017). As the pile rake angle (θ) increased from 0° to 10°, the S ps,max of the inclined pile was decreased by approximately 23% under the low-frequency earthquake excitation but was increased by approximately 16% under the high-frequency earthquake excitation (see Fig. 22b). The S ps,max located at the position near the interface between the liquefied layer and the saturated non-liquefiable layer (see Figs. 23b and 24b). At the top of the pile, the M ps,peak decreased with the increasing θ but the S ps,peak increased as the θ increased. Thus, the inclined pile is more likely to have the shear failure at the head than the vertical pile when undergoing large lateral displacements. This result also explains the shear failure developed at the head of some inclined piles during the earthquake (SEAOC 1991;Gerolymos et al. 2008;NILIM and PWRI 2014). For example, the liquefaction-induced lateral displacement of ground caused the shear failure at the head of inclined piles supporting the abutments of Nishikawa bridge during the 2011 Great East Japan Earthquake (NILIM and PWRI 2014). Moreover, special consideration should be given to the earthquake-induced shear force at deeper location because the shear force could be significantly greater than that at the head of both inclined and vertical piles in certain cases (see Figs. 23b and 24b). In addition, N ps,max generally increased with the increasing θ (see Figs. 22,23c,and 24c), which was consistent with the findings reported in the previous study (Tazoh et al. 2008(Tazoh et al. , 2010Gerolymos et al. 2008;Capatti et al. 2020). Figure 25 shows the residual deformation diagrams of the SPSQ system after the earthquake in two typical cases (i.e., B_P00_L7_0.22_0.37 and B_P10_L7_0.22_0.37). Note that the deformation of pile group with a scale factor of 10 was also given in Fig. 25. It is found that the deformation of the inclined pile group was different from that of the vertical pile group subjected to the lateral displacement of liquefied ground. A schematic diagram was superimposed to further illustrate the deformation modes of the vertical and inclined pile groups. In the case of vertical pile group, both the upstream and downstream piles deflected equally. However, the lateral displacement caused an amount of rotation of the pile cap due to different deflection patterns observed from the inclined upstream and downstream group piles.

Conclusions
A 2D fully coupled dynamic effective-stress analysis model was developed to investigate the dynamic response of the SPSQ system with liquefaction-induced lateral displacement of ground under horizontal and vertical earthquake excitations. After the validation of the proposed numerical model against the centrifuge model tests, various influential factors were investigated based on the benchmark model. Some important implications obtained from this study are listed as follow: (1) At the same PHA, the low-frequency earthquake excitation could induce greater lateral displacements of pile cap and quay wall than the high-frequency earthquake excitation.
(2) The residual lateral displacements of pile cap and quay wall generally increase with the increasing PVA for both vertical and inclined piles. Increasing the fixity of the quay wall could significantly reduce the influence of vertical ground motion on the residual lateral displacements of pile cap and quay wall. The polarity reversal of vertical ground Distribution of earthquake-induced peak internal forces of vertical and inclined upstream piles in two typical cases under the high -frequency earthquake excitation: a M ps,peak ; b S ps,peak ; and c N ps,peak motion has significant effect on the lateral displacements of pile cap and quay wall and should be considered in the analyses. (3) The inclined pile group can significantly reduce the lateral displacement of pile cap with respect to the vertical pile group, and its deformation mode is different from the vertical pile group that undergoes the lateral displacement of liquefied ground. (4) Special consideration should be given to the earthquake-induced shear force at deeper location because the shear force could be significantly greater than that at the head of both inclined and vertical piles in certain cases.