Abstract: Given modification of earthquake motions due to liquefaction-induced soil softening, the effective stress analysis should be conducted to develop the site-specific design ground motion parameters at liquefiable sites. A loosely coupled nonlinear effective stress method for the site response analysis is proposed. In this method, the nonlinear hysteresis model for soils is incorporated with an excess pore water pressure generation model characterized by cyclic shear-volume strain coupling, which establishes the coupling relationship between the degradation of cyclic stiffness and the generation of excess pore water pressure associated with earthquake events. The material subroutine in ABAQUS/Explicit platform is developed. The method is then used to simulate the seismic response of the downhole array non-liquefied and liquefied sites in Japan. The results show a good consistency between the simulations and the recordings at different depths: (1) The difference between the simulated and the recorded peak ground accelerations for the non-liquefied site is only 2.6% (weak motion) and 11.8% (moderate motion). (2) The difference between the simulated and recorded peak accelerations at different depths for the liquefied site is less than 5%. (3) The simulated and recorded spectral acceleration curves of the non-liquefied and liquefied sites are in good agreement. Compared with the responses of a downhole array multilayer liquefied site calculated by DEEPSOIL V6.0, those of the proposed method show a better agreement with the seismic recordings at different depths.