This paper aims at investigating the existence of localized stationary waves in the shallow subsurface whose constitutive behavior is governed by the hyperbolic model, implying non-polynomial nonlinearity and strain-dependent shear modulus. To this end, we derive a novel equation of motion for a nonlinear gradient elasticity model, where the higher-order gradient terms capture the effect of small-scale soil heterogeneity/micro-structure. We also present a novel finite-difference scheme to solve the nonlinear equation of motion in space and time. Simulations of the propagation of arbitrary initial pulses clearly reveal the influence of the nonlinearity: strain-dependent speed in general and, as a result, sharpening of the pulses. Stationary solutions of the equation of motion are obtained by introducing the moving reference frame together with the stationarity assumption. Periodic (with and without a descending trend) as well as localized stationary waves are found by analyzing the obtained ordinary differential equation in the phase portrait and integrating it along the different trajectories. The localized stationary wave is in fact a kink wave and is obtained by integration along a homoclinic orbit. In general, the closer the trajectory lies to a homoclinic orbit, the sharper the edges of the corresponding periodic stationary wave and the larger its period. Finally, we find that the kink wave is in fact not a true soliton as the original shapes of two colliding kink waves are not recovered after interaction. However, it may have high amplitude and reach the surface depending on the damping mechanisms (which have not been considered). Therefore, seismic site response analyses should not a priori exclude the presence of such localized stationary waves.@en