<p>This study employs advanced 3D nonlinear finite difference analysis to investigate the seismic response of pile groups in layered clay deposits to address critical gaps in soil-pile interaction models in layered soils. A validated modified advanced nonlinear constitutive model captures clayey soil dynamic behavior, including strain-dependent stiffness degradation and damping. The entire numerical model and methodology were verified through centrifuge tests. Parametric analyses systematically evaluate the effects of soil layer depth ratio, strength heterogeneity, and seismic excitation properties using Kobe (1995) and Landers (1992) records. Results reveal that maximum bending moments and shear forces consistently occur at or near soil layer interfaces, with magnitude and distribution governed by the relative stiffness and depth of clay layers. Notably, if the stiff clay layer underlies soft clay, peak seismic responses of the pile increase as the soil layer depth ratio approaches 1.0; however, if piles are confined within a thick stiff layer, peak seismic responses of the pile decrease. In addition, if the top clay layer is stiffer and stronger compared to the bottom layer, it induces a sharp curvature and secondary peak in the pile response diagrams. For stiff-over-soft clay profiles, pile seismic response exhibits higher sensitivity to the soil layer depth ratio when the shear strengths of the layers are higher. Moreover, resonance effects dominate when earthquake frequencies align with system natural frequencies and amplify system seismic responses despite lower peak ground acceleration and duration. This study gives valuable insights into practical implications for optimizing pile design in stratified clays.</p>