*2.4. Methods*

#### 2.4.1. Fractional-Order Derivatives Approximation

The critical matter of modeling wave propagation by using the Equation (16) is the method for computing the fractional time derivative. There are different techniques for approximating different fractional derivatives. The spatial derivatives ∇2σ can be calculated by the fast Fourier transform or the finite difference. Carcione et al. [12] computed the fractional time derivative by the original Grunwald–Letnikov finite difference method with second-order accuracy, which is in the following form:

$$\frac{\partial^{\beta}f(t)}{\partial t^{\beta}} = f\_{\hbar}^{(\beta)}(t) \approx \frac{1}{h^{\beta}} \sum\_{j=0}^{l} (-1)^{j} \binom{\beta}{j} f\left(t + \left(\frac{\beta}{2} - j\right)\hbar\right) \tag{18}$$

where β is the order of fractional derivative with β = 2 − 2γ and hence 1 <β< 2, *h* is the time sampling step and *J* = *t*/*h* − 1, ! β*j* " is the fractional binomial coefficients, which can be expressed in terms of Euler's Gamma function as

$$
\Gamma\left(\begin{array}{c}\beta\\j\end{array}\right) = \frac{\beta(\beta-1)(\beta-2)\cdots(\beta-j+1)}{j!} = \frac{\Gamma(\beta+1)}{j!\Gamma(\beta-j+1)} = \frac{\Gamma(\beta+1)}{\Gamma(j+1)\Gamma(\beta-j+1)}.\tag{19}
$$

Calculation of the fractional derivative using Equation (18) requires storing all the previous values of *f*. Therefore, the original method will consume a lot of memory and computation, even though we can truncate the memory length in some situations.

If the function *f*(*t*) has *m* + 1-order continuous derivative on interval [*a*, *t*] and the integer *m* satisfies *m* <β< *m* + 1, the limit (as *h* → 0) of the Equation (18) can be obtained using the R-L fractional derivative:

$$\begin{split} \, \_aD\_t^\beta f(t) &= \lim\_{h \to 0} f\_h^{(\beta)}(t) \\ = \sum\_{k=0}^m \frac{f^{(k)}(a)(t-a)^{-\beta+k}}{\Gamma(-\beta+k+1)} + \frac{1}{\Gamma(-\beta+m+1)} \int\_a^t (t-\tau)^{m-\beta} f^{(m+1)}(\tau) d\tau \end{split} \tag{20}$$

In this study, we focus on the finite difference discretization of the fractional time derivative ∂βσ ∂*t*β for 1 <β< 2. Then the fractional derivative can be defined as following a more unified form:

$$\begin{split} \frac{\partial^{\beta}f(t)}{\partial t^{\beta}} & \triangleq \lim\_{h \to 0} f\_{h}^{(\beta)}(t) \\ &= c \mathbf{1} \frac{(t-a)^{-\beta}f(t\_{0})}{\Gamma(1-\beta)} + c \mathbf{2} \frac{(t-a)^{1-\beta}f'(t\_{0})}{\Gamma(2-\beta)} + \frac{1}{\Gamma(2-\beta)} \int\_{a}^{t} (t-\tau)^{1-\beta} f^{(2)}(\tau) d\tau, \end{split} \tag{21}$$

where the symbol "-"denotes definition, *c*1 and *c*2 are constants of 0 or 1. When *c*1 = *c*2 = 1, (21) is equivalent to R-L fractional derivative and G-L fractional derivative. When *c*1 = *c*2 = 0, (21) is equivalent to the Caputo fractional derivative.

#### 2.4.2. Finite Difference Method for Integer-Order Derivatives

Define the transfer operator *<sup>T</sup>h*, and difference operators ∇*h*, Δ*h* and δ*h* (respectively refers to the backward difference, forward difference and central difference), where *h* ∈ R is the time sampling stepsize. The operator *T<sup>h</sup>* acts on the function *f*(*t*), *t* ∈ R, gives

$$\begin{aligned} T^h f(t) &= f(t+h), \\ \nabla\_h f(t) &= f(t) - f(t-h), \\ \Delta\_h f(t) &= f(t+h) - f(t), \\ \delta\_h f(t) &= f(t + \frac{1}{2}h) - f(t - \frac{1}{2}h). \end{aligned} \tag{22}$$

Assume that *f*(*t*) is sufficiently smooth, we have that the *m*-order derivative of *f*(*t*) can be approximated by *f*(*m*)(*t*) = *Dm f*(*t*), *m* ∈ N, satisfying

$$D^m f(t) = h^{-m} \nabla\_h^m f(t) + O(h) = h^{-m} (I - T^{-h})^m f(t) + O(h)$$

or

$$D^m f(t) = h^{-m} \delta\_h^m f(t) + O(h^2) = h^{-m} (T^{h/2} - T^{-h/2})^m f(t) + O(h^2)$$

for backward and central difference, respectively; in which, *I* is the identity. The power *m* of the two difference operators ∇*h* and δ*h* can be expanded as

$$\nabla\_{\mathbf{h}}^{\mathrm{nr}} = \sum\_{j=0}^{m} (-1)^{j} \binom{m}{j} T^{-j\mathbf{h}}$$

and

$$\delta\_h^m = \sum\_{j=0}^m (-1)^j \binom{m}{j} T^{(m/2-j)h}$$

respectively. The above integer-order derivative forms can be extended to fractional cases, i.e., for any β ∈ R, the fractional-order operator can be related by

$$\begin{aligned} \nabla\_h^\beta &= \sum\_{j=0}^\infty (-1)^j \binom{\beta}{j} T^{-jh} \, \mathrm{d}h \\\ \delta\_h^\beta &= \sum\_{j=0}^\infty (-1)^j \binom{\beta}{j} T^{(\beta/2-j)h} \, \mathrm{d}h \end{aligned} \tag{23}$$

#### 2.4.3. Fractional Differencing Scheme

With the above preparation, we present some details about how to calculate the fractional-order time derivatives. Equation (20) is in the limit case where *h* approaches 0, which is more rigorous and accurate than Equation (18). Jia and Li [34] using the rigorous form to compute the spatial fractional derivative in the fractional advection-dispersion equation. Here, we consider the time fractional derivative discretization. We discretize time domain by (*i* = 0, 1, 2, ··· , *L*) in interval [*<sup>a</sup>*,*<sup>t</sup>*] (*<sup>t</sup>*0 = *a*). The positive integer *L* (*L* = (*ti* − *<sup>t</sup>*0)/*h*) is the memory length. Since 1 <β< 2, so *m* = 1. Using the discretization scheme (23) to calculate the fractional time derivative, the Equation (21) can be discretized as follows:

$$\begin{split} \frac{\partial^{\beta}f(t)}{\partial t^{\beta}}|\_{t=t\_{L}} &= c1\frac{(t\_{L}-t\_{0})^{-\beta}f(t\_{0})}{\Gamma(1-\beta)} + c2\frac{(t\_{L}-t\_{0})^{1-\beta}f'(t\_{0})}{\Gamma(2-\beta)} + \frac{1}{\Gamma(2-\beta)}\int\_{t\_{0}}^{t\_{L}} (t\_{L}-\tau)^{1-\beta}f^{(2)}(\tau)d\tau \\ &= c1\frac{(t\_{L}-t\_{0})^{-\beta}f(t\_{0})}{\Gamma(1-\beta)} + c2\frac{(t\_{L}-t\_{0})^{1-\beta}(f(t\_{1})-f(t\_{0}))}{\Gamma(2-\beta)\text{h}} + \\ &\frac{h^{-\beta}}{\Gamma(3-\beta)}\sum\_{j=0}^{L-1} \left[ (j+1)^{2-\beta} - (j)^{2-\beta} \right] \left( f(t\_{L-j+1}) - 2f(t\_{L-j}) + f(t\_{L-j-1}) \right) + R^{L}\_{h'} \end{split} \tag{24}$$

where *RLh* is the truncation error of magnitude *<sup>O</sup>*(*h*). It is obvious that the fractional derivative at time *t* depends on the past value between [*<sup>a</sup>*,*<sup>t</sup>*]. Recalling that we need to numerically solve the Equation (16). Using above fractional difference scheme, an explicit difference scheme for the stress σ(*t*) can be obtained as follows:

$$\begin{array}{lcl} \sigma(t\_{L+1}) = \frac{\Gamma(3-\beta)}{h^{-\beta}} c\_0^2 w\_0^{-2\gamma} \cos^2(\pi \chi/2) \nabla^2 \sigma(t\_L) + 2\sigma(t\_L) - \sigma(t\_{L-1}) \\ - \sum\_{j=1}^{L-1} \left[ (j+1)^{2-\beta} - (j)^{2-\beta} \right] (\sigma(t\_{L-j+1}) - 2\sigma(t\_{L-j}) + \sigma(t\_{L-j-1})) \\ - c \mathbf{1} \frac{(2-\beta)(1-\beta)(t\_L-t\_0)^{-\beta} \sigma(t\_0)}{h^{-\beta}} - c \mathbf{2} \frac{(2-\beta)(t\_L-t\_0)^{1-\beta} \left[ \sigma(t\_L) - \sigma(t\_0) \right]}{h^{1-\beta}} + \mathcal{S}(t\_L). \end{array} \tag{25}$$

In which, *S*(*tL*) denotes the seismic source at time *tL*, which can be chosen by users, e.g., the Ricker wavelet. Equation (25) does not need to compute the term ! β*j* ", and hence the speed of calculation is accelerated.SpatialderivativesarecalculatedbythefastFouriertransformwithstaggeredgrid.

In the following, we will call the simulation method based on the Equation (18), i.e., Grunwald–Letnikov finite difference method, as the "original method", the simulation method based on the Equation (24) as the "new method", where we choose *c*1 = *c*2 = 0.
