The spectral element method （SEM）-based numerical simulation of seismic wave motion has been widely applied in the study of earthquake source rupture process, large-scale seismic wave propagation, seismic response of regional complex sites without/with engineering structures （systems）, seismic tomography, and so forth. This technique is currently a frontier hotspot of common concern in the fields of earthquake science and technology including earthquake engineering, seismology, geophysics, etc. Spectral element method, which is sometimes also termed as spectral finite element method （SPECFEM）, spectral element method, or hp-type element method, is a combination of spectral method and finite element method （FEM）. Hence, it shares the advantages of both the two methods, i.e., the high accuracy and fast convergence of spectral method, and the regularity and flexibility of FEM.

In early times, the Chebyshev spectral element method （CSEM） and Legendre spectral element method （LSEM） originated from the domain decomposition of spectral methods, and therefore they inherit the complicated formulations of the latter, in which each of the interpolation basis functions is a linear combination of Chebyshev or Legendre orthogonal polynomials. Consequently, both the methods are as accurate as the spectral methods, but their applications are severely limited by the cumbersome and inefficient multi-layer nested computational structure that is resulted from those basis functions. Nowadays, the most frequently-used SEM is a concise form of LSEM developed by Komatitsch et al. In this method, the early-used complicated basis functions are simplified to the Lagrange shape functions that are commonly adopted in FEM, and those orthogonal polynomial-based analytical Gauss-Lobatto-Legendre （GLL） quadrature formulae are replaced by a simple numerical list of the GLL point coordinates and integration weights. Specifically, the non-equally distributed GLL points serve as the element nodes and the GLL high-precision numerical integration formula is applied to calculate the element mass, stiffness matrices, etc. This configuration makes the LSEM enjoy the same solution procedure and computational formulations as that of FEM, but avoid the significant defects of the classical high-order finite element method, including the intrinsic numerical error of the high-order polynomial interpolation based on equally-spaced grid and the lower computational efficiency due to the high-order consistent mass matrix. In a word, this LSEM has actually become a high-performance lumped-mass high-order finite element method. In addition to the above methods, the family of non-conforming spectral element methods has also been broadly studied and successfully applied in many problems, making themselves an important branch of the SEM. By introducing the so-called Lagrange multiplier or interior penalty term as a glue to effectively connect spectral elements with quite different sizes, orders, shapes and so on, the non-conforming SEMs are more flexible and highly efficient in dealing with multi-scale or discontinuous problems, which appear frequently in large-scale or complicated seismic wave simulations.

The great success of SEM is not only due to the high accuracy, regularity and flexibility of the algorithm itself, but also attributed to those well-designed open-source SEM programs represented by SPECFEM2D/3D, SPECFEM_GLOBE, SPEED, etc. These programs have incorporated a variety of key technologies that are required in complex simulations, such as three-dimensional complex models, different seismic source models or plane wave input method, large-scale parallel computing, global simulation, adjoint method, multi-scale or discontinuous modeling and so on. In the field of earthquake engineering, the SEM has been applied to perform physics-based deterministic numerical simulation of strong ground motion and to realize the “end-to-end” seismic response analysis that is from the source rupture to engineering structures or even engineering systems. The massive simulation data is a good supplement to the insufficient strong ground motion records, and the modeling of seismic wave propagation in actual geolocial structures can compensate for the weak physical background of traditional ground motion prediction equations （GMPEs） or stochastic methods. These simulations, which have reached a certain level of reliability, bring new vitality to earthquake engineering research and application. In the fields of seismology or geophysics, the highly-efficient forward simulation of SEM has been combined with the adjoint method, enabling a simultaneous modeling of the seismic wave fields generated from a number of observation stations, thus the number of forward simulations required for an inversion process can be reduced to an acceptable level. In this way, the advanced full wave inversion （FWI） or seismic tomography technique has been practically used to investigate seismic source mechanisms and to reveal regional or even global velocity structures. Finally, the development of SEM in China is elaborated. The SEM was introduced into China around the year of 2000, and the related studies mainly focused on the basic performance of the method and some preliminary applications until early 2010s. In the past decade, the Chinese researchers have been conducting more and more innovative work on the SEM algorithms and various engineering applications, and some of the work has reached the research forefront of the world.