**1. Introduction**

Regarding the numerical modeling of wave propagation in anisotropic media, numerous studies have been implemented using finite difference and spectral methods in the time domain [1–3] to obtain the numerical solution in discrete geological media. These media may contain rocks with various properties, defined as acoustic, elastic, viscoelastic isotropic, and anisotropic. With the flourishing development of computational technology and numerical algorithms, wavefield simulation methods have played a crucial role in imaging subsurface characteristics. Various imaging methods, such as tomography [4,5], reverse time migration [6–8], and full-waveform inversion [9–11], require an efficient wavefield simulation method to generate synthetic waveform data to match the observed data.

Accurately performing wavefield modeling with high efficiency in complex media, which contains anisotropy and attenuation, is critical. In contrast to time-domain wavefield modeling, although it is challenging to solve the sparse linearized matrix, multi-source modeling can be performed efficiently, and attenuation can be easily introduced [12]. In

**Citation:** Yang, J.; He, X.; Chen, H. Processing the Artificial Edge-Effects for Finite-Difference Frequency-Domain in Viscoelastic Anisotropic Formations. *Appl. Sci.* **2022**, *12*, 4719. https://doi.org/ 10.3390/app12094719

Academic Editors: Guofeng Liu, Xiaohong Meng and Zhifu Zhang

Received: 22 March 2022 Accepted: 5 May 2022 Published: 7 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

addition, the physical parameters of the reservoir can be estimated by seismic inversion and imaging methods within a limited number of frequency components [13]. Operto et al. [14] developed frequency-domain finite-difference viscoacoustic wavefield modeling in tilted transversely isotropic (TTI) media. Jeong et al. [15] used the finite element method to simulate the frequency-domain wavefield of vertically transversely isotropic (VTI) media. Zhou et al. [16] implemented seismic wave modeling in elastic anisotropic media using the generalized stiffness reduction method. Yang et al. [17] developed a new version of the generalized stiffness reduction method for anisotropic media and incorporated it into the 2.5D spectral element method. Qiao et al. [3] derived attenuative anisotropic wave equations involving fractional time derivatives and solved them using the pseudospectral method. Reducing artificial reflections is still essential despite numerous frequency-domain modeling methods.

Figure 1 shows a sketch of the absorbing boundary condition (ABC). The middle area is the computational area, and the dashed-dot area is the truncated boundary. When the wavefield propagates to the boundary, the ABC needs to attenuate the reflections inwards. Since Berenger [18] developed a perfectly matched layer (PML) to reduce artificial reflections of the electromagnetic wavefield, many studies have shown success in its applications to acoustic and elastic wave modeling [19,20], and improvements have also been made to the perfectly matched layer method [21]. The PML method has become popular for wave modeling problems because it is applicable to first-order wave equations in the time domain [22] and frequency domain [23]. However, it is not straightforward how to implement it for second-order wave equations [24]. In addition, the most severe problem of PML may exist in absorbing artificial reflections in elastic anisotropic media. This is because PML can only function effectively when all slowness vector components of the waves are parallel to their corresponding group velocity vectors [25]. As is known, an elastic anisotropic wavefield generates multiple modes of waves and becomes more complex when the PML becomes unstable in dealing with the edge-effect problem. Thus, we need to replace PML in wavefield simulations of elastic anisotropic media.

**Figure 1.** Sketch map of the absorbing boundary condition (ABC). The square area in the middle is the computational area. The thick dashed-dot line represents the ABC applied outside the computational area.

A multi-axial perfectly matched layer (M-PML) boundary condition was developed in the time domain [26], and its stability has been verified in time-domain isotropic elastic modeling [27,28]. In this paper, the M-PML boundary condition is applied to frequencydomain finite-difference modeling in several elastic anisotropic media. The key point of this method is based on a more general version of coordinate stretching for PML, where waves are absorbed in all directions because the damping factor is specified in more than one direction. Then, by adjusting the proportion coefficients of the damping factors according to

a specific model, M-PML can be stably performed. Comparing the wavefield generated in two homogenous media shows that PML, frequently used in the first-order wave equation, has the disadvantage of removing artificial reflections. At the same time, the M-PML performs well due to the damping profiles specified in multiple directions. Numerical experiments conducted in a strongly anisotropic medium show that M-PML remains stable in the frequency and time domains for viscoelastic anisotropic wavefield propagation. In addition, the M-PML requires only a modification of the damping factors. Therefore, the additional computational cost remains unchanged for techniques that rely on wavefield simulation, such as wave equation migration and waveform inversion.

### **2. Theoretical Foundations**
