562 The Leading Edge May 2009 SPECIAL SECTION: S e i s m i c m o d e l i n g Seismic modeling New developments in the finite-element method for seismic modeling S everal versions of the finite-element method (FEM) have been proposed for seismic modeling since it was first applied in this field in the late 1960s and early 1970s. Here we bring forward some recently developed FEMs suitable for seismic modeling due to their high-order accuracy. hese methods have been successfully applied to wave-propagation problems and improved accuracy, yielding an improved performance compared to classic FEM. he classical advantages of FEM include the flexibility with which it can accommodate surface topography, discon- tinuities in the subsurface model and boundary conditions, and the ability to approximate the wavefield with polynomi- als of arbitrarily high degree. hese advantages, however, ad- dress accuracy but not necessarily efficiency and performance issues. Recall that, unlike finite-difference methods (FDM) that simply use a Taylor series to approximate the derivatives in the wave equation, FEM uses the so-called weak formula- tion of the wave equation (see appendix) and a representation of the solution as a linear combination of shape functions. his results in a linear system of equations, the solution of which involves costly inversions of large (albeit sparse) matri- ces rendering the implementation slow in comparison with the competing FDMs. Furthermore, earlier implementa- tions considered low-order polynomial approximations of the wavefield, and, therefore, the results were contaminated with grid-dispersion error. As a result, the popularity of the method declined when compared to FDM. Fortunately, several recent developments have helped FEM overcome the limitations and have made the method attrac- tive once more. he main ideas in these developments are: he use of high-degree polynomials to approximate the • wavefields he diagonalization of the matrices to be inverted (called • mass matrices) through mass-lumping techniques he use of high-order time-stepping schemes • he methods that combine these ideas, and that we will review here, are the spectral-element method, mixed FEM, and discontinuous Galerkin method. he spectral-element method he main characteristic of the spectral-element method (SEM) is that it uses shape functions of high polynomial degree. his method, originally developed for fluid dynam- ics, was first implemented for the acoustic wave equation by Seriani and Priolo (1994) using the Gauss-Lobatto-Cheby- shev (GLC) nodes. he method was extended to elastic wave propagation by Komatitsch and Vilotte (1998) using the Gauss-Lobatto-Legendre (GLL) nodes, and it quickly gained credibility within the seismological community. It has been JONAS D. DE BASABE and MRINAL K. SEN, The University of Texas, Austin applied to problems in regional and global seismology, ob- taining results close to the empirical data (Komatitsch et al., 2002). Although the ability to use high-degree polynomials is now perceived to be an advantage of FEM, this has not always been the case. In fact, the early work studying the applicabil- ity of the high-order FEM to wave propagation discouraged this practice because it was believed that it would introduce spurious or nonphysical modes. Furthermore, the bandwidth of the matrices increases proportionally with the polynomial degree; therefore the high-order methods have a bigger com- putational cost if they have to invert these matrices. Fortunately, both concerns have been addressed by SEM. First, about the spurious modes introduced by the high-order methods, it was shown that their contribution is negligible (Mulder, 1999). Furthermore, it has now become clear that using high-degree polynomials practically eliminates the grid- dispersion error (De Basabe and Sen, 2007), and therefore a lower sampling ratio can be used, typically 4–5 nodes per wavelength. Secondly, the penalty associated with the matrix inversions is avoided in SEM by a clever selection of the grid points and quadrature rules that yield a diagonal mass matrix. his is achieved by using the GLL nodes for the distribution of the nodes inside the elements and the GLL quadrature rules for the approximation of the entries in the mass and stiffness matrices (Cohen, 2002). he use of the GLL nodes and quadrature rules is tantamount to mass lumping, with the advantage that it preserves the accuracy. he accuracy of SEM in 2D using the GLC nodes was empirically investigated by Seriani and Priolo for the acous- tic scheme. hey concluded that an average of 4.5 nodes per wavelength with an eighth-order method eliminates the grid dispersion. Although this showed promising results, the use of GLC nodes does not lead to mass lumping. he accuracy of the 1D acoustic SEM scheme, using the GLL, GLC, or equis- paced nodes, was examined by Mulder. He analyzed the grid- dispersion error of the high-order methods and concluded that the error introduced by the spurious modes can be neglected and that SEM using GLL nodes was more accurate than using the GLC or equispaced nodes. Further, Cohen analyzed the grid dispersion of the 1D, 2D, and 3D acoustic SEM schemes using GLL nodes and quadrature rules. He showed dispersion curves for the 1D case using second- or third-order methods and various time-stepping schemes. hese papers consider the acoustic case only. he grid dispersion and stability analysis have been recently extended to the elastic case in De Basabe and Sen. he mixed finite-element method he mixed finite-element method (MFEM) uses the velocity- stress formulation of the wave equation and allows different SPECIAL SECTION: S e i s m i c m o d e l i n g Downloaded 06 May 2009 to 129.116.220.172. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/