Numerical investigation on p–y method of group piles under static and dynamic loads

The present study aims at scrutinizing the static and dynamic behavior of pile-soil interaction in the context of p–y method using finite element analysis. For this purpose, a series of static and dynamic analyses are carried out for single pile, 2 × 2 and 3 × 3 group piles in homogeneous clayey soil. In the analysis, the nonlinear behavior of soil is taken into account using a kinematic hardening model, while piles are modeled as elastic. The soil model parameters are calibrated using experimental modulus degradation and damping curves. The pile-soil interface is modeled considering normal and shear behavior to account for separation and sliding between soil and pile elements. In the static analyses, variation of pile moments and displacements with depth and the back-calculated p–y curves are evaluated. The soil resistance is directly obtained by extracting the shear and normal forces of the nodes in a pile at any depth. The p-multipliers for 2 × 2 and 3 × 3 group piles are calculated and compared with that of the experiments. The variation of p-multipliers with depth and lateral displacement is also evaluated. In dynamic analyses, first, site response analyses are carried out and validated against one-dimensional results under different loading frequencies. Infinite elements are applied at the boundaries to provide non-reflective boundaries. Later, single, 2 × 2 and 3 × 3 group pile analyses are executed. The influence of loading amplitude and frequency on the response are investigated using moment-depth and displacement-depth relationships. Dynamic p–y curves are back-calculated and the results are deeply assessed by comparing with static curves.


Introduction
In parallel with the developments in industry and technology in the last century, the diversity of buildings has increased to meet the needs. New building types (high buildings, viaducts, bridges, etc.) are built in city centers for prestige purposes, usually to meet 1 3 commercial demands. These unique structures today are not only indicators of power and prestige but also the pinnacle of engineering. In the selection of foundations of such structures, piled systems, which are a type of deep foundation, are preferred instead of shallow foundations due to reasons such as increased loads from the structure, weak soil, or material losses. On the other hand, earthquake effects should also be considered in the design of engineering structures built in regions where seismicity is intense.
The p-y method (Matlock 1970), which is proposed to represent the nonlinear interaction between the pile and the soil, is frequently preferred in pile-soil interaction problems. The method is later developed for different soil types (Welch and Reese 1972;Reese et al. 1974Reese et al. , 1975Reese and Welch 1975). In these experimental studies, the behavior of piles in offshore structures under static and cyclic loads is examined and the process steps are suggested to obtain p-y curve sets for different soil types. However, these studies' most important and common deficiency is that they are recommended for single piles under static or reversible (repetitive, cyclic) loads. As a matter of fact, piles are generally designed and analyzed as a group and exposed to earthquake effects.
Firstly, Brown et al. (Brown et al. 1988b) compared the behavior of piles in group piles under static lateral loads with that of a single pile. In this study, it was concluded that the behavior of the piles in the back rows is quite different from that of the single pile due to the shadowing effect. In the same study, P-multiplier values were suggested for each row of piles in the direction of loading. This value is expressed as the ratio of the load on each pile in the group to the single pile. In the following years, extensive studies were carried out on this subject. Similar to the first studies, the behavior under static and cyclic loads was examined in these studies. McVay et al. (1998) investigated pile group interaction effects in loose and medium dense sand. Based on the test results, they found that p-multiplier are independent of soil density and only a function of the pile group geometry. Rollins et al. (1998) performed a full-scale pile group under static lateral loads to determine pilesoil-pile interaction effects. Compared with the single pile behavior, they found out that shadowing/group effects significantly reduced load capacity for all rows. A series of centrifuge tests are conducted to scrutinize laterally loaded pile groups in normally and over consolidated clays by Ilyas et al. (2004). Unlike other studies, they employed the individual performance of piles in the group. They concluded that the simplification of studying with the average performance of piles in the same row rather than the performance of individual piles is acceptable for lead and rear row but center piles in the middle rows carry much less loads than outer piles. Also, different equations have been proposed in the literature to determine the p-multipliers (Rollins et al. 2003;AASHTO 2017;TBEC 2018;Reese et al. 2019) depending on pile interval and geometry/row position.
There are also several 3D numerical studies on the p-y method in the literature (Yang and Jeremic 2003;Kim and Jeong 2011;Lu and Zhang 2020;Xiao-ling et al. 2020;Timurağaoğlu et al. 2021). Rahmani et al. (2018) evaluated static and dynamic p-y curves under lateral loading by comparing them with API (API 2007) curves. They indicated that API springs could not capture the mechanism in soil-pile interaction, although they are more often applied in design. By performing 3D finite element analysis to develop p-multipliers for fixed and free headed 3 × 3 to 5 × 5 vertically squared pile groups, Muhammed et al. (Adeel et al. 2020) found out that, for fixed-head pile groups, p-multipliers are lower because of higher group interactions compared with the free-head pile groups.
However, the literature on the behavior of piles under dynamic loads is limited. Considering two different single pile arrangements in a profile of soft clay overlying dense sand, Boulanger et al. (1999) investigated the dynamic p-y method under various shaking amplitudes. The dynamic p-y analyses were able to model the recorded responses reasonably. Ting et al. (1987) stated that the secant slope of the dynamic p-y curve is highly dependent on the loading frequencies in dynamic pile load tests. El Naggar and Bentley (2000) observed that the soil resistance increases under dynamic loading due to the damping effect. They also stated that the dynamic p-y curves depend on the loading frequency. In two different experimental studies (Yang et al. 2011;Lim and Jeong 2018) on the behavior of single piles under dynamic loads, it is stated that dynamic p-y curves give very different responses than static curves and the parameters affecting the response are determined by the acceleration amplitude and frequency given to the system (Yang et al. 2011;Lim and Jeong 2018), relative stiffness (Lim and Jeong 2018), pile bending stiffness and mass at the pile end (Yang et al. 2011). However, no experimental study focusing on the p-y method has been found in the literature on the behavior of group piles under dynamic loads. Especially in group piles under dynamic loads, the need to determine the parameters affecting the p-y method arises.
This study examines the dynamic p-y curves in group piles and the loading amplitude and frequency affecting these curves using three-dimensional (3D) numerical models. For this purpose, 3 static and 27 dynamic analyses are conducted for single pile, 2 × 2 and 3 × 3 group piles. The results are evaluated in terms of moment-depth, displacement-depth and p-y curves. P-multipliers are assessed by comparing with experimental results available in the literature. The variation of P-multipliers with lateral displacement and depth is also scrutinized. Static and dynamic p-y curves are generated and results are compared to assess the load sharing mechanism between each row and individual pile.

Soil constitutive model
A yield surface, a hardening and a flow rule to overcome limitations of the deformation theory of plasticity generally describe a constitutive model that can most accurately represent complicated soil behavior. Yield surface estimates whether the material should behave elastic or plastic under a load increment. The hardening rule governs the shape of the stress-strain curve due to plastic straining. The flow rule defines the direction of the plastic strain increment. It is well known that models with isotropic hardening rules generally apply to monotonic loading. In contrast, kinematic or mixed hardening models are more appropriate for cyclic or Bauschinger effects (Chen and Mizuno 1990). For the above considerations, an advanced constitutive soil model is selected for utilization in nonlinear dynamic analysis. This model is first proposed to simulate the cyclic response of shallow foundations (Anastasopoulos et al. 2011). Later, it is implemented in several different areas (Gerolymos et al. 2005;Giannakos et al. 2012;Kampitsis et al. 2015;Mucciacciaro and Sica 2018). The model, which is based on the known failure criterion of Von Mises, is originally offered for the cyclic response of metals to account for Bauschinger effects and is modified by inserting nonlinear kinematic and isotropic hardening components. The addition of the nonlinear kinematic hardening component specifies the yield surface translation in stress space throughout the backstress, α. The insertion of an isotropic hardening component characterizes the variation of the equivalent stress delineating the size of the yield surface, σ 0 , as a function of plastic deformation.

3
The evolution of kinematic hardening component is described to be the sum of a purely kinematic term (Ziegler 1959) and a relaxation term called the recall term, which introduces the nonlinearity (Lemaitre and Chaboche 1990) by the expression: where C k is defined as initial kinematic hardening or young's modulus, γ k specifies the reduction rate of kinematic hardening modulus with augmenting plastic deformation and pl stands for the equivalent plastic strain. Simple expressions to specify parameters of the model for clays and sands is suggested in the literature (Anastasopoulos et al. 2011). These parameters are given in Table 1 for clay. In the table, S u defines undrained shear strength of clay.
The yield surface size evolution, σ 0 , is described by isotropic hardening component as a function of the equivalent plastic strain, pl In Eq. (2), σ 0 is the yield stress at zero plastic strain, Q ∞ represents the maximum change in the yield surface size and b specifies the rate of the change at the yield surface size under increasing plastic straining. As highlighted earlier, we desire to apply only kinematic hardening and exclude isotropic hardening behavior for cyclic loadings. Thus, considering Q ∞ = 0, the model drops to a nonlinear kinematic hardening and the size of the yield surface remains constant. It must be noted that the model utilized in this study does not account for strain softening. The detailed information of KH model is available in mentioned papers and program documentation (Abaqus 2020).

Absorbing boundary conditions
Special/artificial boundary conditions should be provided for modeling purposes to absorb the waves at the boundaries, especially in the dynamic analysis of 3D models. These boundaries should be able to perfectly absorb all types of waves (body waves and surface waves) at all angles of incident and frequencies propagating outward. The infinite element available in Abaqus has been chosen as the boundary condition for the soil. These elements are placed at boundaries on four sides of the soil (two sides for 1-directional analysis). For infinite elements, only the elastic material needs to be defined (Abaqus 2020). Infinite elements provide the "silent" (non-reflective) boundary condition to the finite element model with a damping matrix effect, but do not contribute to the stiffness matrix of the finite element. Infinite elements contribute nothing to the intrinsic modes of the system. These elements are used for the far-field region of the static (Zienkiewicz et al. 1983) and dynamic (Lysmer and Kuhlemayer 1969) response. Infinite elements are successfully applied in literature (Van et al. 2017;Fatahi et al. 2018).

Soil-pile interaction
When the surfaces are in contact, they generally transmit shear forces and normal forces between the interfaces. These two force components form the basis of the piles' horizontal and vertical load transfer mechanism. The representation of pile-soil interaction is important in modeling piled systems using the finite element method. In the pile-soil interaction, separation, coupling and friction behavior must be accurately defined. This behavior can be represented numerically by defining normal and frictional behavior at the interface of soil and pile. Abaqus provides its users with a variety of contact formulations. Each formulation is based on a Contact Discretization, a Tracking Approach, and the assignment of "master" and "slave" roles to contact surfaces. In addition, the pressure-penetration relationship and constraint enforcement method needs to be defined.
For the contact discretization, the surface to surface formulation is adopted to provide the contact conditions on an average sense not only at the slave nodes, but also close to these nodes. In the tracking approach, there are two types; One is the finite sliding, which is the most general and allows any random movement of the surfaces, and the other is the small sliding, which assumes that even though the two bodies undergo large movements, there will be relatively little slip of one surface along the other. The finite sliding approach is the most general tracking approach and was chosen because it allows for random relative separation, sliding and rotation of the contacting surfaces. In assigning a master-slave role, when both surfaces in the contact pair are element-based and attached to the deformable masses or the deformable mass is defined as rigid, slave and master surfaces must be chosen. This selection is essential for node-to-surface contact. In general, if a smaller surface is in contact with a larger surface, it is optimal to select the smaller surface as the slave surface. If this distinction cannot be made, the main surface should be chosen as the surface of the more rigid mass, or if the two surfaces have comparable stiffnesses, the mass with the coarser mesh should be selected as the master surface.
When defining the pressure-overclosure relationship, using the hard contact minimizes the penetration of the slave surfaces to the master surface at the constraint regions. It does not allow the transfer of tensile stress across the interface. When the surfaces come into contact, any contact pressure can be transmitted between them. The surfaces are separated when the contact pressure drops to zero. The separated surfaces come into contact again when the gap between them drops to zero (Fig. 1a). However, over-constraining problems have emerged in the analysis when this relationship is applied. To overcome this problem, the penalty method is used in the study. This method is a stiff approximation of hard contact. The method attempts to predict the hard pressure-penetration behavior. In this method, the contact force is proportional to the penetration distance so that some penetration will occur (Fig. 1a).
On the other hand, the relationship known as friction between contacting bodies is usually expressed in terms of stresses at the interface of the bodies. The friction models available in the Abaqus program include the classical Coulomb friction models. In general, the Coulomb model allows the coefficient of friction to be defined by slip rate, contact pressure, the average surface temperature at the contact point, and area variables. The basic friction model assumes that the coefficient of friction is the same in all directions (isotropic friction). For a three-dimensional simulation, there are two orthogonal components of shear stresses τ 1 and τ 2 along with the interface between two bodies. These components move in local tangential directions to the contact surfaces or contact elements. The model combines two shear stress components as "equivalent shear stress" ( ) for stick/slip calculations ( = √ 2 1 + 2 2 ). In addition, the two shear rate components are combined ( ̇e q ) at an equivalent slip rate ( ̇e q = √̇2 1 +̇2 2 ). stick/slip calculations describe a surface where a point in the contact pressure-shear stress plane transforms from stick to slip (Fig. 1b). The friction coefficient was accepted as 0.5 in accordance with the recommended range in the literature (Mat Nor et al. 2015;Belinchon et al. 2016).

Boundary conditions and loadings
The analyses are conducted in two separate steps for static and dynamic analysis. In the static analysis of the first step, the bottom of the system is fixed while sides are constrained in x and y directions and vertical motions are released. In the dynamic analysis of the first step, the bottom is fixed and the horizontal displacement is constrained in the direction perpendicular to the loading direction while infinite elements are placed to the boundaries in loading directions to avoid wave reflections. In the second step, loadings are applied to the system. The load is applied to the pile cap as a force in static analysis. However, in implicit dynamic analysis, the displacements at the bottom in the loading direction is released. The loading was applied from the bottom of the system in the form of a sine wave (asin(2πft)) with different amplitudes (0.1, 0.3 and 0.5 g) and frequencies (1, 2 and 4 Hz). Applied sine input is given in Fig. 2 for 2 Hz and 0.5 g. The time step for dynamic analysis is selected as 0.004 depending on the mesh size and frequency of loading as stated in the literature (Jeremić et al. 2009). The pile diameter was divided into six element. Horizontal mesh size is around 0.17 m, while vertical mesh size is selected as 1 m. 8-node linear hexahedral solid elements with reduced integration (C3D8R) are used for both soil and pile. The reason we use a reduced integration element is that fully integrated elements don't hourglass, but can suffer from "locking" behavior. Also, these type of elements help to reduce the computational cost. When a reduced integration element is utilized, the locking doesn't occur, but the element is not rigid enough in non-critical bending to model the soil. Additionally, since the integration point is located in the middle of this element type, small elements are needed to capture the stress concentration at the interface. Thus, a finer mesh is adopted in the vicinity of the pile, while the coarser mesh is used through the external boundaries. 8-node linear, one-way infinite (CIN3D8) elements are used for absorbing boundaries. Generated model, together with the boundary conditions, are given in Figs. 3, 4 and 5. In both static and dynamic analysis, pile and pile cap are merged together.

Generated 3D FE model
The soil depth is selected as 30 m. The pile diameter is 1 m, which is the commonly used diameter. The pile is assumed to be reinforced concrete and modeled as elastic. The pile length is 15.2 m and the height of the cap is 1 m. The properties of the pile are summarized in Table 2. The pile seems to be flexible according to the equation suggested by Poulos and Hull (1989) who defined the rigidity parameter R = (E p I p /E s ) 0.25 , where E p and I p are Young's modulus and 2nd moment of area of pile and E s is soil Young's modulus. They suggested that if the length of the pile is less than 1.48 R, the pile behaves rigidly, and if the length exceeds 4.44 R, the pile behaves flexibly. According to the equation, the pile is a flexible pile with a rigidity (R) of 2.3. In soil pile analysis, the distance between piles was chosen as three times the pile diameter (3D). In this study, a 0.2 m gap was left between  the soil surface and the cap to neglect the influence of the pile cap. A schematic view of the 3 × 3 group pile model used in the analysis is given in Fig. 3. In single pile analyses, the dimensions do not change; only one pile is used. Homogeneous clayey soil with a PI of 50% and a shear wave velocity of 100 m/s was preferred in the analyses. This type of soil is defined as a ZE class according to the Turkish Building Earthquake Code (TBEC 2018). While smaller element size was chosen in the regions close to the piles, larger elements were preferred towards the boundaries. The generated 3D models are seen in Figs. 4 and 5 for single pile and 3 × 3 group piles systems.

Calibration of constitutive model parameters
Model parameters were numerically calibrated with the help of a strain-controlled simple shear test performed under cyclic loads (Fig. 6). The numerical tests were carried out in a single element. In these tests, calibration was performed for clay by considering the experimental G-γ and ξ-γ curve sets. It has been shown in the study that the plasticity index (PI) is one of the most effective parameters for many different soils (Vucetic and Dobry 1991). For this reason, PI-varying curve sets were used for clayey soils in this study. The results obtained from the cyclic simple shear test model for a plasticity index (PI) ratio of 50% are illustrated in Fig. 6, depending on shear strain magnitude. Numerical tests are implemented in Abaqus (2020) finite element computer program.
The detailed calibration procedure of the model parameters is available in the literature (Anastasopoulos et al. 2011). The hysteretic damping ratio is typically calculated as illustrated in Fig. 7 with the equation given. In the equation, ΔW represents the energy dissipated in a loop (the area inside the curve), while W represents the energy stored during the initial loading (the area of the shaded triangle).   After utilizing cyclic simple shear tests, one can calculate the change in shear modulus (G/G max ) and damping ratios (ξ) from shear stress-shear strain curves for various strain values. It is confirmed that PI is the main parameter governing G-γ and ξ-γ curves for various soils (Vucetic and Dobry 1991). After sensitivity analyses, a reasonable agreement with the G-γ curve is captured. The computed and experimental curves are compared for modulus degradation and damping ratios in Fig. 8. Numerical results suit the experimental curve of PI = 50% clay soil in the modulus degradation curve, whereas a reasonable agreement is encountered in the hysteretic damping curve. The parameters obtained after the calibration procedure are given in Table 3.

Validation of 3D model
The first step of soil-structure interaction is site response analysis. (Wolf 1985). For this reason, site response analyses were carried out under different frequencies (2, 4, 6 and 8 Hz). By changing the length/dimension of the system, the response spectra of 3D model (Abaqus 2020) obtained on the ground surface were compared with the 1D Fig. 9 Comparison of response spectra at the ground surface for different frequencies Deepsoil (Hashash et al. 2017) results and the horizontal dimensions of the model was decided. 1D and 3D elastic analysis results are given in Fig. 9. Here, while compatible responses are obtained for 1D and 3D analyses at some frequencies for 200 m length, it is understood that there are significant differences in the results at frequencies close to the resonance frequency. When the figures are examined, it is seen that the largest reactions occur at a frequency of 2 Hz. Undoubtedly, the resonant frequency is close to this value. When the results obtained by choosing different lengths in the horizontal dimensions of the model for this frequency are examined and considering the dimensions of the 3 × 3 group piles, 300 m has been selected as the horizontal dimensions of the model.

Static analysis results
Static analyses are conducted for single, 2 × 2 and 3 × 3 group piles to investigate the static p-multiplier parameter. Firstly, the shear stress-strain behavior of soil at ground surface of the soil-pile interface is given in Fig. 10. The figure shows that the soil behavior around the pile is nonlinear. When the figure is examined, it is seen that the shear stress values increase as the number of group piles increases at the same strain levels. This increase is due to the shadowing effect (Brown et al. 1988b), as known in the literature.
The displacement-depth relationship of a single pile is compared with 2 × 2 and 3 × 3 group piles in Fig. 11. It can be seen from the figures that the displacement behavior throughout the depth is different than single piles behavior for 2 × 2 and 3 × 3 group piles. The comparison of moment-depth relationships for the corresponding displacement values at the pile cap in Fig. 11 is given in Fig. 12. The results are shown for the same displacement level of the pile cap. According to the result obtained from the figures, as the number of group piles increases, the depth at which the maximum moment and zero lateral displacements at the pile occur increases. According to the results of both 2 × 2 and 3 × 3 group piles, while the piles in the front rows carry more moment, the piles in the back rows carry less moment. The maximum moments for piles in the trailing rows occurred at greater depths, in line with the findings of Brown et al. (1988b). Soil resistance (p) and pile displacement are generally calculated as the second derivative of the moment profile of the pile according to simple beam theory (Yoo et al. 2013;Haouari and Bouafia 2020). However, this method is difficult to implement especially when the Fig. 11 Comparison of displacement-depth relationships Fig. 12 Comparison of moment-depth relationships number of piles is increased and also, the results may alter when derivation is not well utilized. On the other hand, the soil resistance can be obtained numerically by dividing the sum of the shear and normal forces of the nodes in a pile at any depth to the vertical mesh size as below: where h is the vertical mesh size of the pile elements. Lateral displacement of pile (y) is directly obtained from the analysis results. The p-y curves at the first 3 m depths of 2 × 2 group piles are compared with the results of single piles in Fig. 13. In these figures, while pile 2 in the front row behaves closer to the single pile, the behaviors of piles 1 and 3 in the back row diverge from those of the single piles. These results are consistent with existing experimental studies in the literature (Rollins et al. 2003;Ilyas et al. 2004). However, although the behavior of the leading row piles is close to the behavior of the single pile at the ground surface, the behavior of all piles deviates from the behavior of the single pile as the depth increases. As expected, similar behavior is observed in piles 1 and 3 due to symmetrical boundary conditions and loadings.
The p-y curves at the first 3 m depths of 3 × 3 group piles are compared with the results of single piles in Fig. 14. Similar to the results of 2 × 2 group piles, while the piles in the leading row take more load and behave closer to the single pile results, the piles in the Fig. 13 Comparison of p-y curves for different depths of 2 × 2 group piles middle and rear rows take less load and move away from the single pile behavior. Also, complex behavior is observed especially in 3 × 3 group piles at smaller displacement levels (approximately 1 cm for this study). As the displacement increases, the loads on the piles become more pronounced and the behavior becomes easier to understand. Combining the maximum soil resistances across the depth, as shown in Fig. 15, yields results similar to Figs. 13 and 14. When this figure is examined, it is understood that the resistance of the single pile is higher than the resistance of each of the group piles for any depth. Using the p-y curves from Figs. 13 and 14 for 2 × 2 and 3 × 3 group piles, respectively, p-multipliers (Pm) are calculated. The calculated Pm values are compared with the existing literature for three rows as shown in Fig. 16 and given in Table 4. The equations used in Fig. 16 is also provided in Table 5. When the figure is examined, it is understood that the numerical p-multiplier values obtained under static loading in this study are compatible with the experimental studies in the literature. Additionally, within the scope of this study, the variation of the p-multiplier with the first 3 m depth was also examined for 2 × 2 and 3 × 3 group piles as given in Fig. 17. According to the results in the figure, the p-multiplier values generally have a decreasing trend.
Variations of p-multiplier with lateral displacement for 2 × 2 and 3 × 3 group piles at the ground surface are given in Figs. 18 and 19, respectively. At low displacement levels, the p-multiplier value increases in the piles in the front and middle rows, while the p-multiplier

Dynamic analysis results
Dynamic analyzes were performed using three different frequencies (1, 2 and 4 Hz) and three different amplitudes (0.1, 0.3 and 0.5 g) for single, 2 × 2 and 3 × 3 group piles. Thus, a total of 27 nonlinear dynamic analyzes were performed. The analysis is time consuming because the model is large. High-capacity computers are needed for such analyzes. In the analyzes made, each single pile models take about 24 h, while each 3 × 3 group pile models take approximately 48 h with a high-performance computer of 32 cores. Shear stress-strain behavior at the ground surface of the soil-pile interface of static and dynamic analyses (2 Hz and 0.5 g) for a single pile is compared in Fig. 20. The results of tie contact, in which the displacement and rotational degrees of freedom of the pile and the soil are equal, are also given in the figure. Since the system and loading are symmetric, soil nodes at the left and right side of the pile give similar results. Shear stress at the soil for  Rollins et al. (2003) 1st tie contact seems to be much higher than the frictional slip contact (labeled as "Int"). On the other hand, static analysis results are also higher than dynamic results under the same strain levels when the interaction is defined. The relative pile displacement is obtained by subtracting the pile displacement from the free field displacement (Yoo et al. 2013;Lim and Jeong 2018). On the other hand, free field displacement values should be taken from a distance greater than six times the pile diameter (6D) according to literature (see Fig. 16). In this study, the free field displacement values were taken from a distance of approximately 14D to stay on the safe side. Moment-depth and relative displacement-depth relationships under 2 Hz seismic loading for a single pile are compared in Figs. 21 and 22, respectively. As the amplitude increases, the increase at the moment and displacement values are seen. This behavior is in line with the findings in the literature (Rovithis et al. 2009;Yoo et al. 2013). At small amplitudes, the moment changes the sign according to the loading direction, while the sign becomes the same as the amplitude increases. Static results are also given in the figures. Where the maximum displacements at the pile top are the same for the static and dynamic cases, while the displacement changes along the pile are close, there is a difference in the moment values.
Shear stress-strain behavior at the ground surface of the soil-pile interface of soil node near pile 2 in static and dynamic analyses (2 Hz and 0.5 g) for 2 × 2 pile group are Fig. 18 Variation of p-multipliers with lateral displacement for 2 × 2 pile group compared in Fig. 23. It is seen that static and dynamic behavior are compatible with each other.
Moment-depth and relative displacement-depth relationships under 2 Hz seismic loading for 2 × 2 group piles are compared in Figs. 24 and 25, respectively. As indicated for single model results, the increase in amplitude increases the moments and displacements. The front piles take more moments in both loadings (+ or -) while experiencing less displacements. When both figures are examined, it can be seen that pile in the same row gives similar results due to symmetrical boundary conditions, meshes and loadings. Shear stress-strain behavior at the ground surface of the soil-pile interface of soil node near pile 3 in static and dynamic analyses (2 Hz and 0.5 g) for 3 × 3 pile group are compared in Fig. 26. From the figure, it is clear that static results generate greater stresses than dynamic behavior.
Moment-depth and relative displacement-depth relationships under 2 Hz seismic loading for 3 × 3 group piles are compared in Figs. 27 and 28, respectively. Considering the load distribution in the 3 × 3 group piles, it was concluded that the piles in the leading row received the highest load/moment at small amplitudes, but the piles in the middle row carried more load than those in the other rows at larger amplitudes. This is because there is an increase in strength in the soil around the middle row piles under dynamic loads. Here, there is a situation contrary to the load distribution obtained from the static analysis results in the literature and numerical results. According to studies in the literature, the piles in the front row receive the highest load, while the load decreases towards the rear rows. It has been understood that this result is different from the experimental results obtained from static horizontally loaded piles in the literature (Brown et al. 1988b;McVay et al. 1998;Ilyas et al. 2004). When Fig. 29 is examined, it can be seen that pile in the same row (piles 1 and 7) gives similar results due to symmetrical boundary conditions, meshes and loadings.
Soil resistance is expressed as the sum of the normal and shear/frictional resistances. The normal and shear forces occurring at 2 m depth on the contact surface between the pile and the soil are shown in Fig. 29. As expected, the normal contact forces increase in the middle of the pile and decrease towards the edges. Similarly, shear contact forces increase towards the edges and decrease towards the center.

Discussion on dynamic p-y behavior
According to the Korean Society of Civil engineering (KSCE 2008), pile displacement must be at least 1% of the pile diameter to create p-y curves. For this reason, p-y curves were evaluated considering this criterion in the study. As explained in Eq. (6), soil resistance can be calculated from the pile nodes' normal and shear contact forces in contact with the soil elements.
Static and dynamic p-y curves for a single pile under 2 Hz and 0.5 g seismic loading are compared in Fig. 30 for different depths. When the figures are scrutinized, it is seen that Fig. 23 Shear stress-strain curve for 2 × 2 group pile while the piles are moving in one direction, the resistance of soil in that direction increases. When the pile turns in the other direction, there is no resistance until it reaches the soil displaced in the previous cycle, but when the contact with the soil begins, the soil resistance starts to increase. When the pile changes direction, there is no contact and it is understood that the soil resistance approaches zero. This behavior reflects the actual soil-pile interaction in which separation, coupling and frictional behavior are defined. On the other hand, static p-y results are higher than dynamic results for each depth. Static and dynamic p-y curves for pile 2 of 2 × 2 pile group under 2 Hz and 0.5 g seismic loading are compared in Fig. 31 for different depths. Unlike single piles, peaks occur in the cyclic behavior of 2 × 2 group piles. This behavior is thought to be due to the wellknown shadowing effect. Static and dynamic p-y curves for pile 3 of 3 × 3 pile group under 1 Hz and 0.5 g seismic loading are compared in Fig. 32 for different depths. When the results of single, 2 × 2 and 3 × 3 group piles are evaluated together, it can be said that the behavior becomes more complex as the number of piles increases due to the shadowing effect.  (2007) under cyclic loading is also added to the figure for comparison. It is stated in the literature (Yang et al. 2011;Yoo et al. 2013) that dynamic p-y backbone curves are obtained by combining the maximum displacements of increasing amplitude and the corresponding soil resistances. However, the backbone of cyclic p-y curves is used in the study because the load is given in an increasing amplitude. Since the pile has a displacement of less than 1% at the 4 Hz frequency, the results of this frequency are not evaluated in the graphs. The results in Fig. 33 indicate that the dynamic p-y curves are independent of frequency. This result is compatible with the experimental study available in the literature (Yoo et al. 2013). More importantly, although the dynamic curves are compatible with the static curves at small displacements, they are below the static curves at large displacements for any depths. The curves suggested by API seem to underestimate static results. Figure 34 shows the static and dynamic p-y backbone curves of a 2 × 2 pile group under different frequencies for 1 m and 2 m depths in positive directions. As in single pile results, p-y curves are independent of frequency in 2 × 2 group piles. Similarly, dynamic curves fall under static p-y curves. API results in the figure are reduced by p-multipliers obtained from TBEC (2018). Like single piles, API curves underestimate static results in 2 × 2 group piles. Figure 35 shows the static and dynamic p-y backbone curves of a 3 × 3 pile group under different frequencies for 1 m and 2 m depths in positive directions. Similar to other results, p-y curves were independent of frequency in 3 × 3 group piles. Dynamic p-y curves stay under static ones at the leading row, although closer to static ones at the back rows. It is also evident that the API p-y curves underestimate the static results at deeper depths while an overestimation is encountered at shallow depths, as also found out by Li et al. (2017). The poor estimation of API curves is also emphasized in different studies (McGann et al. 2010;Rahmani et al. 2018) Dynamic p-y backbone curves for 2 × 2 and 3 × 3 group piles are compared in Figs. 36 and 37, respectively. According to the experimental and numerical static analysis results, while the piles in the leading row receive more load, the piles in the back row receive less load. According to the dynamic analysis of 2 × 2 group piles, as shown in Fig. 36, p-y curves indicate that the soil resistance of piles at the back row (P1) is higher than that of Fig. 26 Shear stress-strain curve for 3 × 3 group pile the leading pile (P2). Similar to 2 × 2 group piles, soil resistance of piles at the middle row (P2) seems to be greater than leading (P3) and back rows (P1), as given in Fig. 37. This is because there is an increase in soil resistance due to the compression in the soil around the middle row piles under dynamic loads. It has been understood that this result is different from the static experimental and numerical results obtained from horizontally loaded piles.

Conclusions
This paper investigates the static and dynamic behavior of p-y curves and p-multipliers for group piles. For soil, the kinematic hardening model is used and its parameters are calibrated by the experimental study frequently used in the literature. For the dynamic analysis, non-reflecting infinite elements are applied at the boundaries.
The conclusions based on the analysis results and comparisons can be summarized as follows. • Considering static analyses, as the number of piles increases, the depth at which the maximum moment, zero lateral displacements at the pile and the shear stress values at the soil around the pile increase. • Static analyses on p-y curves confirmed that the pile at the front row carries more load than that of trailing rows, although more complex behavior is encountered under low displacement levels for both 2 × 2 and 3 × 3 group piles.  • P-multipliers are generated based on static numerical analyses for 2 × 2 and 3 × 3 group piles. The generated values seem to be compatible compared with the experimentally calculated ones. The p-multipliers are in a decreasing trend throughout pile depth. • At low displacement levels, the p-multiplier value increases in the piles in the front and middle rows, while the p-multiplier starts to remain constant when the displacement increases. On the other hand, in the piles at the trailing rows, there is no significant change in the p-multiplier value while the displacement increases. • Comparing the shear stress-strain behavior of static and dynamic analyses for a single pile, shear stress at the soil for tie contact seems much higher than the frictional slip contact. On the other hand, static analysis results are also higher than dynamic results under the same strain levels when the interaction is defined. • In dynamic analyses of single, 2 × 2 and 3 × 3 group piles, as the amplitude of loading increases, the moment and displacement values increase both positively and negatively. • The moments in the leading row in 2 × 2 group piles under dynamic loading are greater than the trailing row as in static analyses. However, in 3 × 3 group piles, the piles in the leading row received the highest load/moment at small amplitudes, but the piles in the middle row carried more load than those in the other rows at larger amplitudes. This is because there is an increase in strength in the soil around the middle row piles under dynamic loads. This load distribution contradicts the static load distribution of 3 × 3 group piles. • Similar to moments, at lower amplitudes, the displacement profile of piles in the middle row is less than the piles behind and higher than the leading piles. However, the displacements of middle row piles at larger amplitudes seem smaller than other piles in the group. • When static and dynamic p-y curves are compared, it is understood that dynamic p-y curves stay under static results for single pile, 2 × 2 and 3 × 3 group piles. On the other hand, API (2007) suggestions seem to underestimate the static p-y backbone curves for single pile and all individual piles in groups. • The dynamic p-y backbone curves in clay are independent of frequency for all piles.
Although the dynamic curves are compatible with the static curves at small displacements for single and 2 × 2 group piles, they are below the static curves at large displacements for any depths. Dynamic p-y curves for 3 × 3 group piles stay under static ones at the leading row, although they are closer to static ones at the back rows. • Comparing the p-y curves of the individual pile in 2 × 2 group piles, the soil resistance at the leading row piles is slightly greater than back row piles. However, in 3 × 3 group piles, the soil resistance of the piles in the middle row is greater than the leading row piles. This is because there is an increase in soil resistance due to the compression in the soil around the middle row piles under dynamic loads. • The most important contribution of this study to the literature is that the designer/ expert who will carry out dynamic analysis in practical applications of soil-pile interaction can take into account the dynamic results obtained from this study. The findings Fig. 37 Comparison of dynamic p-y backbone curves for 3 × 3 group piles of this paper can provide designers/engineers with a better understanding how seismic soil-pile interaction effects the dynamic behavior of pile groups.
The conclusions drawn in the study are limited to the pile and soil parameters and soil material model used in the analyses. p-y behaviors can also be compared with different material models to expand the scope of future studies.