# SMART layers: a simple and robust alternative to PML approaches for elastodynamics

Josue Tago, Ludovic Métivier, Jean Virieux

<figure><img src="image_1.png" /></figure>

## To cite this version:

Josue Tago, Ludovic Métivier, Jean Virieux. SMART layers: a simple and robust alternative to PML approaches for elastodynamics. Geophysical Journal International, 2014, 199 (2), pp.700-706. 10.1093/gji/ggu298. hal-01072202

HAL Id: hal-01072202

https://hal.science/hal-01072202

Submitted on 6 Jan 2022

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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SMART layers: a simple and robust alternative to PML approaches for elastodynamics

J. Tago,<fnref n="1" />,* L. Métivier<fnref n="1" />,<fnref n="2" /> and J. Virieux<fnref n="1" />

Accepted 2014 July 25. Received 2014 July 22; in original form 2014 May 7

## SUMMARY

For considering elastic seismic wave propagation in large domains, efficient absorbing boundary conditions are required with numerical modelling in finite domains. Since their introduction by Bérenger, the perfectly matched layers (PML) has become the state-of-the art method because of its efficiency and ease of implementation. However, for anisotropic media, theoretical analysis and numerical experiments show that the PML method is amplifying, that is it exhibits numerical instabilities. Numerical experiments can also exhibit numerical instabilities of the PML when dealing with long time simulations even for isotropic media, especially for finite element methods in unstructured grids. Recently, a new method, called SMART layers approach, has been proposed. This method is shown to be stable even for anisotropic media. The drawback is that the SMART layers are not perfectly matched. We have implemented this new approach in a discontinuous Galerkin method and we illustrate that this method does not exhibit numerical instabilities while PML do for an isotropic elastodynamic simulation. We show that this approach is also competitive with respect to the PML method in terms of efficiency and computational cost.

**Key words:** Numerical solutions; Computational seismology; Wave propagation.

## 1 INTRODUCTION

For seismic wave propagation, we often consider the Earth as a semi-infinite medium. However, for the wavefield numerical computation, the domain of computation must be finite and as small as possible. Hence, we need to design absorbing boundary conditions (ABC) or absorbing layers at the edges of our numerical domain to avoid spurious reflections. The first-order ABC were introduced by Clayton & Engquist (1977). For planar waves perpendicular to a rectilinear edge, outgoing waves are fully absorbed. For 2-D and 3-D geometries, the waves arriving at oblique incidences are only partially absorbed and significant spurious reflections are observed at grazing angles. High-order ABC have been proposed to improve the performance of the first-order conditions. For example, Collino (1993) have proposed fractional high-order derivatives, but their implementation is not trivial, and the method is computationally intensive.

The alternative for the absorption of outgoing waves is the introduction of a layer beyond the edge where the wave is damped. In this strategy, the outgoing waves are slowly attenuated as they go deeper into the layer, while preventing any induced incoming wave.

Cerjan et al. (1985) have introduced this strategy where a damping profile is designed inside the sponge layer with a significant efficient performance for grazing angles. Unfortunately, this technique induces also unwanted reflections between the domain of interest and the absorbing layer. A careful design of the damping profile is needed for the mitigation of these reflections, leading to rather thick layers. An improved method has been proposed by Bérenger (1994) for the Maxwell’s equation as the impedance matching is perfect between the domain of interest and the absorbing layer in the continuum media: the reflection coefficient is exactly equal to zero at this interface. Unfortunately, after the domain discretization for numerical simulation, this property disappears and small reflections are generated. Even so, this approach turns out to be more efficient than any precedent strategy and since then they have been continuously developed in different domains of wave simulation.

For elastodynamics it has been used and investigated exhaustively (e.g. Chew & Liu 1996; Hastings et al. 1996; Collino & Tsogka 2001; Komatitsch & Tromp 2003). A successful modification to the perfectly matched layer (PML) strategy is the convolutional PML (C-PML), introduced by Komatitsch & Martin (2007), which improves the absorption at grazing angles and reduces the memory requirements. A strategy proposed to improve the stability of the PMGs is the modified PML (M-PML) method, which introduces additional absorbing terms in the directions tangential to the interface between the domain of interest and the absorbing layer. The

<footnote n="1">¹ISTerre, University of Grenoble Alpes, Grenoble, France. E-mail: josue.tago@gmail.com</footnote>

<footnote n="2">²LJK, CNRS, University of Grenoble Alpes, Grenoble, France</footnote>

<footnote>*Now at: Faculty of Engineering, National Autonomous University of Mexico, Mexico City, Mexico.</footnote>
M-PML strategy loses the perfectly matched property, even at the continuous level, but it has shown to be still efficient for practical applications, at the expense of making the PML thicker (Dmitriev & Lisitsa 2011).

Etienne et al. (2010) applied the C-PML in a discontinuous Galerkin method (DGM) based on an unstructured meshing domain description and reported long-time instabilities, like those reported by Meza-Fajardo & Papageorgiou (2008) using PML. Etienne et al. (2010) introduced additional damping in the directions parallel to the interface mixing the C-PML and the M-PML strategies which allows a time delay in the occurrence of these instabilities. This new strategy was called modified C-PML (M-CPML) and, even though it has proved to be successful in delaying the instabilities, it is unable to fully prevent them.

Many studies on the accuracy of PML approaches have been done (e.g. Martin et al. 2008; Meza-Fajardo & Papageorgiou 2010, 2012). Recently, Zhinan et al. (2014) improved the C-PML stability by using a singularity removal process. However for high-order space approximations, instabilities could still appear in practical situations, even if they can be delayed by using a stretching factor.

As we have observed that these instabilities may dramatically complicate the design of the numerical modelling (different meshes, absorbing layer size and damping profiles should be tried for avoiding these instabilities in the working time window), we investigate an alternative strategy based on the so-called SMART absorbing layers with nice mathematical properties of bounded discrete energy over time introduced by Halpern et al. (2011). Inside the SMART layer, the wavefield is locally decomposed into components propagating inward and outward with respect to the boundary against which waves should be attenuated. A selective damping is applied progressively for outgoing waves as they enter inside the SMART layer. This specific damping of outgoing waves is performed by an additional zero-order operator extracted from the decomposition of the operator related to the first-order hyperbolic system. Unlike PML, SMART layers are intrinsically dissipative as soon as the hyperbolic system on which they are applied satisfies a symmetrizability condition. Métivier et al. (2014a) have demonstrated the robustness of the SMART layers for anisotropic acoustic media. They prove that the associated first-order hyperbolic system is symmetrizable, and they construct the corresponding zero-order dissipative term. Finite-differences numerical examples are provided. Stable results are obtained with the SMART layer while the PMLs exhibit an amplifying behaviour.

We show how to use this strategy with DGM with unstructured meshes. We concentrate on 2-D elastodynamics but the extension to 3-D geometries does not introduce new features. We shall show that the first order hyperbolic system for elastodynamics satisfies also the symmetrizability condition required for the SMART approach. We present how to compute the zero-order term introduced when applying the SMART layers. Then we shall show how they can be introduced in a DGM approach. We shall compare the SMART layers versus the C-PML approach in a specific example where long time instabilities appear.

# 2 SMART LAYERS FOR ELASTODYNAMICS AND THE DGM

## 2.1 SMART equations

The 2-D elastic wave propagation can be modelled with a velocity-stress formulation as a first-order hyperbolic system. Following the

approach of Ben Jemaa et al. (2007), the elastodynamic system can be written in the following pseudo-conservative form

$$ \rho \partial_t \vec{v} = \sum_{\theta \in \{x,y\}} \partial_{\theta} \mathbf{M}_{\theta} \vec{\sigma} + \vec{f} \quad (1) $$

$$ C \partial_t \vec{\sigma} = \sum_{\theta \in \{x,y\}} \partial_{\theta} \mathbf{N}_{\theta} \vec{v}, \quad (2) $$

with the definitions of the velocity and stress vectors as

$$ \vec{v} = [v_x \quad v_y]^T \quad (3) $$

$$ \vec{\sigma} = [\tau \quad \tau' \quad \sigma_{xy},]^T, \quad (4) $$

with

$$ \tau = \frac{1}{2}(\sigma_{xx} + \sigma_{yy}) \quad (5) $$

$$ \tau' = \frac{1}{2}(\sigma_{xx} - \sigma_{yy}). \quad (6) $$

The structure of the elastodynamic system (eqs 1 and 2) involves a change of variables for the stresses (eq. 4) and it also induces a diagonal stiffness matrix. This particular structure allows us to easily compute the inverse of the stiffness matrix, i.e. the compliance matrix C, in terms of the Lamé parameters $\lambda$ and $\mu$

$$ C = \mathrm{diag}(\lambda + \mu)^{-1}, \mu^{-1}, \mu^{-1}). \quad (7) $$

The right-hand side of eq. (2) is left without the physical properties (Burridge 1996) and the matrices related with the differential operators defined are

$$ \mathbf{M}_x = \begin{bmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

$$ \mathbf{M}_y = \begin{bmatrix} 0 & 0 & 1 \\ 1 & -1 & 0 \end{bmatrix} $$

$$ \mathbf{N}_x = \begin{bmatrix} 1 & 0 \\ 1 & 0 \\ 0 & 1 \end{bmatrix} $$

$$ \mathbf{N}_y = \begin{bmatrix} 0 & 1 \\ 0 & -1 \\ 1 & 0 \end{bmatrix}. \quad (8) $$

The elastodynamic system, eqs (1) and (2), can be written in the form of a general hyperbolic system

$$ S \partial_t \vec{u} + \sum_{\theta \in \{x,y\}} \mathbf{A}_{\theta} \partial_{\theta} \vec{u} = \vec{g}, \quad (9) $$

with the velocities and stresses gathered in

$$ \vec{u} = [v_x \quad v_y \quad \tau \quad \tau' \quad \sigma_{xy}]^T, \quad (10) $$

the generalized compliance matrix

$$ S = \mathrm{diag}(\rho, \rho, (\lambda + \mu)^{-1}, \mu^{-1}, \mu^{-1}), \quad (11) $$
the matrices related with the differential operators

$$
\mathbf{A}_x = \begin{bmatrix} 0 & 0 & -1 & -1 & 0 \\ 0 & 0 & 0 & 0 & -1 \\ -1 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 \end{bmatrix}
$$
$$
\mathbf{A}_y = \begin{bmatrix} 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & -1 & 1 & 0 \\ 0 & -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 \end{bmatrix},
$$

(12)

and the source term

$$
\vec{\mathbf{g}} = [ f_x \quad f_y \quad 0 \quad 0 \quad 0 ]^T.
$$

(13)

The system (9) is equivalent to

$$
\partial_t \vec{\mathbf{u}} + \sum_{\theta \in \{x, y\}} \vec{\mathbf{A}}_{\theta} \partial_{\theta} \vec{\mathbf{u}} = \mathbf{S}^{-1} \vec{\mathbf{g}},
$$

(14)

where

$$
\vec{\mathbf{A}}_{\theta} = \mathbf{S}^{-1} \mathbf{A}_{\theta}, \quad \theta \in \{x, y\}
$$

(15)

As S is diagonal with strictly positive entries and the matrices **A**<sub>x</sub>, **A**<sub>y</sub> are symmetric, we see that S is a symmetrizer for the system (14). This ensures the SMART layer will dissipate the energy of the solution, with no possibility of spurious amplification (see Métivier et al. 2014a,b,c, for details).

The SMART layer consists in the addition of a zero-order term, called SMART term, to the general hyperbolic system (14) that will damp the outgoing wavefield. For the smart selection of the damped wavefields an eigendecomposition is done for the space linear operators $\vec{\mathbf{A}}_\theta$, for $\theta \in \{x, y\}$. The SMART term projects the wavefield in the eigenspace span by the eigenvectors associated to the different eigenvalues of $\vec{\mathbf{A}}_\theta$. It only damps the outgoing part of the wavefield. The SMART term construction and addition is explained ahead.

The linear operator $\vec{\mathbf{A}}_x$ has five eigenvalues of multiplicity 1 with four non-zero eigenvalues, related to the P-wave velocity, $\alpha$, and the S-wave velocity, $\beta$. The fifth one is zero. The four non-zero eigenvalues with their corresponding eigenvector are given in Table 1. The linear operator $\vec{\mathbf{A}}_y$ also has five eigenvalues of multiplicity 1, with four non-zero eigenvalues. They are given in Table 2 with their corresponding eigenvectors. These eigenvectors are not identical to those of $\vec{\mathbf{A}}_x$.

Non-zero eigenvalues in Tables 1 and 2 correspond to P- and S-waves travelling in positive and negative direction over the x- and y-coordinates for $\vec{\mathbf{A}}_x$ and $\vec{\mathbf{A}}_y$, respectively. Since each of them has

<table><caption>**Table 1.** Non-zero eigenvalues and their corresponding eigenvector of $\vec{\mathbf{A}}_x$.</caption>
  <thead>
    <tr>
      <th>Eigenvalues</th>
      <th>Eigenvectors</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>e<sup>+</sup><sub>xP</sub> = &alpha;</td>
      <td>v&#773;<sup>+</sup><sub>xP</sub> = [1 0 -(&lambda; + &mu;)/&alpha; -&mu;/&alpha; 0]<sup>T</sup></td>
    </tr>
    <tr>
      <td>e<sup>-</sup><sub>xP</sub> = -&alpha;</td>
      <td>v&#773;<sup>-</sup><sub>xP</sub> = [1 0 (&lambda; + &mu;)/&alpha; &mu;/&alpha; 0]<sup>T</sup></td>
    </tr>
    <tr>
      <td>e<sup>+</sup><sub>xS</sub> = &beta;</td>
      <td>v&#773;<sup>+</sup><sub>xS</sub> = [0 -&beta;/&mu; 0 0 1]<sup>T</sup></td>
    </tr>
    <tr>
      <td>e<sup>-</sup><sub>xS</sub> = -&beta;</td>
      <td>v&#773;<sup>-</sup><sub>xS</sub> = [0 &beta;/&mu; 0 0 1]<sup>T</sup></td>
    </tr>
  </tbody>
</table>

<table><caption>**Table 2.** Non-zero eigenvalues and their corresponding eigenvector of $\vec{\mathbf{A}}_y$.</caption>
  <thead>
    <tr>
      <th>Eigenvalues</th>
      <th>Eigenvectors</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>e<sup>+</sup><sub>yp</sub> = &alpha;</td>
      <td>v&#773;<sup>+</sup><sub>yp</sub> = [0 1 -(&lambda; + &mu;)/&alpha; &mu;/&alpha; 0]<sup>T</sup></td>
    </tr>
    <tr>
      <td>e<sup>-</sup><sub>yp</sub> = -&alpha;</td>
      <td>v&#773;<sup>-</sup><sub>yp</sub> = [0 1 (&lambda; + &mu;)/&alpha; - &mu;/&alpha; 0]<sup>T</sup></td>
    </tr>
    <tr>
      <td>e<sup>+</sup><sub>ys</sub> = &beta;</td>
      <td>v&#773;<sup>+</sup><sub>ys</sub> = [-&beta;/&mu; 0 0 0 1]<sup>T</sup></td>
    </tr>
    <tr>
      <td>e<sup>-</sup><sub>ys</sub> = -&beta;</td>
      <td>v&#773;<sup>-</sup><sub>ys</sub> = [&beta;/&mu; 0 0 0 1]<sup>T</sup></td>
    </tr>
  </tbody>
</table>

multiplicity 1, the projectors for each eigenvalue can be computed
as

$$
\mathbf{B}_{\theta p}^{\pm} = \frac{\vec{\nu}_{\theta p}^{\pm} (\vec{\nu}_{\theta p}^{\pm})^{T}}{\|\vec{\nu}_{\theta p}^{\pm}\|^{2}}, \quad \mathbf{B}_{\theta S}^{\pm} = \frac{\vec{\nu}_{\theta S}^{\pm} (\vec{\nu}_{\theta S}^{\pm})^{T}}{\|\vec{\nu}_{\theta S}^{\pm}\|^{2}}, \quad \forall \theta \in \{x, y\}.
$$

(16)

The SMART term **B**(**x**) is a linear combination of these projectors
weighted by damping coefficients $d_{\theta}^{\pm}$ such that

$$
\mathbf{B}(\vec{\mathbf{x}}) = \sum_{\theta \in \{x, y\}} d_{\theta}^{+}(\vec{\mathbf{x}}) \sum_{y \in \{P, S\}} B_{\theta y}^{+} + \sum_{\theta \in \{x, y\}} d_{\theta}^{-}(\vec{\mathbf{x}}) \sum_{y \in \{P, S\}} B_{\theta y}^{-}.
$$

(17)

The damping coefficients $d_{\theta}^{\pm}(\vec{\mathbf{x}})$ are smooth mono-dimensional functions which are zero in the domain of interest and have a polynomial growth inside the layers.

As well as Etienne et al. (2010), we choose the damping profile
proposed by Collino & Tsogka (2001) for PMLs such that

$$
d_{\theta}^{\pm}(\vec{\mathbf{x}}) = d_{\theta_{\max}}^{\pm} \left( \frac{\delta_{\theta}^{\pm}(\vec{\mathbf{x}})}{L_{\text{smart}}} \right)^{2} \quad \forall \theta \in \{x, y\},
$$

(18)

with

$$
d_{\theta_{\max}}^{\pm} = -3\alpha \frac{\log(R_{\text{coeff}})}{2L_{\text{smart}}} \quad \forall \theta \in \{x, y\},
$$

(19)

where $\delta_{\theta}^{\pm}(\vec{\mathbf{x}})$ is the depth inside the SMART layers in the $\theta$ coordinate, $L_{\text{smart}}$ is the SMART layer thickness and $R_{\text{coeff}}$ is the theoretical reflection coefficient.

Following Halpern et al. (2011), adding the SMART term **B**, eq.
(17), in the general hyperbolic system, eq. (9), leads to the equation

$$
S \partial_t \vec{\mathbf{u}} + \sum_{\theta \in \{x,y\}} \mathbf{A}_{\theta} \partial_{\theta} \vec{\mathbf{u}} + S B \vec{\mathbf{u}} = \vec{\mathbf{g}},
$$

(20)

Since we would rather work with the velocity-stress formulation
given in eqs (1) and (2), we split the damping term **Z** = **S****B** from
eq. (20) into

$$
\mathbf{Z} = \begin{bmatrix} \mathbf{D} & \mathbf{E} \\ \mathbf{G} & \mathbf{H} \end{bmatrix},
$$

(21)

with

$$
\mathbf{D} = \begin{bmatrix} z_{1,1} & z_{1,2} \\ z_{2,1} & z_{2,2} \end{bmatrix} \quad \mathbf{E} = \begin{bmatrix} z_{1,3} & z_{1,4} & z_{1,5} \\ z_{2,3} & z_{2,4} & z_{2,5} \end{bmatrix}
$$

(22)

$$
\mathbf{G} = \begin{bmatrix} z_{3,1} & z_{3,2} \\ z_{4,1} & z_{4,2} \\ z_{5,1} & z_{5,2} \end{bmatrix} \quad \mathbf{H} = \begin{bmatrix} z_{3,3} & z_{3,4} & z_{3,5} \\ z_{4,3} & z_{4,4} & z_{4,5} \\ z_{5,3} & z_{5,4} & z_{5,5} \end{bmatrix}.
$$

We thus incorporate the SMART term in our velocity-stress formulation, eqs (1) and (2), as

$$
\rho \partial_t \vec{\nu} = \sum_{\theta \in \{x, y\}} \partial_{\theta} \mathbf{M}_{\theta} \vec{\sigma} - \mathbf{D} \vec{\nu} - \mathbf{E} \vec{\sigma} + \vec{\mathbf{f}} \quad (23)
$$
$$C \partial_t \vec{\sigma} = \sum_{\theta \in \{x, y\}} \partial_{\theta} \mathbf{N}_{\theta} \vec{\mathbf{v}} - \mathbf{H} \vec{\sigma} - G \vec{\mathbf{v}}. \quad (24)$$

We shall turn our attention how to solve this new system, eqs (23) and (24), by a DGM.

## 2.2 Discontinuous Galerkin method

To solve the hyperbolic system, eqs (23) and (24), we first decompose the domain $\Omega$ into $K$ elements, such that

$$\Omega \approx \Omega_h = \sum_{i=1}^{K} D_i, \quad (25)$$

where each $D_i$ is a straight-sided triangle. The sum of the triangles form the mesh which has to be geometrically conforming. Following Etienne et al. (2010) we apply a nodal DGM to the weak formulation of eqs (23) and (24) to obtain

$$\begin{aligned} \rho_i (\mathbf{I}_2 \otimes \mathbf{K}_i) \frac{\vec{\mathbf{v}}_i^{n+\frac{1}{2}} - \vec{\mathbf{v}}_i^{-\frac{1}{2}}}{\Delta t} & = - \sum_{\theta \in \{x, y\}} (\mathbf{M}_{\theta} \otimes \mathbf{E}_{i\theta}) \vec{\sigma}_i^n + \sum_{k \in N_i} \left[ (\mathbf{P}_{ik} \otimes \mathbf{F}_{ik}) \frac{\vec{\sigma}_i^n + \vec{\sigma}_k^n}{2} \right] \\ & - (\mathbf{D}_i \otimes \mathbf{K}_i) \vec{\mathbf{v}}_i^n - (\mathbf{E}_i \otimes \mathbf{K}_i) \vec{\sigma}_i^n + (\mathbf{I}_2 \otimes \mathbf{K}_i) \vec{\mathbf{r}}_i^{n+\frac{1}{2}} \end{aligned} \quad (26)$$

$$\begin{aligned} (\mathbf{C}_i \otimes \mathbf{K}_i) \frac{\vec{\sigma}_i^{n+1} - \vec{\sigma}_i^n}{\Delta t} & = - \sum_{\theta \in \{x, y\}} (\mathbf{N}_{\theta} \otimes \mathbf{E}_{i\theta}) \vec{\mathbf{v}}_i^{n+\frac{1}{2}} + \sum_{k \in N_i} \left[ (\mathbf{Q}_{ik} \otimes \mathbf{F}_{ik}) \frac{\vec{\mathbf{v}}_i^{n+\frac{1}{2}} + \vec{\mathbf{v}}_k^{n+\frac{1}{2}}}{2} \right] \\ & - (\mathbf{H}_i \otimes \mathbf{K}_i) \vec{\sigma}_i^{n+\frac{1}{2}} - (\mathbf{G}_i \otimes \mathbf{K}_i) \vec{\mathbf{v}}_i^{n+\frac{1}{2}}. \end{aligned} \quad (27)$$

In eqs (26) and (27), $\otimes$ represent the tensor product, $N_i$ is the group of adjacent elements to the i-triangle and the superscript n indicates the time step. The involved matrices are the mass matrix, $\mathbf{K}_i$, the stiffness matrices, $\mathbf{E}_{i\theta}$, for all $\theta \in \{x, y\}$, the flux matrix, $\mathbf{F}_{ik}$ and the auxiliary flux matrices, $\mathbf{P}_{ik}$ and $\mathbf{Q}_{ik}$. The definitions of these matrices can be found, in Etienne et al. (2010) or Tago et al. (2012). The SMART matrices, $\mathbf{D}_i$, $\mathbf{E}_i$, $\mathbf{H}_i$ and $\mathbf{G}_i$ depends on the position and physical properties of the i-triangle.

The proposed DGM, eqs (26) and (27), uses a second-order explicit leap-frog scheme for the time integration and a centred scheme for the fluxes evaluation. The combination of this two makes the DGM energy conservative, efficient and easy to program. However more sophisticated DGM schemes have been successfully proposed, as the ADER-DG proposed by Käser & Dumbser (2006) who uses a modal representation, upwind fluxes and the ADER time integration approach.

The incorporation of the SMART layers makes the right-hand side (RHS) term of the velocity scheme, eq. (26), dependent not only on the stresses but also on the velocities. The leap-frog scheme requires the velocities computed at the n-time, $\vec{\mathbf{v}}_i^n$. A simple approximation is

$$\vec{\mathbf{v}}_i^n = \frac{1}{2} \left( \vec{\mathbf{v}}_i^{n+\frac{1}{2}} + \vec{\mathbf{v}}_i^{-\frac{1}{2}} \right). \quad (28)$$

This simple choice makes the time integration semi-implicit. However, because of the local nature of DGM, we can invert the involved

small local matrices to recover an explicit time integration through the following:

$$\begin{aligned} & \left[ (\mathbf{I}_2 \otimes \mathbf{I}_{\mathrm{dof}}) + \frac{\Delta t}{2\rho_i} (\mathbf{D} \otimes \mathbf{I}_{\mathrm{dof}}) \right] \vec{\mathbf{v}}_i^{n+\frac{1}{2}} \\ & = \vec{\mathbf{v}}_i^{n-\frac{1}{2}} + \frac{\Delta t}{\rho_i} (\mathbf{I}_2 \otimes \mathbf{K}_i)^{-1} \vec{\mathbf{r}}_i^n \\ & - \frac{\Delta t}{2\rho_i} (\mathbf{D}_i \otimes \mathbf{I}_{\mathrm{dof}}) \vec{\mathbf{v}}_i^{n-\frac{1}{2}} - \frac{\Delta t}{\rho_i} (\mathbf{E}_i \otimes \mathbf{I}_{\mathrm{dof}}) \vec{\sigma}_i^n, \end{aligned} \quad (29)$$

where the subscript dof stands for the degrees of freedom inside the i-triangle and

$$\vec{\mathbf{r}}_i^n = - \sum_{\theta \in \{x, y\}} (\mathbf{M}_{\theta} \otimes \mathbf{E}_{i\theta}) \vec{\sigma}_i^n + \sum_{k \in N_i} \left[ (\mathbf{P}_{ik} \otimes \mathbf{F}_{ik}) \frac{\vec{\sigma}_i^n + \vec{\sigma}_k^n}{2} \right]. \quad (30)$$

The treatment of the RHS of the stress scheme, eq. (27), follows the same path as the one done for the velocity scheme, eq. (26). The needed stress at the $n + \frac{1}{2}$-time is approximated with the expression

$$\vec{\sigma}_i^{n+\frac{1}{2}} = \frac{1}{2} (\vec{\sigma}_i^{n+1} + \vec{\sigma}_i^n). \quad (31)$$

The stress scheme, eq. (27), is transformed into the explicit time integration

$$\begin{aligned} & \left[ (\mathbf{I}_3 \otimes \mathbf{I}_{\mathrm{dof}}) + \frac{\Delta t}{2} (\mathbf{C}_i^{-1} \mathbf{H}_i \otimes \mathbf{I}_{\mathrm{dof}}) \right] \vec{\sigma}_i^{n+1} \\ & = \vec{\sigma}_i^n + \Delta t (\mathbf{C}_i \otimes \mathbf{K}_i)^{-1} \vec{\mathbf{o}}_i^{n+\frac{1}{2}} - \frac{\Delta t}{2} (\mathbf{C}_i^{-1} \mathbf{H}_i \otimes \mathbf{I}_{\mathrm{dof}}) \vec{\sigma}_i^n \\ & - \Delta t (\mathbf{C}_i^{-1} \mathbf{G}_i \otimes \mathbf{I}_{\mathrm{dof}}) \vec{\mathbf{v}}_i^{n+\frac{1}{2}}, \end{aligned} \quad (32)$$

with

$$\vec{\mathbf{o}}_i^{n+\frac{1}{2}} = - \sum_{\theta \in \{x, y\}} (\mathbf{N}_{\theta} \otimes \mathbf{E}_{i\theta}) \vec{\mathbf{v}}_i^{n+\frac{1}{2}} + \sum_{k \in N_i} \left[ (\mathbf{Q}_{ik} \otimes \mathbf{F}_{ik}) \frac{\vec{\mathbf{v}}_i^{n+\frac{1}{2}} + \vec{\mathbf{v}}_k^{n+\frac{1}{2}}}{2} \right]. \quad (33)$$

The simple approximations for $\vec{\mathbf{v}}_i^n$ and $\vec{\sigma}_i^{n+\frac{1}{2}}$ degrade the convergence behaviour of the leap-frog scheme. Yet since the SMART matrices are different from zero only inside the SMART layers, the impact will be only in these zones.

## 3 NUMERICAL TEST

We perform the numerical comparison between the C-PML and SMART layers with a very simple model, where long-term C-PML instabilities appear with the DGM scheme. The C-PML implementation is done following the work of Etienne et al. (2010). The ordinary differential equations, related with the memory variables introduced, are solved with the same DGM scheme used to update the velocity-stress hyperbolic system. We consider an homogeneous half-space 10 000 m width and 5000 m depth with the physical parameters $\alpha = 4000 \text{ m s}^{-1}$, $\beta = 2310 \text{ m s}^{-1}$ and $\rho = 2000 \text{ kg m}^{-3}$. The source is a Ricker pulse with a dominant frequency of 3 Hz located near the free surface at the middle of the domain $x_s = 0 \text{ m}$ and $z_s = 10 \text{ m}$. The domain is discretized using a conforming unstructured triangular mesh. We used quadratic interpolation functions inside each triangle, P2 elements, so the mesh ensures three elements per minimum wavelength. The total simulation time is 30 s.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Snapshots of $v_x$ while including 10-element width C-PMLs at both lateral sides and the bottom of the model at: (a) 1 s, (b) 2.4 s, (c) 5 s, (d) 15 s, (e) 25 s and (f) 30 s.</figcaption></figure>

## 3.1 C-PML instabilities and SMART comparison

Absorbing boundary layers are defined at the two lateral sides and at the bottom of the domain. On top we apply a free surface boundary condition. The C-PMLs width is 10 elements, as proposed by Etienne et al. (2010) when using a DGM approach. The damping profile are the ones presented in Section 2. The Fig. 1 presents snapshots of $v_x$ at different times of the simulation. We observe that the C-PMLs absorb effectively the outgoing waves. However, a strong amplification coming from the layer appears later. This can also be observed in Fig. 2(a) where we present seismograms for $v_x$ and $v_y$ obtained using receivers located along the free-surface at the same depth than the source.

We reproduce this test using 10 elements SMART layers, keeping the same mesh and damping profile used for the C-PMLs. We show in Fig. 2(b) that the SMART layers do not suffer from the amplification phenomenon. However, the coupling between the simulation domain and the SMART layers is not as satisfactory as the one obtained with a C-PML. Particularly, for waves propagating at grazing angles, we observe stronger reflections at the interface between the domain of interest and the layer (from 3 to 8 s). In order to reduce the magnitude of the reflections at the order of magnitude of the C-PML reflections, we can increase the width of the SMART layers from 10 to 20 elements. The seismograms for 20 elements SMART layers are given in Fig. 2(c).

To prove the robust dissipative property of the SMART layers, we run the same configuration using 20 elements SMART layers but we increase the total simulation time from 30 to 300 s. In Fig. 3, we present the history of the kinematic energy. The decrease is very fast at the beginning of the simulation. After reaching several orders of magnitude below the magnitude of the wavefield propagated, it seems that the noise level is attained. The kinematic energy keeps being attenuated at a small rate. This emphasizes the dissipative property of the SMART layers.

The increase of the number of elements in the SMART layers makes the computation more expensive. This is necessary to reach the same level of accuracy than the one achieved with the C-PMLs. However, as can be seen in Table 3, an additional feature of the SMART layers method is that for a given width of the layer, it requires less memory and less Floating Point Operations (FLOP's) than C-PMLs. The data presented in Table 3 depend on our specific implementation. However, looking at the system eqs (23) and (24), it is clear that independently of the numerical scheme for a

given absorbing layer width, SMART layers are cheaper in memory and time than C-PMLs because it only consists in the introduction of a zero-order term, and do not require additional memory variables with their respective ordinary differential equations. For our implementation, with 20 elements SMART layers (78.16 per cent SMART elements from the total of elements in the mesh), using Table 3, we can compute that they are approximately 1.1 more time consuming and 1.2 more memory expensive than 10 elements C-PMLs (36.32 per cent CPML elements from the total of elements in the mesh). This is the small price we have to pay to avoid the C-PML instabilities that, as can be seen in Figs 2(a), starts polluting the velocities fields a few seconds after the wavefield has reached the layers. These factors are even smaller when the ratio between the domain elements and absorbing boundary elements increases.

## 4 CONCLUSIONS AND PERSPECTIVES

The SMART layers is a robust alternative to PML approaches for elastodynamic simulations. They are easier to implement and for a given absorbing layer width faster than C-PMLs. We showed that they are robust in a case where C-PMLs are unstable. However, since the SMART layers are not perfectly matched, stronger reflections than those observed in the C-PML approach are generated. In order to increase the efficiency of the SMART layers to the same level of the C-PMLs, we had to double the width of the SMART layers. However, this only implies a small cost of time and memory in comparison with the C-PML approach, as the SMART layers only consist in introducing a zero-order term in the initial system.

For 3-D simulations, the SMART layers width enlargement would increment the ratio of SMART layers elements with respect to regular elements. Nonetheless, the gap of the element memory and FLOP's requirements between SMART and C-PMLs increases too. We believe that the sum of both effects would preserve the 2-D relative cost of the SMART layers for 3-D simulations. Furthermore, for anisotropic elastodynamic simulations, the SMART layers relative cost should be smaller than for isotropic elastodynamic simulations. The anisotropy implies cross-derivatives that would increase the amount of memory variables required by the C-PML strategy.

We believe that SMART layers could be very useful for long time 3-D elastodynamics simulations in regional scales for seismic hazard evaluation where C-PML instabilities could appear. Further
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Seismograms of $v_x$ and $v_y$ for 10-elements C-CPML, (a.1) and (a.2), respectively, 10-elements SMART layers, (b.1) and (b.2), respectively, and 20-elements SMART layers, (c.1) and (c.2), respectively. The scale is saturated two orders below the average of the velocity wavefield to better distinguish the reflections.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** History of the L2 norm of the kinematic energy using 20-elements SMART layers which decreases with time until the noise level.</figcaption></figure>
<table><caption>**Table 3.** Memory and FLOP comparison between the C-PML and the SMART layers using P2 elements in our DGM implementation. The ‘Element’ columns consider the whole memory and FLOP’s per P2 element inside the absorbing layer during an update step. While the ‘Abs. layer’ columns involve the memory and FLOP’s required only by the absorbing layer operations by P2 element inside the absorbing layer during an update step.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Memory (bytes)</th><th colspan="2">FLOP</th></tr><tr><th>Element</th><th>Abs. layer</th><th>Element</th><th>Abs. layer</th></tr></thead><tbody><tr><td>C-PML</td><td>468</td><td>256</td><td>5073</td><td>2262</td></tr><tr><td>SMART</td><td>332</td><td>120</td><td>3141</td><td>330</td></tr><tr><td>C-PML/SMART</td><td>1.409</td><td>2.133</td><td>1.615</td><td>6.854</td></tr></tbody></table>

studies to better match the SMART layers should be done to decrease the layer's width.

## ACKNOWLEDGEMENTS

We thank Andreas Fichtner and Roland Martin for their remarks and suggestions. This work has been supported by the French Agence Nationale de la Recherche under the grant ANR-2011-BS56-017 and partially funded by sponsors of the SEISCPE II consortium (http://seiscpe2.osug.fr) and from the European Community's Seventh Framework Programme [FP7/2007-2013] under the grant agreement PEOPLE-2011-IRSES 295217.

## REFERENCES

Ben Jemaa, M., Glinsky-Olivier, N., Cruz-Atienza, V.M., Virieux, J. & Piperno, S., 2007. Dynamic non-planar crack rupture by a finite volume method, *Geophys. J. Int.*, **171**, 271–285.

Bérénger, J.-P., 1994. A perfectly matched layer for absorption of electromagnetic waves, *J. Comput. Phys.*, **114**, 185–200.

Burridge, R., 1996. *Elastic Waves in Anisotropic Media*, Schlumberger-Doll Research.

Cerjan, C., Kosloff, D., Kosloff, R. & Reshef, M., 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equations, *Geophysics*, **50**(4), 2117–2131.

Chew, W.C. & Liu, Q.H., 1996. Perfectly matched layers for elastodynamics: a new absorbing boundary condition, *J. Comput. Acoust.*, **4**, 341–359.

Clayton, R. & Engquist, B., 1977. Absorbing boundary conditions for acoustic and elastic wave equations, *Bull. seism. Soc. Am.*, **67**, 1529–1540.

Collino, F., 1993. High order absorbing boundary conditions for wave propagation models: straight line boundary and corner cases, in *Proceedings of the Second International Conference on Mathematical and Numerical Aspects of Wave Propagation*, pp. 161–171, SIAM, Delaware.

Collino, F. & Tsogka, C., 2001. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, *Geophysics*, **66**, 294–307.

Dmitriev, M.M. & Lisitsa, V.V., 2011. Application of M-PML reflectionless boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part I: reflectivity, *Numer. Anal. Appl.*, **4**(4), 271–280.

Etienne, V., Chaljub, E., Virieux, J. & Glinsky, N., 2010. An hp-adaptive discontinuous Galerkin finite-element method for 3D elastic wave modelling, *Geophys. J. Int.*, **183**(2), 941–962.

Halpern, L., Petit-Bergez, S. & Rauch, J., 2011. The analysis of matched layers, *Confluentes Mathematici*, **3**(2), 159–236.

Hastings, F.D., Schneider, J.B. & Broschat, S.L., 1996. Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation, *J. acoust. Soc. Am.*, **100**, 3061–3069.

Käser, M. & Dumbser, M., 2006. An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes I: the two-dimensional isotropic case with external source terms, *Geophys. J. Int.*, **166**, 855–877.

Komatitsch, D. & Martin, R., 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation, *Geophysics*, **72**(5), SM155–SM167.

Komatitsch, D. & Tromp, J., 2003. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation, *Geophys. J. Int.*, **154**(1), 146–153.

Martin, R., Komatitsch, D. & Gedney, S.D., 2008. A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation, *Comput. Model. Eng. Sci.*, **37**, 274–304.

Métivier, L., Brossier, R., Labbé, S., Operto, S. & Virieux, J., 2014a. A robust absorbing layer for anisotropic seismic wave modeling, *J. Comput. Phys.*, in press.

Métivier, L., Brossier, R., Operto, S. & Virieux, J., 2014b. A robust absorbing layer method for seismic wave simulation in anisotropic media, in *Proceedings of the 76th Annual EAGE Meeting*, Expanded Abstracts, Amsterdam.

Métivier, L., Brossier, R., Operto, S. & Virieux, J., 2014c. Smart: Robust absorbing layer and s-waves filtering for acoustic anisotropic wave simulation, in *Proceedings of the 84th Annual SEG Meeting*, Expanded Abstracts, Denver.

Meza-Fajardo, K. & Papageorgiou, A., 2008. A nonconvolutional, split-field, perfectly matched layer for wave propagation in isotropic and anisotropic elastic media: stability analysis, *Bull. seism. Soc. Am.*, **98**(4), 1811–1836.

Meza-Fajardo, K.C. & Papageorgiou, A.S., 2010. On the stability of a nonconvolutional perfectly matched layer for isotropic elastic media, *Soil Dyn. Earthq. Eng.*, **30**, 68–81.

Meza-Fajardo, K.C. & Papageorgiou, A.S., 2012. Study of the accuracy of the multiaxial perfectly matched layer for the elastic-wave equation, *Bull. seism. Soc. Am.*, **102**, 2458–2467.

Tago, J., Cruz-Atienza, V., Virieux, J., Etienne, V. & Sánchez-Sesma, F., 2012. A 3D hp-adaptive discontinuous Galerkin method for modelling earthquake dynamics, *J. geophys. Res.*, **117**, B09312, doi:10.1029/2012JB009313.

Zhinan, X., Komatitsch, D., Martin, R. & Matzen, R., 2014. Improved forward wave propagation and adjoint-based sensitivity kernel calculations using a numerically stable finite-element PML, *Geophys. J. Int.*, **198**(3), 1714–1747.