A method is presented for calculating the seismic response of a system of horizontal soil layers. The essential element of the method is a rheological model suggested by Iwan which takes account of the nonlinear hysteretic behavior of soils and has considerable flexibility for incorporating laboratory results on the dynamic behavior of soils. Finite rigidity is allowed in the underlying elastic medium, permitting energy to be radiated back into the underlying medium. Three alternate ways of integrating the equations of motion are compared, an implicit technique, an explicit technique, and integration along characteristics. An example is set up for comparing the different methods of integration and for comparing the nonlinear solution with a solution based on the widely used equivalent linear assumption. The example consists of a 200-m section of firm alluvium excited at its base by the N21E component of the Taft accelerogram multiplied by four to produce a peak acceleration of 0.7 g and a peak velocity of 67 cm/sec. The three techniques of integration give very similar results, but integration along characteristics has the advantage of avoiding spurious high-frequency oscillations in the acceleration time history at the surface. For the chosen example, which has a thick soil column and a strong input motion, the equivalent linear solution underestimates the intensity of surface motion for periods between 0.1 and 0.6 sec by factors exceeding two. The discrepancies, however, would probably be less for input motion of lower intensity. At longer periods the equivalent linear solution is in essential agreement with the nonlinear solution. For the same example both solutions show that, compared to a site with rock at the surface, motion at the surface of the soil is amplified for periods longer than 1.5 sec by as much as a factor of two. At shorter periods the amplitude is reduced.