Next Article in Journal
Optimal Energy Efficiency Tracking in the Energy-Stored Quasi-Z-Source Inverter
Next Article in Special Issue
Factors Affecting Shale Gas Chemistry and Stable Isotope and Noble Gas Isotope Composition and Distribution: A Case Study of Lower Silurian Longmaxi Shale Gas, Sichuan Basin
Previous Article in Journal
Thermo-Mechanical Stress Comparison of a GaN and SiC MOSFET for Photovoltaic Applications
Previous Article in Special Issue
An Analytical Model for Production Analysis of Hydraulically Fractured Shale Gas Reservoirs Considering Irregular Stimulated Regions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Fractional Time Derivative Seismic Wave Equation Modeling for Natural Gas Hydrate

1
Key Laboratory of Petroleum Resources Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
2
College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
3
Innovation Academy of Earth Science, Chinese Academy of Sciences, Beijing 100029, China
*
Author to whom correspondence should be addressed.
Energies 2020, 13(22), 5901; https://doi.org/10.3390/en13225901
Submission received: 1 October 2020 / Revised: 29 October 2020 / Accepted: 3 November 2020 / Published: 12 November 2020
(This article belongs to the Special Issue Development of Unconventional Reservoirs 2020)

Abstract

:
Simulation of the seismic wave propagation in natural gas hydrate (NGH) is of great importance. To finely 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. Differently, 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 finite-difference 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.

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,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.

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 G L D t α f ( t ) = lim h 0 m h = t a 1 h α k = 0 m ω k ( α ) f ( t k h ) ,
where ω k ( α ) = ( 1 ) k ( α k ) = ( 1 ) k Γ ( α + 1 ) Γ ( k + 1 ) Γ ( α k + 1 ) , Γ ( ) 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
Γ ( x ) = 0 t x 1 e t d t ,     x > 0 .
For which it is well-known that Γ ( 1 ) = 1 and Γ ( x + 1 ) = x Γ ( x ) , for any x > 0 .
It is clear that if α = 1 ,   a G L D t α f ( t ) = f ( t ) f ( t h ) h ; and if α = 2 ,   a G L D t α f ( t ) = f ( t ) 2 f ( t h ) + f ( t 2 h ) h 2 . 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., L p ( a , b ) ) is defined by
  a + I t α f ( t ) = 1 Γ ( α ) a t f ( τ ) ( t τ ) 1 α d τ .
The Riemann–Liouville derivative with the lower integration limit a would be
  a + R L D t α f ( t ) = 1 Γ ( 1 α ) d d t a t f ( τ ) ( t τ ) α d τ ,
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 α f ( t ) = 1 Γ ( m α ) a t f ( m ) ( τ ) ( t τ ) α + 1 m d τ ,
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 ) = | t | ( α + 1 ) 2 Γ ( α ) cos ( β ) ,   α > 1 ,
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
D α f ( t ) = f ( t ) k ( t ) = 1 2 Γ ( α ) cos ( π α 2 ) | t τ | α 1 f ( τ ) d τ ,
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
α f ( t ) | t | α = C α ( a D t α f ( t ) + t D b α f ( t ) ) ,
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:
a D t α f ( t ) = 1 Γ ( 1 α ) d d t a t f ( τ ) ( t τ ) α d τ ,
t D b α f ( t ) = 1 Γ ( 1 α ) d d t t b f ( τ ) ( τ t ) α d τ .
In which, we assume that α ( 0 , 1 ) . If α is greater than 1, e.g., α ( 1 , 2 ] , the definitions can be made correspondingly as
a D t α f ( t ) = 1 Γ ( 2 α ) d 2 d t 2 a t f ( τ ) ( t τ ) α 1 d τ ,
t D b α f ( t ) = 1 Γ ( 2 α ) d 2 d t 2 t b f ( τ ) ( τ t ) α 1 d τ .

2.1.5. Relations between the above Fractional-Order Derivatives

For α ( m 1 , m ) , f ( t ) C m 1 ( [ a , t ] ) , and f ( m ) is continuous on [ a , t ] , there holds
a R L D t α f ( t ) =   a + G L D t α f ( t ) = k = 0 m 1 f ( k ) ( a ) ( x a ) k α Γ ( k + 1 α ) + 1 Γ ( m α ) a t f ( m ) ( ξ ) ( x a ) k α ( x ξ ) α m + 1 d ξ ,
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:
σ ( t ) = φ ( t ) ε ˙ ( t ) ,
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]
φ ( t ) = M 0 Γ ( 1 2 γ ) ( t t 0 ) 2 γ H ( t ) ,
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 0 2 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
2 2 γ σ t 2 2 γ = c 0 2 ω 0 2 γ cos 2 ( π γ / 2 ) 2 σ .
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 = 1 ρ 0 ( r ) σ ( r , t ) + S , t ε ( r , t ) = v ( r , t ) , σ ( r , t ) = M 0 ( r ) [ η ( r ) ( 2 ) γ ( r ) + τ ( r ) t ( 2 ) γ ( r ) 1 / 2 ] ε ( r , 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 0 2 γ ( r ) ( r ) ω 0 2 γ ( r ) cos ( π γ ( r ) ) and τ ( r ) = c 0 2 γ ( r ) 1 ( r ) ω 0 2 γ ( r ) 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.

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:
β f ( t ) t β = f h ( β ) ( t ) 1 h β j = 0 J ( 1 ) j ( β j ) f ( t + ( β 2 j ) h ) ,
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
( β j ) = β ( β 1 ) ( β 2 ) ( β j + 1 ) j ! = Γ ( β + 1 ) j ! Γ ( β j + 1 ) = Γ ( β + 1 ) Γ ( j + 1 ) Γ ( β j + 1 ) .
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:
a D t β f ( t ) = lim h 0 f h ( β ) ( t )                           = k = 0 m f ( k ) ( a ) ( t a ) β + k Γ ( β + k + 1 ) + 1 Γ ( β + m + 1 ) a t ( t τ ) m β f ( m + 1 ) ( τ ) d τ
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:
β f ( t ) t β lim h 0 f h ( β ) ( t ) = c 1 ( t a ) β f ( t 0 ) Γ ( 1 β ) + c 2 ( t a ) 1 β f ( t 0 ) Γ ( 2 β ) + 1 Γ ( 2 β ) a t ( t τ ) 1 β f ( 2 ) ( τ ) d τ ,
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 T h , and difference operators h , Δ h and δ h (respectively refers to the backward difference, forward difference and central difference), where h is the time sampling stepsize. The operator T h acts on the function f ( t ) , t , gives
T h f ( t ) = f ( t + h ) , h f ( t ) = f ( t ) f ( t h ) , Δ h f ( t ) = f ( t + h ) f ( t ) , δ h f ( t ) = f ( t + 1 2 h ) f ( t 1 2 h ) .
Assume that f ( t ) is sufficiently smooth, we have that the m -order derivative of f ( t ) can be approximated by f ( m ) ( t ) = D m f ( t ) , m , satisfying
D m f ( t ) = h m h m f ( t ) + O ( h ) = h m ( I T h ) m f ( t ) + O ( h )
or
D m f ( t ) = h m δ 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
h m = j = 0 m ( 1 ) j ( m j ) T j h
and
δ h m = j = 0 m ( 1 ) j ( m j ) T ( m / 2 j ) h
respectively. The above integer-order derivative forms can be extended to fractional cases, i.e., for any β , the fractional-order operator can be related by
h β = j = 0 ( 1 ) j ( β j ) T j h , δ h β = j = 0 ( 1 ) j ( β j ) T ( β / 2 j ) h .

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 [ 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:
β f ( t ) t β | t = t L = c 1 ( t L t 0 ) β f ( t 0 ) Γ ( 1 β ) + c 2 ( t L t 0 ) 1 β f ( t 0 ) Γ ( 2 β ) + 1 Γ ( 2 β ) t 0 t L ( t L τ ) 1 β f ( 2 ) ( τ ) d τ                       = c 1 ( t L t 0 ) β f ( t 0 ) Γ ( 1 β ) + c 2 ( t L t 0 ) 1 β ( f ( t 1 ) f ( t 0 ) ) Γ ( 2 β ) h +                             h β Γ ( 3 β ) j = 0 L 1 [ ( j + 1 ) 2 β ( j ) 2 β ] ( f ( t L j + 1 ) 2 f ( t L j ) + f ( t L j 1 ) ) + R h L ,
where R h L 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:
σ ( t L + 1 ) = Γ ( 3 β ) h β c 0 2 ω 0 2 γ cos 2 ( π γ / 2 ) 2 σ ( t L ) + 2 σ ( t L ) σ ( t L 1 )                         j = 1 L 1 [ ( j + 1 ) 2 β ( j ) 2 β ] ( σ ( t L j + 1 ) 2 σ ( t L j ) + σ ( t L j 1 ) )                         c 1 ( 2 β ) ( 1 β ) ( t L t 0 ) β σ ( t 0 ) h β c 2 ( 2 β ) ( t L t 0 ) 1 β [ σ ( t 1 ) σ ( t 0 ) ] h 1 β + S ( t L ) .
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 c 1 = c 2 = 0 .

2.5. 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 h i related to discrete step h is bounded by the infinitesimal equivalent value of h 2 , and the fact that for any 1 < α < 2 ,
( x + 1 ) α x α = k = 1 α ( α 1 ) ( α k + 1 ) k ! 1 x k α =   α 1 x 1 α + α ( α 1 ) 2 1 x 2 α + k = 3 α ( α 1 ) ( α k + 1 ) k ! 1 x k α .
The above expression is bounded. In next section, we will show how these values of memory length influence the results.

3. Results

3.1. 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. Figure 1, Figure 2, Figure 3 and Figure 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 Figure 1, Figure 2, Figure 3 and Figure 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
e r r rms = i = 1 n t ( d i p d i h ) 2 / i = 1 n t ( d i p ) 2 ,
where d i p denotes the reference solution calculated by the original method with all previous values and d i h 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.
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.

3.2. 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,b show snapshots and seismograms calculated by the original method with all previous values of the length L . Obvious false disturbances can be observed in the Figure 8a,b, which differ greatly from the reference solution Figure 10a,b. Figure 9a,b are almost identical to Figure 10a,b. These numerical examples demonstrate that the solution calculated by the new method is more accurate and stable in such a large contrast velocity and Q model. For the two-layers model, the computation time of the new method is also about 20% of the original method at the same memory length. By the new method, we can choose a smaller memory length and get a better solution, thereby we can further save memory resources and cost of computation.

3.3. 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.

3.4. 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 Figure 13 and Figure 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.
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 Figure 18, Figure 19 and Figure 20, respectively; for the later one, the models, seismic wave propagation and the section of the seismic records are shown in Figure 21, Figure 22 and Figure 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.

3.5. 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 Figure 25 and Figure 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.

4. 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
J ( m ) = 1 2 d obs d 2 2 + α Ω ( m ) = 1 2 d L m 2 2 + α Ω ( m ) min
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.

5. 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.

Author Contributions

Y.W. (Yanfei Wang) designed the study. Y.N. and Y.W. (Yanfei Wang) conducted experiments. Y.W. (Yanfei Wang), Y.N. and Y.W. (Yibo Wang) wrote the paper. All authors contributed to synthetic data interpretation and provided significant input to the final manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Original Innovation Program of CAS (grant no. ZDBS-LY-DQC003), the Key Research Program of the Institute of Geology & Geophysics, CAS (grant no. IGGCAS-201903), and the National Key R & D Program of the Ministry of Science and Technology of China (grant nos. 2018YFC0603500 and 2018YFC1504203).

Acknowledgments

We are grateful to four reviewers’ important questions and suggestions which make an improvement of our paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Guo, P.; Mcmechan, G.A. Evaluation of three first-order isotropic viscoelastic formulations based on the generalized standard linear solid. J. Seism. Explor. 2017, 26, 199–226. [Google Scholar]
  2. Kang, I.B.; McMechan, G.A. Separation of intrinsic and scattering Q based on frequency-dependent amplitude ratios of transmitted waves. J. Geophys. Res. 1994, 99, 23875–23885. [Google Scholar] [CrossRef]
  3. Zhu, T.; Harris, J.M. Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians. Geophysics 2014, 79, T105–T116. [Google Scholar] [CrossRef] [Green Version]
  4. Sun, J.; Fomel, S.; Zhu, T.; Hu, J. Q-compensated least-squares reverse time migration using low-rank one-step wave extrapolation. Geophysics 2016, 81, S271–S279. [Google Scholar] [CrossRef]
  5. Carcione, J.M. Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media; Elsevier: Amsterdam, The Netherlands, 2007. [Google Scholar]
  6. Štekl, I.; Pratt, R.G. Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators. Geophysics 1998, 63, 1779–1794. [Google Scholar] [CrossRef]
  7. Zhu, T.; Carcione, J.M.; Harris, J.M. Approximating constant-Q seismic propagation in the time domain. Geophys. Prospect. 2013, 61, 931–940. [Google Scholar] [CrossRef]
  8. Aki, K.; Richards, P.G. Quantitative Seismology: Theory and Methods; W. H. Freeman: San Franccisco, CA, USA, 1980; Volume 1. [Google Scholar]
  9. Blanch, J.O.; Robertsson, J.O.A.; Symes, W.W. Modeling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique. Geophysics 1995, 60, 176–184. [Google Scholar] [CrossRef] [Green Version]
  10. Blanch, J.O.; Symes, W.W. Efficient iterative viscoacoustic linearized inversion. In SEG Technical Program Expanded Abstracts 1995; Society of Exploration Geophysicists: Houston, TX, USA, 1995; pp. 627–630. [Google Scholar] [CrossRef]
  11. Robertsson, J.O.; Blanch, J.O.; Symes, W.W. Viscoacoustic finite-difference modeling. Geophysics 1994, 59, 1444–1456. [Google Scholar] [CrossRef] [Green Version]
  12. Carcione, J.M.; Cavallini, F.; Mainardi, F.; Hanyga, A. Time-domain seismic modeling of constant-Q wave propagation using fractional derivatives. Pure Appl. Geophys. 2002, 159, 1719–1736. [Google Scholar] [CrossRef]
  13. Kjartansson, E. Constant Q-wave propagation and attenuation. J. Geophys. Res. Space Phys. 1979, 84, 4737–4748. [Google Scholar] [CrossRef] [Green Version]
  14. Caputo, M.; Mainardi, F. A new dissipation model based on memory mechanism. Pure Appl. Geophys. 1971, 91, 134–147. [Google Scholar] [CrossRef]
  15. Carcione, J.M. Theory and modeling of constant-Q P- and S-waves using fractional time derivatives. Geophysics 2009, 74, T1–T11. [Google Scholar] [CrossRef]
  16. Chen, W.; Holm, S. Fractional Laplacian time-space models for linear and nonlinear lossy media exhibiting arbitrary frequency power-law dependency. J. Acoust. Soc. Am. 2004, 115, 1424–1430. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Carcione, J.M. A generalization of the Fourier pseudospectral method. Geophysics 2010, 75, A53–A56. [Google Scholar] [CrossRef]
  18. Song, G.; Zhang, X.; Wang, Z.; Chen, Y.; Chen, P. The asymptotic local finite-difference method of the fractional wave equation and its viscous seismic wavefield simulation. Geophysics 2020, 85, T179–T189. [Google Scholar] [CrossRef]
  19. Sun, J.; Zhu, T.; Fomel, S. Viscoacoustic modeling and imaging using low-rank approximation. Geophysics 2015, 80, A103–A108. [Google Scholar] [CrossRef] [Green Version]
  20. Fomel, S.; Ying, L.; Song, X. Seismic wave extrapolation using low rank symbol approximation. Geophys. Prospect. 2013, 61, 526–536. [Google Scholar] [CrossRef]
  21. Sun, J.; Fomel, S. Low-rank one-step wave extrapolation. SEG Tech. Program Expand. Abstr. 2013, 2013, 3905–3910. [Google Scholar] [CrossRef]
  22. Yao, J.; Zhu, T.; Hussain, F.; Kouri, D.J. Locally solving fractional Laplacian viscoacoustic wave equation using Hermite distributed approximating functional method. Geophysics 2017, 82, T59–T67. [Google Scholar] [CrossRef] [Green Version]
  23. Wang, N.; Zhu, T.; Zhou, H.; Chen, H.; Zhao, X.; Tian, Y. Fractional Laplacians viscoacoustic wavefield modeling with k-space-based time-stepping error compensating scheme. Geophysics 2020, 85, T1–T13. [Google Scholar] [CrossRef]
  24. Xu, Y.; Li, J.-Y.; Chen, X.; Pang, G. Solving fractional Laplacian visco-acoustic wave equations on complex-geometry domains using Grünwald-formula based radial basis collocation method. Comput. Math. Appl. 2020, 79, 2153–2167. [Google Scholar] [CrossRef]
  25. Abd-Elhameed, W.M.; Youssri, Y.H. A Novel Operational Matrix of Caputo Fractional Derivatives of Fibonacci Polynomials: Spectral Solutions of Fractional Differential Equations. Entropy 2016, 18, 345. [Google Scholar] [CrossRef]
  26. Youssri, Y.H.; Abd-Elhameed, W.M. Numerical spectral Legendre-Galerkin algorithm for solving time fractional telegraph equation. Rom. J. Phys. 2018, 63, 3–4. [Google Scholar]
  27. Podlubny, I. Fractional Differential Equations: Mathematics in Science and Engineering; Academic Press: Cambridge, MA, USA, 1999; Volume 198. [Google Scholar]
  28. Samko, S.G.; Kilbas, A.A.; Marichev, D.I. Fractional Integrals and Derivatives: Theory and Applications; Gordon and Breach Science Publishers: Philadelphia, PA, USA, 1993. [Google Scholar]
  29. White, J.E. Computed seismic speeds and attenuation in rocks with partial gas saturation. Geophysics 1975, 40, 224–232. [Google Scholar] [CrossRef]
  30. Shapiro, S.A.; Müller, T.M. Seismic signatures of permeability in heterogeneous porous media. Geophysics 1999, 64, 99–103. [Google Scholar] [CrossRef]
  31. Helgerud, M.B.; Dvorkin, J.; Nur, A.; Sakai, A.; Collett, T. Elastic-wave velocity in marine sediments with gas hydrates: Effective medium modeling. Geophys. Res. Lett. 1999, 26, 2021–2024. [Google Scholar] [CrossRef]
  32. Carcione, J.M.; Helle, H.B.; Pham, N.H. White’s model for wave propagation in partially saturated rocks: Comparison with poroelastic numerical experiments. Geophysics 2003, 68, 1389–1398. [Google Scholar] [CrossRef]
  33. Sumelka, W.; Voyiadjis, G.Z. A hyperelastic fractional damage material model with memory. Int. J. Solids Struct. 2017, 124, 151–160. [Google Scholar] [CrossRef]
  34. Jia, X.; Li, G. Numerical simulation and parameters inversion for non-symmetric two-sided fractional advection-dispersion equations. J. Inverse Ill Posed Probl. 2016, 24, 29–39. [Google Scholar] [CrossRef]
  35. Priest, J.A.; Best, A.I.; Clayton, C.R.I. Attenuation of seismic waves in methane gas hydrate-bearing sand. Geophys. J. Int. 2006, 164, 149–159. [Google Scholar] [CrossRef] [Green Version]
  36. Geng, Z.; Wang, Y. Automated design of a convolutional neural network with multi-scale filters for cost-efficient seismic data classification. Nat. Commun. 2020, 11, 3311. [Google Scholar] [CrossRef] [PubMed]
  37. He, Q.; Wang, Y. Reparameterized full waveform inversion using deep neural networks. Geophysics 2020. [Google Scholar] [CrossRef]
  38. Wang, Y.F.; Volkov, V.T.; Yagola, A.G. Basic Theory of Inverse Problems: Variational Analysis and Geoscientific Applications, 1st ed.; Science Press: Beijing, China, 2020; in press. [Google Scholar]
  39. Zhu, T.; Harris, J.M.; Biondi, B.L. Q-compensated reverse-time migration. Geophysics 2014, 79, S77–S87. [Google Scholar] [CrossRef] [Green Version]
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 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).
Energies 13 05901 g001
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 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).
Energies 13 05901 g002
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 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).
Energies 13 05901 g003
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 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.
Energies 13 05901 g004
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 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.
Energies 13 05901 g005aEnergies 13 05901 g005b
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.
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.
Energies 13 05901 g006
Figure 7. Comparison of seismograms in a two-layers model. 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 7. Comparison of seismograms in a two-layers model. 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.
Energies 13 05901 g007
Figure 8. (a) Snapshot (0.05 s) and (b) seismogram calculated by the original method with memory length L = 14.
Figure 8. (a) Snapshot (0.05 s) and (b) seismogram calculated by the original method with memory length L = 14.
Energies 13 05901 g008
Figure 9. (a) Snapshot (0.05 s) and (b) seismogram calculated by the new method with memory length L = 9.
Figure 9. (a) Snapshot (0.05 s) and (b) seismogram calculated by the new method with memory length L = 9.
Energies 13 05901 g009
Figure 10. (a) Snapshot (0.05 s) and (b) seismogram calculated by the original method with L equaling all previous values.
Figure 10. (a) Snapshot (0.05 s) and (b) seismogram calculated by the original method with L equaling all previous values.
Energies 13 05901 g010
Figure 11. The curve of velocity (a) and inverse quality factor (b) with hydrate saturation.
Figure 11. The curve of velocity (a) and inverse quality factor (b) with hydrate saturation.
Energies 13 05901 g011
Figure 12. (a) Velocity model; (b) Quality factor Q.
Figure 12. (a) Velocity model; (b) Quality factor Q.
Energies 13 05901 g012
Figure 13. (a) Seismic snapshot; (b) Seismic record.
Figure 13. (a) Seismic snapshot; (b) Seismic record.
Energies 13 05901 g013
Figure 14. One section of the seismic record.
Figure 14. One section of the seismic record.
Energies 13 05901 g014
Figure 15. (a) Velocity model; (b) Quality factor Q.
Figure 15. (a) Velocity model; (b) Quality factor Q.
Energies 13 05901 g015
Figure 16. (a) Seismic snapshot; (b) Seismic record.
Figure 16. (a) Seismic snapshot; (b) Seismic record.
Energies 13 05901 g016
Figure 17. One section of the seismic record.
Figure 17. One section of the seismic record.
Energies 13 05901 g017
Figure 18. (a) Velocity model; (b) Quality factor Q.
Figure 18. (a) Velocity model; (b) Quality factor Q.
Energies 13 05901 g018
Figure 19. (a) Seismic snapshot; (b) Seismic record.
Figure 19. (a) Seismic snapshot; (b) Seismic record.
Energies 13 05901 g019
Figure 20. One section of the seismic record.
Figure 20. One section of the seismic record.
Energies 13 05901 g020
Figure 21. (a) Velocity model; (b) Quality factor Q.
Figure 21. (a) Velocity model; (b) Quality factor Q.
Energies 13 05901 g021
Figure 22. (a) Seismic snapshot; (b) Seismic record.
Figure 22. (a) Seismic snapshot; (b) Seismic record.
Energies 13 05901 g022
Figure 23. One section of the seismic record.
Figure 23. One section of the seismic record.
Energies 13 05901 g023
Figure 24. (a) The seismic velocity model; (b) The quality factor model.
Figure 24. (a) The seismic velocity model; (b) The quality factor model.
Energies 13 05901 g024
Figure 25. (a) Seismic snapshots at 1.2 s; (b) Seismic snapshots at 1.5 s.
Figure 25. (a) Seismic snapshots at 1.2 s; (b) Seismic snapshots at 1.5 s.
Energies 13 05901 g025
Figure 26. The seismic record.
Figure 26. The seismic record.
Energies 13 05901 g026
Table 1. Simulation time with different memory lengths of L.
Table 1. Simulation time with different memory lengths of L.
Memory LengthsL = 9L = 14L = 24L = 34L = 44
CPU Time (s)Original method28644076010971413
New method5788152219282
Table 2. Parameters of each rock component.
Table 2. Parameters of each rock component.
Rock ParametersValues
Clay bulk modulus (Pa)20.9 × 109
Clay shear modulus (Pa)6.85 × 109
Clay density (kg/m3)2580
Quartz bulk modulus (Pa)36.6 × 109
Quartz shear modulus (Pa)45 × 109
Quartz density (kg/m3)2650
Sea water bulk modulus (Pa)2.55 × 109
Density of sea water (kg/m3)1050
Seawater viscosity coefficient (Pa·s)0.0018
Methane bulk modulus (Pa)0.01 × 109
Methane gas density (kg/m3)100
Methane viscosity coefficient (Pa·s)0.00002
Hydrate bulk modulus (Pa)5.6 × 109
Hydrate shear modulus (Pa)2.4 × 109
Hydrate density (kg/m3)920
Permeability (m2)100 × 10−14
Interparticle connection coefficient9
Table 3. High hydrate saturation model with gas bearing layer.
Table 3. High hydrate saturation model with gas bearing layer.
Stratum Serial NumberDepth (m)PorosityHydrate SaturationGas Saturation
1 (Seawater)1300---
2 (General sediment)2000.500
3 (Hydrate)500.50.40
4 (Hydrate)500.450.360
5 (Hydrate)500.40.320
6 (Hydrate)500.350.280
7 (Hydrate)500.350.280
8 (Free gas)1000.4200.01
9 (General sediment)2000.4200
Table 4. Low hydrate saturation model with gas bearing layer.
Table 4. Low hydrate saturation model with gas bearing layer.
Stratum Serial NumberDepth (m)PorosityHydrate SaturationGas Saturation
1 (Seawater)1300---
2 (General sediment)2000.500
3 (Hydrate)500.50.14290
4 (Hydrate)500.450.12860
5 (Hydrate)500.40.1140
6 (Hydrate)500.350.10
7 (Hydrate)500.350.10
8 (Free gas)1000.4200.01
9 (General sediment)2000.4200
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wang, Y.; Ning, Y.; Wang, Y. Fractional Time Derivative Seismic Wave Equation Modeling for Natural Gas Hydrate. Energies 2020, 13, 5901. https://doi.org/10.3390/en13225901

AMA Style

Wang Y, Ning Y, Wang Y. Fractional Time Derivative Seismic Wave Equation Modeling for Natural Gas Hydrate. Energies. 2020; 13(22):5901. https://doi.org/10.3390/en13225901

Chicago/Turabian Style

Wang, Yanfei, Yaxin Ning, and Yibo Wang. 2020. "Fractional Time Derivative Seismic Wave Equation Modeling for Natural Gas Hydrate" Energies 13, no. 22: 5901. https://doi.org/10.3390/en13225901

APA Style

Wang, Y., Ning, Y., & Wang, Y. (2020). Fractional Time Derivative Seismic Wave Equation Modeling for Natural Gas Hydrate. Energies, 13(22), 5901. https://doi.org/10.3390/en13225901

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop