[
Free and smooth boundaries for elastic waves]
Free and smooth boundaries in 2D finitedifference schemes for transient elastic waves
B. Lombard, J. Piraux, C. Gélis, J. Virieux]
B. Lombard, J. Piraux, C. Gélis, J. Virieux
CNRS, Laboratoire de Mécanique et d’Acoustique, Marseille, France
CNRS, Géosciences Azur, SophiaAntipolis, France.
\pagerange[–5
\volume?
200?
A method is proposed for accurately describing arbitraryshaped free boundaries in finitedifference schemes for elastodynamics, in a timedomain velocitystress framework. The basic idea is as follows: fictitious values of the solution are built in vacuum, and injected into the numerical integration scheme near boundaries. The most original feature of this method is the way in which these fictitious values are calculated. They are based on boundary conditions and compatibility conditions satisfied by the successive spatial derivatives of the solution, up to a given order that depends on the spatial accuracy of the integration scheme adopted. Since the work is mostly done during the preprocessing step, the extra computational cost is negligible. Stressfree conditions can be designed at any arbitrary order without any numerical instability, as numerically checked. Using 10 grid nodes per minimal Swavelength with a propagation distance of 50 wavelengths yields highly accurate results. With 5 grid nodes per minimal Swavelength, the solution is less accurate but still acceptable. A subcell resolution of the boundary inside the Cartesian meshing is obtained, and the spurious diffractions induced by staircase descriptions of boundaries are avoided. Contrary to what occurs with the vacuum method, the quality of the numerical solution obtained with this method is almost independent of the angle between the free boundary and the Cartesian meshing.
ree surface, Seismic modeling, Velocitystress formulation, Numerical methods, Finitedifference methods, ADER schemes, Boundary conditions, Compatibility conditions.
1 Introduction
Various approaches have been proposed for simulating the propagation of elastic waves with free boundaries. The first approach is based on variational methods, as done in finite elements (DAY77), spectral finite elements (KOMATITSCH98) and discontinuous Galerkin (BENJEMAA07). These methods provide a fine geometrical description of boundaries by adapting the mesh to the boundaries. Boundary conditions are accounted for weakly by the underlying variational formulation. However, a gridgenerating tool is required, and small time steps may result from the smallest geometrical elements and from the stability condition. The SAT methods based on energy estimates (CARPENTER94) avoid these limitations by introducing Cartesian grids and give timestable highorder schemes with interfaces. However and up to our knowledge, these methods have not been applied so far to elastodynamics with free boundaries.
The second approach used in this context is based on the strong form of elastodynamics, as done in finite differences and spectral methods (TESSMER94). In seismology, finite differences are usually implemented on staggered Cartesian grids, either with completely staggered stencils (CSS) or with the recently developed partially staggered stencils (PSS). With CSS, the velocity and stress components are distributed between different node positions (VIRIEUX86). With PSS, all the velocity components are computed at a single node, as are the stress components, although the latter are shifted by half a node in two separate grids. Secondorder (SAENGER00; SAENGER04) and fourthorder (BOHLEN03; CRUZ04) spatiallyaccurate PSS have been developed; for further discussion, we denote them PSS2 and PSS4, respectively. Unlike variational methods, finite differences require special care to incorporate the free boundary conditions strongly. There exist two main strategies for this purpose:

First, the boundaries can be taken into account implicitly by adjusting the physical parameters locally (KELLY76; VIRIEUX86; MUIR92). The bestknown implicit approach is the socalled vacuum method (ZAHRADNIK95; GRAVES96; MOCZO02; GELIS05). For instance, the vacuum method applied to PSS involves setting the elastic parameters in the vacuum to zero, and using a small density value in the first velocity node in the vacuum to avoid a division by zero. However, this easytoimplement method gives at best secondorder spatial accuracy. In addition, a systematic numerical study has shown that the accuracy of the solution decreases dramatically when the angle between the boundary and the meshing increases (BOHLEN06). Lastly, applying the vacuum method sometimes gives rise to instabilities: see for instance PSS4 (BOHLEN03).

A second idea is to explicitly change the scheme near the boundaries (KELLY76). The bestknown explicit approach is the socalled image method, which was developed for dealing with flat boundaries to fourthorder accuracy (LEVANDER88) and then extended to variable topographies (JIH88; ROBERTSSON96; ZHANGW06). However image methods require a fine grid to reduce the spurious diffractions up to an acceptable level. To avoid this spatial oversampling, various techniques have been proposed, such as grid refinement in the vicinity of the boundary (RODRIGUES93) or adjusted finitedifference approximations: see (MOCZO07) for a review on these subjects.
The aim of this paper is to present a finitedifference approach accounting for free boundaries without introducing the aforementioned drawbacks of the vacuum and image methods. The requirements are as follows: smooth arbitraryshaped boundaries must be treated as simply as straight boundaries; the accuracy of the method must not depend on the position of the boundary relative to the meshing; and lastly, the computations must be stable even with very long integration times. We establish that these requirements can be met by applying an explicit approach involving fictitious values of the solution in the vacuum. In previous studies, interface problems in the context of elastodynamics were investigated in a similar way (PIRAUX01; LOMBARD04; LOMBARD06). The fictitious values are highorder Taylor expansions of the boundary values of the solution. Estimating these boundary values involves some mathematical background, in order to compute the highorder boundary conditions; to determine a minimal set of independent boundary values; lastly, to perform a leastsquare numerical estimation of this minimal set. To help the reader, subroutines in FORTRAN are proposed freely at the web page http://w3lma.cnrsmrs.fr/~MI/Software/. These subroutines enable a straightforward implementation of the algorithms detailed in the present paper.
The disadvantage here is that the above requirements cannot be fully satisfied if staggeredgrid schemes are used. Singlegrid finitedifference schemes are therefore chosen, where all the unknowns are computed at the same grid nodes. Our numerical experiments are based on the highorder ADER schemes which are widely used in aeroacoustics (MUNZ05). Although these schemes are not yet widely used in the field of seismology (DUMBSER06), they have also great qualities because of their accuracy and their stability properties: using 10 grid nodes per minimal Swavelength with a propagation distance of 50 wavelengths gives highly accurate results. Moreover, on Cartesian grids, these methods do not require much more computational memory than staggeredgrid schemes.
This paper is organized as follows. Section 2 deals with the continuous problem: the highorder boundary conditions and compatibility conditions are stated. These conditions are useful for handling the discrete problem presented in section 3, where the focus is on obtaining fictitious values of the solution in the vacuum. In section 4, numerical experiments confirm the efficiency of this method in the case of various topographies. In section 5, conclusions are drawn and some prospects suggested.
2 The continuous problem
2.1 Framework
Let us consider a solid separated from the vacuum by a boundary (Figure 1). The configuration is inplane and twodimensional, with a horizontal axis and a vertical axis pointing respectively rightward and downward. The boundary is described by a parametric expression where the parameter describes the sampling of the boundary. The tangential vector and the normal vector are and , respectively, with , , and refers to the transposed vector. We assume the spatial derivatives at any point of the boundary to be available, as specified below.
The solid is assumed to be linearly elastic, isotropic, and to have the following constant physical parameters: the density and the Lamé coefficients , . The P and Swave velocities are and . A velocitystress formulation is adopted, hence the unknowns are the horizontal and vertical components of the elastic velocity , and the independent components of the elastic stress tensor . Setting
the solution satisfies the firstorder hyperbolic system (VIRIEUX86)
(1) 
2.2 Highorder boundary conditions
At any point on the free surface (Figure 1), the stress tensor satisfies the homogeneous Dirichlet conditions . These zeroth order boundary conditions are written compactly
(2) 
where is the limit value of at and is the matrix
From now on, the dependence on is generally omitted. To determine the boundary conditions satisfied by the firstorder spatial derivatives of , two tasks are performed. First, the zerothorder boundary conditions (2) are differentiated in terms of . The time derivative is replaced by spatial derivatives using the conservation laws (1), which gives
(3) 
Secondly, the zerothorder boundary conditions (2) are differentiated in terms of the parameter describing . The chainrule gives
(4) 
Since the matrix in (4) involves and , it accounts for the curvature of at . Setting the block matrix
equations (2), (3) and (4) give the boundary conditions up to the firstorder
with
Let be an integer whose value will be discussed in section 3. To get the boundary conditions up to the th order, one deduces from (2)
(5) 
The derivatives are replaced by spatial derivatives by applying times the chain rule. The derivatives are replaced by spatial derivatives by injecting times the conservation laws (1). The boundary conditions soobtained up to the th order can be written compactly
(6) 
with
(7) 
where and . The vector has components. is a matrix, with . This matrix involves the successive derivatives of the curvature of at . Computing with is a tedious task, which can be greatly simplified by using computer algebra tools.
2.3 Compatibility conditions
The second spatial derivatives of stress components are linked together by the compatibility condition of Barréde Saint Venant (LOVE)
(8) 
with
This compatibility condition is a necessary and sufficient condition for the strain tensor to be symmetrical. If , it can be differentiated times in terms of and . With , one obtains relations; with , . Unlike the boundary conditions, these compatibility conditions are satisfied everywhere in : in particular, they are satisfied at on . The vector of boundary values can therefore be expressed in terms of a shorter vector with independent components
(9) 
An algorithm for building the matrix is given in LOMBARD06.
3 The discrete problem
3.1 Numerical scheme
To integrate the hyperbolic system (1), we introduce a single Cartesian lattice of grid points: , where is the mesh spacing and is the time step. Unlike with staggered grids, all the unknowns are computed at the same grid nodes. The approximation of is computed using any explicit, twostep, and spatiallycentred finitedifference scheme. A review of the huge body of literature on finitedifferences is given in LEV90 and MOCZO07.
Here we propose to use ADER schemes, that allow to reach easily arbitrary highorder of time and space accuracy (MUNZ05). On Cartesian grids, these finitevolume integration schemes originally developed for aeroacoustic applications are equivalent to finitedifference LaxWendrofftype integration schemes (LORCHER06). In the numerical experiments described in section 4, we use a fourthorder ADER integration scheme. This scheme is stable under the CourantFriedrichsLewy (CFL) condition in 2D; as usually with singlegrid schemes, it is slightly dissipative (MUNZ05).
Many other singlegrid schemes can be used in this context. In particular, the method described in the next subsections has been successfully combined with fluxlimiter schemes (LEV90) and with the standard secondorder LaxWendroff scheme. Difficulties have been encountered with dissipativefree schemes based on centred staggeredgrid finitedifference schemes, as we will see in section 3.6.
3.2 Use of fictitious values
Timemarching at gridpoints where the stencil crosses requires fictitious values of the solution in the vacuum, which have to be determined. The question arises as to how to compute, for instance, the fictitious value at the grid point in the vacuum, as sketched in Figure 2. Let be the orthogonal projection of on , with coordinates . At any grid point , we denote
the matrix containing the coefficients of th order Taylor expansions in space at , where is the identity matrix, , and . The fictitious value is defined as the Taylorlike extrapolation
(10) 
where defined by (7) still remains to be estimated.
3.3 Reduced vector of boundary values
Before determining in (10), we first reduce the number of independent components it contains. The expressions obtained in section 2 are used for this purpose. The linear homogeneous system following from (6) and (9) is
(11) 
This system has fewer equations () than unknowns (). It therefore has an infinite number of possible solutions that constitute a space with the dimension . Let be a matrix containing the basis vectors of the kernel of . The general solution of (11) is therefore
(12) 
where the components of are real numbers. Injecting (12) into (9) gives
(13) 
The computation of is a key point. For this purpose, we use a classical linear algebra tool: singular value decomposition of . Technical details can be found in the Appendix A of LOMBARD04.
3.4 Computation of fictitious values
Let us now consider the grid points of in the circle centred at with a radius ; for instance, in Figure 2. At these points, we write the th order Taylor expansion in space of the solution at , and then we use the expression (13). This gives
(14) 
The set of equations (14) is written compactly via a matrix
(15) 
where is the vector containing the exact values of the solution at the grid nodes of inside . These exact values are replaced by the known numerical values , and Taylor rests are removed. From now on, numerical values and exact values of the fields are used indiscriminately. The discrete system thus obtained is overdetermined (see the remark (i) about and typical values of in subsection 3.5). We now compute its leastsquares solution
(16) 
where the matrix denotes the pseudoinverse of . From (10), (13) and (16), the fictitious value in the vacuum at is
(17) 
The matrix is called the extrapolator at . The fictitious values have no clear physical meaning. They only allow, by interpolation with numerical values inside , to recover the highorder Dirichlet conditions (7).
3.5 Comments and practical details
The extrapolation method described in section 3.4 has to be applied at each grid point in the vacuum where a fictitious value is required for the timemarching procedure. Useful comments are proposed about this method:

The radius of must ensure that the number of equations in (15) is greater than the number of unknowns:
(18) No theoretical results are available about the optimal value of . However, numerical studies have shown that a definite overestimation ensures longterm stability: typically, . Various strategies can be used to ensure (18), such as an adaptative choice of depending on the local geometry of at . Here we adopt a simpler strategy consisting in using a constant radius . With , numerical experiments have shown that is a good candidate for this purpose. In this case, one typically obtains .

Since the boundary conditions do not vary with time, the extrapolators in (17) can be computed and stored during a preprocessing step. At each time step, only small matrixvector products are required.

The extrapolators account for the local geometry of at the projection points on via (section 2.2). Moreover, they incorporate the position of relative to the Cartesian meshing, via (14) and (17). The set of extrapolators therefore provides a subcell resolution of in the meshing, avoiding the spurious diffractions induced by a naive description of the boundaries.

The stability of the method has not been proved. However, numerical experiments clearly indicate that the CFL condition of stability is not modified compared with the case of a homogeneous medium. The solution does not grow with time, even in the case of longtime simulations (see section 4.5).

In a previous onedimensional study (PIRAUX01), the local truncation error of the method has been rigorously analysed, leading to the following result: using the fictitious values (17) ensures a local th order spatial accuracy if , where is the order of spatial accuracy of the scheme. In 2D configrations with material interfaces (LOMBARD04; LOMBARD06), no proof has been conducted, but numerical experiments have shown that the th order overall accuracy is also maintained by taking . Note that a slightly smaller order of extrapolation can be used: suffices to provide th order overall accuracy (GUSTAFSSON75). The value is therefore used for the fourthorder ADER scheme.

The extrapolators do not depend on the numerical scheme adopted. They depend only on and on physical and geometrical features. Standard subroutines for computing the extrapolators can therefore be developed and adapted to a wide range of schemes. Subroutines of this kind are freely available in FORTRAN at the web page http://w3lma.cnrsmrs.fr/~MI/Software/.
3.6 Case of staggeredgrid schemes
(a)  (b) 

Instead of using a singlegrid scheme as proposed in section 3.1, readers may be interested in adapting our approach to staggeredgrid schemes such as CSS or PSS (see section 1 for the definition of these terms). However, in the case of some of the boundary positions relative to the meshing, computational instabilities occur, especially when longtime integration is considered.
To understand why this is so, let us consider PSS2. Taking a simple flat boundary to exist between the medium and the vacuum leads to two typical geometrical configurations. At one position of the free surface, the boundary discretization will require only the stress field to be extrapolated (Figure 3(a)). Our procedure works satisfactorily with this type of discretization at any order . It also yields stable and accurate solutions when dealing with PSS4, contrary to the vacuum method. Using 10 grid nodes per minimal Swavelength gives similar performance in this case to those of our numerical experiments based on the ADER scheme, which are shown in section 4.
At another position of the free surface where only extrapolated velocities are required within a wide zone (Figure 3(b)), our procedure results in instabilities. The reason for this problem is as follows: fictitious velocities involve firstorder boundary conditions (3) and higherorder conditions (see section 2.2), but they do not involve the fundamental zerothorder Dirichlet conditions (2). Since the latter conditions are never enforced, an increasing oscillating drift occurs near the boundary, which invalidates the computations. Similar behavior is observed with PSS4, but after a longer time: the numerical solution generally works well during a few thousand time steps, before growing in a unstable manner.
The extrapolation method presented here is therefore not recommended for use with staggeredgrid schemes, especially PSS2, except in the trivial case sketched in Figure 3(a).
(a) 
(b) 
(c) 
(a) 
(b) 
(c) 
(d) 
3.7 Case of nonsmooth geometries
Up to know, we have assumed that the boundary was sufficiently smooth at the projection points, being at least at each , where is the order of differentiation defined in section 3.5. Let us assume now that is only at a point , with . Then, the components of in (6) involving the derivatives and () of the parametric representation are discontinuous, invalidating locally the method proposed. In our software, we have implemented the following rough treatment:

If , the boundary owns a corner and the solution has an integrable singularity. The corner is replaced by an arc of circle centred at with radius (figure 4), leading to a curve.

If , as in the previous case at and with , the values of and () are taken indiscriminately on one side or the other of the point considered.
No numerical instabilities were observed if (in case (i)) or the radius of curvature (in case (ii)) is much greater than . It is agreed that the accuracy of computations is no more controlled in the cases (i) and (ii), especially the convergence towards the exact solution.
More sophisticated treatments of geometrical singularities, such as spacetime mesh refinement (BERGER98), require further investigation, which is out of the scope of the present paper. New studies are also needed in the case of merging boundaries, occuring for instance when an internal material interface reaches the free surface (MOCZO04).
4 Numerical experiments
4.1 Configurations
The time evolution of the source is a Ricker wavelet
(19) 
where is the central frequency, and . The maximal frequency defined by (the tilde designates the Fourier transform) is . We will adopt this frequency for our rule of thumb about the number of grid nodes per Swavelength. The following values of the physical parameters will be used in all the following tests: , m/s, and m/s. Lastly, the mesh size and the time step satisfy .
The simulations are performed on a PC Pentium at 3 GHz with 2 GB of RAM. The results of tests 1 and 2, with constant and null curvature of , compare quantitatively with analytical solutions denoted by a solid line. Test 3, with a variable curvature, is purely qualitative. Test 4 shows the slow decrease in the mechanical energy which occurs during very long integration times, which confirms the stability of the method.
4.2 Test 1: circular boundary
(a) 
(b) 
Computational efficiency. Let us consider a circular cavity containing vacuum, with radius 1 km, at the center of a 18 km 18 km domain. In a first part, the mesh spacing is m. The source is a rightwardmoving plane wave, with Hz, ensuring 22 grid nodes per minimal Pwavelength and 10 grid nodes per minimal Swavelength at that frequency.
During the preprocessing step, the program finds the 616 grid nodes where fictitious values are required; it also computes and stores the 616 extrapolators defined by the expression (17). Time integration is performed in 550 time steps, which corresponds to a propagation time of 2.75 s and a propagation distance of 22 minimal wavelengths. The preprocessing step takes 21 s of CPU time. The time integration takes 1100 s of CPU time, including 28.6 s induced by the computation and by the use of fictitious values, which amounts to an extra time cost of only 2.6 %. Figure 5 shows snapshots of at the initial instant (a), after 275 time steps (b) and after 550 time steps (c).
Quantitative study. In a second part, three discretizations are considered: m, 50 m, 60 m, corresponding respectively to 10, 5 and 4 grid nodes per minimal Swavelength. A receiver is set just above the cavity, at the position (9 km, 10.2 km), that can be seen on Figure 5; it mainly records the waves propagating along the boundary. Numerically, these waves are highly sensitive to the quality of the fictitious values defined in section 3.
Figure 6(a) shows the time history of at the receiver. In this time window, three main wave packets are generated; with the scale of Figure 6(a), the third packet cannot be seen. The amplitude is divided by a factor of approximately 30 from one packet to the following one. For the sake of clarity, zooms around each wave packet are shown in Figure 6(b,c,d). These solutions are compared with an exact solution computed thanks to inverse Fourier transforms on 4096 frequencies, with Hz as the sampling frequency; each harmonic component is expanded into 60 Bessel modes. The agreement between the numerical and the analytical values is very good when 10 grid nodes per wavelength are used, even at very small amplitudes (d). For 5 and 4 grid nodes per wavelength, the solution is slightly less accurate, but it is still acceptable.
4.3 Test 2: inclined straight boundary
(a) 
(b) 
(c) 
(a) 
(b) 
The Garvin’s problem. As a second test, we take a plane boundary inclined against the Cartesian mesh. The domain under investigation is 18 km wide and 12 km high, with the origin of the coordinates on the top and left. The mesh spacing is m. Four receivers at (10 km, 1.8 km), (11 km, 1.6 km), (12 km, 1.4 km) and (13 km, 1.2 km) belong to the free boundary which is inclined at an angle of relative to .
An explosive source S is buried at (9 km, 2.1 km), with Hz. The distance between the source and the free surface is roughly 100 m , where is the wavelength of the compressional waves at frequency , and hence large Rayleigh waves are generated, with a velocity m/s. To prevent the occurence of spurious oscillations, the source is spread out numerically over a radius of m. The source is weighted by a gaussian law with a standard deviation m. The spatial discretization ensures the sampling of roughly 18 grid nodes per minimal Pwavelength, and 9 grid nodes per minimal Swavelength at the frequency .
Figure 7(a) shows a snapshot of after 1200 time steps, corresponding to a propagation time of 2.25 s and a propagation distance of 55 minimal wavelengths. Direct cylindrical waves are observed, together with converted PP waves, converted PS waves (with an almost linear wavefront), and Rayleigh waves. In Figure 7(b), the time history of recorded at the receivers can be favourably compared with an exact solution. The latter is obtained by convolving the Green’s function obtained by the wellknown Cagniardde Hoop method (GARVIN55; SANCHEZ06) with the source wavelet (19) and with the discrete source spreading.
Influence of the slope. To quantify the effects of the angle between the boundary and the Cartesian meshing on the numerical solution, we perform a parametric study of the error in terms of . Ten angles are considered, from to in steps of . In each configuration, the waves are measured at the free boundary after propagating for 65 minimal wavelengths. The error of is measured in norm , and then it is normalized by the norm of the exact time history of .
The results of this study are shown in Figure 8, with various discretizations: 5, 10, 20 grid nodes per minimal Swavelength. With a given , the error is almost constant and independent of . This constitutes a crucial advantage of our method over the vacuum method, where the error at is much greater than that at : it means that an extremely fine discretization is required to obtain accurate results with the vacuum method when arbitraryshaped boundaries are encountered (BOHLEN06).
4.4 Test 3: sinusoidal boundary
Since boundaries not related to the finitedifference grid can be included, the third test is performed on a sinusoidal free boundary, with a peaktopeak amplitude of 800 m and various wavelengths: 0.5 km, 1 km and 2 km. The sinusoidal curve is centred around km. The source S is located at (9 km, 1 km). The other parameters are the same as in test 2. Figure 9 shows snapshots of at the final instant 2.25 s. One can clearly see how the wavelength of the sinusoidal boundary influences the diffracted fields.
Convergence studies (not shown here) were performed in these three cases, by comparing solutions computed on finer grids. We again concluded that accurate solutions can be obtained when the simulations involve approximately 10 grid nodes per minimal Swavelength at the frequency of the source wavelet, even in the case of complex topographies with variable curvatures.
4.5 Test 4: longterm stability
The fourth test focuses on longterm stability (STACEY94; HESTHOLM03). For this purpose, we consider a circular elastic domain with a radius of 150 m, surrounded by vacuum. The source S is located inside the circle, at (320 m, 200 m). This configuration is obviously not realistic, but it enlights the influence of the boundary on the numerical solution after many reflections, and especially on the possible excitation of numerical spurious modes leading to longterm instability. The mesh size is m. Time integration is performed during time steps, with Hz.
Figure 10(a) shows a snapshot of at the final instant: no instability is observed, and the antisymmetry of is satisfied. Once the source is extincted (), the mechanical energy is theoretically maintained. It can be written in terms of and
(20) 
At each time step, the integral in (20) is estimated by a basic trapezoidal rule at the grid nodes inside . Figure 10(b) shows the time history of this mechanical energy soobtained. It slightly decreases, due to the numerical diffusion of the scheme, which confirms that the method is stable.
5 Conclusion
Here we have presented a method of incorporating free boundaries into timedomain singlegrid finitedifference schemes for elastic wave simulations. This method is based on fictitious values of the solution in the vacuum, which are used by the numerical integration scheme near boundaries. These highorder fictitious values accurately describe both the boundary conditions and the geometrical features of the boundaries. The method is robust, involving negligible extra computational costs.
Unlike the vacuum method, the quality of the numerical solution thus obtained is almost independent of the angle between the free boundaries and the Cartesian meshing. Since the free boundaries do not introduce any additional artefacts, one can use the same discretization as in homogeneous media. Typically, when a fourthorder ADER scheme is used on a propagation distance of 50 minimal wavelengths, 10 grid nodes per minimal Swavelength yield to a very good level of accuracy. With 5 grid nodes per minimal Swavelength, the solution is less accurate but still acceptable.
For the sake of simplicity, we have dealt here with academic cases, considering twodimensional geometries, constant physical parameters, and simple elastic media. Let us examine briefly the generalization of our approach to more realistic configurations:

Extending the method to 3D topographies a priori does not require new tools. The main challenge will concern the computational efficiency of parallelization. A key point is that the determination of each fictitious value is local, using numerical values only at neighboring grid nodes. Particular care will however be required for fictitious values near frontiers between computational subdomains, in order to minimize the exchanges of data.

Near free boundaries, the domains of propagation are usually smoothly heterogeneous. To generalize our method to continuously variable media, the main novelty expected concerns the highorder boundary conditions detailed in section 2.2. With variable matrices and indeed and , the procedure (5) will involve the following quantities, to be estimated numerically:

Realistic modeling of wave propagation requires to incorporate attenuation. The only rheological viscoelastic models able to approximate constant quality factor over a frequency range are the generalized Maxwell body (EMMERICH84) and the generalized Zener body (CARCIONE01). These two equivalent models (MOCZO05) yield to additional unknowns called memory variables. In the time domain, the whole set of unknowns satisfies a linear hyperbolic system with source term
(21) where is a definite positive matrix. Compared with the elastic case (1) examined in the present paper, the main difference expected concerns the time differentiation of the boundary condition (2). Indeed, equation (3) has to be modified accordingly to (21). Similar modifications are also foreseen in the case of poroelasticity in the lowfrequency range (DAI95), where the evolution equations can be put in the form (21).