To study the effects of nonlinearity in the seismic response of soils, a numerical simulation of the propagation of vertically incident seismic waves in horizontal soil layers were performed. Shear noiselike and monochromatic seismic waves of various intensities were used as input signals. The behavior of soils was described by a nonlinear hysteretic model. To extract and study nonlinear components in the ground response, the nonlinear system identification method and analysis of higher-order spectra of oscillations on the surface were applied. Even for weak input signals, the response of the simulated soils contained a noticeable nonlinear component. An increase in the intensity of input signals led to increasing distortions of propagating signals, due to the generation and growth of combination-frequency harmonics. The results show that odd types of nonlinearity are most typical for soils, such as cubic and fifth-order nonlinearities, causing generation of the third and fifth higher harmonics of main frequencies of input signals. Nonlinearities of even types, such as quadratic, fourth-order, and sixth-order, concerned with asymmetry, or skewness, of oscillations (i.e., quasi-static deformations of the surface) are usually weak, except some special cases, in which a stress-strain relationship of a soil can be represented by functions with noticeable even components. A weak nonlinearity results in an increase in high-frequency components, due to the generation of higher harmonics. In cases of strong nonlinearity, in which a decrease in amplification and in shear moduli become noticeable, changes in spectra of propagating signals achieve their maximum. As a result, input signals with arbitrary spectra are transformed into output signals with spectra of the type of E(f) ∼ fk, where k depends on the properties of the medium.