*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 di fferent 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 di ffer 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 ge<sup>t</sup> a better solution, thereby we can further save memory resources and cost of computation.

**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 8.** (**a**) Snapshot (0.05 s) and (**b**) seismogram calculated by the original method with memory length *L* = 14.

**Figure 9.** (**a**) Snapshot (0.05 s) and (**b**) seismogram calculated by the new method with memory length *L* = 9.

**Figure 10.** (**a**) Snapshot (0.05 s) and (**b**) seismogram calculated by the original method with *L* equaling all previous values.

#### *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.


**Table 2.** Parameters of each rock component.

**Figure 11.** The curve of velocity (**a**) and inverse quality factor (**b**) with hydrate saturation.

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


**Table 3.** High hydrate saturation model with gas bearing layer.

**Figure 12.** (**a**) Velocity model; (**b**) Quality factor *Q*.

**Figure 13.** (**a**) Seismic snapshot; (**b**) Seismic record.

**Figure 14.** One section of the seismic record.

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.

**Table 4.** Low hydrate saturation model with gas bearing layer.


**Figure 15.** (**a**) Velocity model; (**b**) Quality factor *Q*.

**Figure 16.** (**a**) Seismic snapshot; (**b**) Seismic record.

**Figure 17.** One section of the seismic record.

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 Figures 18–20, 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.

**Figure 18.** (**a**) Velocity model; (**b**) Quality factor *Q*.

**Figure 19.** (**a**) Seismic snapshot; (**b**) Seismic record.

**Figure 20.** One section of the seismic record.

**Figure 21.** (**a**) Velocity model; (**b**) Quality factor *Q*.

**Figure 22.** (**a**) Seismic snapshot; (**b**) Seismic record.

**Figure 23.** One section of the seismic record.

#### *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.

**Figure 24.** (**a**) The seismic velocity model; (**b**) The quality factor model.

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.

**Figure 25.** (**a**) Seismic snapshots at 1.2 s; (**b**) Seismic snapshots at 1.5 s.

**Figure 26.** The seismic record.

#### **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 grea<sup>t</sup> 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) = \frac{1}{2} \|d\_{\text{obs}} - d\|\_2^2 + \alpha \cdot \Omega(m) = \frac{1}{2} \|d - Lm\|\_2^2 + \alpha \cdot \Omega(m) \to \min \tag{27}$$

where <sup>Ω</sup>(*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–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.
