**3. Results**

#### *3.1. Di*ff*erent Q Media*

First, we investigate the accuracy of the solution of fractional time wave equation using different discrete algorithms in uniform media. The media is discretized on a grid of 151 × 151 points with the same spatial sampling interval 2 m in x- and z-directions and the time step is 0.05 ms. The phase velocity *c*0 is 2200 m/s at a reference frequency *f*0 = ω0/(<sup>2</sup>π) = 1500 Hz. The source is a Ricker wavelet with the central frequency of 100 Hz located at the center of the model. Figures 1–4 show the seismograms calculated by the original method and the new method with different values of *Q* = 5, 10, 30, 100 and different memory lengths of *L* = 9, 14, 24, 34, 44. When *Q* and *L* are both small, the deviation of the seismograms calculated by the original method is large. In order not to affect the display effect of other curves, these seismograms are not displayed. The results are also partially magnified for clear presentation of details. The solution calculated by the original method with a long memory length is excellently consistent with the analytical solution [12]. However, when the memory length is smaller, a degraded numerical solution will be obtained. We choose the solution calculated by the original method with all previous values as the reference solution. From Figures 1–4, when *L* is relatively small, it can be seen that the seismograms calculated by the original method has a large deviation from the reference curve, and the amplitude cannot return to 0 after the propagation of wavelet, which will seriously disturb the wavefield. In particular, when *Q* is relatively small, this error will be more remarkable. However, the new method is relatively stable when *L* is small and is closer to the reference curve. Figure 5 shows the comparison of the snapshots and seismograms calculated by different methods at *Q* = 10 and *L* = 34. Receivers are at the same depth as the source. The snapshot and seismogram calculated by the original method have obvious false disturbances, and pseudo hyperbolic fluctuations occur in the seismogram. The original method is inaccurate at a small memory length *L*. We access the accuracy of numerical solutions by the root-mean-square (rms) error, which is defined by

$$err\_{\text{rms}} = \sum\_{i=1}^{nt} \left(d\_i^p - d\_i^{\text{lr}}\right)^2 / \sum\_{i=1}^{nt} \left(d\_i^p\right)^2.$$

where *d p i* denotes the reference solution calculated by the original method with all previous values and *dh i* denotes the calculated value using the di fference method with limited memory length *L*. The relationship between the rms error and the number of memory length for di fferent values of *Q* is shown in Figure 6a. The error decreases with the increasing value of *L*, and the smaller value of *Q* results in a larger error. The new method has higher accuracy and stability at a small value of *L*. Table 1 shows the simulation time that is used to solve the fractional time derivative wave equation by di fferent methods. The new method is faster at the same memory length, whose calculation time is about 20% of the original method (round numbers of CPU time is used). Under the same error requirements, we can choose a smaller memory length using the new method, which can further save a lot of calculation time.

**Figure 1.** (**a**) Comparisons of seismograms at *Q* = 100 and 30 m away from the source. The first number '1' represents the original method, and the first number '2' represents the new method. The second number represents the memory length *L*. (**b**) Zoomed part of (**a**)**.**

**Figure 2.** (**a**) Comparisons of seismograms at *Q* = 30 and 30 m away from the source. The first number '1' represents the original method, and the first number '2' represents the new method. The second number represents the memory length *L*. (**b**) Zoomed part of (**a**)**.**

**Figure 3.** (**a**) Comparisons of seismograms at *Q* = 10 and 30 m away from the source. The first number '1' represents the original method, and the first number '2' represents the new method. The second number represents the memory length *L*. (**b**) Zoomed part of (**a**)**.**

**Figure 4.** Comparisons of seismograms at *Q* = 5 and 30 m away from the source. The first number '1' represents the original method, and the first number '2' represents the new method. The second number represents the memory length *L*.

**Figure 5.** *Cont.*

**Figure 5.** Calculation results of *Q* = 10, *L* = 34: (**a**) and (**b**) are the snapshots (0.05 s) and seismogram calculated by the original method, respectively; (**c**) and (**d**) are the snapshots (0.05 s) and seismogram calculated by the new method, respectively.

**Figure 6.** (**a**) The relationship between the error and the number of memory length for different values of *Q*: solid lines represent the original method; dotted lines represent the new method. (**b**) Four snapshot parts calculated by the new methods: (1) *Q* = 100; (2) *Q* = 30; (3) *Q* = 10; (4) *Q* = 5; Snapshots are recorded at 0.05 s.


**Table 1.** Simulation time with different memory lengths of *L*.

Figure 6b (1)–(4) shows the snapshots calculated by the new method at *Q* = 100, 30, 10, 5, respectively. As the value of *Q* decreases, we can see that the amplitude gradually decreases, and the phase gradually delays.
