A spectral element method for modelling surface wave dispersion of Love and Rayleigh waves is presented. This method uses standard Gauss–Lobatto–Legendre polynomials coupled with a Gauss–Laguerre–Legendre element to represent a half-space, and in doing so, improves both the efficiency and accuracy of the calculation of phase and group velocities, particularly at lower frequencies. It is demonstrated that this method is able to directly represent 1-D earth models with smoothly varying structure, layered structure or combinations thereof. The method is both efficient and accurate, while the solution error can be tuned across a range of frequencies to balance the trade-off with computational cost. In addition, the adjoint technique is used to develop an efficient algorithm for the calculation of the gradient of arbitrary misfit functions with respect to earth model parameters, which is of prime interest in the inversion of surface waves. It is demonstrated that accurate gradients can be computed for misfit functions based on phase velocity, group velocity, and Rayleigh wave ellipticity. To demonstrate both the spectral element method and the adjoint method developed in this paper, two simulated inverse problems are presented using Love wave phase velocity observations and Rayleigh wave ellipticity observations.