Free and smooth boundaries for elastic waves] Free and smooth boundaries in 2-D finite-difference 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, Sophia-Antipolis, France.
\pagerange[5 \volume? 200?


A method is proposed for accurately describing arbitrary-shaped free boundaries in finite-difference schemes for elastodynamics, in a time-domain velocity-stress 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. Stress-free conditions can be designed at any arbitrary order without any numerical instability, as numerically checked. Using 10 grid nodes per minimal S-wavelength with a propagation distance of 50 wavelengths yields highly accurate results. With 5 grid nodes per minimal S-wavelength, 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, Velocity-stress formulation, Numerical methods, Finite-difference 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 grid-generating 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 time-stable high-order 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. Second-order (SAENGER00; SAENGER04) and fourth-order (BOHLEN03; CRUZ04) spatially-accurate PSS have been developed; for further discussion, we denote them PSS-2 and PSS-4, respectively. Unlike variational methods, finite differences require special care to incorporate the free boundary conditions strongly. There exist two main strategies for this purpose:

  1. First, the boundaries can be taken into account implicitly by adjusting the physical parameters locally (KELLY76; VIRIEUX86; MUIR92). The best-known implicit approach is the so-called 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 easy-to-implement method gives at best second-order 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 PSS-4 (BOHLEN03).

  2. A second idea is to explicitly change the scheme near the boundaries (KELLY76). The best-known explicit approach is the so-called image method, which was developed for dealing with flat boundaries to fourth-order 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 finite-difference approximations: see (MOCZO07) for a review on these subjects.

The aim of this paper is to present a finite-difference approach accounting for free boundaries without introducing the aforementioned drawbacks of the vacuum and image methods. The requirements are as follows: smooth arbitrary-shaped 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 high-order Taylor expansions of the boundary values of the solution. Estimating these boundary values involves some mathematical background, in order to compute the high-order boundary conditions; to determine a minimal set of independent boundary values; lastly, to perform a least-square numerical estimation of this minimal set. To help the reader, subroutines in FORTRAN are proposed freely at the web page http://w3lma.cnrs-mrs.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 staggered-grid schemes are used. Single-grid finite-difference schemes are therefore chosen, where all the unknowns are computed at the same grid nodes. Our numerical experiments are based on the high-order 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 S-wavelength 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 staggered-grid schemes.

This paper is organized as follows. Section 2 deals with the continuous problem: the high-order 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 in-plane and two-dimensional, 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.

Figure 1: Boundary between a solid and vacuum.

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 S-wave velocities are and . A velocity-stress 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 first-order hyperbolic system (VIRIEUX86)


2.2 High-order boundary conditions

At any point on the free surface (Figure 1), the stress tensor satisfies the homogeneous Dirichlet conditions . These zero-th order boundary conditions are written compactly


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 first-order spatial derivatives of , two tasks are performed. First, the zeroth-order boundary conditions (2) are differentiated in terms of . The time derivative is replaced by spatial derivatives using the conservation laws (1), which gives


Secondly, the zeroth-order boundary conditions (2) are differentiated in terms of the parameter describing . The chain-rule gives


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 first-order


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)


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 so-obtained up to the -th order can be written compactly




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)



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


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, two-step, and spatially-centred finite-difference scheme. A review of the huge body of literature on finite-differences is given in LEV90 and MOCZO07.

Here we propose to use ADER schemes, that allow to reach easily arbitrary high-order of time and space accuracy (MUNZ05). On Cartesian grids, these finite-volume integration schemes originally developed for aeroacoustic applications are equivalent to finite-difference Lax-Wendroff-type integration schemes (LORCHER06). In the numerical experiments described in section 4, we use a fourth-order ADER integration scheme. This scheme is stable under the Courant-Friedrichs-Lewy (CFL) condition in 2D; as usually with single-grid schemes, it is slightly dissipative (MUNZ05).

Many other single-grid schemes can be used in this context. In particular, the method described in the next subsections has been successfully combined with flux-limiter schemes (LEV90) and with the standard second-order Lax-Wendroff scheme. Difficulties have been encountered with dissipative-free schemes based on centred staggered-grid finite-difference schemes, as we will see in section 3.6.

3.2 Use of fictitious values

Figure 2: Determination of the fictitious value required for time-marching at neighboring grid nodes. P is the orthogonal projection of on . The grid nodes in and inside the circle centred at with a radius are denoted by +.

Time-marching at grid-points 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 Taylor-like extrapolation


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


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


where the components of are real numbers. Injecting (12) into (9) gives


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


The set of equations (14) is written compactly via a matrix


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 least-squares solution


where the matrix denotes the pseudo-inverse of . From (10), (13) and (16), the fictitious value in the vacuum at is


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 high-order 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 time-marching procedure. Useful comments are proposed about this method:

  1. The radius of must ensure that the number of equations in (15) is greater than the number of unknowns:


    No theoretical results are available about the optimal value of . However, numerical studies have shown that a definite overestimation ensures long-term 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 .

  2. Since the boundary conditions do not vary with time, the extrapolators in (17) can be computed and stored during a pre-processing step. At each time step, only small matrix-vector products are required.

  3. 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.

  4. 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 long-time simulations (see section 4.5).

  5. In a previous one-dimensional 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 fourth-order ADER scheme.

  6. 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.cnrs-mrs.fr/~MI/Software/.

3.6 Case of staggered-grid schemes

(a) (b)
Figure 3: Staggered-grid schemes with a plane boundary parallel to the meshing: two cases can be distinguished, depending on the position of relative to the meshing. Case (a), where the fictitious stress is estimated, works well, while case (b), where the fictitious velocity is estimated, leads to long-term instabilities.

Instead of using a single-grid scheme as proposed in section 3.1, readers may be interested in adapting our approach to staggered-grid 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 long-time integration is considered.

To understand why this is so, let us consider PSS-2. 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 PSS-4, contrary to the vacuum method. Using 10 grid nodes per minimal S-wavelength 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 first-order boundary conditions (3) and higher-order conditions (see section 2.2), but they do not involve the fundamental zeroth-order 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 PSS-4, 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 staggered-grid schemes, especially PSS-2, except in the trivial case sketched in Figure 3-(a).

Figure 4: Boundary with a corner at , replaced locally by an arc of circle with a radius between and .
Figure 5: Test 1: snapshots of at the initial instant (a), at mid-term (b) and at the final instant (c).
Figure 6: Test 1: time history of (a). Zooms on successive time windows, with various discretizations (b,c,d): the number after # denotes the number of grid nodes per minimal S-wavelength.

3.7 Case of non-smooth 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:

  1. 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.

  2. 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 space-time 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


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 S-wavelength. 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

Figure 7: Test 2: snapshot of at the final instant (a). Numerical and exact time histories of (b).
Figure 8: Test 2: parametric study of the relative error in terms of the boundary’s angle , with various discretizations. The number after # denotes the number of grid nodes per minimal S-wavelength.

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 rightward-moving plane wave, with Hz, ensuring 22 grid nodes per minimal P-wavelength and 10 grid nodes per minimal S-wavelength at that frequency.

During the pre-processing 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 S-wavelength. 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

Figure 9: Test 3: snapshots of , for various sinusoidal topographies.
Figure 10: Test 4: snapshot of (a) and time history of the mechanical energy (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 P-wavelength, and 9 grid nodes per minimal S-wavelength 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 well-known Cagniard-de 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 S-wavelength. 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 arbitrary-shaped boundaries are encountered (BOHLEN06).

4.4 Test 3: sinusoidal boundary

Since boundaries not related to the finite-difference grid can be included, the third test is performed on a sinusoidal free boundary, with a peak-to-peak 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 S-wavelength at the frequency of the source wavelet, even in the case of complex topographies with variable curvatures.

4.5 Test 4: long-term stability

The fourth test focuses on long-term 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 long-term 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


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 so-obtained. 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 time-domain single-grid finite-difference 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 high-order 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 fourth-order ADER scheme is used on a propagation distance of 50 minimal wavelengths, 10 grid nodes per minimal S-wavelength yield to a very good level of accuracy. With 5 grid nodes per minimal S-wavelength, the solution is less accurate but still acceptable.

For the sake of simplicity, we have dealt here with academic cases, considering two-dimensional geometries, constant physical parameters, and simple elastic media. Let us examine briefly the generalization of our approach to more realistic configurations:

  1. Extending the method to 3-D 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.

  2. 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 high-order 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:

  3. 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


    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 low-frequency range (DAI95), where the evolution equations can be put in the form (21).

The authors thank the reviewers P. Moczo and I. Oprsal for their instructive comments and bibliographic insights.


This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.