Fractional Time Derivative Seismic Wave Equation Modeling for Natural Gas Hydrate

: Simulation of the seismic wave propagation in natural gas hydrate (NGH) is of great importance. To ﬁnely portray the propagation of seismic wave in NGH, attenuation properties of the earth’s medium which causes reduced amplitude and dispersion need to be considered. The traditional viscoacoustic wave equations described by integer-order derivatives can only nearly describe the seismic attenuation. Di ﬀ erently, the fractional time derivative seismic wave-equation, which was rigorously derived from the Kjartansson’s constant- Q model, could be used to accurately describe the attenuation behavior in realistic media. We propose a new fractional ﬁnite-di ﬀ erence method, which is more accurate and faster with the short memory length. Numerical experiments are performed to show the feasibility of the proposed simulation scheme for NGH, which will be useful for next stage of seismic imaging of NGH.


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][2][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][6][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][10][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.

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] Energies 2020, 13 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.

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 () ft with respect to t and the terminal value a is given by [12,27,28] where () , ()  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 For which it is well-known that (1) 1  and ( 1) where ω , Γ(·) 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 Energies 2020, 13, 5901

of 24
For which it is well-known that Γ(1) = 1 and Γ(x + 1) = xΓ(x), for any . They correspond to the first-order difference quotient and second-order difference quotient, respectively.

Riemann-Liouville Fractional-Order Derivative
The Riemann-Liouville (R-L) integral of non-integer order α > 0 for f (t) in suitable space (e.g., L p (a, b)) is defined by The Riemann-Liouville derivative with the lower integration limit a would be which is called the Riemann-Liouville fractional derivative of order α.

Caputo's Fractional-Order Derivative
The Caputo's derivative of fractional order α > 0 is defined by where m = α is the smallest integer such that m > α, f (m) denotes the classical derivative of integer order m.

The Riesz Fractional-Order Derivative
The kernel k(t) of the Riesz fractional-order derivative is given by where β = (π/2)α. Clearly, if α = 1, the kernel would be k(t) = − |t| −2 π . With the above kernel, the α-th order Riesz fractional derivative of function f (t) is defined by 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 where C α = 1 2 cos(πα/2) , α 1, a D α t and t D t D α b denotes the left-and right-side Riemann-Liouville derivatives given by (9) and (10), respectively: Energies 2020, 13, 5901 4 of 24 In which, we assume that α ∈ (0, 1). If α is greater than 1, e.g., α ∈ (1, 2], the definitions can be made correspondingly as 2.1.5. Relations between the above Fractional-Order Derivatives 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.

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.

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: where the symbol " * " refers to the convolution operator.
Energies 2020, 13, 5901 5 of 24 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] 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 2 0 cos 2 (πγ/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/ω 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 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 ∂v(r,t) ∂t where S denotes the source, σ(r, t) and ε(r, t) are the stress and strain fields at position r at time t, η(r) and τ(r) are two coefficients vary in space and are in the form of η(r) = −c 2γ(r) 0 (r)ω −2γ(r) 0 cos(πγ(r)) and 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.

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: 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 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: (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: where the symbol " "denotes definition, c1 and c2 are constants of 0 or 1. When c1 = c2 = 1, (21) is equivalent to R-L fractional derivative and G-L fractional derivative. When c1 = c2 = 0, (21) is equivalent to the Caputo fractional derivative.

Finite Difference Method for Integer-Order Derivatives
Define the transfer operator T 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 h acts on the function f (t), t ∈ R, gives Energies 2020, 13, 5901 Assume that f (t) is sufficiently smooth, we have that the m-order derivative of f (t) can be approximated by 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 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

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 [a, t] (t 0 = a). The positive integer L (L = (t i − t 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: where R L h is the truncation error of magnitude O(h). It is obvious that the fractional derivative at time t depends on the past value between [a, t]. 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: Energies 2020, 13, 5901 8 of 24 In which, S(t L ) denotes the seismic source at time t L , 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. Spatial derivatives are calculated by the fast Fourier transform with staggered grid. 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 c1 = c2 = 0.

Proposition
In the above fractional differencing scheme, the memory length L is variable and is a fixed number in each differencing calculation. This is particularly important for efficient simulation. Otherwise, the memory length t − a will be a function of the variable t. According to Equations (18) and (24), when L is large, the calculation time with difference scheme (24) will be less than that of difference scheme of (18). However, as the value of L increases, its influence on the derivative will decrease. Therefore, proper choice of the memory length L will balance the accuracy of calculation and the time cost. In addition, the fractional differencing scheme (24) is stable, since the residual R i h related to discrete step h is bounded by the infinitesimal equivalent value of h 2 , and the fact that for any 1 < α < 2, The above expression is bounded. In next section, we will show how these values of memory length influence the results.

Different 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 /(2π) = 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 Energies 2020, 13, 5901 9 of 24 length L. We access the accuracy of numerical solutions by the root-mean-square (rms) error, which is defined by where d p i denotes the reference solution calculated by the original method with all previous values and d h i denotes the calculated value using the difference method with limited memory length L. The relationship between the rms error and the number of memory length for different 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 different 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.
length L . The relationship between the rms error and the number of memory length for different 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 different 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.  length L . The relationship between the rms error and the number of memory length for different 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 different 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.

Layered Model
To

Layered Model
To    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.

Layered Model
To verify the accuracy and stability of the new method in the media with large contrast in velocity and Q, we construct a two-layers model. The velocity and the Q of the top layer are 1200 m/s and 30, respectively. In the bottom layer, the velocity and the Q are 2200 m/s and 300, respectively. The interface is at a depth of 120 m. The media is discretized on a grid of 151 × 151 points with the same spatial sampling interval 2 m in x-and z-directions. The time step is 0.05 ms. A Ricker wavelet with the central frequency of 100 Hz is located at the center of the model (150, 150 m). Receivers are at the same depth as the source. Figure 7 shows the traces at x = 150 m extracted from seismograms calculated by the original method and the new method with different memory length L = 9, 14, 24. Similarly, we chose the solution calculated by the original method with all previous values as the reference solution. Using the original method with the memory length L = 14 leads to a significant deviation between the numerical solution and the reference curve. However, the solution calculated by the new method with a smaller memory length L = 9 is excellently consistent with the reference curve. Figure 8a,b show the snapshot and seismogram calculated by the original method with memory length L = 14. Figure 9a,b show the snapshot and seismogram calculated by the new method with memory length L = 9. Figure 10a

Simulations on Velocity and Quality Factor Model of Hydrate Layer
The elastic modulus of each mineral component and fluid we used is shown in Table 2. We assume that the rock solid phase is composed of 60% clay and 40% quartz, and the pores can contain different proportions of water, natural gas hydrate, and methane gas. The curve of velocity and inverse quality factor with hydrate saturation is calculated by white theory as shown in Figure  11. The seismic frequency used in the calculation is 35 Hz.
The calculation result of White theory shows that longitudinal wave velocity increases with the increase of hydrate saturation. When the hydrate saturation is small, the longitudinal wave attenuation increases rapidly to the maximum with the increase of saturation, and then the longitudinal wave attenuation decreases rapidly. Priest et al. [35] used hydrate resonance column (Gas hydrate resonant column) device to measure the attenuation value of 13 sandstone samples containing hydrate. The measurement results show that in the range of hydrate saturation less than 3-5%, the longitudinal wave attenuation increases rapidly to the maximum value with the increase of saturation, and then the longitudinal wave attenuation decreases rapidly, and then shows a slow increase trend. Therefore, it shows that there is a certain consistency between the White theoretical model and the experimental measurement results. In all subsequent seismic wave simulations of hydrate stratum, White theory is used to establish the velocity and quality factor models of hydrate stratum.

Simulations on Velocity and Quality Factor Model of Hydrate Layer
The elastic modulus of each mineral component and fluid we used is shown in Table 2. We assume that the rock solid phase is composed of 60% clay and 40% quartz, and the pores can contain different proportions of water, natural gas hydrate, and methane gas. The curve of velocity and inverse quality factor with hydrate saturation is calculated by white theory as shown in Figure  11. The seismic frequency used in the calculation is 35 Hz.
The calculation result of White theory shows that longitudinal wave velocity increases with the increase of hydrate saturation. When the hydrate saturation is small, the longitudinal wave attenuation increases rapidly to the maximum with the increase of saturation, and then the longitudinal wave attenuation decreases rapidly. Priest et al. [35] used hydrate resonance column (Gas hydrate resonant column) device to measure the attenuation value of 13 sandstone samples containing hydrate. The measurement results show that in the range of hydrate saturation less than 3-5%, the longitudinal wave attenuation increases rapidly to the maximum value with the increase of saturation, and then the longitudinal wave attenuation decreases rapidly, and then shows a slow increase trend. Therefore, it shows that there is a certain consistency between the White theoretical model and the experimental measurement results. In all subsequent seismic wave simulations of hydrate stratum, White theory is used to establish the velocity and quality factor models of hydrate stratum.

Simulations on Velocity and Quality Factor Model of Hydrate Layer
The elastic modulus of each mineral component and fluid we used is shown in Table 2. We assume that the rock solid phase is composed of 60% clay and 40% quartz, and the pores can contain different proportions of water, natural gas hydrate, and methane gas. The curve of velocity and inverse quality factor with hydrate saturation is calculated by white theory as shown in Figure 11. The seismic frequency used in the calculation is 35 Hz.  (a) (b) Figure 11. The curve of velocity (a) and inverse quality factor (b) with hydrate saturation.

Layered NGH Model
Now we consider seismic simulations of NGH model using fractional-order differencing of the constant-Q viscoacoustic wave equation. The amplitude of seismic wave propagation in hydrate layer always is very weak, which is called the amplitude blank zone. In order to simulate this phenomenon, we assume that the saturation of hydrate varies with porosity, and the ratio of hydrate saturation to porosity is constant. We also assume that when the depth increases, the porosity of the rock generally decreases as the pressure increases, and the hydrate saturation also decreases as the porosity decreases.
The first layered NGH model is a high hydrate saturation model and the pores of the underlying layer are filled with free gas. The model parameters are given in Table 3. The velocity model and the quality factor Q are shown in Figure 12a,b, respectively. From which we can see that there are high-speed and high-Q anomalies in the hydrate layer. The seismic simulations with fractional-order differencing is shown in Figures 13 and 14. It can be seen that bottom-simulating reflector (BSR) features are obvious, and its polarity is reversed, and the amplitude is significantly The calculation result of White theory shows that longitudinal wave velocity increases with the increase of hydrate saturation. When the hydrate saturation is small, the longitudinal wave attenuation increases rapidly to the maximum with the increase of saturation, and then the longitudinal wave attenuation decreases rapidly. Priest et al. [35] used hydrate resonance column (Gas hydrate resonant column) device to measure the attenuation value of 13 sandstone samples containing hydrate. The measurement results show that in the range of hydrate saturation less than 3-5%, the longitudinal wave attenuation increases rapidly to the maximum value with the increase of saturation, and then the longitudinal wave attenuation decreases rapidly, and then shows a slow increase trend. Therefore, it shows that there is a certain consistency between the White theoretical model and the experimental measurement results. In all subsequent seismic wave simulations of hydrate stratum, White theory is used to establish the velocity and quality factor models of hydrate stratum.

Layered NGH Model
Now we consider seismic simulations of NGH model using fractional-order differencing of the constant-Q viscoacoustic wave equation. The amplitude of seismic wave propagation in hydrate layer always is very weak, which is called the amplitude blank zone. In order to simulate this phenomenon, we assume that the saturation of hydrate varies with porosity, and the ratio of hydrate saturation to porosity is constant. We also assume that when the depth increases, the porosity of the rock generally decreases as the pressure increases, and the hydrate saturation also decreases as the porosity decreases.
The first layered NGH model is a high hydrate saturation model and the pores of the underlying layer are filled with free gas. The model parameters are given in Table 3. The velocity model and the quality factor Q are shown in Figure 12a,b, respectively. From which we can see that there are high-speed and high-Q anomalies in the hydrate layer. The seismic simulations with fractional-order differencing is shown in Figures 13 and 14. It can be seen that bottom-simulating reflector (BSR) features are obvious, and its polarity is reversed, and the amplitude is significantly increased. Above the BSR, an amplitude blank zone can be seen in the hydrate layer. The reflection of the upper layer of gas hydrate is relatively weak. increased. Above the BSR, an amplitude blank zone can be seen in the hydrate layer. The reflection of the upper layer of gas hydrate is relatively weak. Table 3. High hydrate saturation model with gas bearing layer.   The second layered NGH model is a low hydrate saturation model and the pores of the underlying layer are also filled with free gas. The model parameters are given in Table 4. The velocity model and the quality factor Q are shown in Figure 15a,b, respectively. From which we can see that there are high-speed and low-Q anomalies in the hydrate layer. The seismic simulations with fractional-order differencing is shown in Figures 16 and 17. The amplitude of BSR is weaker than that of the first NGH model, but it is still obvious, and its polarity is reversed. There also exists an amplitude blank zone above the BSR and the relatively weak reflection of the upper layer of gas hydrate.  The second layered NGH model is a low hydrate saturation model and the pores of the underlying layer are also filled with free gas. The model parameters are given in Table 4. The velocity model and the quality factor Q are shown in Figure 15a,b, respectively. From which we can see that there are high-speed and low-Q anomalies in the hydrate layer. The seismic simulations with fractional-order differencing is shown in Figure 16, Figure 17. The amplitude of BSR is weaker than that of the first NGH model, but it is still obvious, and its polarity is reversed. There also exists an amplitude blank zone above the BSR and the relatively weak reflection of the upper layer of gas hydrate.        We also perform numerical experiments on the layered NGH model without gas bearing layer for high hydrate saturation model and low hydrate saturation model, respectively. For the former one, the models, seismic wave propagation and the section of the seismic records are shown in respectively; for the later one, the models, seismic wave propagation and the section of the seismic records are shown in Figures 21-23, respectively. Weak BSR can be seen in high hydrate saturation model, but it is almost invisible in low hydrate saturation model. Through the above layered NGH model test, when the underlying layer contains free gas, the BSR phenomenon is more obvious. When the underlying layer does not contain free gas, the BSR phenomenon is relatively weak. However, as the hydrate saturation increases, the BSR phenomenon becomes more obvious.

Stratum Serial Number Depth (m) Porosity Hydrate Saturation Gas
Energies 2020, 13, x FOR PEER REVIEW 19 of 25 We also perform numerical experiments on the layered NGH model without gas bearing layer for high hydrate saturation model and low hydrate saturation model, respectively. For the former one, the models, seismic wave propagation and the section of the seismic records are shown in respectively; for the later one, the models, seismic wave propagation and the section of the seismic records are shown in Figures 21-23, respectively. Weak BSR can be seen in high hydrate saturation model, but it is almost invisible in low hydrate saturation model. Through the above layered NGH model test, when the underlying layer contains free gas, the BSR phenomenon is more obvious. When the underlying layer does not contain free gas, the BSR phenomenon is relatively weak. However, as the hydrate saturation increases, the BSR phenomenon becomes more obvious.   We also perform numerical experiments on the layered NGH model without gas bearing layer for high hydrate saturation model and low hydrate saturation model, respectively. For the former one, the models, seismic wave propagation and the section of the seismic records are shown in respectively; for the later one, the models, seismic wave propagation and the section of the seismic records are shown in Figures 21-23, respectively. Weak BSR can be seen in high hydrate saturation model, but it is almost invisible in low hydrate saturation model. Through the above layered NGH model test, when the underlying layer contains free gas, the BSR phenomenon is more obvious. When the underlying layer does not contain free gas, the BSR phenomenon is relatively weak. However, as the hydrate saturation increases, the BSR phenomenon becomes more obvious.

Complex Seabed NGH Model
We consider using the White theory to model the hydrate formation, and then the fractional wave equation method is used to simulate the seismic wave of the established hydrate mode. The relationship between the content and elastic parameters of each rock component and the seismic wave velocity, attenuation characteristics, and seismic wave propagation law of the rock as a whole, lays a theoretical foundation for the use of seismic exploration to identify natural gas hydrates and estimate the hydrate content. The geological structure of the seabed is usually complex, and a large number of gas hydrates are distributed on the sea floor. The main identification mark of natural gas hydrate is bottom-simulating reflector (BSR). The position of BSR is usually consistent with the bottom interface of gas hydrate, which reflects that the fluctuation of bottom interface of gas hydrate is similar to that of seabed. We have established a complex seabed model containing natural gas hydrate, and there are a lot of "gas chimney" under the hydrate. The velocity and quality factor models are shown in the Figure 24a,b, respectively. Using the fractional Laplace operator wave Equation (17) to simulate wave field in hydrate formation, meanwhile the time derivative adopts finite difference, and the space derivative adopts the staggered grid fast Fourier transform method. In calculation, the first two first-order conservation equations in Equation (17) apply PML absorbing boundary conditions, and do not perform absorbing boundary processing on the equations containing Laplace operator. The snapshots of wave propagation and the seismic record are shown in Figures 25 and 26, respectively. As shown in Figure 25, BSR phenomenon can be clearly observed. Its main characteristics are the

Complex Seabed NGH Model
We consider using the White theory to model the hydrate formation, and then the fractional wave equation method is used to simulate the seismic wave of the established hydrate mode. The relationship between the content and elastic parameters of each rock component and the seismic wave velocity, attenuation characteristics, and seismic wave propagation law of the rock as a whole, lays a theoretical foundation for the use of seismic exploration to identify natural gas hydrates and estimate the hydrate content. The geological structure of the seabed is usually complex, and a large number of gas hydrates are distributed on the sea floor. The main identification mark of natural gas hydrate is bottom-simulating reflector (BSR). The position of BSR is usually consistent with the bottom interface of gas hydrate, which reflects that the fluctuation of bottom interface of gas hydrate is similar to that of seabed. We have established a complex seabed model containing natural gas hydrate, and there are a lot of "gas chimney" under the hydrate. The velocity and quality factor models are shown in the Figure 24a,b, respectively.

Complex Seabed NGH Model
We consider using the White theory to model the hydrate formation, and then the fractional wave equation method is used to simulate the seismic wave of the established hydrate mode. The relationship between the content and elastic parameters of each rock component and the seismic wave velocity, attenuation characteristics, and seismic wave propagation law of the rock as a whole, lays a theoretical foundation for the use of seismic exploration to identify natural gas hydrates and estimate the hydrate content. The geological structure of the seabed is usually complex, and a large number of gas hydrates are distributed on the sea floor. The main identification mark of natural gas hydrate is bottom-simulating reflector (BSR). The position of BSR is usually consistent with the bottom interface of gas hydrate, which reflects that the fluctuation of bottom interface of gas hydrate is similar to that of seabed. We have established a complex seabed model containing natural gas hydrate, and there are a lot of "gas chimney" under the hydrate. The velocity and quality factor models are shown in the Figure 24a  Using the fractional Laplace operator wave Equation (17) to simulate wave field in hydrate formation, meanwhile the time derivative adopts finite difference, and the space derivative adopts the staggered grid fast Fourier transform method. In calculation, the first two first-order conservation equations in Equation (17) apply PML absorbing boundary conditions, and do not perform absorbing boundary processing on the equations containing Laplace operator. The snapshots of wave propagation and the seismic record are shown in Figures 25 and 26, respectively. As shown in Figure 25, BSR phenomenon can be clearly observed. Its main characteristics are the Using the fractional Laplace operator wave Equation (17) to simulate wave field in hydrate formation, meanwhile the time derivative adopts finite difference, and the space derivative adopts the staggered grid fast Fourier transform method. In calculation, the first two first-order conservation equations in Equation (17) apply PML absorbing boundary conditions, and do not perform absorbing boundary processing on the equations containing Laplace operator. The snapshots of wave propagation and the seismic record are shown in Figures 25 and 26, respectively. As shown in Figure 25, BSR phenomenon can be clearly observed. Its main characteristics are the polarity reversal and the significantly increase in amplitude, especially at the junction of the hydrate bottom interface and the gas chimney. The same phenomenon can be seen in the seismic record Figure 26.

Discussion and Future Research
In order to better find gas hydrate accumulation areas through seismic exploration and estimate the scale of gas hydrate resources and hydrate content, it is necessary to in-depth study of the law and characteristics of seismic wave propagation in gas hydrate-bearing formations and establish the relationship of hydrate formation between seismic attributes and seismic response. For complex oceanary geological conditions, the mathematical analytical solution of seismic wave propagation is difficult to find. Numerical simulation based on the assumption of known underground physical parameters and medium models, uses computers and various numerical simulation methods to obtain seismic waves. The approximate solution of the field is a relatively simple and efficient method.
When the underlying stratum contains free gas, the BSR phenomenon is more obvious, with a significant increase in amplitude and polarity reversal. At the same time, a reflective blank zone in the hydrate layer can be seen above the BSR. When the underlying formation does not contain free gas, the BSR phenomenon is relatively weak. As the hydrate saturation increases, the BSR

Discussion and Future Research
In order to better find gas hydrate accumulation areas through seismic exploration and estimate the scale of gas hydrate resources and hydrate content, it is necessary to in-depth study of the law and characteristics of seismic wave propagation in gas hydrate-bearing formations and establish the relationship of hydrate formation between seismic attributes and seismic response. For complex oceanary geological conditions, the mathematical analytical solution of seismic wave propagation is difficult to find. Numerical simulation based on the assumption of known underground physical parameters and medium models, uses computers and various numerical simulation methods to obtain seismic waves. The approximate solution of the field is a relatively simple and efficient method.
When the underlying stratum contains free gas, the BSR phenomenon is more obvious, with a significant increase in amplitude and polarity reversal. At the same time, a reflective blank zone in the hydrate layer can be seen above the BSR. When the underlying formation does not contain free gas, the BSR phenomenon is relatively weak. As the hydrate saturation increases, the BSR

Discussion and Future Research
In order to better find gas hydrate accumulation areas through seismic exploration and estimate the scale of gas hydrate resources and hydrate content, it is necessary to in-depth study of the law and characteristics of seismic wave propagation in gas hydrate-bearing formations and establish the relationship of hydrate formation between seismic attributes and seismic response. For complex oceanary geological conditions, the mathematical analytical solution of seismic wave propagation is difficult to find. Numerical simulation based on the assumption of known underground physical parameters and medium models, uses computers and various numerical simulation methods to obtain seismic waves. The approximate solution of the field is a relatively simple and efficient method.
When the underlying stratum contains free gas, the BSR phenomenon is more obvious, with a significant increase in amplitude and polarity reversal. At the same time, a reflective blank zone in the hydrate layer can be seen above the BSR. When the underlying formation does not contain free gas, the BSR phenomenon is relatively weak. As the hydrate saturation increases, the BSR phenomenon becomes more obvious. On the basis of phenomena, such as BSR and amplitude blank zone, hydrates can be further identified, and the hydrate content can be estimated based on the characteristics of speed and quality factor. The fractional time derivative wave equation can accurately describe the constant-Q behavior, but the fractional time derivative will introduce huge memory consumption. We proposed a new finite-difference method, which can save memory resources and cost of computation. But for large-scale models, the amount of calculation is still relatively large. It is necessary to further study the fast algorithm of fractional time derivative. In addition, the current gas hydrate model is still relatively rough. Establish a more sophisticated model hydrate is still of great importance.
In this study, we only consider the forward simulation of the seismic propagation behavior in NGH. Next stage, we will consider the seismic imaging. The basic idea is performing a least-squares error minimization problem, i.e., if we denote d the simulated data by aforementioned seismic simulation steps with fractional differencing, d obs the observed data, m the velocity model, and L the forward operator that acts on m to generate d, then we will solve the following regularized minimization problem where Ω(m) denotes the a priori knowledge of the model, α > 0 is the regularization parameter.
The above optimization problem requires forward simulation and gradient or Hessian information calculation, e.g., least-squares migration (LSM) or full waveform inversion (FWI) can be developed. Of course, some advanced regularization techniques either physics-based or learned with artificial neural network (ANN) should be fully employed [36][37][38]. It is well-understanding that the easy way of seismic imaging is using reverse time migration (RTM) [39], however performing RTM requires a big deal of computation. We will continue to report our research in this field in future works.

Conclusions
We have introduced a more rigorous form of definition of fractional derivative to calculate the fractional time derivative of the viscoasoustic wave equation, which is more accurate and simpler than the other nearly constant-Q methods for describing the propagation of viscoasoustic wave. The rigorous form is the limit of the original definition when time step approaching zero. We discretize the rigorous form and bring it into the fractional time derivative wave equation to obtain the explicit expression of viscoelastic wave extrapolation. The computation time of the new method is merely about 20% of the original method at the same memory length. Some numerical examples demonstrate that the solution calculated by the new method is more accurate and stable in the uniform model or the large contrast velocity and Q model. By the new method, we can choose a smaller memory length at an acceptable accuracy; thereby, we can save memory resources and computation time. For the wavefield simulation of gas hydrate layer, White theory is selected to establish velocity and quality factor models, and then the fractional wave equation is used to simulate the propagation of seismic waves in the hydrate model. Finally, the connection between the elastic parameters and content of each component of the rock in the hydrate layer, the seismic wave speed, attenuation characteristics, and the law of seismic wave propagation are established. The above contents provide theoretical foundation for identifying gas hydrate and estimating hydrate content by seismic exploration.