## 3 Finite Element (FE) Model  
The 3D nonlinear FE model of the pile group is shown in Figure 1. The size of the computational domain was determined after a sensitivity analysis such that the calculated responses were not affected by the boundaries. The length and width of the numerical model were set to \(50D\) and \(33D\) from the center of the foundation, where \(D\) is the pile diameter. The pile configurations considered in this study are 3 \(\times\) 3, 4 \(\times\) 4, and 5 \(\times\) 5 pile groups. The size of the computational domain is 30, 20, and 15 m in length, width, and height, respectively. The convergence analysis for the finite element mesh was also performed to determine optimum element sizes to obtain accurate results. The mesh was generated in such a way that it was finer near the piles and coarser towards the boundaries of the computational domain. The width of the smallest element was \(0.15D\). Eight-node brick elements (C3D8) were used to model both the piles and soil. The interface between the piles and soil was modeled using a surface-to-surface contact model that allows for both slipping and normal separation (gapping). The Coulomb model was used to simulate the tangential slip, where the friction coefficient was set to \(\tan(2/3\varphi)\), as used in the study of [PERSON] et al. (2018).  
The pile group was placed at the center of the computational domain. The length of the piles was fixed to 12 m. Pile bottoms were tied to the soil elements. The bottom of the computational model was fixed in the horizontal and vertical directions. The horizontal displacement constraints were applied at the lateral boundaries. No constraint was applied at the surface of the soil domain. To simulate the fixed pile head condition, pile heads were tied with a pile cap. The pile cap was fixed in the vertical direction, whereas lateral movement was allowed. The piles were modeled using the linear elastic model. The properties of the hollow steel pipe pile are listed in Table 4. The piles were modeled as solid rods with an outer diameter of 0.3 m. To achieve identical flexural rigidity as the reference steel pipe pile, the modulus of elasticity of the solid piles was adjusted such that it was equivalent to that of the steel pipe pile, as summarized in Table 4.  
Figure 1: Finite element (FE) model for the \(5\times 5\) pile group: (**a**) free-head and (**b**) fixed-head condition.  
The nonlinear soil stress-strain behavior is simulated using the bounding surface plasticity model of [PERSON] [39]. The model was selected because it has been widely used in the simulation of the seismic response of soils. It is reported to have no purely elastic region, and therefore effective in modeling even the small-strain nonlinearity of soil. The shear modulus reduction curve derived from the plasticity model is as follows:  
\[\frac{G}{G_{max}}=\ 1-\frac{3}{2\gamma_{0}}\int\limits_{0}^{2\tau_{0}}\!\!\! \left[h\left(\frac{R/\sqrt{2}+\tau_{0}-\tau}{\tau}\right)^{m}+H_{0}\right]^{-1 }\!\!d\tau \tag{1}\]  
where \(\frac{G}{G_{max}}\) is the secant shear modulus normalized to maximum shear modulus; \(\gamma_{0}\) is the shear strain; \(\tau_{0}\) is the shear stress; \(R\) is radius of the bounding surface; \(h\) and \(m\), which are the coefficient and exponent of the exponential hardening function, respectively; and \(H_{0}\) is the kinematic hardening parameter. This plasticity model was implemented in ABAQUS using the UMAT subroutine code developed by [PERSON] et al. [40]. The numerical models of single and \(3\times 3\) group piles were validated against the field test measurements of [PERSON] et al. [25]. The details of the numerical model and validation results are reported in [PERSON] et al. [41].