SUMMARYA spectral element method for the approximate solution of linear elastodynamic equations, set in a weak form, is shown to provide an e cient tool for simulating elastic wave propagation in realistic geological structures in two-and three-dimensional geometries. The computational domain is discretized into quadrangles, or hexahedra, deÿned with respect to a reference unit domain by an invertible local mapping. Inside each reference element, the numerical integration is based on the tensor-product of a Gauss -Lobatto -Legendre 1-D quadrature and the solution is expanded onto a discrete polynomial basis using Lagrange interpolants. As a result, the mass matrix is always diagonal, which drastically reduces the computational cost and allows an e cient parallel implementation. Absorbing boundary conditions are introduced in variational form to simulate unbounded physical domains. The time discretization is based on an energy-momentum conserving scheme that can be put into a classical explicit-implicit predictor=multicorrector format. Long term energy conservation and stability properties are illustrated as well as the e ciency of the absorbing conditions. The accuracy of the method is shown by comparing the spectral element results to numerical solutions of some classical two-dimensional problems obtained by other methods. The potentiality of the method is then illustrated by studying a simple three-dimensional model. Very accurate modelling of Rayleigh wave propagation and surface di raction is obtained at a low computational cost. The method is shown to provide an e cient tool to study the di raction of elastic waves and the large ampliÿcation of ground motion caused by three-dimensional surface topographies.