**1. Introduction**

Seismic exploration is a main technique in natural gas hydrate (NGH) survey. To gain a better estimation of the NGH content, we need to understand the characteristics of the seismic wave propagation in NGH. Traditional seismic modeling and inversion techniques are usually built on perfectly elastic medium model. However, the real underground medium usually has attenuation properties, which will cause amplitude loss and phase distortion of seismic waves. In particular, gas hydrate layer often shows abnormal velocity and attenuation characteristics due to hydrate filling. Ignoring the attenuation of the medium will make the numerical simulation and inversion results different from the real situation [1–3], which will result in the inability to obtain the true and accurate structural features of the underground. Studying the numerical simulation method of viscoacoustic seismic wave will help to obtain the actual propagation of underground seismic wave. Moreover, using the viscoacoustic wave equation for migration and inversion can effectively compensate the loss of amplitude, accelerate the convergence rate, and make the inversion results closer to the actual underground geological characteristics [4].

At present, many methods have been proposed for viscoacoustic seismic wave simulation [5–7]. Aki et al. [8] achieves viscoacoustic wave simulation by introducing complex velocity in the frequency domain, but the calculation is huge. Another method to construct viscoacoustic models is by combining mechanical elements [5,9–11], such as standard linear solid model (SLS) and Maxwell model, which are approximately constant-*Q* model and require a large amount of memory and computational time [5,7]. In addition to the above integer-order method, the fractional-order wave equation can better describe the amplitude loss and phase distortion of seismic waves. Carcione et al. [12] proposed a fractional

time derivative wave equation using the stress-strain relationship proposed by Kjartansson [13]. The equation can accurately describe the constant-*Q* behavior, but the fractional time derivative will introduce huge memory consumption [14,15]. In order to reduce memory consumption, Chen and Holm [16] proposed using fractional Laplacian to simulate irregular attenuation behavior. Carcione [17] uses fractional Laplacian to give a new wave equation, which effectively reduces memory and describes amplitude attenuation and velocity dispersion in a single term. Song et al. [18] proposed an asymptotic local finite difference based on the truncated difference stencil and applied the method to the fractional wave equation pointed out by Carcione [17]. Zhu and Harris [3] proposed a nearly constant-*Q* (NCQ) wave equation by approximating the fractional time derivative wave equation. The amplitude attenuation and velocity dispersion are described by two decoupled fractional Laplacians in the NCQ equation, which are helpful to compensate the attenuation loss in the inverse problem. Sun et al. [4,19] further numerically simulated the equation proposed by Zhu and Harris using low-rank one-step wave approximation [20,21]. Yao et al. [22] developed a local-spectral approach to implement the fractional Laplacian in the viscoacoustic wave equation. Wang et al. [23] improved the numerical discretization of the temporal derivatives in the NCQ equation. Xu et al. [24] adopted radial basis collocation method to solve the NCQ equation. The NCQ equation avoids solving the fractional time derivative and therefore saves the amount of calculation and memory, but it is actually an approximate modification of the fractional time derivative wave equation and has a certain loss of accuracy. In this paper we use the fractional time derivative wave equation to simulate the wave field. In recent years, many scholars have proposed different methods for solving fractional time derivatives [25,26], which have promoted its development.

In the previous researches, the realization of the fractional time derivative was approximated by the original Grunwald–Letnikov finite difference method. In this paper, a more unified definition of fractional derivatives is proposed to calculate the fractional time derivative. The rigorous form is the limit of the original finite difference when time step approaches zero. This method is more stable at different memory lengths and has higher accuracy when the memory length is short. We first give a review of the derivation of the fractional time derivative wave equation. Then, we introduce the more unified definition of fractional derivatives with finite difference method and apply this method to calculate the fractional time derivative. Finally, some numerical examples on NGH model illustrate the accuracy and stability of this new method.

#### **2. Materials and Methods**

#### *2.1. Definitions of Fractional Derivatives*

#### 2.1.1. Grunwald–Letnikov Fractional-Order Derivative

The definition of Grunwald–Letnikov (G-L) fractional-order derivative is an extension of the definition of the limit of integer-order derivative to real-order. Let α > 0, the Grunwald–Letnikov α-th order fractional derivative of the function *f*(*t*) with respect to *t* and the terminal value *a* is given by [12,27,28]

$$\, \, \_a^{GL} D\_i^a \, f(t) = \lim\_{h \to 0 \atop mh \to t \to a} \frac{1}{h^a} \sum\_{k=0}^m \alpha\_k^{(a)} \, f(t - kh) \, , \tag{1}$$

where ω(α) *k* = (−<sup>1</sup>)*<sup>k</sup>*! α*k* " = (−<sup>1</sup>)*<sup>k</sup>*Γ(α+<sup>1</sup>) <sup>Γ</sup>(*k*+<sup>1</sup>)Γ(<sup>α</sup>−*k*+<sup>1</sup>), <sup>Γ</sup>(·) denotes gamma function given below and α represents the order of the fractional operator. Binomial coefficients can be generalized to real arguments by means of the gamma function

$$\Gamma(\mathbf{x}) = \int\_0^\infty t^{\mathbf{x}-1} e^{-t} dt,\ \mathbf{x} > 0. \tag{2}$$

For which it is well-known that Γ(1) = 1 and <sup>Γ</sup>(*x* + 1) = *<sup>x</sup>*<sup>Γ</sup>(*x*), for any *x* > 0.

It is clear that if α = 1, *GL a <sup>D</sup>*<sup>α</sup>*t f*(*t*) = *f*(*t*)−*f*(*<sup>t</sup>*−*h*) *h* ; and if α = 2, *GL a <sup>D</sup>*<sup>α</sup>*t f*(*t*) = *f*(*t*)−<sup>2</sup> *f*(*<sup>t</sup>*−*h*)+ *f*(*<sup>t</sup>*−2*h*) *h*<sup>2</sup> . They correspond to the first-order difference quotient and second-order difference quotient, respectively.

#### 2.1.2. Riemann–Liouville Fractional-Order Derivative

The Riemann–Liouville (R-L) integral of non-integer order α > 0 for *f*(*t*) in suitable space (e.g., *Lp*(*<sup>a</sup>*, *b*)) is defined by

$$I\_{a^{+}}I\_{t}^{a}f(t) = \frac{1}{\Gamma(a)} \int\_{a}^{t} \frac{f(\tau)}{\left(t-\tau\right)^{1-a}}d\tau. \tag{3}$$

The Riemann–Liouville derivative with the lower integration limit *a* would be

$${}^{RL}\_{a+}D\_t^a f(t) = \frac{1}{\Gamma(1-a)} \frac{d}{dt} \int\_a^t \frac{f(\tau)}{\left(t-\tau\right)^a} d\tau,\tag{4}$$

which is called the Riemann–Liouville fractional derivative of order α.

2.1.3. Caputo's Fractional-Order Derivative

> The Caputo's derivative of fractional order α > 0 is defined by

$$\,\_{a+}^{C}D\_{t}^{a}f(t) = \frac{1}{\Gamma(m-\alpha)} \int\_{a}^{t} \frac{f^{(m)}(\tau)}{\left(t-\tau\right)^{\alpha+1-m}}d\tau,\tag{5}$$

where *m* = α is the smallest integer such that *m* > α, *f*(*m*) denotes the classical derivative of integer order *m*.

#### 2.1.4. The Riesz Fractional-Order Derivative

The kernel *k*(*t*) of the Riesz fractional-order derivative is given by

$$k(t) = \frac{|t|^{-(a+1)}}{2\Gamma(-\alpha)\cos(\beta)}, \alpha > -1,\tag{6}$$

where β = (π/2)<sup>α</sup>. Clearly, if α = 1, the kernel would be *k*(*t*) = −|*t*|−<sup>2</sup> π . With the above kernel, the α-th order Riesz fractional derivative of function *f*(*t*) is defined by

$$D^a f(t) = f(t) \* k(t) = \frac{1}{2\Gamma(a)\cos(\frac{\pi a}{2})} \int\_{-\infty}^{\infty} |t - \tau|^{-a-1} f(\tau) d\tau,\tag{7}$$

where the notation "∗" denotes the convolution operator. Equation (7) can be regarded as the two-sided fractional derivatives. With above definitions, suppose *t* ∈ [*a*, *b*], for the order of α > 0, the Riesz fractional derivative can be written as

$$\frac{\partial^a f(t)}{\partial |t|^\alpha} = -\mathbb{C}\_a(\_a D\_t^a f(t) + \_t D\_b^a f(t)),\tag{8}$$

where *C*α = 1 2 cos(πα/2), α - 1, *<sup>a</sup>D*<sup>α</sup>*t* and *<sup>t</sup>DtD*<sup>α</sup>*b* denotes the left- and right-side Riemann–Liouville derivatives given by (9) and (10), respectively:

$$\, \_aD\_t^a f(t) = \frac{1}{\Gamma(1-a)} \frac{d}{dt} \int\_a^t \frac{f(\tau)}{\left(t-\tau\right)^a} d\tau,\tag{9}$$

$$\,\_1D\_b^\alpha f(t) = -\frac{1}{\Gamma(1-\alpha)} \frac{d}{dt} \int\_t^b \frac{f(\tau)}{\left(\tau - t\right)^\alpha} d\tau. \tag{10}$$

In which, we assume that α ∈ (0, <sup>1</sup>). If α is greater than 1, e.g., α ∈ (1, 2], the definitions can be made correspondingly as

$$\_aD\_t^a f(t) = \frac{1}{\Gamma(2-\alpha)} \frac{d^2}{dt^2} \int\_a^t \frac{f(\tau)}{\left(t-\tau\right)^{a-1}} d\tau,\tag{11}$$

$$\_dD\_b^\alpha f(t) = \frac{1}{\Gamma(2-\alpha)} \frac{d^2}{dt^2} \int\_t^b \frac{f(\tau)}{\left(\tau - t\right)^{\alpha - 1}} d\tau. \tag{12}$$

#### 2.1.5. Relations between the above Fractional-Order Derivatives

For α ∈ (*m* − 1, *<sup>m</sup>*), *f*(*t*) ∈ *Cm*−<sup>1</sup>([*<sup>a</sup>*, *<sup>t</sup>*]), and *f*(*m*) is continuous on [*a*, *t*], there holds

$${}^{RL}\_{a}D\_{t}^{a}f(t) = {}^{GL}\_{a+}D\_{t}^{a}f(t) = \sum\_{k=0}^{m-1} \frac{f^{(k)}(a)(x-a)^{k-a}}{\Gamma(k+1-a)} + \frac{1}{\Gamma(m-a)} \int\_{a}^{t} \frac{f^{(m)}(\xi)(x-a)^{k-a}}{(x-\xi)^{a-m+1}}d\xi,\tag{13}$$

i.e., the definition of R-L is equal to that of G-L for smooth functions.

Clearly, we can see from (13) that the three definitions can be identical if *f*(*t*) is smooth enough and satisfies zero initial conditions of *f*(*k*)(*a*) = 0 (*k* = 0, 1, ··· , *m* − 1).

It also indicates from the three definitions that R-L is suitable for general functions, G-L is easy to be utilized for discretization, and Caputo's is often related with continuous processes depending on time.

With above preparations, we now turn to the gas hydrate model establishment and the related seismic wave modeling in fractional-order form.

#### *2.2. Establishing Velocity and Quality Factor Model of Hydrate Layer*

Gas hydrate layer often shows abnormal velocity and attenuation characteristics due to hydrate filling. Accurately establishing the velocity and quality factor model is very important for the simulation of the seismic wave field of gas hydrate. White [29] first proposed the concept of a mesoscopic scale theoretical model. The so-called mesoscopic scale is an intermediate scale that is much larger than rock particles but much smaller than the wavelength. The object of White's theoretical research is the non-uniform infiltration phenomenon of immiscible multiphase fluids infiltrating into the porous medium. Many scholars have found through analysis that the mesoscopic scale porous model can better explain the attenuation of the seismic frequency band [30]. Therefore, the White theory is chosen to model the velocity and quality factor of the hydrate layer.

The calculation steps for the velocity and quality factor of the hydrate formation are as follows. First, through the effective medium model [31], the elastic modulus of the solid phase and dry rock skeleton can be calculated from the elastic modulus of each mineral component of the rock, and then the velocity and quality factor of the saturated fluid rock can be calculated by the White theory [29,32]. Through the above method, the velocity and quality factor of stratums with different mineral components and hydrate saturation can be calculated. Therefore, seismic wave field simulations can be performed on different hydrate stratums, which is very helpful for studying the law of seismic wave propagation in hydrate stratums.

#### *2.3. Viscoacoustic Wave Equation*

To simulate seismic wave propagation of the gas hydrate stratums, in anelastic media, the stress σ(*t*) can be expressed as a convolution of the variation of the strain ε(*t*) and the relaxation function ϕ(*t*) in the following form:

$$
\sigma(t) = \varphi(t) \* \dot{\varepsilon}(t),
\tag{14}
$$

where the symbol "∗" refers to the convolution operator.

Kjartansson [13] gives a relaxation function, which exactly describes the constant *Q* characteristic and is widely used in many seismic applications. The relaxation function is written as [3,12]

$$\varphi(t) = \frac{M\_0}{\Gamma(1-2\gamma)} \left(\frac{t}{t\_0}\right)^{-2\gamma} H(t),\tag{15}$$

where γ = arctan(1/*Q*)/π is a dimensionless parameter and we can know that 0 <γ< 1/2 by the range of *Q*(*Q* > 0); *M*0 = ρ*c*<sup>2</sup> 0 cos<sup>2</sup>(πγ/2) is a bulk modulus, Γ is the Euler's Gamma function, *H*(*t*) is the Heaviside step function and *t*0 is a reference time, e.g., *t*0 = 1/<sup>ω</sup>0, and *c*0 is the phase velocity at reference frequency ω0.

Combining the first-order conservation equation ∂ ∂*t v* = 1 ρ0 ∇σ (*v* is the phase velocity, ∇ denotes the gradient operator) with the strain-velocity equation ∂ ∂*t* ε = ∇ · *v* (here ∇· denotes the divergence operator), along with Equations (14) and (15), the fractional time derivative wave equation with uniform density ρ can be derived as

$$\frac{\partial^{2-2\gamma}\sigma}{\partial t^{2-2\gamma}} = \varepsilon\_0 {}^2 \omega\_0^{-2\gamma} \cos^2(\pi \gamma/2) \nabla^2 \sigma. \tag{16}$$

The above Equation (16) was rigorously derived from Kjartansson's constant *Q* model [13], hence it could accurately describe the attenuation behavior in realistic media. The Equation (16) reduces to the classical acoustic equation when γ → 0 (corresponding to *Q* → ∞). When γ → 1 2 (corresponding to *Q* → 0), this equation describes the propagation of the seismic wave in infinite attenuation medium. When 0 <γ< 1 2 , this equation describes the propagation of the seismic wave in attenuating medium. Moreover, in discrete calculation, a short-term memory principle shown in [33] is applied. The fractional time wave equation describes the loss mechanisms only by two parameters, i.e., the phase velocity *c*0 and the quality *Q*. Therefore, the method based on the above equation is more accurate and simpler than the other nearly constant *Q* methods, such as the methods based on SLS model or complex velocity [3,17].

In Zhu and Harris [3], the authors approximated the time fractional wave equation, and expressed the amplitude attenuation and velocity dispersion with two independent fractional Laplacians to obtain an approximately constant *Q* wave equation. The separated amplitude attenuation and frequency dispersion term are beneficial to compensate for attenuation loss in the inverse problem. The fractional Laplacian form of the wave equation reads as

$$\begin{cases} \frac{\partial \boldsymbol{v}(\boldsymbol{r},t)}{\partial t} = \frac{1}{\rho\_0(\boldsymbol{r})} \nabla \boldsymbol{\sigma}(\boldsymbol{r},t) + \mathcal{S},\\ \frac{\partial}{\partial t} \boldsymbol{\varepsilon}(\boldsymbol{r},t) = \boldsymbol{\nabla} \cdot \boldsymbol{v}(\boldsymbol{r},t),\\ \boldsymbol{\sigma}(\boldsymbol{r},t) = M\_0(\boldsymbol{r}) \Big[ \boldsymbol{\eta}(\boldsymbol{r}) \big( -\nabla^2 \boldsymbol{\rangle} \big)^{\boldsymbol{\gamma}(\boldsymbol{r})} + \boldsymbol{\tau}(\boldsymbol{r}) \frac{\partial}{\partial t} \big( -\nabla^2 \big)^{\boldsymbol{\gamma}(\boldsymbol{r}) - 1/2} \Big] \boldsymbol{\varepsilon}(\boldsymbol{r},t), \end{cases} \tag{17}$$

where *S* denotes the source, <sup>σ</sup>(*<sup>r</sup>*, *t*) and <sup>ε</sup>(*<sup>r</sup>*, *t*) are the stress and strain fields at position *r* at time *t*, η(*r*) and τ(*r*) are two coe fficients vary in space and are in the form of η(*r*) = −*c* <sup>2</sup>γ(*r*) 0 (*r*)ω<sup>−</sup>2γ(*r*) 0 cos(πγ(*r*)) and τ(*r*) = −*c* <sup>2</sup>γ(*r*)−<sup>1</sup> 0 (*r*)ω<sup>−</sup>2γ(*r*) 0 cos(πγ(*r*)), respectively, and all other parameters are defined the same as above. In calculation, the spatial fractional Laplace operator adopts fast Fourier transform method. Separating the amplitude attenuation term and dispersion term are beneficial to compensate attenuation loss in the inverse problem. However, γ(*r*) changes with space and is taken as the average of the entire space in the Fourier transform, which does not conform to the actual situation.
