*2.1. Frequency Domain Elastic VTI Equation*

The 2D frequency-domain anisotropic elastic Equation [10] is given as:

$$\begin{cases} \rho \omega^2 u\_x + \mathbb{C}\_{11} \frac{\partial^2 u\_x}{\partial x^2} + \mathbb{C}\_{44} \frac{\partial^2 u\_x}{\partial z^2} + (\mathbb{C}\_{13} + \mathbb{C}\_{44}) \frac{\partial^2 u\_x}{\partial x \partial z} + f(\omega) = 0, \\\\ \rho \omega^2 u\_z + \mathbb{C}\_{44} \frac{\partial^2 u\_z}{\partial x^2} + \mathbb{C}\_{33} \frac{\partial^2 u\_z}{\partial z^2} + (\mathbb{C}\_{13} + \mathbb{C}\_{44}) \frac{\partial^2 u\_x}{\partial x \partial z} + g(\omega) = 0, \end{cases} \tag{1}$$

where *ρ* is the density, *f*(*ω*) and *g*(*ω*) are Fourier components of horizontal and vertical body forces, respectively; *ux* and *uz* are the horizontal and vertical displacements, respectively; and *ω* is the angular frequency. The elastic stiffness is written as:

$$\begin{aligned} \text{C}\_{11} &= (1 + 2\varepsilon)\rho v\_{pr}^2\\ \text{C}\_{33} &= \rho v\_{pr}^2\\ \text{C}\_{13} &= \sqrt{2\delta \text{C}\_{33} (\text{C}\_{33} - \text{C}\_{44}) + (\text{C}\_{33} - \text{C}\_{44})^2} - \text{C}\_{44} \\\\ \text{C}\_{44} &= \text{C}\_{55} = \rho v\_{sr}^2 \end{aligned} \tag{2}$$

where *vpr* and *vsr* are the P-wave and the SV-wave velocity, respectively; *ε* and *δ* are Thomsen's anisotropic parameters [29].

The frequency-domain finite-difference method can easily introduce attenuation [30]. The constant-Q model is adopted here, and the complex velocity is given as:

$$\begin{split} \frac{1}{v\_{p}(\omega)} &= \frac{1}{v\_{pr}} \left( 1 - \frac{1}{\pi Q\_{p}} \ln \left| \frac{\omega}{\omega\_{r}} \right| \right) \cdot \left( 1 - \frac{\text{i}}{2Q\_{p}} \right) , \\\\ \frac{1}{v\_{s}(\omega)} &= \frac{1}{v\_{sr}} \left( 1 - \frac{1}{\pi Q\_{s}} \ln \left| \frac{\omega}{\omega\_{r}} \right| \right) \cdot \left( 1 - \frac{\text{i}}{2Q\_{s}} \right) , \end{split} \tag{3}$$

where *vp*, *vs* is the complex p-wave and SV-wave velocity, respectively; and *Qp*, *Qs* is the quality factor of the P-wave and SV-wave; *ω<sup>r</sup>* is the dominant frequency of the source wavelet; *i* represents the imaginary unit. Substituting Equation (3) into Equation (2) yields the viscoelastic anisotropic wave equation.

### *2.2. Absorbing Boundary Condition*

In order to introduce the M-PML absorbing condition to process unwanted reflections, the frequency-domain anisotropic elastic wave equation can be written as:

$$\begin{aligned} \rho\omega^{2}u\_{x} + \frac{1}{\varepsilon\_{x}} \Big[ \mathcal{C}\_{11} \frac{1}{\varepsilon\_{x}} \frac{\partial^{2}u\_{x}}{\partial x^{2}} + \mathcal{C}\_{11} \frac{\partial}{\partial x} \Big( \frac{1}{\varepsilon\_{x}} \right) \frac{\partial u\_{x}}{\partial x} + \mathcal{C}\_{13} \frac{1}{\varepsilon\_{x}} \frac{\partial^{2}u\_{x}}{\partial x \partial z} + \mathcal{C}\_{13} \frac{\partial}{\partial x} \Big( \frac{1}{\varepsilon\_{x}} \right) \frac{\partial u\_{x}}{\partial z} \Big] \\ + \frac{1}{\varepsilon\_{z}} \Big[ \mathcal{C}\_{44} \frac{1}{\varepsilon\_{x}} \frac{\partial^{2}u\_{x}}{\partial x \partial x} + \mathcal{C}\_{44} \frac{\partial}{\partial z} \Big( \frac{1}{\varepsilon\_{x}} \right) \frac{\partial u\_{x}}{\partial x} + \mathcal{C}\_{44} \frac{1}{\varepsilon\_{z}} \frac{\partial^{2}u\_{x}}{\partial z^{2}} + \mathcal{C}\_{44} \frac{\partial}{\partial z} \Big( \frac{1}{\varepsilon\_{z}} \right) \frac{\partial u\_{x}}{\partial z} \Big] + f(\omega) = 0, \\\\ \rho\omega^{2}u\_{z} + \frac{1}{\varepsilon\_{x}} \Big[ \mathcal{C}\_{44} \frac{1}{\varepsilon\_{x}} \frac{\partial^{2}u\_{x}}{\partial x^{2}} + \mathcal{C}\_{44} \frac{\partial}{\partial x} \Big( \frac{1}{\varepsilon\_{x}} \right) \frac{\partial u\_{x}}{\partial x} + \mathcal{C}\_{44} \frac{1}{\varepsilon\_{z}} \frac{\partial^{2}u\_{z}}{\partial z} + \mathcal{C}\_{44} \frac{\partial}{\partial$$

where *ex* = 1 − i *dx*(*x*) *<sup>ω</sup>* and *ez* = 1 − i *dz* (*z*) *<sup>ω</sup>* , *x* and *z* are respectively the distance from the inner boundary of the horizontal and vertical direction; *dx*(*x*) is the damping profiles, given as:

$$d\_x(\mathbf{x}) = 2\pi a\_0 f\_0 \left(\frac{\mathbf{x}}{L}\right),\tag{5}$$

where *x* is the distance from the inner boundary of the horizontal direction; *α*<sup>0</sup> is the optimized parameter, and here *α*<sup>0</sup> is given as 1.79 [30]; *f*<sup>0</sup> is the dominant frequency of the source wavelet and *L* is the thickness of the PML absorbing boundary condition. The new vertical coordinate has a similar form to the horizontal one:

$$d\_z(z) = 2\pi a\_0 f\_0\left(\frac{z}{L}\right),\tag{6}$$

where *z* is the distance from the inner boundary of the vertical direction.

Artificial boundary reflection can seriously affect the effect of numerical modeling. Therefore, a frequency-domain multi-axial perfectly matched layer (M-PML) is exploited to absorb the unwanted reflections.

In classical PML, unwanted waves are only attenuated in one direction (uniaxial). Within the left and right boundaries, only the damping function along the x-direction is nonzero and can be formulated as:

$$d\_x = d\_x^{\mathbf{x}}(\mathbf{x}), \ d\_z = 0. \tag{7}$$

Similarly, within the bottom and left boundaries, only the damping function along the z-direction takes effect and can be given as:

$$d\_x = 0, \; d\_z = d\_z^z(z). \tag{8}$$

To attenuate unwanted reflections efficiently, M-PML was developed [26]. The basic idea of M-PML is that the damping functions are proportional to each other. The damping function along the x-direction can be defined as:

$$d\_x = d\_x^x(\mathbf{x}), \ d\_z = p^{(z/x)} d\_x^x(\mathbf{x}), \tag{9}$$

where *p*(*z*/*x*) is the proportion coefficient in either the right or left PML boundary.

Similarly, the damping function along the z-direction can be written as:

$$d\_x = p^{(x/z)} d\_z^z(z), \; d\_z = d\_z^z(z), \tag{10}$$

where *p*(*x*/*z*) is the proportion coefficient in either the bottom or top PML boundary. During implementation, the proportions *p*(*x*/*z*) are between 0 and 1 [27].
