Modelling Seismic Wave Propagation for Geophysical Imaging

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Modelling Seismic Wave Propagation for Geophysical Imaging J. Virieux, V. Etienne, V. Cruz-Atienza, R. Brossier, E. Chaljub, O. Coutant, S. Garambois, D. Mercerat, V. Prieux, S. Operto, et al.


Introduction
The Earth is an heterogeneous complex media from the mineral composition scale ( 10 −6 m) to the global scale ( 10 6 m). The reconstruction of its structure is a quite challenging problem because sampling methodologies are mainly indirect as potential methods Rücker et al., 2006), diffusive methods (Cognon, 1971;Druskin & Knizhnerman, 1988;Goldman & Stover, 1983;Hohmann, 1988;Kuo & Cho, 1980;Oristaglio & Hohmann, 1984) or propagation methods (Alterman & Karal, 1968;Bolt & Smith, 1976;Dablain, 1986;Kelly et al., 1976;Levander, 1988;Marfurt, 1984;Virieux, 1986). Seismic waves belong to the last category. We shall concentrate in this chapter on the forward problem which will be at the heart of any inverse problem for imaging the Earth. The forward problem is dedicated to the estimation of seismic wavefields when one knows the medium properties while the inverse problem is devoted to the estimation of medium properties from recorded seismic wavefields.
The Earth is a translucid structure for seismic waves. As we mainly record seismic signals at the free surface, we need to consider effects of this free surface which may have a complex topography. High heterogeneities in the upper crust must be considered as well and essentially in the weathering layer zone which complicates dramatically the waveform and makes the focusing of the image more challenging.
Among the main methods for the computation of seismic wavefields, we shall describe some of them which are able to estime the entire seismic signal considering different approximations as acoustic or elastic, isotropic or anisotropic, and attenuating effects. Because we are interested in seismic imaging, one has to consider methods which should be efficient especially for the many-sources problem as thousands of sources are required for imaging. These sources could be active sources as explosions or earthquakes. We assume that their distribution are known spatially as punctual sources and that the source time function is the signal we need to reconstruct aside the medium properties.
Asymptotic methods based on the high frequency ansatz (see (Virieux & Lambaré, 2007) for references or textbooks (Červený, 2001;Chapman, 2004)) and spectral methods based on a spatial and time Fourier transformations (Aki & Richards, 2002;Cagniard, 1962;de Hoop, 1960;Wheeler & Sternberg, 1968) are efficient methods which are difficult to control: whispering galeries for flat layers are efficiently considered using spectral methods. These two methods may be used either for local interpretation of specific phases or as efficient alternatives when media is expected to be simple. They could be used as well for scattering inverse problems. In the general heterogeneous case, we have to deal with volumetric methods where the medium properties are described through a volume while seismic wave fields satisfy locally partial differential equations. Although one may consider boundaries as the free surface or the water/solid interface, we may consider that variations of the medium properties are continuous at the scale of the wavelength which we want to reconstruct: the best resolution we could expect is half the wavelength (Williamson & Worthington, 1993). Therefore a volumetric grid discretization is preferred where numerical expressions of boundary conditions should be mostly implicit through properties variations.
A quite popular method because of this apparent simplicity is the finite difference method where partial derivatives are transformed into finite difference expressions as soon as the medium has been discretized into nodes: discrete equations should be exactly verified. We shall consider first this method as it is an excellent introduction to numerical methods and related specific features. We will consider both time and frequency approaches as they have quite different behaviours when considering seismic imaging strategies.
Applications will enhance the different properties of this numerical tool and the caveats we must avoid for the various types of propagation we need.
Another well-known approach is the finite element method where partial differential equations are asked to be fulfilled in a average way (to be defined) inside elements paving the entire medium. We shall concentrate into the discontinuous Galerkin method as it allows to mix acoustic and elastic wave propagation into a same formalism: this particular method shares many features of finite element formalism when describing an element, but differs by the way these elements interact each other. We avoid the description of the continuous finite element method for compactness and differences will be pointed out when necessary. Again, we shall discuss both time-domain and frequency-domain approaches.
Applications will illustrate the different capabilities of this technique and we shall illustrate what are advantages and drawbacks compared to finite difference methods while specific features will be identified compared to continuous finite element methods.
We shall conclude on the strategy for seismic imaging when comparing precision of solutions and numerical efforts for both volumetric methods.

Equations of seismic wave propagation
In a heterogeneous continuum medium, seismic waves verify partial differential equations locally. Integral equations may provide an alternative for the evolution of seismic fields either in the entire domain or at the scale of an elementary element of a given mesh describing the medium structure.
Fundamental laws of dynamics require the conservation of linear and angular momentum in a Galilean reference frame. In the continuum, a force applied across a surface, oriented by the unit normal n at a given position x = (x, y, z) in a Cartesian coordinate system (O, x, y, z), by one side of the material on the other side defines the traction vector t i = σ ij n j where the second-rank stress tensor σ has been introduced. The conservation of the angular momentum makes the stress tensor symmetrical σ ij = σ ji . We shall introduce as well a volumetric force per unit mass at the current position denoted as f = ( f x , f y , f z ). The conservation of linear momentum allows one to write the acceleration of the displacement motion u(x) of a given particle at the current position as where the density is denoted by ρ(x).
Aside the translation and the rotation transformations preserving the distances inside the body we consider, the deformation of the continuum body is described by defining a strain tensor expressed as The symmetrical definition of the deformation ensures that no rigid-body rotations are included. The particle motion is decomposed into a translation, a rotation and a deformation: the two formers transformations preserve distances inside the solid body while the third one does not preserve distances, inducing stress variations inside the solid body. In the framework of linear elasticity, there is a general linear relation between the strain and stress tensors by introducing fourth-rank tensor c ijkl defined as follows Because of symmetry properties of stress and strain tensors, we have only 36 independent parameters among the 81 elastic coefficients while the positive strain energy leads to a further reduction to 21 independent parameters for a general anisotropic medium. For the particular case of isotropic media, we end up with two coefficients which can be the Lamé coefficients λ and μ. The second one is known also as the rigidity coefficient as it characterizes the mechanical shear mode of deformation. The following expression of elastic coefficients, with the Kronecker convention for δ ij gives the simplified expression linking the stress tensor to the deformation tensor for isotropic media as One may prefer the inverse of the relation (5) where the deformation tensor is expressed from the stress tensor by the introduction of the Young modulus E. Still, we have two independent coefficients. By injecting the relation (5) into the fundamental relation of dynamics (1), we end up with the so-called elastic wave propagation system, which is an hyperbolic system of second order, where only the displacement u has to be found. This system can be written as where we have neglected spatial variations of Lamé coefficients. Therefore, we must reconstruct over time the three components of the displacement or equivalently of the velocity or the acceleration. Choosing the stress is a matter of mechanical behaviour in a similar way for seismic instruments which record one of these fields.
For heterogeneous media, spatial differential rules for Lamé coefficients have to be designed. We shall see how to avoid this definition in the continuum by first considering hyperbolic system of first-order equations, keeping stress field. More generally, any hyperbolic equation with n-order derivatives could be transformed in a hyperbolic system with only first derivatives by adding additional unknown fields. This mathematical transformation comes naturally for the elastodynamic case by selecting the velocity field v and the stress field σ as fields we want to reconstruct. In a compact form, this first-order system in particle velocity and stresses is the following with i, j = x, y, z. We may consider other dual quantities as (displacement, integrated stress) or (acceleration, stress rate) as long as the medium is at rest before the dynamic evolution. Let us underline that time partial derivatives are on the left-hand side and that spatial variations and derivations are on the right-and side.
Using simple linear algebra manipulations, alternative equivalent expressions may deserve investigation: the three components σ ii could be linearly combined for three alternative components considering the trace σ 1 = trace(σ)/3, the x-coordinate deviatoric stress σ 2 = (2σ xx − σ yy − σ zz )/3 and the y-coordinate deviatoric stress σ 3 = (−σ xx + 2σ yy − σ zz )/3 which allows to separate partial spatial derivatives in the right hand side and material properties in the left hand side. The system (7) becomes 3 2μ which could be useful when we move from differential formulation to integral formulation over elementary volumes. Partial differential operators only in the right-hand side of the system (8) are separated from spatial variations of model parameters on the left-hand side as a diagonal matrix Λ = ( 3 3λ+2μ , 3 2μ , 3 2μ , 1 μ , 1 μ , 1 μ ). Similar strategies could be applied for 2D geometries.
Finally, for easing discussions on the numerical implementation, let us write both the 1D scalar second-order acoustic wave equation in the time domain as or, in frequency domain, away from sources where one can see the importance of the mixed operator ∂ x E(x)∂ x .
We have introduced the Young modulus E related to unidirectional compression/delation motion. The 1D vectorial first-order acoustic wave equation can be written as from which one can deduce immediatly the system of equations in the frequency domain Please note that the mixed operator does not appear explicitly. By discretizing this system and by eliminating the stress discrete values, one can go back to an equation involving only the velocity: a natural and systematic procedure for discretizing the mixed operator as proposed by Luo & Schuster (1990).
For an isotropic medium, two types of waves -compressional and shear waves -are propagating at two different velocities v p and v s . These velocities can be expressed as except for the 1D medium where only compression/dilatation motion could take place. The displacement induced by these two different waves is such that compressive waves u P verify ∇ × u P i = 0 and shear waves u S verify ∇ · u S i = 0. Applying these operators to the numerical displacement will separate it into these two wavefields.

Time-domain or frequency-domain approaches
These systems of equations could be solved numerically in the time domain or in the frequency domain depending on applications. For seismic imaging, the forward problem has to be solved for each source and at each iteration of the optimisation problem. The time approach has a computational complexity increasing linearly with the number of sources while precomputation could be achieved in the frequency domain before modelling the propagation of each source. Let us write a compact form in order to emphasize the time/frequency domains approaches. The elastodynamic equations are expressed as the following system of second-order hyperbolic equations, where M and S are the mass and the stiffness matrices (Marfurt, 1984). The source term is denoted by s and the seismic wavefield by w. In the acoustic approximation, w generally represents pressure, while in the elastic case, w generally represents horizontal and vertical particle displacements. The time is denoted by t and the spatial coordinates by x. Equation (14) is generally solved with an explicit time-marching algorithm: the value of the wavefield at a time step (n + 1) at a spatial position x is inferred from the value of the wavefields at previous time steps (Dablain, 1986;Tal-Ezer et al., 1990). Implicit time-marching algorithms are avoided as they require solving a linear system (Marfurt, 1984;Mufti, 1985). If both velocity and stress wavefields are helpful, the system of second-order equations can be recast as a first-order hyperbolic velocity-stress system by incorporating the necessary auxiliary variables (Virieux, 1986). The time-marching approach could gain in efficiency if one consider local time steps related to the coarsity of the spatial grid (Titarev & Toro, 2002): this leads to a quite challenging load balancing program between processors when doing parallel programming as most processors are waiting for the one which is doing the maximum of number crunching as illustrated for the ADER scheme (Dumbser & Käser, 2006). Adapting the distribution of the number of nodes to each processor depending on the expected complexity of mathematical operations is still an open problem. Other integration schemes as the Runge-Kutta scheme or the Stormer/Verlet symplectic scheme (Hairer et al., 2002) could be used as well.
Seismic imaging requires the cross-correlation in time domain or the product in frequency domain of the incidents field of one source and the backpropagated residues from the receivers for this source. In order to do so, one has to save at each point of the medium the incident field from the source which could be a time series or one complex number. The storage when considering a time-domain approach could be an issue: a possible strategy is storing only few time snapshots for recomputing the incident field on the fly (Symes, 2007) at intermediate times.
An additional advantage is that the attenuation effect could be introduced as well. in the time-domain approach, the complexity increases linearly with the number of sources.
In the frequency domain, the wave equation reduces to a system of linear equations, the right-hand side of which is the source, and the solution of which is the seismic wavefield. This system can be written compactly as where B is the so-called impedance matrix (Marfurt, 1984). The sparse complex-valued matrix B has a symmetric pattern, although is not symmetric because of absorbing boundary conditions (Hustedt et al., 2004;Operto et al., 2007). The fourier transform is defined with the following convention Solving the system of equations (15) can be performed through a decomposition of the matrix B, such as lower and upper (LU) triangular decomposition, which leads to the so-called direct-solver techniques. The advantage of the direct-solver approach is that, once the decomposition is performed, equation (15) is efficiently solved for multiple sources using forward and backward substitutions (Marfurt, 1984). This approach has been shown to be efficient for 2D forward problems (Hustedt et al., 2004;Jo et al., 1996;Stekl & Pratt, 1998). However, the time and memory complexities of the LU factorization, and its limited scalability on large-scale distributed memory platforms, prevents the use of the direct-solver approach for large-scale 3D problems (i.e., problems involving more than ten millions of unknowns) (Operto et al., 2007).
Iterative solvers provide an alternative approach for solving the time-harmonic wave equation (Erlangga & Herrmann, 2008;Plessix, 2007;Riyanti et al., 2006;. Iterative solvers are currently implemented with Krylov-subspace methods (Saad, 2003) that are preconditioned by the solution of the dampened time-harmonic wave equation. The solution of the dampened wave equation is computed with one cycle of a multigrid. The main advantage of the iterative approach is the low memory requirement, while the main drawback results from the difficulty to design an efficient preconditioner, because the impedance matrix is indefinite. To our knowledge, the extension to elastic wave equations still needs to be investigated. As for the time-domain approach, the time complexity of the iterative approach increases linearly with the number of sources or, equivalently, of right-hand sides, in contrast to the direct-solver approach.
An intermediate approach between the direct and the iterative methods consists of a hybrid direct-iterative approach that is based on a domain decomposition method and the Schur complement system (Saad, 2003;Sourbier et al., 2011): the iterative solver is used to solve the reduced Schur complement system, the solution of which is the wavefield at interface nodes between subdomains. The direct solver is used to factorize local impedance matrices that are assembled on each subdomain. Briefly, the hybrid approach provides a compromise in terms of memory saving and multi-source-simulation efficiency between the direct and the iterative approaches.
The last possible approach to compute monochromatic wavefields is to perform the modeling in the time domain and extract the frequency-domain solution, either by discrete Fourier transform in the loop over the time steps (Sirgue et al., 2008) or by phase-sensitivity detection once the steady-state regime has been reached (Nihei & Li, 2007). An arbitrary number of frequencies can be extracted within the loop over time steps at a minimal extra cost. Time windowing can easily be applied, which is not the case when the modeling is performed in the frequency domain. Time windowing allows the extraction of specific arrivals (early arrivals, reflections, PS converted waves) for the full waveform inversion (FWI), which is often useful to mitigate the nonlinearity of the inversion by judicious data preconditioning Sears et al., 2008).
Among all of these possible approaches, the iterative-solver approach has theoretically the best time complexity (here, complexity denotes how the computational cost of an algorithm grows with the size of the computational domain) if the number of iterations is independent of the frequency (Erlangga & Herrmann, 2008). In practice, the number of iterations generally increases linearly with frequency. In this case, the time complexity of the time-domain approach and the iterative-solver approach are equivalent (Plessix, 2007).
For one-frequency modeling, the reader is referred to those articles (Plessix, 2007;Virieux et al., 2009) for more detailed complexity analysis of seismic modeling based on different numerical approaches. A discussion on the pros and cons of time-domain versus frequency-domain seismic modeling relating to what it is required for full waveform inversion is also provided in Vigh & Starr (2008) and Warner et al. (2008).

Boundary conditions
In seismic exploration, two boundary conditions are implemented for wave modeling: absorbing boundary conditions to mimic an infinite medium and free surface conditions on the top side of the computational domain to represent the air-solid or air-water interfaces which have the highest impedance contrast. For internal boundaries, we assume that effects are well described by variations of the physical properties of the medium: the so-called implicit formulation (Kelly et al., 1976;Kummer & Behle, 1982).

PML absorbing boundary conditions
For simulations in an infinite medium, an absorbing boundary condition needs to be applied at the edges of the numerical model. An efficient way to mimic such an infinite medium can be achieved with Perfectly-Matched Layers (PML), which has been initially developed by Berenger (1994) for electromagnetics, and adapted for elastodynamics by Chew & Liu (1996); Festa & Vilotte (2005). PMLs are anisotropic absorbing layers that are added at the periphery of the numerical model. The classical PML formulation is based on splitting of the elastodynamic equations. A new kind of PML, known as CPML, does not require split terms. The CPML originated from Roden & Gedney (2000) for electromagnetics was applied by Komatitsch & Martin (2007) and Drossaert & Giannopoulos (2007) to the elastodynamic system. CPML is based on an idea of Kuzuoglu & Mittra (1996), who has obtained a strictly causal form of PML by adding some parameters in the standard damping function of Berenger (1994), which enhanced the absorption of waves arriving at the boundaries of the model with grazing incidence angles.
In the frequency domain, the implementation of PMLs consists of expressing the wave equation in a new system of complex-valued coordinatesx defined by (e.g., Chew & Weedon, 1994): In the PML layers, the damped 1D acoustic wave equation could be deduced from the equation (10) as where ξ x (x) = 1 + iγ x (x)/ω and γ x (x) is a 1D damping function which defines the PML damping behavior in the PML layers. In the CPML layers, the damping function ξ x (x) becomes with angular frequency ω and coefficients κ x ≥ 1 and α x ≥ 0. The damping profile d x varies from 0 at the entrance of the layer, up to a maximum real value d θ max at the end (Collino & Tsogka, 2001) such that and with δ x as the depth of the element barycentre inside the CPML, L cpml the thickness of the absorbing layer, and R coe f f the theoretical reflection coefficient. Suitable expressions for κ x , d x and α x are discussed in Collino & Monk (1998); Collino & Tsogka (2001); Drossaert & Giannopoulos (2007); Komatitsch & Martin (2007); Kuzuoglu & Mittra (1996); Roden & Gedney (2000). We often choose R coe f f = 0.1% and the variation of the coefficient α x goes from a maximum value (α x max = π f 0 ) at the entrance of the CPML, to zero at its end. If κ x = 1 and α x = 0, the classical PML formulation is obtained.
One can use directly these frequency-dependent expressions when considering the frequency approach. The formulation in the time domain is slightly more involved. The spatial derivatives are replaced by with where H(t) denotes the Heaviside distribution. Roden & Gedney (2000) have demonstrated that the time convolution in equation (21) can be performed in a recursive way using memory variables defined by The function ψ x represents a memory variable in the sense that it is updated at each time step. Komatitsch & Martin (2007) have shown that the term κ x has a negligible effect on the absorbing abilities, and it can be set to 1. If we take κ x = 1, we derive the equation (23) using the equation (22) as One equation is generated for each spatial derivative involved in the elastodynamic system, which can be a memory-demanding task. Once they are computed at each time step, we can introduce the memory variables into the initial elastodynamic system which requires two additional variables for the 1D equations (11) with the definition of ψ x (v) and ψ x (σ) leading to the following system At the outer edge of the PML zone, one could apply any conditions as simple absorbing conditions (Clayton & Engquist, 1977) or free surface conditions  as fields go to zero nearby the outer edge.
We must underline that the extension to 2D and 3D geometries is straightforward both in the frequency domain  and in the time domain Komatitsch & Martin, 2007).

Free surface
Planar free surface boundary conditions can be simply implemented using a strong formulation or a weak formulation.
In the first case which is often met in the finite-difference methods, one requires that the stress is zero at the free surface. The free surface matches the top side of the FD grid and the stress is forced to zero on the free surface (Gottschamer & Olsen, 2001). Alternatively, the method of image can be used to implement the free surface along a virtual plane located half a grid interval above the topside of the FD grid (Virieux, 1986). The stress is forced to vanish at the free surface by using a virtual plane located half a grid interval above the free surface where the stress is forced to have opposite values to that located just below the free surface. In case of more complex topographies, one strategy is to adapt the topography to the grid structure at the expense of numerical dispersion effect (Robertsson, 1996) or to deform the underlying meshing used in the numerical method to the topography (Hestholm, 1999;Hestholm & Ruud, 1998;Tessmer et al., 1992). In the first case, because of stair-case approximation, a local fine sampling is required (Hayashi et al., 2001).
Owing to the weak formulation used in finite-element methods, the free surface boundary condition are naturally implemented by zeroing test functions on these boundaries which follow edges of grid elements (Zienkiewicz & Taylor, 1967). This approach could be used as well for finite-difference methods through the summation-by-parts (SBP) operators based on energy minimization combined with Simultaneous Approximation Term (SAT) formulation based on a boundary value penalty method (Taflove & Hagness, 2000). In this case, boundaries on which stress should be zero are not requested to follow any grid discretisation.
Finally, one may consider an immersed boundary approach where the free surface boundary is not related to the discretisation of the medium as promoted by LeVeque (1986). Extensive applications have been proposed by Lombard & Piraux (2004); Lombard et al. (2008) where grid discretisation does not influence the application of boundary conditions. This approach might be seen as an extension of the method of images following extrapolation techniques above the free surface at any a priori order of precision.

Source implementation
There are different ways of exciting the numerical grid by the source term. The simplest one is the direct contribution of the source term in the discrete partial differential equations: for example, we may just increment by the source term after each time step or we may consider the right-hand side source term for solving the linear system in the frequency domain. Depending on the numerical approach, it is necessary to consider specific influences coming from the discretisation as we shall see for numerical methods we consider.
In order to avoid singularities of solutions nearby the source, on can use the injection technique as proposed in the pioneering work of Alterman & Karal (1968) where a specific box around the source is defined. Inside the box, only the scattering field is computed. The incident field is estimated at the edges of the box and it is substracted when propagations are estimated inside the box and added when propagations are estimated outside the box. A more general framework is proposed by Opršal et al. (2009) related to boundary integral approaches.

Finite-difference methods: solving the equations through a strong formulation
We shall first consider the discretization based on simple and intuitive approaches as finite-difference methods for solving these partial differential equations while focusing on techniques useful for seismic imaging which means a significant number of forward problems for many sources in the same medium. We shall identify features which might be interesting for seismic imaging. If these approaches are intuitive for solving differential equations, the numerical implementation of boundary conditions and source excitation is less obvious and requires specific strategies as we shall see.

Spatial-domain finite-difference approximations
Whatever is our strategy for the reconstruction of the wave field u, one has to discretize it. We may be very satisfied by considering a set of discrete values (u 1 , u 2 , ..., u I−1 , u I ) along one direction at a given specific time t n which can be discretized as well. Therefore, a simple way of solving this first-order differential system is by making finite difference approximations of spatial derivatives.
Still considering a 1D geometry, the partial operator (∂/∂x) could be deduced from a Taylor expansion using Lagrange polynomial. A quite fashionable symmetrical estimation using a centered finite difference approximation is expressed as which is a three-nodes stencil as three nodes are involved: two for the derivative estimation and one for the updating. Let us mention that the discrete derivative is shifted with respect to discrete values of the field. Because of the very specific antisymmetrical structure of our first-order hyperbolic system where time evolution of velocity requires only stress derivatives (and vice versa), we may consider centered approximations both in space and in time. This will lead to a leap-frog structure or a red/black pattern. Of course, we have truncation errors expressed by the function O(Δx n ) which depends on the power n of the spatial stepping and by the function O(Δt k ) on the power k of the time stepping.
We may require a greater precision of the derivative operator by using more points for this partial derivative approximation and a very popular centered finite difference approximation of the first derivative is the following expression where c 1 = 9/8 et c 2 = −1/24 (Levander, 1988). This fourth-order stencil is compact enough (few discrete points inside the stencil) for numerical efficiency while having a small local truncation error. This stencil is a five-nodes stencil. Let us underline that centered approximations lead to have field quantities not at the same position in the numerical grid as derivative approximations (figure 1). In other words, stress and velocity components should be specified on different positions of the spatial grid. If we still consider a full grid where stress and velocities are known at the same position, this stencil could be recast as a seven-nodes where Δx = Δx/2 and where we have following specific coefficients d 1 = c 1 , d 2 = 0. and d 3 = c 2 . The fourth-order scheme would require the following theoretical coefficients d 1 = 15/20, d 2 = −3/20 and d 3 = 1/60. For fourth-order stencils, the two sub-grids are not entirely decoupled and are weakly coupled leading to a dispersion behaviour as if the discretization is Δx. Let us remind that these sub-grids are completely decoupled when considering second-order stencils, leading to the staggered structure. Therefore, solving partial differential equations in the staggered grid structure has a less accurate resolution but improves significantly the efficiency of the method than solving equations in the full grid even with dispersion-relation-preserving stencils (Tam & Webb, 1993). The memory saving can be easily seen when comparing nodes for staggered grid and nodes for full grid (figure 1) When dealing with 2D and 3D geometries, we may exploit the extra freedom and estimate derivatives along the direction x from nodes shifted by half the grid step in x but also by half the grid step in y (and eventually in z). This leads to another compact stencil as shown in the figure 1 where all components of the velocity are discretized in one location while all components of the stress field are discretized half the diagonal of the grid as proposed by Saenger et al. (2000). This grid is still partially staggered and could be named as a partial grid.
These standard and partial staggered structures are sub-grids of the full grid as shown in the figure 1 which is used in aeroacoustics (Tam & Webb, 1993).
These different discretizations related to various stencils may lead to preferential directions of propagation. Numerical anisotropy effect is observed even when considering isotropic wave propagation. The figure 2 shows error variations in velocities with respect to angles of propagation for the standard grid and the partial grid: one can see that the anisotropy behavior is completely different with a rotation shift of 45 0 . In 2D, the two stencils provide the same anisotropic error while the partial grid has a slightly improved numerical anisotropic performance (percentage differences go from 3 % down to 2 % in 3D geometries). Of course, the spatial sampling is such that the error should be negligible and few percentages is considered to be acceptable except nearby the source. Other spatial interpolations are possible. Previous discrete expressions are based on Lagrange interpolations while other interpolations are possible such as Chebychev or Laguerre polynomial or Fourier interpolations (Kosloff & Baysal, 1982;Kosloff et al., 1990;Mikhailenko et al., 2003). Interpolation basis could be local (Lagrange) or global(Fourier) ones based on equally spaced nodes or judiciously distributed nodes for keeping interpolation errors as small as possible: this will have a dramatic impact on the accuracy of the numerical estimation of the derivative and, therefore, on the resolution of partial differential equations. We should stress that local stencils should be preferred for seismic imaging for efficiency in the computation of the forward modeling.

Time-domain finite-difference approximations
Similarly, one may consider finite difference approximation for time derivatives which can be illustrated on the simple scalar wave equation. A widely used strategy is again the centered differences through the expression For understanding how the procedure of computing new values in time is performing, let us consider the simple 1D second-order scalar wave equation for displacement u. This equation could be discretized through these finite difference approximations The next value at the discrete time n + 1 comes from older values known at time n and time n − 1 through the expression A more compact notation of this equation as shows the quantity known as the Courant number in the literature. This quantity is quite important for understanding the numerical dispersion and stability of finite difference schemes. The related stencil on the spatio-temporal grid as shown in the left panel of the figure 3 clearly illustrates that the value at time n + 1 is explicitly computed from values at time n − 1 and time n. In this explicit formulation, the selection of the time step Δt should verify that the Courant number is lower than 1 for any point of the medium. On the contrary, we may consider spatial derivatives at time n + 1. This leads us to an implicit scheme where more than one value at time n + 1 is present in the discretisation. The equation which can be described by the Courant number S through the equation The right panel of the figure 3 illustrates the structure of the stencil and that three unknowns have to be estimated through the single equation (35). By considering different spatial nodes, we may find these three unknowns by solving a linear system. The Courant number could take any value for time integration as long as discrete sampling is correctly performed.
Other implicit stencils might be designed by averaging the spatial derivatives over the three times n − 1, n and n + 1. We may as well average the time derivative over the three positions i − 1, i and i + 1. This lead to another linear system to be solved. These weighting strategies could be designed for reducing numerical noise as numerical dispersion and/or anisotropy: a road for further improvements.
As discretisation in space and time goes to zero, one expects the solution to be more precise but cumulative rounding errors should prevent to have too small values. In expressions (26) and (29), truncation error O[Δx 2 ] goes to zero as the square of the discrete increment. We shall say that this is a second-order precision scheme both in space and in time. One consider often the stencil with the fourth-order precision in space and second-order precision in time, denoted as O Δx 4 , Δt 2 , as an optimal one for finite-differences simulations.

Frequency-domain finite-difference approximations
The second-order acoustic equation (10) provides a generalization of the Helmholtz equation.
In exploration seismology, the source is generally a local point source corresponding to an explosion or a vertical force.
Attenuation effects of arbitrary complexity can be easily implemented in equations (10) and (12) using complex-valued wave speeds in the expression of the bulk modulus, thanks to the correspondence principle transforming time convolution into products in the frequency domain: in the frequency domain, one has to replace elastic coefficients by corresponding viscoelastic complex moduli for considering visco-elastic behaviors (Bland, 1960). For example, according to the Kolsky-Futterman model (Futterman, 1962;Kolsky, 1956), the complex wave speedc is given bȳ where the P wave speed is denoted by c = E/ρ, the attenuation factor by Q and a reference frequency by ω r . The function sgn gives the sign of the function.
Since the relationship between the wavefields and the source terms is linear in the first-order and second-order wave equations, one can explicitely expressed the matrix structure of equations (10) and (12) through the compact expression, where M is the mass matrix, S is the complex stiffness/damping matrix. This expression holds as well in 2D and 3D geometries. The dimension of the square matrix B is the number of nodes in the computational domain multiplied by the number of wavefield components. System (37) could be solved using a sparse direct solver. A direct solver performs first a LU decomposition of B followed by forward and backward substitutions for the solutions (Duff et al., 1986) as shown by the following equations: Exploration seismology requires to perform seismic modeling for a large number of sources, typically, up to few thousands for 3D acquisition. The use of direct solver is the efficient computation of the solutions of the system (37) for multiple sources. Combining different stencils for constructing a compact and accurate stencil can follow strategies developped for acoustic and elastic wave propagation (Jo et al., 1996;Operto et al., 2007;Stekl & Pratt, 1998). The numerical anisotropy is dramatically reduced The mass matrix M is a diagonal matrix although never explicitly constructed when considering explicit time integration. In the frequency domain formulation, we may spread out the distribution of mass matrix over neighboring nodes in order to increase the precision without increasing the computer cost as we have to solve a linear system in all cases. This strategy is opposite to the finite element approach where often the mass matrix is lumped into a diagonal matrix for explicit time integration (Marfurt, 1984). For a frequency formulation, considering the mass matrix as a non-diagonal matrix does not harm the solver. The weights of distribution are obtained through minimization of the phase velocity dispersion in an infinite homogeneous medium Jo et al., 1996): the numerical dispersion is dramatically reduced.

PML absorbing boundary condition implementation
Implementation of PML conditions in the frequency domain is straightforward using unsplit variables while, in the time domain, we need to introduce additional variables for handling the convolution through memory variables or to use split unphysical field variables (Cruz-Atienza, 2006). These additional variables are only necessary in the boundary layers following the figure 4 We first consider an infinite homogeneous medium which is embedded into a cubic box of a 16 km size and a grid stepping of h = 100 m. The thickness of the PML layer is 1 km leading to nsp = ten nodes inside the PML zone. The P-wave velocity is 4000 m/s while the S-wave velocity is 2300 m/s and the density 2500 kg/m 3 . The figure 5 shows various time sections of the 3D volume for the vertical particle velocity where one can see that the explosive wavefront is entering the PML zone at the time 2.8 s. The last two snapshots shows the vanishing of the wavefront with completely negligible residues at the final time (the decrease of the elastic energy is better than 0.2 % for ten nodes and could reach 0.03 % for twenty nodes).
When we have discontinuous interfaces crossing the PML zone, we may expect difficulties coming from various angles of propagation waves (Chew & Liu, 1996;Festa & Nielsen, 2003;Marcinkovich & Olsen, 2003). Therefore, a simple heterogeneous medium is considered   figure 5 shows that, in spite of the complexities of waves generated at the horizontal flat interface, the PML layer succeeds to absorb seismic energy with a residual energy of 0.3 % in this case, still far better than standard paraxial absorbing boundary conditions (Clayton & Engquist, 1977).

Source and receiver implementation on coarse grids
Seismic imaging by full waveform inversion is initiated at an initial frequency as small as possible to mitigate the non linearity of the inverse problem resulting from the use of local optimization approach such as gradient methods. The starting frequency for modeling in exploration seismics can be as small as 2 Hz which can lead to grid intervals as large as 200 m.
In this framework, accurate implementation of point source at arbitrary position in a coarse grid is critical. One method has been proposed by Hicks (2002) where the point source is approximated by a windowed Sinc function. The Sinc function is defined by where x = (x g − x s ), x g denotes the position of the grid nodes and x s denotes the position of the source. The Sinc function is tapered with a Kaiser function to limit its spatial support (Hicks, 2002) . For multidimensional simulations, the interpolation function is built by tensor product construction of 1D windowed Sinc functions. If the source positions matches the position of one grid node, the Sinc function reduces to a Dirac function at the source position and no approximation is used for the source positioning. If the spatial support of the Sinc function intersects a free surface, part of the Sinc function located above the free surface is mirrored into the computational domain with a reverse sign following the method of image. Vertical force can be implemented in a straightforward way by replacing the Sinc function by its vertical derivative. The same interpolation function can be used for the extraction of the pressure wavefield at arbitrary receiver positions. The accuracy of the method of Hicks (2002) is illustrated in Figure 6a which shows a 3.75 Hz monochromatic wavefield computed in a homogeneous half space. The wave speed is 1500 m/s and the density is 1000 kg/m 3 . The grid interval is 100 m. The free surface is half a grid interval above the top of the FD grid and the method of image is used to implement the free surface boundary condition. The source is in the middle of the FD cell at 2 km depth. The receiver line is oriented in the Y direction. Receivers are in the middle of the FD cell in the horizontal plane and at a depth of 6 m just below the free surface. Comparison between the numerical and the analytical solutions at the receiver positions are first shown when the source is positioned at the closest grid point and the numerical solutions are extracted at the closest grid point (Figure 6b). The amplitude of the numerical solution is strongly overestimated because the numerical solution is extracted at a depth of 50 m below free surface (where the pressure vanishes) instead of 6 m. Second, a significant phase shift between numerical and analytical solutions results from the approximate positioning of the sources and receivers. In contrast, a good agreement between the numerical and analytical solutions both in terms of amplitude and phase is shown in Figure 6c where the source and receiver positioning is implemented with the windowed Sinc interpolation.

Realistic examples for acoustic and elastic propagations using FD formulations
We shall provide two simple examples of seismic modeling using finite-differences methods both in the frequency and time approaches. The first example concerns seismic exploration problem where the acoustic approximation is often used while the second one is related to seismic risk mitigation where free surface effects including elastic propagation are quite important.

3D EAGE/SEG salt model
The  (2000) and a small number of processors (48) are used (Table 1). For a number of processors equal to the number of RHS members, the two approaches have the same cost. Of note, in the latter configuration (N P =N rhs ), the cost of the two methods is almost equal in the case of the salt model (0.94 h versus 0.816 h).
Over the last decades, simulations of wave propagation in complex media have been efficiently tackled with finite-difference methods (FDMs) and applied with success to numerous physical problems (Graves, 1996;Moczo et al., 2007). Nevertheless, FDMs suffer from some critical issues that are inherent to the underlying Cartesian grid, such as parasite diffractions in cases where the boundaries have a complex topography. To reduce these artefacts, the discretisation should be fine enough to reduce the 'stair-case' effect at the free surface. For instance, a second-order rotated FDM requires up to 60 grid points per wavelength to compute an accurate seismic wavefield in elastic media with a complex topography (Bohlen & Saenger, 2006). Such constraints on the discretisation drastically restrict the possible field of realistic applications. Some interesting combinations of FDMs and finite-element methods (FEMs) might overcome these limitations (Galis et al., 2008). The idea is to use an unstructured FEM scheme to represent both the topography and the shallow part of the medium, and to adopt for the rest of the model a classical FDM regular grid. For the same reasons as the issues related to the topography, uniform grids are not suitable for highly  (right) processors. The number of sources is 2000. Pre. denotes the elapsed time for the source-independent task during seismic modeling (i.e., the LU factorization in the frequency-domain approach). Sol. denotes the elapsed time for multi-RHS solutions during seismic modeling (i.e., the substitutions in the frequency-domain approach). Fig. 8. On the left, the French Riviera medium with complex topography and bathymetry: an hypothetical earthquake of magnitude 4.5 is at a depth of 10 km below the epicenter shown by a red ball. The simulation medium is 20 km by 20 km by 15 km. On the right, the related peak ground acceleration (PGA). Please note that the acceleration is always lower than one tenth of the Earth acceleration g heterogeneous media, since the grid size is determined by the shortest wavelength. Except in some circumstances, like mixing grids (Aoi & Fujiwara, 1999) or using non uniform Cartesian grids (Pitarka, 1999) in the case of a low velocity layer, it is almost impossible to locally adapt the grid size to the medium properties in the general case. From this point of view, FEMs are appealing, since they can use unstructured grids or meshes. Due to ever-increasing computational power, these kinds of methods have been the focus of a lot of interest and have been used intensively in seismology (Aagaard et al., 2001;Akcelik et al., 2003;Ichimura et al., 2007).

PGA estimation in the French Riviera
Peak ground acceleration (PGA) are estimated using empirical attenuation laws calibrated through databases of seismic records of various areas: these laws should be adapted to each area around the world and European moderate earthquakes require a specific calibration (Berge-Thierry et al., 2003). Aside these attenuation laws, numerical tools as finite-differences time-domain methods allows the deterministic estimation of the peak ground acceleration (PGA) in specific areas of interest once the medium is known and the source specified.
Small areas as the French Riviera where a complex topography and bathymetric makes the simulation difficult. We would like to illustrate the procedure of time-domain simulation on this specific example ). The Figure 8 shows a very simple model surrounding the city of Nice: the box is 20 km by 20 km by 15 km in depth. The P-wave velocity is 5700 m/s while the S-wave velocity is 3300 m/s and the density 2600 km/m 3 . The water is characterized by a P-wave velocity of 1530 m/s while the density is about 1030 km/m 3 . The grid step is 50 m and the time integration step is 0.005 s.
The numerical simulation of an hypothetical earthquake of magnitude 4.5 at a depth of 10 km in the Mediterranean Sea provides us a deterministic estimation of the PGA as shown in the Figure 8. This small source is characterized upto a frequency of 3 Hz and we select a source time function with this expected spectral content.
Successful applications have been proposed in the Los Angeles basin and is improved as we increase our knowledge about the medium of propagation and about the source location and its characterization. The PGA is estimated everywhere and one can see that increase of the PGA is observed at the sea bottom and nearby the coast. One can show that the amplification of PGA is decreased when considering the water layer at the expense of a longer duration of the seismic signal.
Of course, various simulations should be performed using different models of the medium and for various source scenarii. These simulations could help to assess the variability of the acceleration for possible potential earthquakes and may be used for the mitigation of seismic risks. The importance of constraining the model structure should be emphasized and we can accumulate this knowledge through various and different initiatives performed for a more accurate reconstruction of the velocity structure (Rollet et al., 2002). One tool is the seismic imaging procedure we have underlined in this chapter.

Finite-elements discontinuous Galerkin methods: a weak formulation
Finite element methods, often more intensive in computer resources, introduce naturally boundary conditions in an explicit manner. Therefore, we expect improved accurate solutions with this numerical approach at the expense of computer requirements. The system of equations (14) in time has now a non-diagonal mass matrix while the system of equations (15) has a impedance matrix particularly ill-conditioned in 3D geometry taking into account its dimensionality. Therefore, for 2D geometries, the frequency formulation is still a quite feasible option while time domain approaches are there appealing when considering 3D geometries. Due to ever-increasing computational power, finite element methods using unstructured meshes have been the focus of increased interest and have been used extensively in seismology (Aagaard et al., 2001;Akcelik et al., 2003;Ichimura et al., 2007).
Usually, the approximation order remains low, due to the prohibitive computational cost related to a non-diagonal mass matrix. However, this high computational cost can be avoided by mass lumping, a standard technique that replaces the large linear system by a diagonal matrix (Chin-Joe- Kong et al., 1999;Marfurt, 1984) and leads to an explicit time integration. Another class of FEMs that relies on the Gauss-Lobatto-Legendre quadrature points has removed these limitations, and allows for spectral convergence with high approximation orders. This high-order FEM, called the spectral element method (SEM) (Komatitsch & Vilotte, 1998;Seriani & Priolo, 1994) has been applied to large-scale geological models up to the global scale (Chaljub et al., 2007;Komatitsch et al., 2008). The major limitation of SEM is the exclusive use of hexahedral meshes, which makes the design of an optimal mesh cumbersome in contrast to the flexibility offered by tetrahedral meshes. With tetrahedral meshes (Frey & George, 2008), it is possible to fit almost perfectly complex topographies or geological discontinuities and the mesh width can be adapted locally to the medium properties (h-adaptivity). The extension of the SEM to tetrahedral elements represents ongoing work, while some studies have been done in two dimensions on triangular meshes (Mercerat et al., 2006;Pasquetti & Rapetti, 2006). On the other hand, another kind of FEM has been proven to give accurate results on tetrahedral meshes: the Discontinuous Galerkin finite-element method (DG-FEM) in combination with the arbitrary high-order derivatives (ADER) time integration (Dumbser & Käser, 2006). Originally, DG-FEM has been developed for the neutron transport equation (Reed & Hill, 1973). It has been applied to a wide range of applications such as electromagnetics (Cockburn et al., 2004), aeroacoustics (Toulopoulos & Ekaterinaris, 2006) and plasma physics (Jacobs & Hesthaven, 2006), just to cite a few examples. This method relies on the exchange of numerical fluxes between adjacent elements. Contrary to classical FEMs, no continuity of the basis functions is imposed between elements and, therefore, the method supports discontinuities in the seismic wavefield as in the case of a fluid/solid interface. In such cases, the DG-FEM allows the same equation to be used for both the elastic and the acoustic media, and it does not require any explicit conditions on the interface (Käser & Dumbser, 2008), which is, on the contrary, mandatory for continuous formulations, like the SEM (Chaljub et al., 2003). Moreover, the DG-FEM is completely local, which means that elements do not share their nodal values, contrary to conventional continuous FEM. Local operators make the method suitable for parallelisation and allow for the mixing of different approximation orders (p-adaptivity).

3D finite-element discontinuous Galerkin method in the time domain
Time domain approaches are quite attractive when considering explicit time integration. However, in most studies, the DG-FEM is generally used with high approximation orders. We present a low-order DG-FEM formulation with the convolutional perfectly matched layer (CPML) absorbing boundary condition (Komatitsch & Martin, 2007;Roden & Gedney, 2000) that is suitable for large-scale three-dimensional (3D) seismic wave simulations. In this context, the DG-FEM provides major benefits.
The p-adaptivity is crucial for efficient simulations, in order to mitigate the effects of the very small elements that are generally encountered in refined tetrahedral meshes. Indeed, the p-adaptivity allows an optimised time stepping to be achieved, by adapting the approximation order according to the size of the elements and the properties of the medium. The benefit of such a numerical scheme is particularly important with strongly heterogeneous media. Due to the mathematical formulation we consider, the medium properties are assumed to be constant per element. Therefore, meshes have to be designed in such a way that this assumption is compatible with the expected accuracy. The discretization must be able to represent the geological structures fairly well, without over-sampling, while the spatial resolution of the imaging process puts constraints on the coarsest parameterisation of the medium. If we consider full waveform inversion (FWI) applications, the expected imaging resolution reaches half a wavelength, as shown by Sirgue & Pratt (2004). Therefore, following the Shannon theorem, a minimum number of four points per wavelength is required to obtain such accuracy. These reasons have motivated the development of DG-FEM with low orders. We focus on the quadratic interpolation, which yields a good compromise between accuracy, discretisation and computational cost.

3D time-domain elastodynamics
It is worth to provide notations for specific manipulation of equations for DG-FEM approaches. The first-order hyperbolic system (8) under the so-called pseudo-conservative form can be written following the approach of Ben Jemaa et al. (2007) as with the definitions of the velocity and stress vectors as v t = (v x v y v z ) t and σ = (σ 1 σ 2 σ 3 σ xy σ xz σ yz ) t . Under this pseudo-conservative form, the RHS of (41) does not include any term that relates to the physical properties. The diagonal matrix Λ has been introduced in the system (8) and its inverse is required for the computation of the stress components (equation (41)). Matrices M θ and N θ are constant real matrices . The extension of the pseudo-conservative form for the visco-elastic cases could be considered with the inclusion of memory variables while the anisotropic case should be further analysed since the change of variable may depend on the physical parameters. Finally, in the equation (41), the medium density is denoted by ρ, while f and σ 0 are the external forces and the initial stresses, respectively.

Spatial discretisation
Following standard strategies of finite-element methods (Zienkiewicz et al., 2005), we want to approximate the solution of the equation (41) by means of polynomial basis functions defined in volume elements. The spatial discretisation is carried out with non-overlapping and conforming tetrahedra. We adopt the nodal form of the DG-FEM formulation (Hesthaven & Warburton, 2008), assuming that the stress and velocity vectors are approximated in the tetrahedral elements as follows where i is the index of the element, x is the spatial coordinates inside the element, and t is the time. d i is the number of nodes or degrees of freedom (DOF) associated with the interpolating Lagrangian polynomial basis function ϕ ij relative to the j-th node located at position x j . Vectors v ij and σ ij are the velocity and stress vectors, respectively, evaluated at the j-th node of the element. Although it is not an intrinsic limitation, we have adopted here the same set of basis functions for the interpolation of the velocity and the stress components. In the following, the notation P k refers to a spatial discretisation based on polynomial basis functions of degree k, and a P k element is a tetrahedron in which a P k scheme is applied. The number of DOF in a tetrahedral element is given by d i = (k + 1)(k + 2)(k + 3)/6. For instance, in a P 0 element (Figure 9.a), there is only one DOF (the stress and velocity are constant per element), while in a P 1 element (Figure 9.b), there are four DOF located at the four vertices of the tetrahedron (the stress and velocity are linearly interpolated). It is worth noting that the P 0 scheme corresponds to the case of the finite-volume method (Ben Jemaa et al., 2009;Brossier et al., 2008). For the quadratic approximation order P 2 , one node is added at the middle of each edge of the tetrahedron, leading to a total of 10 DOF per element (Figure 9.c). The first step in the finite-element formulation is to obtain the weak form of the elastodynamic system. To do so, we multiply the equation (41) by a test function ϕ ir and integrate the system over the volume of the element i. For the test function, we adopt the same kind of function as used for the approximation of the solution. This case corresponds to the standard Galerkin method and can be written as where the volume of the tetrahedral element i is denoted by V i . For the purpose of clarity, we have omitted the external forces and stresses in the equation (43). Standard manipulations of finite-elements methods (integration by parts, Green theorem for fluxes along boundary surfaces) are performed as well as an evaluation of centered flux scheme for its non-dissipative property (Ben Jemaa et al., 2007;Delcourte et al., 2009;Remaki, 2000). Moreover, we assume constant physical properties per element. We define the tensorial product ⊗ as the Kronecker product of two matrices A and B given by where (n × m) denotes the dimensions of the matrix A. We obtain the expression (45) where I 3 represents the identity matrix. In the system (45), the vectors v i and σ i should be red as the collection of all nodal values of the velocity and stress components in the element i. The system (45) indicates that the computations of the stress and velocity wavefields in one element require information from the directly neighbouring elements. This illustrates clearly the local nature of DG-FEM. The flux-related matrices P and Q are defined as follows where the component along the θ axis of the unit vector n ik of the face S ik that points from element i to element k is denoted by n ikθ , while we also introduce the mass matrix, the stiffness matrix and the flux matrices with θ ∈ {x, y, z} respectively, It is worth noting that, in the last equation of the system (46), the DOF of elements i and k appear (d i and d k , respectively) indicating that the approximation orders are totally decoupled from one element to another. Therefore, the DG-FEM allows for varying approximation orders in the numerical scheme. This feature is referred to as p-adaptivity. Moreover, given an approximation order, these matrices are unique for all elements (with a normalisation according to the volume or surface of the elements) and they can be computed before hand with appropriate integration quadrature rules. The memory requirement is therefore low, since only a collection of small matrices is needed according to the possible combinations of where the superscript n indicates the time step. We chose to apply the definition of the time step as given by , which links the mesh width and time step as follows where r i is the radius of the sphere inscribed in the element indexed by i, V Pi is the P-wave velocity in the element, and k i is the polynomial degree used in the element. Equation (49) is a heuristic stability criterion that usually works well. However, there is no mathematical proof for unstructured meshes that guarantees numerical stability.
28 WillbesetbyINTECH approximation orders. The maximum size of these matrices is (d max × d max ) where d max is the maximum number of DOF per element and the number of matrices to store is given by the square of the number of approximation orders mixed in the numerical domain. The four matrices K i , E i , F ik and G ik are computed by numerical integration using Hammer quadrature (Hammer & Stroud, 1958) and explicit forms of these matrices could be found in Etienne et al. (2010) for P 0 , P 1 and P 2 orders.
It should be mentioned that, in order to retrieve both the velocity and the stress components, the system (45) requires the computation of K −1 i , which can also be performed before hand. Note that, if we want to consider variations in the physical properties inside the elements, the pseudo-conservative form makes the computation of flux much easier and computationally more efficient than in the classical elastodynamic system. These properties come from the fact that, in the pseudo-conservative form, the physical properties are located in the left-hand side of the system (41). Therefore, no modification of the stiffness and flux matrices nor additional terms are needed in the system (45) to take into account the variation of properties. Only the mass matrix needs to be evaluated for each element and for each physical property according to the expression where χ i ( x) represents the physical property (ρ i or one of the Λ i components) varying inside the element.

Time discretisation
The time integration of the system (45) relies on the second-order explicit leap-frog scheme that allows to compute alternatively the velocity and the stress components between a half time step. The system (45) can be written as

Computational aspects
The DG-FEM is a local method, and therefore it is naturally suitable for parallel computing. In our implementation, the parallelism relies on a domain-partitioning strategy, assigning one subdomain to one CPU. This corresponds to the single program mutiple data (SPMD) architecture, which means that there is only one program and each CPU uses the same executable to work on different parts of the 3D mesh. Communication between the subdomains is performed with the message passing interface (MPI) parallel environment (Aoyama & Nakano, 1999), which allows for applications to run on distributed memory machines. For efficient load balancing among the CPUs, the mesh is divided with the partitioner METIS (Karypis & Kumar, 1998), to balance the number of elements in the subdomains, and to minimise the number of adjacent elements between the subdomains. These two criteria are crucial for the efficiency of the parallelism on large-scale numerical simulations. Figure 10 shows the observed speed-up (i.e. the ratio between the computation time with one CPU, and the computation time with N CPUs) when the number of MPI processes is increased from 1 to 256, for strong scaling calculations on a fixed mesh of 1.8 million P 2 elements. This figure shows good efficiency of the parallelism, of around 80%. In our formulation, another key point is the time step, which is common for all of the subdomains. The time step should satisfy the stability condition given in equation (49) for every element. Consequently, the element with the smallest time step imposes its time step on all of the subdomains. We should mention here a more elaborate approach with local time stepping (Dumbser et al., 2007) that allows for elements to have their own time step independent of the others. Nevertheless, the p-adaptivity offered by DG-FEM allows mitigation of the computational burden resulting from the common time step as we shall see.

Source excitation
We proceed with the addition of the excitation to incremental increase of each involved field component. The excitation of a point source is projected onto the nodes of the element that contains the source as follows with s n i the nodal values vector associated to the excited component, t = nΔt, x s the position of the point source and s(t) the source function. Equation (50) gives the source term that should be added to the right-hand side of equation (48) for the required components. It should be noticed that this term is only applied to the element containing the source. Depending on the approximation order, the spatial support of the source varies. Figure 11.a shows that the support of a P 0 element is actually the whole volume of the element (represented on the cross-section with a homogeneous white area). In this case, no precise localisation of the source inside the element is possible due to the constant piece-wise interpolation approximation. On the other hand, in a P 1 element (Figure 11.b), the spatial support of the source is linear and allows for a rough localisation of the source. In a P 2 element (Figure 11.c), the quadratic spatial support tends to resemble the expected Dirac in space close to the source position. It should be noted that the limitations concerning source localisation also apply to the solution extraction at the receivers, according to the approximation order of the elements containing the receivers.

Free surface condition
Among the various approaches presented previously, we proceed by considering that the free surface follows the mesh elements. For the element faces located on the free surface, we use an explicit condition by changing the flux expression locally. This is carried out with the concept of virtual elements, which are exactly symmetric to the elements located on the free surface. Inside the virtual elements, we impose a velocity wavefield that is identical to the wavefield of the corresponding inner elements, and we impose an opposite stress wavefield on this virtual element. Thanks to the nodal formulation, the velocity is seen as continuous across the free surface, while the stress is equal to zero on the faces related to the free surface. This is a quite natural approach similar to the one used in continuous finite-element methods where the test function is set to zero on the free surface boundary.

Absorbing boundary condition
We proceed through some simulations of wave propagation in a homogeneous, isotropic and purely elastic medium for an illustration of CPML conditions. The model size is 8 km × 8 km × 8 km, and the medium properties are: V P = 4000 m/s, V S = 2310 m/s and ρ = 2000 kg/m 3 . An explosive source is placed at coordinates (x s = 2000 m, y s = 2000 m, z s = 4000 m) and a line of receivers is located at coordinates (3000 m ≤ x r ≤ 6000 m, y r = 2000 m, z r = 4000 m) with 500 m between receivers. The conditions of the tests are particularly severe, since the source and the receivers are located close to the CPMLs (at a distance of 250 m), thus favouring grazing waves. The source signature is a Ricker wavelet with a dominant frequency of 3 Hz and a maximum frequency of about 7.5 Hz. Due to the explosive source, only P-wave is generated and the minimum wavelength is about 533 m. The mesh contains 945,477 tedrahedra with an average edge of 175 m, making a discretisation of about 3 elements per λ min . Figures 12.c and 12.d show the results obtained with the P 2 interpolation and CPMLs of 10-elements width (L cpml = 1750 m) at all edges of the model. With the standard scale, no reflection can be seen from the CPMLs. When the amplitude is magnified by a factor of 100, some spurious reflections are visible. This observation is in agreement with the theoretical reflection coefficient (R coe f f = 0.1%) in equation (20). As shown by Collino & Tsogka (2001), the thickness of the absorbing layer plays an important role in the absorption efficiency. In Figures 12.a and 12.b, the same test was performed with CPMLs of 5-elements width (L cpml = 875 m) at all edges of the model. Compared to Figures 12.c and 12.d, the amplitude of the reflections have the same order of magnitude. Nevertheless, in the upper and left parts of the model, some areas with a strong amplitude appear close to the edges. These numerical instabilities arise at the outer edges of the CPMLs, and they expand over the complete model during the simulations.
Instabilities of PML in long time simulations have been studied in electromagnetics (Abarbanel et al., 2002;Bécache et al., 2004). For elastodynamics, remedies have been proposed by Meza-Fajardo & Papageorgiou (2008) for an isotropic medium with standard PML. These authors proposed the application of an additional damping in the PML, onto the directions parallel to the layer, leading to a multiaxial PML (M-PML) which does not follow strictly the matching property of PML in the continuum and which has a less efficient absorption power. Through various numerical tests, Etienne et al. (2010) has shown that instabilities could be delayed outside the time window of simulation when considering extended M-PML from CPML. Table 2 gives the computation times for updating the velocity and stress wavefields in one element for one time step, for different approximation orders, without or with the update of the CPML memory variables (i.e. elements located outside or inside the CPMLs). These computation times illustrate the significant increase with respect to the approximation order, and they allow an evaluation of the additional costs of the CPML memory variables computation from 40% to 60%. The effects of this additionnal cost have to be analysed in Approximation order Element outside CPML Element inside CPML P 0 2.6 μs 3.6 μs P 1 5.0 μs 8.3 μs P 2 21.1 μs 29.9 μs Table 2. Computation times for updating the velocity and stress wavefields in one element for one time step. These values correspond to average computation times for a computing platform with bi-processor quad core Opteron 2.3 GHz CPUs interconnected with Infiniband 20 at Gb/s. the context of a domain-partitioning strategy. The mesh is divided into subdomains, using a partitioner. Figure 13.a shows the layout of the subdomains that were obtained with the partitioner METIS (Karypis & Kumar, 1998) along the xy plane used in the previous validation tests. The mesh was divided into 32 partitions, although only a few of these are visible on the cross-section in Figure 13.a. We used an unweighted partitioning, meaning that each partition contains approximately the same number of elements.

Saving computation time and memory
The subdomains, partially located in the CPMLs, contain different numbers of CPML elements. In large simulations, some subdomains are totally located inside the CPMLs, and some others outside the CPMLs. In such a case, the extra computation costs of the subdomains located in the absorbing layers penalise the whole simulation. Indeed, most of the subdomains spend 40% to 60% of the time just waiting for the subdomains located in the CPMLs to  (Karypis & Kumar, 1998) along the xy plane that contains the source location. Grey lines, the limits of the CPMLs. The mesh was divided into 32 partitions, although only a few of these are visible on this cross-section. (b) View of the approximation order per element along the same plane. Black, the P 2 elements; white, the P 1 elements.
complete the computations at each time step. For a better load balancing, we propose to benefit from the p-adaptivity of DG-FEM, using lower approximation orders in the CPMLs. Indeed, inside the absorbing layers, we do not need a specific accuracy, and consequently the approximation order can be decreased. Table 2 indicates that such a mixed numerical scheme is advantageous, since the computation time required for a P 0 or P 1 element located in the CPML is shorter than the computation time of a standard P 2 element. Figure 13.b shows the approximation order per element when P 1 is used in the CPMLs and P 2 in the rest of the medium. We should note here that the interface between these two areas is not strictly aligned to a cartesian axis, and has some irregularities due to the shape of the tetrahedra. Although it is possible to constrain the alignment of the element faces parallel to the CPML limits, we did not observe significant differences in the absorption efficiency whether the faces are aligned or not. Figure 14.a shows the seismograms computed when the modelling was carried out with P 2 inside the medium and P 1 in the CPMLs. Absorbing layers of 10-elements width are applied at all edges of the model. For comparison, Figure 14.b shows the results obtained with P 2 inside the medium and P 0 in the CPMLs. In this case, the spurious reflections have significant amplitudes, preventing any use of these seismograms. On the other hand, the seismograms computed with the mixed scheme P 2 /P 1 show weak artefacts, and are reasonably comparable with the seismograms obtained with complete P 2 modelling. Therefore, taking into account that the computation time and the memory consumption of the P 2 /P 1 simulation are nearly half of those required with the full P 2 modelling, we can conclude that this mixed numerical scheme is of interest. It should be noticed that it is possible to adopt a weighted partitioning approach to overcome partly load balancing issues We should also stress that the saving in CPU time and memory provided with this kind of low-cost absorbing boundary condition is crucial for large 3D simulations, and this becomes a must in the context of 3D seismic imaging applications that require a lot of forward problems, such as FWI.

Accuracy of DG-FEM with tetrahedral meshes
There are a variety of studies in the literature concerning the dispersive and dissipative properties of DG-FEM with reference to wave-propagation problems. Let us quote few examples: Ainsworth et al. (2006) provided a theoretical study for the 1D case; Basabe et al. (2008) analysed the effects of basis functions on 2D periodic and regular quadrilateral meshes; and  discussed the convergence of the DG-FEM combined with ADER time integration and 3D tetrahedral meshes. More related to our particular concern here, Delcourte et al. (2009) provided a convergence analysis of the DG-FEM with a centred flux scheme and tetrahedral meshes for elastodynamics. They demonstrated the sensitivity of the DG-FEM to the mesh quality, and they proved that the convergence is limited by the second-order time integration we have used in the present study, despite the order of the basis function. Specific analysis of the convergence in the scheme we have presented could be found in Etienne et al. (2010).

2D finite-element discontinuous Galerkin method in the frequency domain
On land exploration seismology, there is a need to perform elastic wave modeling in area of complex topography such as foothills and thrust belts (Figure 15) in the frequency domain. Moreover, onshore targets often exhibit weathered layers with very low wave speeds in the near surface which require a locally-refined discretisation for accurate modeling. In shallow water environment, a mesh refinement is also often required near the sea floor for accurate modeling of guided and interface waves near the sea floor. Accurate modeling of acoustic and elastic waves in presence of complex boundaries of arbitrary shape and the local adaptation of the discretisation to local features such as weathered near surface layers or sea floor were two of our motivations behind the development of a discontinuous element method on unstructured meshes for acoustic and elastic wave modeling.

hp-adaptive discontinuous Galerkin discretisation
Similarly to the time formulation we adopt the nodal form of the DG formulation, assuming that the wavefield vector is approximated in triangular elements for 2D geometry which leads to the following expression, where u is the wavefield vector of components such as the following vector u = (p, v x , v y , v z ) for acoustic propagation. The index of the element in an unstructured mesh is denoted by i. The expression u i (ω, x, y, z) denotes the wavefield vector in the element i and (x, y, z) are the coordinates inside the element i. In the framework of the nodal form of the DG method, ϕ ij denotes Lagrange polynomial and d i is the number of nodes in the element i. The position of the node j in the element i is denoted by the local coordinates (x j , y j , z j ).
In the frequency domain, the pseudo-conservative form (41) could be written in a 2D geometry as where N θ u are linear fluxes and the source vector is denoted by s. Expressions of matrices M and N could be found in .
The weak form of the system (52) is similar in the frequency domain and proceed by selecting a test function ϕ ir and then an integration over the element volume V i which gives where the quantity r ∈ [1, d i ]. In the framework of Galerkin methods, we used the same function for the test function and the shape function. Similar procedures as for the 3D case and related to standard steps of the finite-element method lead to the discrete expression, where the mass matrix K i , the stiffness matrix E i and the flux matrices F i and G i are similar to those defined for the 3D case (equation (46)). The matrix Q is also defined as for the 3D case (equation (46)) It is worth repeting that, in the equation (46), arbitrary polynomial order of the shape functions can be used in elements i and k indicating that the approximation orders are totally decoupled from one element to another. Therefore, the DG allows for varying approximation orders in the numerical scheme, leading to the p-adaptivity.
The equation (54) can be recast in matrix form as B u = s. (55)

Which interpolation orders to choose?
For the shape and test functions, we used low-order Lagrangian polynomials of orders 0, 1 and 2, referred to as P k , k ∈ 0, 1, 2 in the following (Brossier, 2009;Etienne et al., 2009). Let us remind that our motivation behind seismic modeling is to perform seismic imaging of the subsurface by full waveform inversion, the spatial resolution of which is half the propagated wavelength and that the physical properties of the medium are piecewise constant per element in our implementation of the DG method. The spatial resolution of the FWI and the piecewise constant representation of the medium direct us towards low-interpolation orders to achieve the best compromise between computational efficiency, solution accuracy and suitable discretisation of the computational domain. The P 0 interpolation (or finite volume scheme) was shown to provide sufficiently-accurate solution on 2D equilateral triangular mesh when ten cells per minimum propagated wavelength are used (Brossier et al., 2008), while 10 cells and 3 cells per propagated wavelengths provide sufficiently-accurate solutions on unstructured triangular meshes with the P 1 and the P 2 interpolation orders, respectively (Brossier, 2011). Of note, the P 0 scheme is not convergent on unstructured meshes when centered fluxes are used (Brossier et al., 2008). This prevents the use of the P 0 scheme in 3D medium where uniform tetrahedral meshes do not exist . A second remark is that the finite volume scheme on square cells is equivalent to second-order accurate n d 1 1 3 6 n z 9 5-9 13-25 24-48 Table 3. Number of nodes per element (n d ) and number of non-zero coefficients per row of the impedance matrix (n z ) for the FD and DG methods. The number n z depends on the number of wavefield components involved in the r.h.s of the first-order wave equation n der .
FD stencil (Brossier et al., 2008) which is consistent with a discretisation criterion of 10 grid points per wavelength (Virieux, 1986). Use of interpolation orders greater than 2 would allow us to use coarser meshes for the same accuracy but these coarser meshes would lead to an undersampling of the subsurface model during imaging. On the other hand, use of high interpolation orders on mesh built using a criterion of 4 cells par wavelength would provide an unnecessary accuracy level for seismic imaging at the expense of the computational cost resulting from the dramatic increase of the number of unknowns in the equation (55).
The computational cost of the LU decomposition depends on the numerical bandwidth of the matrix, the dimension of the matrix (i.e., the number of rows/columns) and the number of non-zero coefficients per row (n z ). The dimension of the matrix depends in turn of the number of cell (n cell ), of the number of nodes per cell (n d ) and the number of wavefield components (n wave ) (ranging from 3 to 5 in 2D geometry). The number of nodes in a 2D triangular element is given by Hesthaven & Warburton (2008) and leads to the following expression n d = (k + 1)(k + 2)/2 where k denotes the interpolation order similar to what is done in the 3D geometry.
The numerical bandwidth is not significantly impacted by the interpolation order. The dimension of the matrix and the number of non-zero elements per row of the impedance matrix are respectively given by n wave × n d × n cell and (1 + n neigh ) × n d × n der + 1, where n neigh is the number of neighbor cell (3 in 2D geometry) and n der is the number of wavefield components involved in the r.h.s of the velocity-pressure wave equation, equation (52). Table 3 outlines the number of non-zero coefficients per row for the mixed-grid FD and DG methods. Increasing the interpolation order will lead to an increase of the number of non-zero coefficients per row, a decrease of the number of cells in the mesh and an increase of the number of nodes in each element. The combined impact of the 3 parameters n z , n cell , n d on the computational cost of the DG method makes difficult the definition of the optimal discretisation of the frequency-domain DG method. The medium properties should rather drive us towards the choice of a suitable discretisation.
One must underline that the LU factorization is quite demanding in computer memory and has also some drawbacks for scalability, suggesting that nodes with high memory should be preferred at the expense of the CPU numbers.

Boundary conditions and source implementation
Absorbing boundary conditions are implemented with unsplitted PML in the frequency-domain DG method (Brossier, 2011) following the same approach than for the FD method: one can see that the PML implementation in the frequency is straightforward.
We have found that constraining the meshing to have edges of elements in the PML zone parallel to the direction of dissipation of the waves improves the efficiency.
Free surface boundary condition is implemented with the method of image. A virtual cell is considered above the free surface with the same velocity and the opposite pressure components to those below the free surface. This allows us to fulfill the zero pressure condition at the free surface while keeping the correct numerical estimation of the particle velocity at the free surface. Using these particle velocities and pressures in the virtual cell, the pressure flux across the free surface interface vanishes, while the velocity flux is twice the value that would have been obtained by neglecting the flux contribution above the free surface. As in the FD method, this boundary condition has been implemented by modifying the impedance matrix accordingly without introducing explicitely the virtual element in the mesh. The rigid boundary condition is implemented following the same principle except that the same pressure and the opposite velocity are considered in the virtual cell.
Concerning the source excitation, the point source at arbitrary positions in the mesh is implemented by means of the Lagrange interpolation polynomials for k ≥ 1. This means that the source excitation is performed at the nodes of the cell containing the source with appropriate weights corresponding to the projection of the physical position of the source on the polynomial basis. When the source is located in the close vicinity of a node of a triangular cell, all the weights are almost zero except that located near the source. In the case of the P 2 interpolation, a source close to the vertex of the triangular cell is problematic because the integral of the P 2 basis function over the volume of the cell is zero for nodes located at the vertex of the triangle. In this case, no source excitation will be performed (see equation (54)).
To overcome this problem specific to the P 2 interpolation, one can use locally a P 1 interpolation in the element containing the source at the expense of the accuracy or distribute the source excitation over several elements or express the solution in the form of local polynomials (i.e., the so-called modal form) rather than through nodes and interpolating Lagrange polynomials (i.e., the so-called nodal form). Another issue is the implementation of the source in P 0 equilateral mesh. If the source is excited only within the element containing the source, a checker-board pattern is superimposed on the wavefield solution. This pattern results from the fact that one cell out of two is excited in the DG formulation because the DG stencil does not embed a staggered-grid structure (the unexcited grid is not stored in staggered-grid FD methods; see Hustedt et al. (2004) for an illustration). To overcome this problem, the source can be distributed over several elements of the mesh or P 1 interpolation can be used in the area containing the sources and the receivers, while keeping P 0 interpolation in the other parts of the model . Of note, use of unstructured meshes together with the source excitation at the different nodes of the element contribute to mitigate the checker-board pattern in the in P 1 and P 2 schemes. The same procedure as for the source is used to extract the wavefield solution at arbitrary receiver positions.

Realistic examples for highly contrasted and strongly heterogeneous media using finite-elements methods
We shall consider two examples for the illustration of the Discontinuous Galerkin approach. The first one is related to the problem of 3D wave propagation inside an active volcano using the time-domain approach while the second one deals with the problem of 2D wave propagation above a oil reservoir using the frequency-domain approach.

Characteristics of the model
La Soufrière of Guadeloupe (France) is one of nine active volcanoes of Lesser Antilles. It belongs to a recent volcanic system situated in the south part of the Basse-Terre. A P-wave velocity model of the volcano has been obtained by first arrival time tomography (Coutant et al., 2010). Figure

Construction of the tetraedral mesh
The mesh has been built with the mesher TETGEN (Si, 2006) combined with an iterative h-refinement procedure to obtain a locally adapted mesh to the velocity field(with an average of 3 elements per minimum wavelength λ min ): a cross-section is shown in the left panel of the Figure 17. For building this mesh, we have started our iterative reconstruction with a uniform mesh shown in the right panel of the Figure 17. After the sixth refinement iteration, the discretization criteria are met. Areas of high velocities are correlated with the parts of the  mesh where the elements are the largest ones. On the contrary, near the free surface, we find the finest elements.

Numerical result
We have performed 3D simulations with the Discontinuous Galerkin Finite-Element Method in the time domain. The computations have been performed on a Blue Gene machine with 512 processors. The statistics for these computations are given in Table 4.
The configuration of the seismic acquisition is given in Figure 16. This is a quasi-2D system with a profile according to the East-West direction, which includes 100 single-component receivers (v z ) with 10 m between receivers. The source is a shot of dynamite. For the numerical simulations, we used an explosive source with a Ricker function of dominant frequency of 10 Hz (maximum frequency 25 Hz). We present in Figure 18 a comparison between the observed and computed data. Despite significant uncertainties and approximations (source function, Poisson ratio, density, absence of attenuation, low signal to noise ratio), there are striking similarities in the data. In particular, the seismic traces exhibit well marked discontinuities related to the strong velocity contrasts and the complex topography of the volcano La Soufrière.

Application 2D in the frequency-domain: the synthetic Valhall application
This 2D application is based on a synthetic representation of the Valhall zone in the North Sea, Norway. This model is representative of oil and gas fields in shallow water environments of the North Sea (Munns, 1985). The model is described as an heterogeneous P-and S-wave velocity model (Figure 19a-b). The water layer is only 70 m depth. The main targets are a gas cloud in the large sediment layer, and in a deeper part of the model, the trapped oil underneath the cap rock, which is formed of chalk. Gas clouds are easily identified by the low P-wave velocities, whereas their signature is much weaker in the V S model, as gas does not affect S-waves propagation. In order to investigate seismic imaging in such environment, the selected acquisition mimics a four-component ocean-bottom cable survey (Kommedal et al., 2004), as is deployed on the field. A line of 315 explosive sources is positioned 5 m below water surface to simulated air-gun sources and a cable of 315 3-components sensors is located on the sea floor (1 hydrophone and 2 geophones). This geological setting, which is composed of a significant soft sea-bed with high Poisson'ratio due to soft and unconsolidated sediments leads to a particularly ill-posed problem for S-wave velocity reconstruction, due to the relatively small shear-wave velocity contrast at the sea bed, which prevents recording of significant P-to-S converted waves.
For the meshing of the model, the narrow velocity range in most parts of the model requires the use of a regular mesh as much as possible for computational efficiency. However, to correctly discretize the shallow-water layer and liquid-solid interface, a p-adapted mesh implemented with a mixed P 0 -P 1 interpolation is chosen ( Figure 19c): a refined unstructured P 1 layer of cells is used for the first 130 m of the subsurface for accurate modeling of the interface waves at the liquid/solid interface, and for accurate positioning of the sources, located 5 m below the surface, and of the receivers, located on the sea floor. Below this 130 m depth zone, a regular equilateral mesh is used in combination of a P 0 interpolation. Figure 20 illustrates an example of frequency-domain data (real-part of the complex-value wavefield) at 4Hz. These data are plotted in the source/receiver domain for the full acquisition survey. The diagonal part of the figure represents the collocation of the source and the closest receiver, as the source moves in the acquisition. The Figure 21 illustrates a time-domain shot-gather for the 3 components of the sensors and a source located at position 4 km : the frequency-domain solutions at all the receiver positions have been computed for the single source at all the frequencies of the source spectrum, between 0 and 13 Hz. These frequency-domain complex-values data have then been Fourier transform to the time-domain. The time-domain show specific properties of propagation in such environments : the hydrophone and the vertical geophones are mainly sensitive to P-wave arrivals that dominate the elastic propagation in soft-sediment zones. The horizontal geophone allows however to record some late P-to-S conversions, which could be used to image the V S model from seismic imaging methods.

Conclusion
We have presented mainly two families of techniques for solving partial differential equations for elastodynamics: some finite-differences formulations in both time domain and frequency domain and some finite-element methods also in both time domain and frequency domain. Both approaches have appealing features, especially when considering seismic imaging where numerous forward problems should be performed. Such classification helps to understand the advantages and limitations of each particular method to model a specific physical phenomenon The discretization of the strong formulation of the partial differential equations has been presented through finite-difference techniques. These approaches are easy to implement and quite flexible. They are currently the methods of choice for large-scale modelling and inversion in exploration geophysics, especially in the marine environment. They may however demand a very fine discretization when the earth model contains large contrasts; and accurately modelling the responses around a sharp interface is quite challenging. We have introduced various perspectives as summation-by-parts formulation or the immersed boundary approach as well as simple mesh deformation might broaden the use of finite-differences techniques by avoiding the stair-case approximation.
The weak formulation expressed in the finite-element methods has been considered under the specific family of Discontinuous Galerkin approaches. The use of test functions gives us more freedom and the integral form provides us flexibility in the meshing. However, they lead to numerical challenges: they are more difficult to implement than the finite-difference method, they are often more expensive in computational time and memory, and they are more complicated to use because the accuracy of the response depends on the quality of the meshing. Therefore, they are not intensively used for seismic imaging and are until now more oriented to seismic modeling in the final reconstructed model.
It should be noticed that attempts exist to combine the advantages of these methods in one approach for computing elastic fields, at least for specific applications. Even, one can think that decoupling the inverse problem procedure and the forward problem is possible: we can flip-flop between the two forward problem formulations inside iterations of the inverse problem.
When the modelling method serves as the kernel of an inversion algorithm, additional constrains generally appear because the gradient of the misfit functional needs to be evaluated. The choice of the modelling approach notably depends (1) on the needed accuracy, (2) the efficiency in evaluating the solution and the gradient of the misfit functional in an inversion algorithm, and (3) the simplicity of use.
Finally, the practical implementation shall probably be adapted to the data acquisition. Densely sampled acquisition in exploration geophysics with or without blending, or in lithospheric investigation with the recent deployment of sensors such as the USarray experiment challenges our modelling choice. This seems to indicate that development in modelling and associated inversion approaches remain crucial to improve our knowledge of the subsurface, notably by extracting more information from the, ever larger, recorded data sets.