The present paper deals with the multiple scattering by randomly distributed elastodynamic systems at the surface of a horizontally layered elastic halfspace due to an incident plane wave. Instead of solving this problem for a particular configuration of the system, multiple scattering theory is used to compute the ensemble response statistics. The Dyson equation is used to calculate the mean field, while the nonstationary second order statistics are obtained by means of the Bethe-Salpeter equation. This allows for the determination of the mean square response of the system in the time and frequency domains. This model is used to study multiple scattering between buildings under seismic excitation. The influence of multiple scattering on the seismic site response is verified. Furthermore, the influence of the footprint and the damping of the buildings are investigated. The results are compared to results of a coupled finite element/boundary element solution for a group of buildings.