*3.1. Comparative Analysis of the Stability of PML and M-PML*

In this section, a medium of strong anisotropy was employed to test the stability and performance of PML and M-PML. The material properties of the strongly anisotropic medium are given in Table 1: medium 3. We performed wavefield simulations to test the instability of the boundary conditions and confirmed that M-PML outperformed PML in this strongly anisotropic medium. A square model was considered embedded in an absorbing boundary of four sides. The computational model of medium 3, similar to the model used in Meza-Fajardo [25], has a size of 720 m × 720 m, and it is discretized by 480 × 480 grids with a grid size of 1.5 m. A boundary with a thickness of 120 grids is considered to observe the phenomenon of wave propagating inside truncated boundary conditions. A vertical source dominated by 15 Hz is located at the grid point (230 m, 230 m), which is close to the artificial edge. The real parts of the monochromatic wavefields using classical PML and M-PML, respectively, are given in Figure 3. In Figure 3, the monochromatic wavefield

at 15 Hz was calculated using classical PML and M-PML, respectively. The asterisk in Figures 3 and 4 is the vertical source. The black dashed line in the figure shows the onset of the absorbing zone. From the results, it can be seen that the PML appears unstable at the boundary intersection in the upper right corner. This phenomenon is due to the complexity of the wavefield in the absorption boundary region, which leads to numerical instability of the PML. Since the PML is a uniaxially defined damping profile, absorption incompleteness is also observed in the other three corners. In the M-PML, however, since the damping profile is multi-axial, the instability that appears in the upper right corner of the PML is well resolved, and the attenuation of the waves is better in the corner regions.


**Table 1.** Model parameters used for numerical experiments.

**Figure 3.** Real part of the horizontal component from the monochromatic wavefield at 15 Hz with (**a**) classical PML and (**b**) M-PML.

Snapshots of the time-domain wavefield corresponding to Figure 3 are given in Figure 4. The black dashed line in Figure 4 is the onset of the absorbing boundary condition, and the source is indicated by an asterisk. Figure 4a,c,e,g are the wavefields calculated by PML, and Figure 4b,d,f,h are the wavefields calculated by M-PML. The results show that the PML-calculated wavefields display a significant instability because at 0.0818 s, the wavefield has not yet propagated to the boundary, but strong energy reflections have already occurred, and these reflections result from the instability remain almost constant as the wavefield propagates. Therefore, PML is unstable in strongly anisotropic media. In contrast, the results of M-PML in the right panel have no unstable reflections, and the waves are successfully absorbed in the truncated boundary regions, so the M-PML with multi-axial damping profiles solves the problem of instability of PML in strongly anisotropic media.

**Figure 4.** Snapshots of the horizontal component simulated with classical PML (**a**), (**c**), (**e**), (**g**) (left panel) and M-PML (**b**), (**d**), (**f**), (**h**) (right panel) at 0.0818 s, 0.2636 s, 0.4455 s, and 1.3046 s.

### *3.2. Comparative Analysis in Homogeneous Elastic Anisotropic Media*

In the previous section, we verified that M-PML is superior to PML in terms of stability. However, the absorption performance advantages of M-PML and PML were not compared. Therefore, in this section, wavefield simulations are performed using M-PML and PML in different anisotropic media to compare the absorption of the two boundary conditions by observing the time-domain wavefield snapshots and comparing the energy curves. The energy of the wavefield snapshot is the sum of the squares of the amplitude values of each point in the wavefield. Then, the energy decay curve is obtained for each moment of the wavefield snapshot and plotted as a curve with the horizontal axis of time and the vertical axis of amplitude.

In the first test, the physical properties are given in Table 1: medium 1. We set the model size to be 1.0 km × 1.0 km with a grid space of 2.5 m, and employed an explosive Ricker wavelet source with a dominant frequency of 20 Hz. Figure 5 exhibits snapshots of horizontal displacement and its corresponding analytical group velocity calculated from medium 1. These snapshots were simulated using PML and M-PML, respectively, at different times. Figure 5a,b are the snapshots computed inside the PML and M-PML at 0.1667 s. The consistency between the wavefront calculated from the analytical group velocity and the wavefront of snapshots can be observed in Figure 5c. Similar to the figures above, Figure 5d–f are recorded at 0.2444 s. The wavefront of snapshots is almost the same as the wavefront calculated from the analytical group velocity, and the absorptive capacity of M-PML is superior to that of PML at four corners. Figure 5g,h are recorded at 0.7 s, and the absorptive advantage of M-PML over PML can be observed. Figure 5i shows the wavefield energy decay curves in the computational domain. These results show that the absorption performance of M-PML is better than that of PML when both ABCs remain stable, which can be verified by the comparative observations in Figure 5g,h. Because at 0.7 s, the effective wavefield has all propagated within the boundary, while the boundary reflections observed in PML are stronger at this time, and that of M-PML is much weaker. To further quantitatively compare the performance of the two ABCs, we gave the energy curves of them, where the horizontal axis is time, and the vertical axis is the logarithm of energy, and we can see that the energy of the waves of both ABCs becomes weaker as time passes, but the energy of M-PML decays faster than that of PML. In addition, in order to verify the absorptive capacity and reliability of M-PML in viscous media, the wavefield snapshots of M-PML in viscoelastic anisotropic media are given in Figure 6. From these results, we know that M-PML is still effective in viscous anisotropic media.

Regarding the verification of the results, we calculated the wavefront using the analytical group velocity in anisotropic media [31], and the results are displayed in Figure 5c,f, where the fast propagating wavefront is the P-wave, and the slow one is the S-wave. It can be seen that the wavefield snapshots at the corresponding moment also contain Pwaves and S-waves. The wavefront in snapshots can correspond well with the analytical wavefront in Figure 5c,f, which illustrates the reasonableness and applicability of the numerical results.

We further tested our algorithm with medium 2 in this test. The physical properties are given in Table 1. We also set the model size to be 1.0 km × 1.0 km with a grid space of 2.5 m and employed an explosive Ricker wavelet source with a dominant frequency of 20 Hz. Figure 7 shows snapshots of horizontal displacement and its analytical wavefront calculated from medium 2. These snapshots are calculated within the PML and M-PML boundary conditions, respectively. Figure 7a,b are the snapshots simulated inside the PML and M-PML at 0.1222 s, respectively; Figure 7d,e are simulated inside the PML and M-PML at 0.1889 s, respectively; and the wavefront calculated by the analytical group velocity in Figure 7c,f are consistent with the wavefront of snapshots at 0.1222 s and 0.1889 s, respectively. Figure 7g,h are recorded at 0.6 s, and significant reflections can be observed at the four corners in the snapshot of the PML; however, the M-PML does not. Figure 7i shows the energy decay curves of the wavefield calculated from medium 2. This numerical experiment was conducted to confirm that M-PML has superior absorption

performance to PML in arbitrarily different anisotropic media. It can still be seen from the results that the wavefront in the wavefield snapshots Figure 7a,b coincides with the analytical wavefront [31] Figure 7c,f at the corresponding moment, verifying the rationality and applicability of the numerical results. Similar to the previous section, viscosity is still introduced in this medium, and its wavefield snapshots are displayed in Figure 8 to confirm that M-PML can still maintain good absorptive properties and stability in different viscous anisotropic media.

**Figure 5.** Comparison between snapshots of the horizontal displacement in medium 1 at 0.1667 s, 0.2444 s, 0.7 s (**a**), (**d**), (**g**) calculated by M-PML and (**b**), (**e**), (**h**) PML. Wavefront calculated from analytical group velocity at (**c**) 0.1667 s and (**f**) 0.2444 s. (**i**) is the energy decay curve of PML and M-PML.

**Figure 6.** Viscoelastic anisotropic wavefield of medium 1 calculated by M-PML at (**a**) 0.1667 s, (**b**) 0.2444 s, and (**c**) 0.7 s.

**Figure 7.** Comparison between snapshots of the horizontal displacement in medium 2 at 0.1222 s, 0.1889 s, 0.6 s. (**a**), (**d**), (**g**) calculated by M-PML, and (**b**), (**e**), (**h**) PML. Wavefront calculated from analytical group velocity at (**c**) 0.1222 s and (**f**) 0.1889 s. (**i**) is the energy decay curve of PML and M-PML.

**Figure 8.** Viscoelastic anisotropic wavefield of medium 2 calculated by M-PML at (**a**) 0.1222 s, (**b**) 0.1889 s, and (**c**) 0.6 s.

Regarding absorption performance, observing the wavefield snapshots of PML and M-PML shows that both ABCs are stable in this homogenous anisotropic medium, so it is possible to compare their absorption performance visually. At 0.6 s, the effective wavefield wholly left the working area, and we can see that the boundary reflections of PML are strong at the four corners in Figure 7h, but those of M-PML are well attenuated. Further, in Figure 7i, it is illustrated that the energy decay curve of M-PML is lower than that of PML, which quantitatively displays the absorption advantage of M-PML over PML.

### *3.3. Wave Propagation in a Complex Anisotropic Medium*

Finally, we used the frequency domain method based on M-PML to perform wavefield simulation for a more complicated cross-well model. The cross-well model is shown in Figure 9. The corresponding parameters selected are as follows: computing area, 0 ≤ *x* ≤ 108 m, 0 ≤ *z* ≤ 289 m; spatial step, *h* = Δ*x* = Δ*z* = 0.5 m; grid parameter, 578 × 216; the number of M-PML layers, 30; dominant frequency of the source, 130 Hz, the location of the source: (135 m, 0 m); receiver alignment positions: *x* = 108 m, 0.5 ≤ *z* ≤ 288.5 m; receiver interval, 2 m. The cross-well model has eight layers. There are numerical constants in red in Figure 8, which label the layers in the figure. The red constants 1–8 in the figure correspond to layers 1–8 in Table 2, respectively. The asterisk in the figure represents the source. These anisotropic medium parameters are taken from the true sedimentary and are available in the study of Thomsen [28].

**Figure 9.** Multi-layer cross-well model with single-excited source well on the left and multiplereceived receiver well on the right, with red constants for medium numbers.


**Table 2.** Physical parameters used for the cross-well model.

Figure 10 shows the mono-frequency wavefield snapshots of elastic and viscoelastic anisotropy obtained by the proposed method at 163.2 Hz, 329.8 Hz, and 429.5 Hz. Furthermore, the time-domain snapshots at the time point of 0.0199 s, 0.0299 s, and 0.0399 s computed by the proposed method are shown in Figure 11, respectively. It is easy to see that there exist no artificial edge effects near the truncated boundary in either the elastic or viscoelastic anisotropic wavefield. Then, we give the single-shot seismograms of the crosswell model, which are shown in Figure 12, and we also see that there exist no numerical dispersion or artificial edge effects in either elastic or viscoelastic anisotropic seismograms. Thus, we may conclude that the M-PML for the finite-difference frequency domain is suitable for wavefield simulation in a viscoelastic anisotropic and complicated medium.

**Figure 10.** Mono-frequency horizontal elastic anisotropic wavefield snapshots of different frequencies for the cross-well model at (**a**) 163.2 Hz, (**b**) 329.8 Hz, (**c**) 429.5 Hz, and viscoelastic anisotropic wavefield at (**d**) 163.2 Hz, (**e**) 329.8 Hz, and (**f**) 429.5 Hz.

**Figure 11.** Time-domain horizontal elastic anisotropic wavefield snapshots for the cross-well model at (**a**) 0.02 s, (**b**) 0.03 s, (**c**) 0.04 s, and viscoelastic anisotropic wavefield at (**d**) 0.02 s, (**e**) 0.03 s, and (**f**) 0.04 s.

**Figure 12.** Time-domain horizontal (**a**) elastic anisotropic and (**b**) viscoelastic anisotropic seismogram corresponding to Figure 11.
