The goal of this study is to find out how nonlinear soil behavior affects the dynamic impedance of a rigid surface foundation that is put under a harmonic dynamic load. Researchers can use a computational code that combines the boundary element method (BEM) and the thin layer method (TLM) to calculate the nonlinear dynamic impedances (stiffness and damping) of foundations in the frequency domain. This method integrates nonlinear soil behavior into a linear framework. The soil's nonlinear behavior is evident in the changes observed in its dynamic properties, specifically the decrease in the normalized shear modulus (G/Gmax), the increase in the normalized hysteretic damping coefficient (ξ/ξmax), and the variation in the equivalent shear wave velocity (Vs/Vsmax) as a result of the seismic excitation applied to the soil, which is dependent on the unit shear strain (γ/γr). The parameters are derived using the equivalent linear approach employed in the Caldynasoil computational code for various levels of seismic loading. We conducted a comprehensive analysis to assess how soil non-linearity affects the behavior of a soil-foundation system under various conditions. This included considering factors such as the soil's homogeneity or heterogeneity, the shape of the foundation, and whether the foundation was supported by a semi-infinite soil or constrained by bedrock. The study also considered different levels of seismic deformation, using five seismic accelerations ranging from 0.01g to 0.4g as reference points. The linear and nonlinear dynamic impedance coefficients of the massless foundation were determined for all vibration modes (translation, rocking, and torsion) as a function of the dimensionless excitation frequency a₀