Numerical modeling of mechanical wave propagation

The numerical modeling of mechanical waves is currently a fundamental tool for the study and investigation of their propagation in media with heterogeneous physical properties and/or complex geometry, as, in these cases, analytical methods are usually not applicable. These techniques are used in geophysics (geophysical interpretation, subsoil imaging, development of new methods of exploration), seismology (study of earthquakes, regional and global seismology, accurate calculation of synthetic seismograms), in the development of new methods for ultrasonic diagnostics in materials science (non-destructive methods) and medicine (acoustic tomography). In this paper we present a review of numerical methods that have been developed and are currently used. In particular we review the key concepts and pioneering ideas behind finite-difference methods, pseudospectral methods, finite-volume methods, Galerkin continuous and discontinuous finite-element methods (classical or based on spectral interpolation), and still others such as physics-compatible, and multiscale methods. We focus on their formulations in time domain along with the main temporal discretization schemes. We present the theory and implementation for some of these methods. Moreover, their computational characteristics are evaluated in order to aid the choice of the method for each practical situation.


Introduction
Each field of applied sciences has particular requirements for computational modeling and often develops its own suite of numerical techniques. The numerical modeling of mechanical waves in some applications involve two somewhat conflicting requirements: • Complex, heterogeneous structures must be correctly modeled. In particular, interfaces and shapes of geological structures must be taken into account during space discretization. Moreover, high accuracy is needed for avoiding numerical anisotropy, attenuation, or dispersion that mislead interpretation. As a result, the numerical solution requires a huge computational effort, both in memory storage (from Giga to Tera nodes) and CPU time (from hours to weeks); • The medium where waves propagate is iteratively updated to fit recorded data.
Wave simulation is one step of an imaging/inversion algorithm that may be repeated several times, thus it must be fast enough for not compromising the entire process.
These requirements are typical in the analysis of material response and imaging, for instance in • Exploration geophysics, reservoir scale seismics; • Geotechnical and engineering seismology; • Local, regional, and global seismology; • Planetary seismology; • Earth's interior imaging; • Ground-shaking risk analysis -strong ground motion; • Monitoring of volcanic processes; • Earthquake and Tsunami early warning systems; • Global monitoring of nuclear tests.
Traditionally, these demands have been met with high-order schemes [1,2]. They are highly accurate methods that require a low number of grid points per wavelength, thus reducing storage and CPU time requirements. Regardless of the chosen method, an efficient implementation is needed for reducing the total cost of the simulations. Some alternatives are resorting to vector/parallel platforms (massively parallel, clusters, GRID [2]), efficient subroutines and libraries (FFT, Lapack, MPI [3]), and seeking a low count of operations and of primary storage.
The numerical methods that have been developed for the above-mentioned purposes constitute a multidisciplinary field named computational seismology, the numerical simulation of seismic wave propagation in arbitrary 3D models [2]. Its scope is naturally beyond global-scale seismology, reaching other topics of geosciences (such as rock physics, exploration geophysics, volcanology, and geotechnical engineering) and beyond (computational mechanics, materials science, underwater acoustics, and medicine).
The recent literature provides detailed reviews oriented towards specific methods [4][5][6][7] or communities [8][9][10]. The purpose of this paper is to share an overview of computational seismology methods with a broader audience, starting from the mathematical models, visiting general aspects of spatial and temporal discretization and then arriving at the theoretical and computational aspects of the main numerical methods currently in use.

Scalar wave equation
The most elementary mathematical model of wave propagation is the scalar wave equationü where the unknown function u(x, t) may denote, for instance, the acoustic pressure, c is the wave velocity in a homogeneous medium, and f (x, t) is the divergence of an external body force. The dots denote time differentiation. A related model is the one-dimensional shear-wave propagation equation where ρ is the density and μ is the shear modulus. When the latter is constant, Eq. (2) can be written similarly to (1) by defining c = (μ/ρ) 1/2 . In the absence of source functions, Eq. (1) has the general solution where the functions F(·) and G(·) are arbitrary. Some important particular solutions are the d'Alembert solution which satisfies the initial conditions u(x, 0) = u 0 (x) andu(x, 0) = u 00 (x), and the plane-wave solution where the angular frequency ω and the wave number κ satisfy the dispersion relation ω = ±cκ. Equation (1) can be generalized to 2D or 3D media as follows: which admit plane-wave solutions of the form u(x, t) = exp[−i(ωt − κ · x)], with ω = ±c|κ|, when f = 0. Both Eqs. (1) and (6) may be considered in the more general case of a homogeneous velocity field c(x). Moreover, the more general acoustic wave equation ∂ ∂t accounts for variable density. We recall that Eq. (7) arises from the linearized mass and conservation equations 1 Ku ρv + ∇u = F.

Elastic wave equation
The standard model for wave propagation in solids is given by the conservation of linear momentum (Newton's law) where u is the displacement and σ is the stress tensor. In particular, elastic media are described by the linear constitutive relation (Hooke's law) σ (u) = C : (u), (u) = 1 2 ∇u + ∇u T (10) between the stress σ and the linearized strain (u), where C is the elasticity tensor. Due to the symmetry provided by the conservation of angular momentum, it is convenient to use Voigt notation σ = {σ x x , σ yy , σ zz , σ xy , σ xz , σ yz } , (11) or σ = {σ x x , σ yy , σ xy } in 2D, for which the governing equations may be written as where the differential operator D is The elasticity tensor C has up to 21 free parameters, but there can be significantly fewer ones depending on the symmetry assumptions [12]. When the medium is isotropic, we have C = λItr( ) + 2μ , where I is the identity operator and λ, μ > 0 are the Lamé coefficients. In experimental studies, other elastic parameters are more typical, such as the elastic modulus and the Poisson's ratio, which respectively are E:=λ + 2μ and ν:=λ/2(λ + μ) for isotropic media. Under Voigt notation, tensor C has the following matrix representation: In vector notation, Eq. (9) assumes the following standard form: When f = 0 and the elastic parameters are constant, the plane wave u(x, t):= R exp[−i(ωt − κ · x)] is a solution to (15) if where c P := λ + 2μ ρ , c S := μ ρ . (17) are the compressional and shear-wave velocities. The vector Eq. (16) admits the solutions (ω P , R P ) and (ω S , R S ), where the angular frequencies are ω P = ±c P |κ| and ω S = ±c S |κ|, while the propagation directions R P and R S (R SV and R S H in 3D) are parallel and perpendicular to κ, respectively [13]. A relevant class of anisotropic elastic stress-strain relations in computational seismology is that of transversely isotropic media with vertical symmetry axis (VTI): where C 66 = (C 11 − C 12 )/2. Following [14], a plane-wave solution u(x, t):= R exp[−i(ωt − κ · x)] in the three-dimensional case with κ = κ{sin θ, 0, cos θ } (without loss of generality, given the cylindrical symmetry [15]) yields the following dispersion relations: ω P = ±c P 1 + ε sin 2 θ + (θ ) ω S H = ±c S 1 + 2γ sin 2 θ where c P = (C 33 /ρ) 1/2 , c S = (C 44 /ρ) 1/2 , ε = (C 11 − C 33 )/(2C 33 with r C P = 1 − c 2 S /c 2 P . The fact that ω SV and ω S H may not coincide leads to the phenomenon of shear-wave splitting. The dimensionless numbers ε, δ, and γ are known as Thomsen parameters and serve as measures of anisotropy. Under additional assumptions, the elasticity tensor C can be entirely written in terms of c P , c S , and Thomsen parameters [16,17].

Viscoelasticity
The model of viscoelastic media is based on stress-strain relations that account for not only the instantaneous strain, but all its history. This is accomplished by a convolution in time between the strain rate of change and a time-dependent tensor G known as the relaxation tensor [18]: This equation may be recast in differential form through fractional derivatives [12, Sec. 2.5.2] and can be numerically handled with the aid of auxiliary memory variables [19]. If G(t) = H (t)G 0 , where H is the Heaviside function, then we recover the elastic model (10) Important aspects of wave propagation in viscoelastic media are present in the Kelvin-Voigt model, which combines both cases above, in the one-dimensional case: This constitutive relation yields the following wave equation: which serves as a preliminary site-response model for small strains [20]. The dispersion relation for this equation is ω 2 = κ 2 (G 0 − iωη 0 ), hence phase (ω/κ) and group (dω/dκ) are complex and frequency-dependent, highlighting two relevant aspects of viscoelastic wave propagation: attenuation and physical dispersion, respectively.

Poroelasticity
Wave propagation in fluid-saturated porous media had been studied by M. A. Biot for a variety of cases [21][22][23]. For simplicity, we consider an isotropic solid matrix with constant porosity φ. Let us denote the densities and bulk moduli of the constituent solid and the saturating fluid as ρ s , ρ f and K s , K f , respectively. Moreover, the bulk and shear modulus of the dry matrix will be denoted as K d and μ d . For convenience, let us also introduce the Lamé parameter Though the governing equations may be written in terms of the displacement vectors u s , u f at the solid and fluid phases [21], it is convenient to substitute fluid displacement by w = φ(u f − u s ), which represents the flow of the fluid relative to the solid but measured in terms of volume per unit area of the bulk medium [23]. In this case, the equations of motion are where ρ = (1 − φ)ρ s + φρ f is the density of the saturated matrix, T is the tortuosity, η is the fluid viscosity, and k is the permeability. The constitutive relations for the total stress σ b = σ + I(−φ p) and the fluid pressure p can be written as [12]. A plane-wave analysis can be performed by decoupling the P− and S− modes of propagation [12]. Namely, we can apply the divergence operator to Eq. (24) to find a system of equations for ∇ · u s and ∇ · w, and apply the curl operator to the same equations and find another system for curl u s and curl w. Proceeding to the analysis of propagation in a single direction, we find two compressional velocities rather than a single one; the propagation mode associated with the velocity of lower magnitude is known as the Biot's slow wave.

Velocity-stress formulation
The mathematical models reviewed above involve partial differential equations of second order in time for the displacement field. If we seek instead the velocity field (and introduce the stress as an additional variable), we are led to a system of first-order equations, for which a large variety of numerical methods is available.
Let us consider for instance the elastic Eq. (9). By taking time derivatives of both sides of (10), we find where the unknowns are the stress σ and the velocity v =u. In matrix form, we havė A related first-order system is given by the velocity-displacement formulation From the computational point of view, formulation (28) has the advantage of involving less unknowns than in (26). For instance, in the three-dimensional case we have six unknowns rather than nine.

Boundary conditions
Modeling of wave-propagation problems may involve not only physical but also computational (or artificial) boundaries depending on the region of interest, denoted by (the problem domain). Physical boundaries are usually modeled by transmission conditions of the form σ · n = g, where n is the unit vector normal to the boundary and pointing outwards. The case g = 0 is referred to as a free-surface boundary condition.
Ideally, a computational boundary should not interfere with the waves, which makes it very challenging to model. One of the classical approaches is to use absorbing (or non-reflecting) boundary conditions [7,24,25]. As pointed out in [26], absorbing boundary conditions are related to the Sommerfeld radiation condition lim r →0 for Helmholtz equation u + k 2 u = 0, in the sense that the analogue of (30) for the scalar wave equation (6) is Engquist and co-authors [24,27] obtained the following boundary condition at r = a in the two-dimensional case: as well as higher-order conditions, based on paraxial approximations of the wave equation (see also [28]). In Cartesian coordinates, the lowest-order boundary conditions at x = 0 (a typical lateral border in a two-dimensional simulation) are Conditions (33) were also derived by annihilating the reflection coefficient of planewave solutions [29][30][31]. A more general approach was later proposed by Higdon [25,32], who considered boundary conditions of order p in the form which reduce to Engquist-Majda conditions when α j = 0, 0 ≤ j ≤ p. These coefficients may be chosen to minimize the reflection coefficient of plane waves traveling with angle of incidence θ . Bamberger et al. [33] proposed modified conditions that account for Rayleigh waves. Another relevant progress on absorbing boundary conditions is handling corner points [34][35][36][37].
As pointed out in [7], the error of a high-order absorbing boundary condition does not necessarily converge to zero as the order tends to infinity. When the error due to the boundary condition does converge to zero, it is referred to as an exact non-reflecting boundary condition [38][39][40][41].
Moreover, higher-order approximations involve high-order spatial and temporal derivatives, which must be appropriately represented in the numerical discretization and usually incur a higher computational cost. Such a constraint has motivated the study of high-order local non-reflecting boundary conditions (high-order local NRBCs, [7,36]), which introduce auxiliary variables that avoid the need of high-order derivatives. As outlined in [7], one of the first approaches of this class, due to Collino [34], can be written as where θ j = jπ/(2 p +1). Thanks to the auxiliary variables φ 1 , . . . , φ p , the derivatives in (36) have order no greater than two.

Variational formulation
The variational, or weak formulation is a convenient representation of the mathematical model that allows to seek the approximate solution in a functional space with lower regularity requirements, for instance when the material properties are discontinuous.
For conciseness, let us focus on system (9)- (10). Its weak form is obtained by taking the scalar product of both sides of (9) by a test function w and integrating over the domain . In the case of a homogeneous Dirichlet condition, where d = 2 or 3 and H 1 0 ( ) = {u ∈ L 2 ( ) | ∇u ∈ L 2 ( ) and u | ∂ = 0}, where L 2 ( ) is the space of functions that are square-integrable with respect to the Lebesgue measure in and u | ∂ denotes the trace of u over the boundary of [42]. Analogously, the variational formulation of system (26) is where X = {τ ∈ L 2 ( ) d×d ; τ = τ }. A similar variational formulation can be obtained for the velocity-displacement formulation (28). One may also consider a formulation of system (26) with lower regularity on velocities and higher regularity on stresses, which can be discretized with mixed finite elements (see, e.g., [43]), and has been adapted to include perfectly matched layers [44].

Model discretization
The differential or variational formulations presented in the previous section, when complemented with proper initial and/or boundary conditions, provide a unique wave field that is complete in the sense that it can be determined for any point x of the domain and any time t ≥ t 0 , where t 0 is the initial time of observation (for convenience, we consider t 0 = 0 from here on). Since those initial-boundary value problems rarely have an analytical solution, we must resort to approximate solutions.
For some methods such as finite-difference and mimetic methods, the approximate solution corresponds to an array of coefficients U n j such that U n j ≈ u(x, t n ) for any point x in a subset X j that usually consists of a single point x j , but could also be an edge, face, or a three-dimensional shape. The differential/integral operators are approximated or replaced with discrete operators defined at the sets X j for any j and at t = t n , n = 0, 1, . . ..
For another group of methods that is considered in this survey, we have a finite expansion of functions in the form where φ 0 (x), . . . , φ N (x) are previously chosen functions, whileû 0 (t), . . . ,û N (t) are time-dependent coefficients to be determined. Pseudospectral and finite-element methods fall into this category. Through expansion (39), time and space are separately handled, in analogy with the method of separation of variables for partial differential equations. The spatial operators are applied to functions φ j (x) in original form, i.e., they are not discretized, though the operator evaluation may depend on interpolation or numerical integration.
To determine the arrays U n j or the coefficientsû j (t), each method relies on a particular approximation principle, but most methods need to discretize the independent variables, space and time. Time and space sampling could be required in the approximation principle, in auxiliary calculations, or in the output generation.

Spatial discretization
The main requirements for the discretization in space are that the waves and the medium heterogeneities must be sufficiently sampled. While the former is strongly dependent on the chosen numerical method, the latter is similar for most methods.
The geological model is usually discretized in a finite number of cells (also referred to as elements or blocks) where, ideally, • there are no restrictions on the material variability; • interfaces and structures can be honored; • different rheologies can be easily handled.
Following [2], we classify the ensemble of all cells of a domain, called the mesh or grid, as structured or non-structured. A structured mesh is obtained by the mapping from the integer set {0, 1, . . . , N 1 }×· · ·×{0, 1, . . . , N d } to an d-dimensional domain . For example, a cubic mesh with uniform spacing h can defined by the mapping of each {i, j, k} ∈ {0, . . . , N } 3 to the vertex (x 0 + ih, y 0 + jh, z 0 + kh). As data structures, these meshes are uniquely defined by arrays with the coordinates of the vertices. The cells are understood to be formed by points whose coordinates have adjacent indices in a Cartesian fashion.
On the other hand, non-structured meshes do not have a Cartesian structure, hence the coordinates are not sufficient to identify the mesh. Usually there is one mapping from {0, 1} d , or another reference cell, for each cell in the mesh. Despite their complexity, non-structured meshes conform more easily to interfaces and structures, while structured meshes usually need a large amount of grid nodes to avoid the staircase effect [45].

Temporal discretization
The discretization of the time variable can be done in a straightforward manner: a uniform partition of the interval [0, T ] into N intervals with time step t = T /N : where T is the final time of the simulation. The time step will depend on the smallest grid length in space and the chosen approximation method for the time derivatives, and in general it can be dynamically selected by means of error estimators [46][47][48].
One of the greatest difficulties with temporal discretization is the situation where some elements in the spatial grid are extremely small (such as slivers generated from 3D unstructured mesh-generation codes [49]), forcing an equally small time step and undermining the efficiency of the wave-propagation simulation. This has motivated the development of local time-stepping schemes [50][51][52][53][54], where the time step is smaller where the grid is refined and larger where the grid is coarse.

Temporal discretization methods
In the following we review the methods for approximation of time derivatives that are most frequently used in seismic wave propagation. These methods are typically applied to the following second-order linear system of ordinary differential equations (ODEs): Mü where the matrices M, K, and C are known as the mass, stiffness, and damping matrices, respectively. The latter arises from the discretization of viscous terms such as in Eqs. (22) and (24b), but also absorbing boundary conditions, such as (33). One can also seek the block vector v = {u,u} by solving a first-order system of the forṁ which is also the kind of system of ODEs that arises from velocity-stress and velocitydisplacement formulations (27)- (28). The convenience of system (42) resides on the fact that it is a particular case of the classical equationṡ v = f (t, u), (43) for which a vast literature is available (see, e.g., [55]). On the other hand, system (41) as well as and its generalization to time-dependent coefficient matrices [56] and non-linear internal forces [57], among others, has been thoroughly studied by the community of structural dynamics. For simplicity, temporal discretization methods will be described for the case of a uniform partition (40). Perhaps the most popular method is the leapfrog scheme, where the approximation u n ≈ u(t n ) over the time partition is determined from the centered finite-difference approximation which is present in some of the pioneering works on computational seismology [58,59] and has been employed in several different approaches (e.g., [60][61][62]).
Some schemes such as leapfrog can be written as multi-step methods of the form If B 0 = I in (45), the scheme is called explicit, otherwise, it is implicit. In the literature of finite-element methods, the concept of explicit is extended to the cases where B 0 is diagonal [63,64]. In general, we may refer as explicit a method that does not require solving a linear system in any step of the computation of u n+1 . For instance, when M = I and C = 0, then leapfrog method (44) becomes explicit:

Newmark methods
The leapfrog method is one member of a family of time-integration methods known as Newmark methods [65]. As described, e.g., in [64], the Newmark method for (41) with parameters β and γ can be written as follows: where a n ≈ü(t n ) and v n ≈u(t n ). The Newmark scheme naturally provides approximations for not only displacement, but also velocity and acceleration, which are useful to the inversion of three-component data [66]. The leapfrog method corresponds to β = 0 and γ = 1/2. Another well-known method from the Newmark family is the average-acceleration method [67], which corresponds to β = 1/4 and γ = 1/2. This method has been applied to acoustic and elastic wave propagation [68,69]. The Newmark scheme has also been interpreted as a time-staggered velocity-stress algorithm with the purpose of implementing absorbing boundary conditions [70]. Newmark methods are at most second-order accurate [67], i.e., u n − u(t n ) = O( t 2 ), and the use of higher-order methods may increase the computational efficiency despite their higher cost [71]. In the following we review several high-order temporal discretization methods.

Lax-Wendroff methods
One of the most traditional high-order approaches are the Lax-Wendroff methods [72], which use spatial derivatives to replace high-order time derivatives. This principle is also present in the arbitrary high-order derivatives (ADER) method [73,74] and nearly analytical discrete methods [75,76].
This scheme has been widely employed to second-order equations in the form which can be derived from (41) in the absence of damping (C = 0) as, e.g., in [54]. If the spatial discretization is performed with finite differences [71,77], then M = I and there is no need of transformations from (41) to (48). Most authors have considered the acoustic wave equation [54,71,77,78], but the Lax-Wendroff approximation has also been applied to the elastic wave equation [79,80]. Following [77], let us derive the fourth-order Lax-Wendroff method for (48). We refer to [78] for higher-order approximations. A standard Taylor expansion of the second-derivative term in (48) yields On the other hand,v(t n ) = f (t n ) − Av(t n ) and, by taking a second-order derivative of this expression, we find By combining (49) and (50), we arrive at the following explicit scheme: The vectorf n can be approximated by a second-order scheme [81]. The use of the Lax-Wendroff scheme has also been investigated in the first-order scheme (42), not only for acoustic waves [82], but also for structural dynamics [83] and viscoacoustic waves [84].

Runge-Kutta and symplectic methods
The Runge-Kutta (RK) methods [55] can be readily applied to system (42). For instance, the classical fourth-order RK method is given as follows: The RK discretization for first-order systems in the form (42) has been used in conjunction with pseudospectral [85], finite-element [86], and discontinuous Galerkin [87,88] methods. The extension of RK methods to the second-order system (41) can be done by 2-step schemes known as Runge-Kutta-Nyström methods [89,90], which have been thoroughly developed in the context of Hamiltonian systems [91], Let Φ : ( p 0 , q 0 ) → ( p(t), q(t)) be the flow map defined by the solution ( p(t), q(t)) to (53) with initial conditions p(t) = p 0 and q(t) = q 0 . Its Jacobian Φ satisfies (Φ ) JdΦ = J, and if the same property holds for the flow map produced by a numerical method for (53), then this method is called symplectic [92,93]. Symplectic schemes have slower error growth and preserve conservative quantities, which makes them an attractive choice for long-time simulations [94,95]. These schemes have also been proposed for the more general Birkhoffian systems [96], with application to poroelastic wave propagation [97].
The implicit, fourth-order Störmer-Numerov method, which has been applied to acoustic and elastic waves [98][99][100], is an example of a symplectic method [101]. The Störmer-Numerov approximation for system (48) is the following: (54) Makridakis [98] has extended this method to solve system (41), where the damping term accounts for absorbing boundary conditions.

Approximation of evolution operators
Another family of high-order methods is based on approximations of matrix power series that are present in the analytical solutions of systems of ODEs. As shown for instance in [102], the analytical solution of system (48) If a is a scalar, then C(a, t) = cos(a 1/2 t) and S(a, t) = sin(a 1/2 t)/a 1/2 . Taking into account the temporal discretization (40), the solution (55) may be represented as follows: (56) In order to derive numerical schemes from the recursive form (56), we may determine approximations to C(A, t) and S(A, t), and select quadrature schemes to deal with the source term [102]. For instance, when f = 0 and the rational approximation then (56) leads to the following implicit scheme: which corresponds to the Störmer-Numerov method (54) when β = 1/12 and has been revisited in [81,103]. Baker et al. [104] studied sufficient conditions for the convergence of rational-approximation methods in the homogeneous case. One of the high-order schemes that satisfies these conditions is defined by the following approximation: which leads to 6th-order accuracy. Some classical methods are related with Taylor approximation of evolution operators [8,105]. Indeed, a Taylor approximation of degree m ≥ 1 for (56) in the case which are the leapfrog and Lax-Wendroff scheme when m = 1 and m = 2, respectively. Besides rational (Padé) and Taylor expansions, Chebyshev expansions can also be employed [106]. Such an approach is known in the seismic exploration literature as the rapid expansion method [107]. One can also consider the approximation of exponential operators associated with the analytical solution of the first-order system (42): Similarly as in (55), several approximations for exp(At) may be employed, such as Taylor, Fourier [105], Chebyshev [105,107,108], and rational [109,110] expansions. Besides using a truncated expansion of the matrix operator, one may also split the matrices in diagonal and non-diagonal parts, leading to a fixed-point iteration method [111].

Spatial discretization methods
The numerical methods described next provide approximations of the mathematical models described in Sect. 2 as the systems of ODEs (41) or (42).

Finite-difference methods
Most time-discretization schemes described in Sect. 4 are based on finite-difference approximations of time derivatives. These formulas are based on Taylor-series expansions that are combined to reach the desirable accuracy. For instance, the expansion where t n−1 ≤t 1 ≤ t n and t n ≤t 2 ≤ t n+1 , leads to the finite-difference formulä over a space-time grid defined by the points (x j , t n ) = (x 0 + j x, n t). This scheme is second-order accurate if f has continuous fourth-order derivatives. By applying the same approximation to the second partial derivative in space, we arrive at the following finite-difference approximation of the scalar wave equation (1): At the time step n ≥ 1, Eq. (64) for 0 < i < N complemented with boundary conditions lead to the fully discrete system (44). For instance, under boundary conditions In the simpler situation where c = 1, f = 0, and x = t, Eq. (64) reduces to which is considered in the seminal paper by Courant et al. [112]. The same formula along with the centered finite-difference approximation of first-order derivatives was used in the approximation of the elastic wave equation in cylindrical coordinates [58] and in 2D Cartesian coordinates [113]. Alford and co-authors proposed schemes with fourth-order accuracy in space for both acoustic [59] and elastic [114] wave equations, with the purpose of improving accuracy (in particular, reducing numerical dispersion).
For the velocity-stress formulation (27) of the elastic wave equation, the use of staggered grids [115] is a standard practice [116][117][118][119]. A related approach seeks second-order accurate finite-difference formulas with a lower number of grid points, leading to simplicial grids [120]. Igel [2] illustrates staggered grids with the 1D elastic wave equation (2), whose velocity-stress form is The space and time derivatives may approximated by centered finite differences with spacing x/2 and t/2, respectively ( Fig. 1): (68) In this manner, v and σ are computed on disjoint spatial grids with spacing x rather than x/2, reducing computer memory requirements. For higher dimension, one may need to interpolate stress components from separate grids to evaluate stressstrain relations [8].
As illustrated in scheme (64), the mass matrix M is typically the identity matrix in finite-difference methods, thus it is natural to consider explicit temporal discretization schemes. On the other hand, implicit schemes require an efficient implementation to handle the additional cost of solving linear systems, and a common approach is to use sequential splitting techniques [140].
Splitting methods have their roots in the classical alternating-direction [141,142] and fractional-step [143] methods, which have been gathered as locally onedimensional (LOD) methods [144]. LOD has been extended from parabolic to hyperbolic problems [145][146][147], including high-order time discretizations [148]. Perhaps the first high-order ADI method for the acoustic wave equation is due to Ciment and Leventhal [148], though other works have previously considered extending ADI methods from parabolic to hyperbolic problems [145][146][147]. More recently, these methods have employed to seismic wave-propagation problems [81,103,149]. Wave equations with viscous terms have been considered in [150].
Another splitting approach due to Strang [151] and Marchuk [152] has also been adopted in wave-propagation problems [153]. A related technique has been employed to separate the stiff from the nonstiff part of the velocity-pressure poroacoustic wave equations, so that Biot's slow wave can be numerically modeled [154]. A theoretical framework to analyze splitting methods for general second-order systems of ODEs has been proposed in [155].
For instance, let us review the splitting of scheme (58) for systemv + Av = 0 following [81], in the case where A is the two-dimensional, second-order finitedifference operator where u n = A x v n and w n = A y v n are given as follows: Note that A x and A y involve approximations of the x− and y− directions, respectively. Scheme (58) is approximated as the following three-step formula: Note that step (71a) is explicit whereas steps (71b)-(71c) involve tridiagonal linear systems. The error of (71) with respect to (58) is O( t 4 ) [81].
As pointed out by Emerman et al. [156], schemes such as (71) may allow larger time steps than explicit schemes but have very low accuracy. Thus implicit methods must employ high-order finite-difference approximations (see, e.g., [103]) in order to become competitive.
A well-known limitation of traditional finite-difference methods is that they poorly handle irregular interfaces, topography, or boundary conditions, and several approaches have been proposed to circumvent these difficulties [157]. Alterman and Nathaniel have addressed the case of a constant slope by means of a change of coordinates [158]. Ilan has adapted the work in [158] to standard Cartesian coordinates, and by allowing a non-uniform grid, has extended it to polygonal topography. It is worth noting that several schemes have been proposed for non-uniform finite-difference grids [159][160][161]. Jih et al. [162] revisited the same issue and proposed local changes of coordinates, which allow the use of a uniform grid. This technique evolved to boundary-conforming grids defined by curvilinear coordinates [163][164][165][166][167], motivated by the work by Fornberg [168]. Additional approaches are hybrid finite-difference and finite-element/discrete-wavenumber methods [169], the vacuum method [170], the interface method [171,172], and the use of non-matching grids [173].

Pseudospectral methods
As mentioned in Sect. 3, another approach to discretize the initial-boundary value problems arising from wave-propagation models is to employ the finite expansion (39). The schemes that follow this approach are known as spectral methods [174][175][176], and are essentially characterized by the choice of the basis functions φ j (x) and the way to determine the expansion coefficientsû The classical choices for the approximation space are orthogonal trigonometric or polynomial functions, while the approaches to determine the expansion coefficients are classically divided into tau [177], Galerkin [178], and collocation [179] methods.
Both tau and Galerkin methods choose the coefficients such that the solution satisfies the variational formulation in the approximation space, and they differ on how boundary conditions are handled. Currently, the tau method is seldom used [180]. The Galerkin method is better known in the form of finite-and spectral-element methods, which use piecewise-polynomial interpolation basis functions and are discussed later on. The Galerkin technique has also been proposed with wavelet basis functions [181][182][183][184][185]. Global orthogonal polynomials are rarer, and are mostly applied to waves in fluids [186]. On the other hand, spectral methods with the collocation technique, which became known as pseudospectral methods [187], have achieved great popularity thanks to the Fast Fourier Transform (FFT) algorithm [188].
In the one-dimensional case, a pseudospectral method is essentially defined from a set of functions φ 0 , . . . , φ N that are orthogonal with respect to some inner product ·, · and collocation points x 0 , . . . , x N that are chosen such that the orthogonal projection which corresponds to the best approximation of a function u(x) in the vector space In seismic wave propagation, the earlier works are due to Gazdag [45], Kosloff and Baysal [60], which was generalized to account for anisotropy [189] and viscosity [19,190]. They considered Cartesian grids with uniformly spaced collocation points and complex Fourier basis functions, and the resulting method be interpreted as the limit of finite differences with infinite order of accuracy [191].
In the Fourier pseudospectral method, the basis functions are φ j (x) = exp(iκ( j)x), with κ( j) = 2π j/(N x). The corresponding collocation points are x j = a + j x, which have uniform spacing x = (b−a)/N over the interval [a, b]. For convenience, the indices run from 0 to N − 1 so that N coincides with the number of subintervals. The relationship between u N (x j ) andû j can be written in terms of the discrete Fourier transform pairv The expansion coefficients are determined from a system of algebraic equations that is obtained by evaluating the differential equations at the collocation points, which requires evaluating derivatives of the expansion u N . The Fourier method is more naturally derived by approximating the computation of spatial derivatives in the frequency domain. For instance, the approximate solution to (6) with leapfrog time discretization in two dimensions [45] can be written as where with κ x (ĵ) = 2πĵ/(N x x) and κ z (k) = 2πk/(N z z), whileû nĵ ,k is defined by the 2D discrete Fourier transform The calculations from (75) and (76) can be efficiently carried out with the FFT algorithm. This algorithm requires O (N log N ) operations, which is significantly lower than the O(N 2 ) operations of matrix-vector multiplication for large N . On the other hand, this method assumes periodic boundary conditions, demanding additional strategies to implement realistic ones (see, e.g. [192]).
To circumvent this problem, Raggio [193] proposed the Chebyshev pseudospectral method. Chebyshev basis functions have also been used in only one spatial direction, either on polar [85] or Cartesian [68] coordinates in order to better handle free-surface boundary conditions. Its extension to three-dimensional problems can be found in [194].
The Chebyshev basis functions in the interval [−1, 1] are T j (x) = cos( j cos −1 (x)), while the collocation points are x j = cos( jπ/N ), 0 ≤ j ≤ N . Herein, N denotes the maximal polynomial degree. The coefficientsû j in (72) arê These coefficients can be written aŝ which is related to the real part of the discrete Fourier transform. The coefficients of the derivatives of u N can be computed fromû j through recursive relations [195].
The fact that the collocation points x j are clustered at the boundary implies that the distance between grid points can be very small, which in turn leads to a small time step as well. For this reason, a stretching transformation should be employed [196,197]. An alternative approach uses the tau method with Legendre polynomials [198], but the Chebyshev pseudospectral method has become more popular as it can be implemented with the FFT algorithm.
As finite-difference methods, pseudospectral methods benefit from staggered grids [199][200][201]. Their classical implementations were also limited to regular grids, and one alternative to avoid such a restriction is to employ curvilinear coordinates [168,202]. Another approach is to resort to domain decomposition, which was proposed for the elastic wave equation initially for isotropic media [203,204] and later to viscoelastic [205,206] and poroelastic [207] media. Another benefit of using domain decomposition in pseudospectral methods is that the resulting wave operator is not entirely global, avoiding non-causal interactions of the propagating wavefield with parameter discontinuities in the model [200].
The Fourier pseudospectral method has been recently generalized to handle fractional derivatives, which are useful to model attenuating media without the need of memory variables [208][209][210].

Finite-element methods
As mentioned in the previous section, finite-element methods [63,64] belong to the family of Galerkin methods, and typically use continuous Lagrange interpolation basis functions, which are associated with the spatial grid. The earliest applications of finite elements to seismic wave-propagation problems are due to Lysmer and Drake [211] in the frequency domain, and to Smith [212] in the time domain. Later works addressed several relevant aspects of seismic modeling with finite elements [5,213,214].
In the following we describe a finite-element approximation of the variational problem (37). The first step is to decompose the spatial domain into n e non-overlapping elements e such that = ∪ e . Physical elements e are mapped through a transformation x = Λ e (ζ ) onto a reference elementˆ where computations are actually performed. The approximate solution may be written as where the union operator denotes thatũ is defined in andũ(x, ·) =ũ e (x, ·) for any x ∈ e . Moreover, n dof is the number of element degrees of freedom. The basis functions satisfy Φ e j (Λ e (ζ )) =Φ j (ζ ), whereΦ j is the Lagrange interpolation vector function associated with the j-th degree of freedom in a space of polynomial, vectorvalued functions inˆ . For instance, if the polynomial degree is one and the spatial dimension is two withˆ = where ϕ 1 (ζ ) = (1 − ζ )/2 and ϕ 2 (ζ ) = (1 + ζ )/2, for ζ ∈ [−1, 1]. The coefficients u e j (t) are determined from the following Galerkin approximation of (37): whereṼ is the subspace of H 1 0 ( ) d of continuous piecewise-polynomial functions built from local functionsΦ e j (x). After algebraic manipulations, Galerkin equation (81) are written as the system of ordinary differential equations whereM e ,K e , andF e are the elemental mass matrix, stiffness matrix, and load vector in sparse global form, that is, only non-zero entries are used and they are mapped into appropriate global locations by a connectivity map from local to global nodes [215]. The dense elemental arrays M e , K e , and F e are defined by the contributions from element e to the integrals in (81), for 1 ≤ i, j ≤ n dof . These integrals are computed in the reference elementˆ through standard changes of variables [64].
In general, finite-element basis functions are defined from linear, quadratic and cubic polynomials over triangular or quadrilateral elements (tetrahedral or hexahedral in 3D, though other elements such as pyramids and wedges can be used [216]). Low-order finite-element methods are comparable to centered finite differences [213] and thus have low accuracy, but can be useful when the problem geometry leads to highly refined and irregular meshes [217,218]. One alternative to improve the performance of low-order finite elements is to enrich the approximation space [219][220][221][222]. It is worth noticing that isogeometric analysis, a finite-element approach integrated with computer-aided design [223], has been shown to provide a better geometric representation than traditional finite elements. Several authors have studied this technique in wave-propagation problems [224][225][226].

Spectral-element methods
Unlike pseudospectral methods, the use of high-order finite elements was not common in the literature of computational seismology at least until the 1990s, mostly because there were concerns about their accuracy [227].
Standard high-order finite-element bases are based on equally spaced polynomial interpolation, which is an ill-conditioned problem [228,229]. This can be noticed from the behavior of Lagrange basis functions of degree N at the reference element [−1, 1] with equally spaced points j/N , 0 ≤ j ≤ N . As shown in Fig. 2, the Lagrange functions associated with the element midpoint begin to oscillate as N increases, similarly to the Runge phenomenon.
A classical approach to circumvent this problem is to use the Chebyshev collocation points rather than equally spaced points, as shown in Fig. 3. The use of Chebyshev points as well as other orthogonal polynomial roots was initially regarded as unnecessary [230], but in the work of Patera and co-workers [231][232][233], these points provided the link between finite-element and pseudospectral methods, with the appeal of having the geometric flexibility of the former and the rapid convergence properties of the latter.
Spectral elements are usually quadrangular or hexahedral, so that the Lagrange shape functionsΦ j (ζ ) = Φ e j (Λ e (ζ )), defined in the reference elementˆ = [−1, 1] 2 or [−1, 1] 3 , are built from tensor products of one-dimensional Lagrange shape functions ϕ i (ζ ) of degree N such that ϕ i (ζ j ) = δ i, j , as in as in (80). The collocation points ζ i (0 ≤ i ≤ N ) are given by (85) for Chebyshev elements, while for Legendre elements these points are the solutions to where L N (ζ ) is the N th degree Legendre polynomial [1]. The standard implementations of spectral-element methods with Chebyshev (SEM-GLC) and Legendre (SEM-GLL) collocation points may be classified as consistent and lumped finite elements, respectively. In other words, the mass matrix in (83) for Chebyshev elements is calculated without approximations (assuming piecewise constant density) and is non-diagonal, while for Legendre elements this matrix is approximated through reduced integration by a diagonal matrix. SEM-GLL has been adapted to triangular meshes [245][246][247][248], but the selection of collocation points is more complex and may lead to non-diagonal mass matrices.
The computation of elemental matrices in SEM-GLC is based on the properties of Chebyshev polynomials T j (cos θ) = cos( jθ). In particular, the Lagrange shape functions ϕ i (0 ≤ i ≤ N ) are obtained by choosing u in (72) such that u(x j ) = δ i, j . It follows from (72) and (77) that so that the entries of the elemental mass matrix in [−1, 1] arê The integral in (88) [215,231]. These formulas have been generalized to take into account variable material properties, which are represented by expansions in basis functions that do not necessarily coincide with those from the wave field [215,249].
Because fully discrete SEM-GLC schemes are implicit in time, they need efficient linear-system solvers. Some useful strategies are the element-by-element formulation and suitable factorizations of matrix-vector products [250,251]. Moreover, unconditionally stable time-integration schemes should be chosen to allow the use of large time steps.
For SEM-GLL, the standard practice is to employ the GLL quadrature with ζ j satisfying (86). This formula is exact if f is a polynomial of degree ≤ 2N − 1.
In particular, the calculation is not exact, since ϕ i ϕ j has degree 2N . On the other hand, since ϕ i (ζ j ) = δ i, j , we haveM G L L i, j = w j δ i, j , i.e.,M G L L is diagonal. The diagonality of the mass matrix was crucial in the success of spectral elements in computational seismology.
The spectral-element method has been implemented for anisotropic, visco-and poroelastic wave propagation [61,252,253] and the advent of collaborative codes and platforms [254, Appendix A] has encouraged its application in several studies of regional and global seismology [6,[255][256][257]. This method has also been widely used in conjunction with adjoint methods [258,259].

Finite-volume methods
The finite-volume method for elastic wave propagation was initially proposed by Dormy and Tarantola [260] for the velocity-displacement formulation (28) and later by Tadi [261] for the second-order formulation (15).
In [260], the main motivation to introduce the finite-volume method was to generalize minimum grid, second-order finite differences [120] to unstructured grids and irregular boundaries. The key idea was to use the divergence theorem to obtain derivative estimates of a field from its values at surrounding grid points, rather than Taylor series expansions. On the other hand, the main concern in [261] was the enforcement of traction boundary conditions. By using the ADER method [73,74] and a reconstruction algorithm to introduce high-order numerical fluxes, Dumbser et al. [262] have obtained higher accuracy both in space and time. Later on, Zhang and co-authors [263,264] proposed a high-order finite-volume method that combines the reconstruction algorithm from [262] and the element subdivision algorithm from the spectral volume method [265]. The finitevolume method has been used in studies of dynamic rupture [266], hypoelastic [267], and poroelastic media [268].
As in [269], let us introduce the finite-volume methods through the one-dimensional conservation lawu Given a uniform spatial grid with spacing x, let C j = [x j−1/2 , x j+1/2 ] be a grid cell (or finite volume) centered at node x j . By integrating both sides of (91) over C j , we find A subsequent integration from t n to t n+1 yields which can be written as where is the average approximation of the unknown field u(x, t n ) at C j , while (96) are the average fluxes at x j±1/2 . Assuming that F n j±1/2 approximately depends only on the adjacent average values U n j and U n j±1 , we arrive at the finite-volume approximation of (94): where the numerical-flux function F(U − , U + ) approximates the average fluxes (96). For instance, let us consider the linear case f (u) = au with a > 0, which is a one-way wave equation. Taking into account that the information propagates from the left to the right for a > 0, an effective choice is the upwind flux F( which coincides with the classical upwind finite-difference method. In the case where a may change sign, we choose F(U − , U + ) = max{a, 0}U − + min{a, 0}U + , so that The jumps U j−1/2 = U n j −U n j−1 and U j+1/2 = U n j+1 −U n j can be interpreted as waves moving across cells C j and C j+1 with opposite senses; in general, the numerical flux is driven by the solution of a Riemann problem [269].
Let us now consider the 1D elastic wave equation (67) in the homogeneous case, which can be written in a matrix form similar to (91); From the eigenvalue decomposition A = RΛR −1 , the vector w = R −1 u satisfies We can apply the upwind scheme (99) to either w + or w − : The matrix form of (102) yields an upwind scheme for w: By multiplying both sides of (103) by R on the left and replacing the jumps W n j±1/2 with the jumps R −1 U n j±1/2 in terms of the average approximations of u(x, t n ), we arrive at the scheme for the original field u, where The extension of the finite-volume method to 2D and 3D can be obtained through the divergence theorem [2,260].

Discontinuous Galerkin methods
The discontinuous Galerkin method (DG) incorporates the concept of numerical fluxes across element interfaces from the finite-volume method into the finite-element framework. In particular, computations are done on reference element to increase computational efficiency. It is well-suited for parallelization due to the local character of the scheme and low amount of communication. Because continuity between elements is not required, the choice of spatial meshes is more flexible. On the other hand, these methods require more degrees of freedom than continuous methods. In the acoustic case, for instance, each node from an internal edge is associated with two or more local degrees of freedom, rather than a single global degree of freedom.
As finite-volume methods, DG was initially designed for transport problems [270,271]. Interior penalty methods, another class of discontinuous Galerkin approximations, were independently developed for elliptic and parabolic problems [272]. While the former is based on suitable approximation of fluxes across elements, the latter concerns weakly imposing continuity between them. We refer to Arnold et al. [273] for a unified framework for these approaches as well as other discontinuous Galerkin formulations.
The earliest DG methods for seismic wave propagation followed the interior-penalty approach [274][275][276][277][278]. They have been developed for second-order formulations. However, the most popular DG method in computational seismology is derived from first-order, conservation law formulations [74,279]. It employs piecewise high-order polynomial approximation, as spectral-element methods, together with the ADER time-integration approach. The resulting method achieves arbitrary high approximation order in space and time. In contrast with spectral elements, which rely on Lagrange interpolation functions, ADER-DG methods follow a modal rather than nodal approach [2], i.e., they rely on orthogonal polynomial basis functions, in particular for triangular and tetrahedral elements [280,281]. Hence, these methods benefit from automatic mesh generators of unstructured triangular and tetrahedral meshes, which are usually higher developed than those for quadrilateral and hexahedral elements, for example with optimized mesh partitioning techniques based on graph theory [262]. This allows, for instance, precise digital elevation models for the topography of the free surface.
A more recent version of DG, termed hybridizable discontinuous Galerkin (HDG), employs Lagrange multipliers over element boundaries and can provide a higher convergence rate [88,100,282]. As shown in [283], HDG is related to the earlier staggered discontinuous Galerkin methods [284][285][286].
There have been several applications of DG beyond the isotropic wave equation, such as anisotropic [287], viscoelastic [288], and poroelastic [289,290] waves. A unified Riemann solver has been recently proposed to couple these media [291]. Wilcox et al. [87] consider a velocity-strain (rather than velocity-stress) formulation, which allows coupling elastic and acoustic media. As pointed out by [2], DG is very well suited to dynamic-rupture problems [292][293][294][295] and has also benefited from the availability of open-source codes [126,296,297].
In order to compare spectral-element and discontinuous Galerkin methods [2], let us derive the elemental equations for a nodal upwind DG method for the 1D elastic wave system (100). Let e be an interior subinterval of the domain and Φ e 1 , . . . , Φ e n dof as in (79). We perform the scalar product of both sides of (100) with a test function Φ e i and integrate by parts in e , so that where x e L and x e R are the left and right endpoints of e . The latter terms are not present on standard finite-and spectral-element methods as they cancel out with contributions from adjacent intervals. In the discontinuous case, adjacent elements are no longer connected through continuity conditions, which are replaced with numerical fluxes. Similarly to (104), the undefined boundary terms of (106) can be approximated as follows (Fig. 4): By introducing numerical fluxes (107) and the expansion (79) into (106) , we find where the indices i(R) and i(L) satisfy Φ e+1 i (x e L ), respectively. Equation (108) is the same for the spectral-element method, except that the terms in the right-hand side are not present. Fig. 4 Upwind contributions of a scalar, discontinuous piecewise-polynomial fieldũ at the endpoints of element e , when the velocity is positive (a + ) or negative (a − ). Adapted from [2, p. 259]

Other methods
We close this section reviewing some families of numerical methods that are conceptually diverse from the ones presented above though their computational implementation may be developed from traditional spatial discretizations.

Physics-compatible numerical methods
Physics-compatible (also called mimetic or conservative) numerical methods are techniques that try to preserve (mimic) the fundamental physical and mathematical properties of continuous physics models in their finite-dimensional algebraic representations.
The numerical methods presented above, such as finite differences (FD), finite volumes (FV) and finite elements (FE), evolved separately until recently, but in recent years the need to develop more complex algorithms for solving new challenging real problems prompted the search for better and more robust schemes. Investigations and experience on the computational behavior of standard methods (stability, convergence, numerical errors, and efficiency) demonstrated that solving a physical problem by discrete models reproducing fundamental properties of the original continuum model allows for the best results. Among the important properties to be preserved are topology, conservation of energy, monotonicity, stability, maximum principles, symmetries, and involutions of continuum models. For this purpose, differential geometry, external calculus, and algebraic topology are the main mathematical tools for developing compatible discretizations.
Examples are compatible methods for spatial discretizations, variational and geometric integrators, or conservative finite-volume, finite-element and spectral-element methods, etc. The design principles for the development of mimetic discretization methods are described in books [298][299][300] and the references therein, while a general introduction and overview of spatial and temporal mimetic/geometric methods can be found in [301][302][303][304][305][306][307].
The general approach in developing compatible numerical schemes is to formulate the PDEs, which describe the continuum models, using invariant first-order differential operators, such as the divergence of vectors and tensors, the gradient of scalars and vectors, and the curl of vectors. The next step is to work out the compatible discretizations by using equivalent discrete forms of these invariant operators. The divergence, gradient, and curl differential operators satisfy certain integral identities (such as Green', Gauss', and Stokes' theorems) that are closely related to the conservation laws of the continuum models.
Therefore, the equivalent discrete forms of these integral identities are used in building the compatible discrete divergence, gradient, and curl operators since they must satisfy such discrete integral identities. Furthermore, other approaches have also been used, for example, based on algebraic topology, variational principles, or discrete vector calculus as well as for extending the mimetic approach to more general grids including polygonal, polyhedral, locally refined, and non-matching meshes.
For the sake of clarity, we show below the application of the basic principles of mimetic discretizations using the scalar wave equation (6) as described by Solano-Feo et al. [308].
Let us initially consider a one-dimensional grid with points x 0 , . . . , x N and uniform spacing h. We denote the mimetic approximations of a scalar function u at the grid points and their midpoints as u = u 0 , . . . , u N andū = u 1/2 , . . . , u N−1/2 , respectively. We may define the discrete divergence and gradient operators through central finite differences as which can be written in matrix form asv = Du and v = Gū. The discrete divergence and operators yield grid functions defined at midpoints and node values, respectively (Fig. 5). Left-and right-sided approximations should be employed to define the gradient operator at the grid endpoints, and two-and three-dimensional operators can be constructed with the aid of Kronecker products [299]. Solano-Feo et al. [308] pointed out that these operators satisfy the discrete integral identity where B is a boundary operator, u, v A = u Av denotes a discrete inner product with weighting matrix A, d = 1, 2, or 3 is the spatial dimension, Q (P) is the diagonal matrix containing quadrature weights of the compound midpoint (3/8 Newton-Cotes) rule, and I is the identity matrix.
Let us now proceed to the mimetic approximation to (6), considering the leapfrog scheme in time. By writing u = div(∇u), the Laplacian operator can be approximated by the compound discrete operator DG, leading to the explicit schemē whereā = c 2 (x 1/2 ), . . . , c 2 (x N−1/2 ). More details can be found in [299,308] and references therein. Mimetic principles have been applied to modeling wave-propagation problems by many authors [309][310][311][312]. Mimetic finite differences are particularly effective to handle topography and boundary conditions [313][314][315].

Cell method
It has been observed that many physical theories have a very similar formal structure from a geometric, algebraic, and analytical point of view. This principle led to the Tonti diagrams [316], a classification scheme of the physical quantities and the physical theories in which they are involved, such as equations of equilibrium, continuity and motion. These equations can be reformulated in a finite grid using basic concepts of algebraic topology such as completely discrete functions defined on a combination of elements of the grid rather than functions in the continuum. It is therefore possible to directly establish a set of algebraic relations between physical variables associated with the geometric elements of the problem and which are suitable for numerical simulations.
In practice, the CM at first accepts the idea of an approximate solution and focuses on single limited parts of the analyzed domain: the cells. After dividing the domain into cells (named as primal cell complex) a second subdivision is made, coupling a piece of each cell to each of its nodes. With this last subdivision, a domain area is attributed to each node of the primal cell complex, thus creating a second cell system, called the dual system (Fig. 6). In fact, there is full reciprocity (duality) between the geometric elements of the two systems of cells: a cell of the dual system (considered as "tributary region") remains connected to each node of the primal system; vice versa: to each node of the dual system correspond cells of the primal system. Therefore the geometric elements of the primal system (points P, lines L, surfaces S, and volumes V ) correspond to the geometric elements of the dual system (respectively volumesṼ , surfacesS, linesL, and pointsP).
The cell method, in addition to its simplicity, has close compatibility with the physical and experimental reality as a consequence of connecting the physical quantities to the geometric elements of the two cell systems with the same logic with which the quan- tities are investigated experimentally. The association between quantities involved in a physical problem and the geometric elements of the two cell complexes is illustrated and effectively summarized in the Tonti diagrams.
In summary, spatial and temporal quantities are represented by sets of topological entities (cells) of multiple dimensions called primal and dual cell complexes, and a system of inner and outer orientations are assigned to such cell complexes. The physical variables are associated with spatial and temporal elements according to the following classification: • Configuration variables: geometric and kinematic variables that describe the configuration of the (wave) field, such as displacement; • Source variables: static and dynamic variables that describe the sources of the field, such as force and mass flow; • Energy variables: variables that are obtained from the product of configuration and source variables, such as work and kinetic energy.
Configuration variables are associated with elements of the primal cell complex, while source/energy variables are associated with elements of the dual complex. The constitutive and balance laws are then imposed, leading to algebraic equations for the variables of interest.

Homogenization and multiscale methods
Two approaches have been proposed to handle medium heterogeneities that must be taken into account in the wave simulation, but would need a mesh refinement that is impractical to implement. These methods essentially convert the material properties in a fine scale, where relevant variability occurs, into equivalent ones in a coarse scale corresponding to the target wavelength. In general, physical laws may be different in fine and coarse scales [326].
One of these approaches can be seen as a generalization of averaging techniques [327] and lead to the homogenization methods [328][329][330][331][332]. Through asymptotic theory, they obtain effective equations at the macroscopic level that qualitatively account for the fine scales. The asymptotic expansion of order zero usually corresponds to the classical averaging techniques. Though these methods are not restricted to periodic media, they have been mostly developed for rectangular and cuboid grids [333][334][335].
The second approach obtains the effective medium by numerically modeling the fine scales [336]. There are several methods that follow this approach, such as numerical upscaling [337,338], heterogeneous multiscale method [339,340], multiscale finite elements [341], multiscale coupling methods [342], and Fast Fourier homogenization [343]. These techniques are especially useful to finely layered, randomly oriented, and fractured media [344,345].
Homogenization methods have the advantage of a lower computational cost. On the other hand, methods that numerically evaluate the contribution of fine scales to the macroscopic model can be more flexible with respect to the mesh geometry [346].

Numerical approximation of boundary conditions
This section concerns the implementation of the free-surface and computational boundary conditions mentioned in Sect. 2.4. They can also be handled directly at the discrete level, rather than being based on the discretization of analytical boundary conditions [29,347]. An important class of such methods use an artificial layer surrounding the domain to attenuate reflected waves.

Free-surface boundary conditions
Free-surface conditions can be easily imposed in numerical methods based on variational formulations, such as finite/spectral elements and discontinuous Galerkin methods. As an illustration, let us consider the variational formulation of the elastic wave equations (9) with boundary conditions where 1 ∪ 2 = ∂ and 1 denotes the surface boundary. in the space We have from divergence theorem that, for any w ∈ W , The first and second boundary integrals in the right-hand side of (114) vanish due to (112a) and (113), respectively, so that we arrive at the same variational formulation as in (37), except that w ∈ W . A similar result applies to the velocity-stress and velocity-displacement formulations.
Other methods such as finite differences need some strategy to handle z-derivatives present on condition (112a). This condition becomes, for instance in the 2D isotropic case, For a free-surface condition over z = 0, the x-derivatives can be approximated using grid points over this line. For the z-derivatives, one can extend the grid above z = 0 and impose the skew-symmetry of the stress components to evaluate the variables over these additional points [119,157], or to employ one-sided finite-difference expansions to avoid extending the grid [348]. Another approach is to use extrapolation based on characteristic variables [118], which has also been used on pseudospectral methods [68,107,349]. In the finite-volume method, the free-surface boundary condition may be imposed by solving an inverse Riemann problem [262].

Absorbing boundary conditions
The classical absorbing boundary conditions were initially implemented on finitedifference methods for the scalar wave equation, with one-sided difference formulas at the boundary [24,28,32]. Later on, these conditions were implemented on finitedifference methods for the two-and three-dimensional elastic wave equation [350,351].
High-order local absorbing boundary conditions have been mostly implemented with finite-difference and finite-element methods [357][358][359], but have also been considered in other methods [360,361].
It is worth noting that absorbing boundary conditions often involve first-order time derivatives, leading to second-order linear systems of ODEs (41) where the damping matrix C is present [98]. We refer to [36] for conditions for first-order hyperbolic problems.

Absorbing layers and PML
An alternative to designing non-reflecting boundary conditions is to extend the computational domain by surrounding with a layer where the wave field is subject to some form of filtering that attenuates the waves generated by reflection at the outer layer boundary (Fig. 7). This technique can be traced back to the works of Petschek and Hanson [362,363] and became popular in exploration geophysics after the method of Cerjan et al. [192]. The latter attenuates the numerical solution at the end of each time step by multiplication of a factor that tapers gradually towards the center of the grid [364], as suggested by the shading pattern in Fg. 7.
Rather than post-processing the wave field, wave attenuation may be obtained by adopting a modified governing equation in the absorbing layer, as proposed in  Sketch of an absorbing layer (shaded) surrounding a rectangular domain with a homogeneous grid. In the upper right, the idealized effect of the absorbing layer on waves reflected by the outer boundary 366]. For instance, the acoustic wave equation (7) can be modified in the absorbing layer as follows: where the parameter γ (x) is chosen to achieve the best amplitude elimination [364]. Sarma et al. [367] developed modified equations for finite-element methods under the framework of Rayleigh damping. A disadvantage of absorbing layers is that, while waves going through them may be effectively damped, spurious reflections occur at the interface between the domain and the absorbing layer. This limitation motivated the development of perfectly matched layers (PML), which were originally proposed to electromagnetic waves [368] and later extended to acoustic and elastic waves [44,[369][370][371].
Following [372], let us illustrate a PML for the scalar wave equation (6). Firstly, this equation is rewritten as a first-order system and is Laplace transformed: For simplicity, consider the layer portion adjacent to the boundary x = 0. In this part of the layer, system (117) is modified as follows: whereû j ( j = 1, 2, 3) are auxiliary variables such thatû =û 1 +û 2 +û 3 , while ω 1 is a function that vanishes along with its derivative at the interface (for instance, ω 1 (x 1 ) = Ax 1 ). Finally, the equations in time domain are obtained by applying the inverse Laplace transform to (118): In general, we split the derivative operators and the unknowns on components that are normal and tangential to the boundary and apply a complex change of variables in the normal direction [44,373]. Moreover, the modified equations may be obtained on second-order formulations of the wave equation [373].
Later on, the convolutional perfectly matched layer (CPML) was developed to avoid spurious reflections at grazing incidence. This method was originally proposed to the elastic wave equation in the velocity-stress formulation [374] and later extended to the displacement formulation [375] and to poroelastic [376] and viscoelastic [377] media.
Another approach that uses an absorbing layer applies high-order local NRBCs on two parallel artificial boundaries, and it is known as double-absorbing-boundary method [379]. This method has been evaluated in 2D and 3D seismic wave-propagation benchmark problems [380,381].

Numerical errors
Convergence analyses have been proposed for most of the fully discrete methods outlined in the previous section. In the following we list some of these works: • Finite-difference methods: [382,383].
However, convergence analysis usually does not guide the practitioner in the choice of discretization parameters for a wave-propagation simulation. Such an information is mostly provided by stability analysis (often a constituent part of convergence proofs [392]) and dispersion analysis.

Stability
Numerical stability, or the sensitivity of the numerical solution to perturbations, is an essential feature on numerical wave simulations, where such perturbations should not grow over time. The analysis of numerical stability of time-dependent problems is usually done through Von Neumann analysis [393], the matrix method [394], and the energy method [395].
Let us illustrate these analyses with the explicit finite-difference scheme (64) in the absence of source terms. In Von Neumann (or discrete Fourier) analysis, we represent the numerical solution in the form u n j ←û 0 e αt n e ikx j =û 0 ξ n e ikx j , and refer to the method as stable if the amplification factor ξ = u n+1 j /u n j = exp(α t) satisfies |ξ | ≤ 1 for any k [393], and unstable otherwise. By substituting (120) into (64), we obtain the following equation for the amplification factor [393,396]: The solutions ξ = A ± √ A 2 − 1 satisfy |ξ | ≤ 1 if the space and time steps x and t are chosen such that coinciding with the CFL stability criterion [112]. Similarly as in [396], the amplification factor of the implicit version of scheme (64) in the absence of sources, Therefore, the implicit scheme (123) is unconditionally stable, i.e., it is stable for any combination of the grid parameters x and t. For this reason, implicit schemes may be an attractive choice despite their higher computational cost. However, one must keep in mind that numerical dispersion, presented in the next section, constrains the choices of the grid parameters for both explicit and implicit methods.
A general form of this procedure can be found in many textbooks (e.g., [397]). The simplicity of Von Neumann analysis has made it the most frequently used tool for stability analysis of finite-difference methods as well as other techniques [80,277,398,399]. However, since the free-space solution in the form (120) is essentially limited to unbounded or periodic domains, boundary conditions are not taken into account. Moreover, the simplicity of Von Neumann analysis is limited to equations with constant coefficients.
While equations for heterogeneous media can be handled by considering the "frozen-coefficient" equations [400,401], the analysis of boundary conditions require alternative techniques. Even though Von Neumann conditions are sufficient in particular wave-propagation problems, in general the boundary conditions are stable only for a certain range of elastic parameters [402].
One of these alternatives is the matrix method, whose mostly well known source is the book by Mitchell [394] as well as other textbooks [401,403], though its main idea is present in earlier works [404][405][406]. The matrix method analysis is carried out in the physical domain (writing the equations in matrix form) rather than the wave number domain, and has been considered in [153,156,240,261].
Recalling that the matrix form of (64) with homogeneous boundary conditions is u n+1 − 2u n + u n−1 + t 2 Ku n = 0, with K given in (65), we have the following single where Q = 2I − t 2 K. It then follows that v n = A n v 0 ≤ A n v 0 , where · denotes a vector norm and its induced matrix norm (i.e., A = sup{ Av ; v = 1}. A necessary (but in general not sufficient) condition for the boundedness of v n is ρ(A) ≤ 1, where ρ(A) is the spectral radius of A.
A drawback of the matrix method is the need to analyze large matrices, and there are some approaches that alleviate the underlying computational cost. Ilan and Loewenthal [402] restricted the analysis to a portion of the domain close to the boundary. On the other hand, Kamel [408] proposed to seek the largest eigenvalue through power method [228], which is interpreted as updating an initial data through successive time steps.
As pointed out in [401,409], the condition ρ(A) ≤ 1 assures that A n remains bounded as n increases but, under this condition, A n may initially increase before decreasing. Griffiths et al. [409] point out that the condition A ≤ 1 is a sufficient one, and suggest an intermediate condition of the form A k ≤ C.
Let us proceed to the energy method, which seeks a discrete quadratic form that does not grow or moderately grows with time and in the same time bounds from above the discrete L 2 norm [410]. In the context of the previous examples, the discrete L 2 norm and inner product are Let us consider Eq. (1) in the absence of sources with boundary conditions u(a) = u(b) = 0. By multiplying both sides of (1) byu and integrating by parts over [a, b], It then follows that i.e., the quadratic form E(t) (the total energy) is constant over time. Note by writing c 2 = μ/ρ in analogy with (2) that the first and second terms in E(t) are associated with kinetic and potential energy, respectively. In general, E(t) may not be associated with a physical energy [410,411]. Analogously to (127), we multiply (64) by the centered-difference approximation (u n+1 j − u n−1 j )/(2 t) and sum from j = 1 to j = N − 1 arriving at the discrete energy conservation  (129), so does D + u n+1/2 h . The details of (129)- (130) are available in Sec. 9.2 of [412], where the general heterogeneous case is considered.
The energy method has been successfully used to analyze problems with freesurface [413], PML [414,415], and absorbing boundary conditions [28,37,379,411,416]. Besides finite differences, the analysis with energy method is present in finite [35] and spectral [417] elements, finite volumes [266], and discontinuous Galerkin methods [418]. On the other hand, since this method indirectly bounds some norm of the numerical solution, the stability condition is sufficient, but not necessary [410].
Finally, let us point out that other approaches are very well suited to the stability analysis of boundary conditions, such as the normal-mode (also known as GKS or GKSO) analysis [32,419,420] and the geometric stability condition [421,422].

Dispersion and numerical anisotropy
Dispersion analysis is an important tool for assessing the quality of approximation of numerical methods, providing an estimate of the minimum number of grid points per wavelength required to prevent waves from traveling with incorrect speed. A continuous or discrete wave model is dispersive if the wave speed depends on its wavelength.
Let us initially remain on the same 1D problem as in the previous section, recalling that the plane-wave solution (5) satisfies the dispersion relation ω = ±cκ, thus the phase and group velocities, coincide and are constant. On the other hand, if u n j = exp(−i(ω h t n − κ x j )) then scheme (64) in the absence of sources yields sin with r defined as in (121). It follows from (132) where H = κ x/(2π). Noting that exp(κ x) has period (wavelength) 2π/κ, we have that G = 2π/(κ x) = 1/H is the number of grid points per wavelength. If r = 1, then c ph h = c ph and c gr h = c gr , as long as H ≤ 1/2 (since sin x is not invertible for 0 ≤ x ≤ L if L > π/2). The bound H ≤ 1/2, or G ≥ 2, is known as the Nyquist limit. The result for r = 1 is exceptional and is not necessarily observed on less trivial problems [401]. When r < 1, numerical and exact phase/groups velocities do not coincide. The dispersion error is illustrated with r = 1/2 in Fig. 8. For instance, when H = 0.2 (i.e., G = 5 grid points per wavelength), the relative errors of phase and group velocity are approximately 5% and 15%, respectively.
A comprehensive review of the dispersion analysis of finite-difference schemes of higher order and spatial dimension is available in [412]. Trefethen presents in [423] a comprehensive review of group-velocity analysis of finite-difference schemes for the acoustic wave equation, plus a relationship between group velocity and GKS stability for first-order hyperbolic systems. Numerical dispersion has also been studied beyond the acoustic case [122,125,128,213].
Kosloff and Baysal [60] presented the numerical dispersion relations of 1D and 2D Fourier pseudospectral method using a similar procedure as above, while Fornberg [191] focused on the dispersion of the spatial discretization. Spa et al. [424] considered fully discrete schemes with Lax-Wendroff and rapid expansion methods. The numerical dispersion in the case of Chebyshev collocation points has not been studied in the classical works, though an analysis of its multidomain version has been recently proposed [425].
The numerical dispersion of finite-element methods of degree one can be done exactly as finite-difference methods; that is, to plug the discrete plane wave into the finite-element stencil assuming an infinite, periodic mesh (see, e.g., [213]). For 1D quadratic meshes and certain triangular meshes we must separate the nodes into sets which share the same degrees of freedom and are located at the same cyclically repeating location in the mesh pattern [426]. In this case, the numerical dispersion relation is expressed by an eigenvalue problem, whose solutions are analogous to the acoustic and optical branches from the theory of wave propagation into crystal structures [427,428]. Finite elements of higher degree lead to a larger number of solutions, and the classical interpretation is that only one eigenvalue is physically meaningful (in the case of the acoustic wave equation), while the others are regarded as spurious modes [227,429]. For this reason, the use of high-order finite element methods had been discouraged in numerical wave propagation.
Priolo and Seriani [69,234] performed a dispersion analysis of the 1D spectralelement method with Chebyshev collocation points by solving the discrete problem for a large final time, taking a wavelet as the initial condition and periodic boundary conditions. The final approximate and exact solutions are transformed into the Fourier space and the amplitude and the phase of their ratio is found for several wave numbers and degrees of polynomial approximation. The results are similar to the theoretical estimates presented in [191].
Mulder [430] applied the discrete Fourier transform sampled in the mesh nodes to the spatial operator and matched its eigenpairs with the transformed plane waves and their normalized wave numbers. Under this setting, the spurious modes provide reasonable approximations of particular eigenvectors of the exact operator. On the other hand, the spatial operator must be properly ordered to assure eigenpair matching. It is not trivial to find such an ordering for 2D or 3D problems.
A common practice in the dispersion analysis of spectral-element and discontinuous Galerkin methods is to select the eigenvalue mode that approximates the dispersion relation of the continuous wave equation [235,277,412,431], and locating these modes is also not trivial in general. Cohen et al. [432] use a Taylor series expansion. Abboud and Pinsky [433] writes the amplitude-variable discrete plane wave as a linear combination of discrete plane waves and classify the modes with the dominating coefficient of the combination (see also [434]). Seriani and Oliveira [435] identify these modes by a Rayleigh quotient approximation of the constant-amplitude mode. A similar analysis was done for the elastic wave equation [99,436]. The Rayleigh-quotient technique has also been employed in other Galerkin-type methods [225,437].
Another related form of error is numerical anisotropy, which is present when the speed of the approximate wave solution depends on the propagation direction in a different fashion than the exact solution's speed [397,438,439].
Let us illustrate numerical anisotropy with the two-dimensional version of the previous example. Let κ = κ {cos θ, sin θ } be a wave vector with magnitude κ and propagation direction given by the angle θ . It follows from the dispersion relation ω = ±cκ of the scalar wave equation (6) that the phase and group velocity do not depend on θ . On the other hand, by substituting u n j,k = exp(−i(ω h t n −κ(x j cos θ + y k sin θ ))) into the 2D version of (64) we find sin ω h t 2 = ± r 2 x sin 2 κ x cos θ 2 + r 2 y sin 2 κ y sin θ 2 , where r x = c t/ x and r y = c t/ y. Thus the numerical phase and group velocities will depend on θ , even if r x = r y = 1. The detailed study of numerical anisotropy of this scheme is available in [59], and the three-dimensional case readily follows by considering κ = κ {cos θ sin φ, sin θ sin φ, cos φ} . A convenient way to represent an angle-dependent dispersion relation is a polar diagram [397]. Similarly to [2], Fig. 9 shows percent phase-velocity errors of the 2D and 3D versions of the finite-difference scheme (64) in polar form.

Mass lumping and blending
Finite-element methods for second-order wave equations lead to systems of ODEs in the form (41) where the mass matrix M is usually non-diagonal. The mass-lumping technique approximates M by a diagonal matrixM, allowing the use of explicit-time stepping schemes. The classical approach is to row-lump the mass matrix [63,64,428], i.e.,M This concept has also been proposed to discretize the same integral formulations that lead to finite-volume methods [448,449].
In addition to the algebraic form (135), the approximate diagonal mass matrix may be obtained through reduced integration [450]. For high-order finite elements, the natural choice is to employ a Gauss-Lobatto-Legendre quadrature and shift the degrees of freedom so that they coincide with the quadrature nodes [239]. This procedure is followed in GLL spectral-element methods, as mentioned in Sect. 5.4. Mass lumping significantly affects the numerical dispersion of finite elements. As illustrated in Fig. 10, the numerical phase velocity is higher than the actual phase  Rayleigh-quotient dispersion analysis of optimal blended GLC spectral-element methods of degree N = 1, 2, 4, 8 using the leapfrog time integration CFL parameter is r = 0.5/N [452] velocity when a consistent matrix is used but lower when row lumping is employed. Usually the consistent mass matrix produces leading phase and group error, while the lumped mass matrix produces lagging phase and group error [451].
It is natural to seek a combination of lumped and consistent mass matrices that balances the over-and undershoots of these approaches, reducing numerical dispersion [213,428,452,453]. For instance, Fig. 11 shows the numerical phase velocity of the optimal blended operators [452] for the same example provided in Fig. 10. The dispersion remains lower than consistent and lumped elements until nearly the limit of π grid points per wavelength for Chebyshev collocation points [454].
An alternative approach is to seek the coefficients of the mass and stiffness matrices that minimize dispersion [254]. The search for optimal operators is well known in the context of finite-difference methods [455][456][457][458] and can be performed in a framework that is valid for most numerical methods [459].
should push for further improvement of these techniques, and community coding shortens the gap between theoretical advances and practical applications.
It is worth noting that some ideas developed for a method have been transferred to others. The concept of staggered grids from finite differences has been useful to pseudospectral methods, which in turn contributed back through curvilinear coordinates. Spectral elements have inspired discontinuous Galerkin methods to seek higher accuracy by using orthogonal polynomials, which is also a contribution from pseudospectral methods. Such an exchange corroborates the relevance of each family of methods to the overall progress of numerical modeling.