## 2 Theoretical Foundations
### 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:\[\rho\omega^{2}u_{x}+\frac{1}{e_{x}}\Big{[}C_{11}\frac{1}{e_{x}}\frac{ \partial^{2}u_{x}}{\partial x^{2}}+C_{11}\frac{\partial}{\partial x}\Big{(} \frac{1}{e_{x}}\Big{)}\frac{\partial u_{x}}{\partial x}+C_{13}\frac{\partial^{2 }u_{x}}{\partial x\partial x}+C_{13}\frac{\partial}{\partial x}\Big{(}\frac{1}{ e_{z}}\Big{)}\frac{\partial u_{z}}{\partial z}\Big{]} \tag{4}\] \[+\frac{1}{e_{z}}\Big{[}C_{44}\frac{1}{e_{x}}\frac{\partial^{2}u_{ x}}{\partial x\partial x}+C_{44}\frac{\partial}{\partial z}\Big{(}\frac{1}{e_{x}} \Big{)}\frac{\partial u_{z}}{\partial x}+C_{44}\frac{1}{e_{z}}\frac{\partial^ {2}u_{x}}{\partial z^{2}}+C_{44}\frac{\partial}{\partial z}\Big{(}\frac{1}{e_ {z}}\Big{)}\frac{\partial u_{x}}{\partial z}\Big{]}+f(\omega)=0,\] \[\rho\omega^{2}u_{z}+\frac{1}{e_{x}}\Big{[}C_{44}\frac{1}{e_{x}} \frac{\partial^{2}u_{x}}{\partial x^{2}}+C_{44}\frac{\partial}{\partial x} \Big{(}\frac{1}{e_{x}}\Big{)}\frac{\partial u_{z}}{\partial x}+C_{44}\frac{ \partial^{2}u_{x}}{\partial x\partial z}+C_{44}\frac{\partial}{\partial x} \Big{(}\frac{1}{e_{z}}\Big{)}\frac{\partial u_{z}}{\partial z}\Big{]}\] \[+\frac{1}{e_{z}}\Big{[}C_{13}\frac{1}{e_{x}}\frac{\partial^{2}u_{ x}}{\partial z\partial x}+C_{13}\frac{\partial}{\partial z}\Big{(}\frac{1}{e_{x}} \Big{)}\frac{\partial u_{x}}{\partial x}+C_{33}\frac{1}{e_{z}}\frac{\partial^{ 2}u_{x}}{\partial z^{2}}+C_{33}\frac{\partial}{\partial z}\Big{(}\frac{1}{e_ {z}}\Big{)}\frac{\partial u_{z}}{\partial z}\Big{]}+g(\omega)=0,\]  
where \(e_{x}=1-i\frac{d_{x}(x)}{\omega}\) and \(e_{z}=1-i\frac{d_{z}(z)}{\omega}\), \(x\) and \(z\) are respectively the distance from the inner boundary of the horizontal and vertical direction; \(d_{x}(x)\) is the damping profiles, given as:  
\[d_{x}(x)=2\pi\alpha_{0}f_{0}\Big{(}\frac{x}{L}\Big{)}, \tag{5}\]  
where \(x\) is the distance from the inner boundary of the horizontal direction; \(\alpha_{0}\) is the optimized parameter, and here \(\alpha_{0}\) is given as 1.79 [30]; \(f_{0}\) 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\alpha_{0}f_{0}\Big{(}\frac{z}{L}\Big{)}, \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}^{x}(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}(x),\ d_{z}=p^{(z/x)}d_{x}^{x}(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].