The core observations made by Berge et al. (1999) showed that when the gas hydrate saturation exceeds a certain critical value, gas hydrate starts to cement around the sandstone particles. Corresponding laboratory core measurements indicated that when the gas hydrate saturation is greater than a certain critical value (i.e., 0.35), the acoustic velocity of the gas hydrate-bearing sandstone increases significantly with increasing gas hydrate saturation. While for the gas hydrate saturation less than this critical value, the acoustic velocity of the sandstone is not sensitive to changes in the hydrate saturation. The laboratory measurements made by Yun et al. (2005) also show that for high gas hydrate saturation, the shear modulus of the gas hydrate-bearing sandstone increases significantly. Therefore, it is inferred that the variation in the gas hydrate saturation affects depositional modes of the hydrate in sediments, and thus represents significant impact on elastic behaviors of the gas hydrate-bearing rocks.

In this study, we introduced the critical hydrate saturation parameter Sc in the rock physics model to describe various patterns of gas hydrate concentration. Figure 1A shows the case when the gas hydrate saturation S is less than the critical hydrate saturation Sc (S < Sc), where the hydrate exists only as pore filling. Figure 1B shows the case for S ≥ Sc, where the gas hydrate begins to attach to the solid mineral particles, becoming a part of the solid skeleton. At this time, the gas hydrate exists both as pore filling and particle cementation. Priest et al. (2006) showed that when existing as solid skeleton, the gas hydrate will have impact on porosity and permeability and can lead to squirt flow of the pore fluids, which in turn causes P-wave velocity dispersion and attenuation.

Figure 2 shows the rock physical modeling workflow corresponding to Figure 1. The critical gas hydrate saturation parameter Sc was introduced into the modeling process in order to describe different concentration patterns of hydrate depositing in sediments. When S < Sc, the gas hydrates only exist as pore filling that is mixed with brine. Wood’s theory (Wood, 1955) is used to calculate the properties of the fluid mixture. When S ≥ Sc, the gas hydrate saturation of the fluid mixture remains unchanged, and the gas hydrate begins to attach to the mineral particles. At this time, the concentration mode of the gas hydrate represents as the co-existence of pore filling and cementation. The elastic properties of the solid matrix are calculated using the Hill average theory (Hill, 1952) and the elastic modulus of the hydrate-bearing rock skeleton are calculated using the Krief’s theory (Krief et al., 1990). Finally, frequency-dependent P-wave velocity Vp and the attenuation factor 1/Qp caused by the squirt flow of pore fluids in the gas hydrate-bearing sandstone are calculated by the BISQ theory.

Based on the research of Qadrouh et al. (2015), elastic modulus of the solid skeleton of the gas hydrate-bearing sediments are calculated by the Krief model (Krief et al., 1990), where the bulk modulus K_dry and shear modulus G_dry of the solid skeleton are expressed as follows:

where the empirical parameter A is set to 3 to describe the case of unconsolidated sediment; φ is the porosity; G_s and K_s are the bulk modulus and shear modulus of the rock matrix, which are calculated using the Hill average theory (Hill, 1952):

where m is the number of solid components; f_i is the volume fraction of the i-th component; and K_i and G_i are the bulk modulus and shear modulus of the i-th component. In addition, the density of the solid matrix is calculated by averaging volume density of each individual component.

The bulk moduli of the fluid mixture composed of gas hydrate and brine is calculated using the Wood’s theory (Wood, 1955):

where K_h is the bulk modulus of hydrate; K_w is the bulk modulus of brine; S₁ is the saturation of hydrate as pore filling. When the hydrate represents both as the pore filling and particle cementation, S₁ = S_c.

For the pore filling mode, the gas hydrate saturation only affects elastic properties of the fluid mixture. For the model of hydrate co-existing as pore filling and particle cementation, the hydrate that is attached to solid particles reduces the porosity and permeability of the rock skeleton, and thus has the impact on squirt flow of pore fluids. The BISQ theory proposed by Dvorkin and Nur (1993) is used to calculate the velocity dispersion and attenuation of the gas hydrate-bearing reservoir. P-wave velocity V_p, attenuation coefficient α, and inverse quality factor Q⁻¹ are expressed as follows:

where

and

where K_dry and G_dry are the bulk and shear moduli of dry rock; φ is the porosity; ρ_s and ρ_f are the densities of the solid and fluid phases, respectively; ρ_a is the additional density of the inertial coupling between the solid and fluid phases; R is the squirt flow length; K is the bulk modulus of the solid phase; K_f is the bulk modulus of the fluid phase; η is the viscosity of the fluid; k is the permeability.

In rock physics modeling, it is considered that for the co-existence mode of hydrate both as pore filling and particle cementation, the cementation reduces the porosity of the formation by φ = φ₀[1−(S−S_c)], where S is the total saturation of the gas hydrate in the sandstone, and φ₀ is the porosity of the gas hydrate-bearing sandstone when the hydrate has not attached to the solid minerals yet. According to the Kozeny–Carman theory used by Santos et al. (1992), the permeability k of the gas hydrate formation varies with porosity φ and φ₀:

Figure 7 shows a geological model of a gas hydrate-bearing layer, which overlies above a layer containing free gas, and is overlaid by a brine-saturated layer. The P-wave velocity dispersion and attenuation of the gas hydrate-bearing layer are shown in Figure 3.

The P-wave velocity of the brine-bearing formation was 2628 m/s, and the density was 2140 kg/m³. The velocity dispersion and attenuation of the free gas layer are calculated using the patchy saturation theory (White, 1975), where the bulk modulus and shear modulus of the rock matrix are calculated using the Hill average theory (Hill, 1952). Properties used in modeling are shown in Table 1. The elastic modulus of the solid skeleton are calculated using the Krief model (Krief et al., 1990). Porosity and permeability are given in Table 2. Fluid properties used in the model are given in Table 3. Calculation results of the P-wave velocity dispersion and attenuation of the free gas layer are shown in Figure 8.

Seismic responses for the model in Figure 7 are calculated by integrating the built rock physics model and the propagator matrix method.

In the propagator matrix theory, the reflection and transmission coefficient vector r = [R_PP,R_PS,T_PP,T_PS]^T for the P-wave incidence is given by (Carcione, 2001):

where A₁ and A₂ are the propagator matrices related to the elastic parameters of the upper and lower media; B_α = T(0)T⁻¹(h_α)(α = 1,…,N) is the propagator matrix of the intermediate layer with an N layer structure; h_α is the thickness of each layer; and for a single-layered medium, N = 1. i_p is the P-wave incident vector, which is related to the physical parameters of the incident medium. Appendix A shows details of the propagator matrix theory.

The frequency-dependent reflection coefficient R_f of the PP wave (that is, R_pp) was calculated using the propagator matrix theory. It is multiplied by the seismic wavelet W_f in the frequency domain to obtain the corresponding amplitude spectrum U_f:

By implementing the inverse Fourier transform to U_f, the waveform U of reflection wave in the time domain is obtained:

where f is the angular frequency; i is an imaginary unit; t is time.

In Figure 7, the thickness of the gas hydrate reservoir is set as 50 m, and the dominant frequency of the Ricker wavelet is set as 40 Hz. The P-wave velocity dispersion and attenuation of the gas hydrate reservoir and the free gas are shown in Figures 3, 8, respectively.

Figure 9 shows the seismic AVO responses of the PP reflections for gas hydrate saturations of 0.2, 0.35, and 0.5. Among them, and the second type of AVO response is observed at the top interface of the gas hydrate reservoir. The interface between the gas hydrate reservoir and the free gas layer has a negative reflection amplitude, known as the BSR. As shown in Figure 9, a higher gas hydrate saturation corresponds to a higher BSR amplitude.

Figure 10 shows the mode-converted PS reflections for gas hydrate saturations of 0.2, 0.35, and 0.5. The interface between the gas hydrate-bearing reservoir and the free gas layer generates a mode-converted PS wave that has a positive phase of amplitude.

Figure 11A shows the post-stacked PP reflections obtained from AVO responses in Figure 9. Figure 11B shows the post-stacked mode-converted PS reflection obtained from Figure 10. For both the PP and PS reflections in Figure 11, the difference is very small for the two cases where gas hydrate saturations are 0.2 and 0.35, respectively. However, by comparing gas hydrate saturations of 0.35 and 0.5, the differences between the stacked reflection amplitudes are quite large. This can be interpreted by the analysis on Figure 4 where the Vp shows little variation with hydrate saturation for the case that the saturation is lower than the critical value 0.35, and represents significant variations for the hydrate saturation higher than the critical value 0.35.

Berge et al. (1999) and Dvorkin et al. (1999) showed that when free gas migrates upward to gas hydrate-bearing sediments, it may lead to heterogeneous hydrate saturation in the sediments. Thus, we design a model as shown in Figure 12A, where the hydrate saturation gradually decreases from about 0.55 at the bottom interface to the value just above zero at the top of the gas hydrate-bearing layer, which has a layer thickness about 28 m. According to the rock physics model built in this study, the P-wave velocity Vp dispersion and attenuation factor 1/Qp for varied hydrate saturation and frequency are calculated and shown in Figures 12B,C.

Based on the model of heterogeneous hydrate saturation in Figure 12A, we design a 2D model where the gas hydrate saturation varies vertically and layer thickness of the reservoir varies laterally, as shown in Figure 13. The thickness of the gas hydrate-bearing reservoir gradually decreases from 28 to 2.8 m from the center to both sides in the model. We calculate the AVO responses for each model in Figure 13, and then obtain the stacked waveforms as shown in Figure 14. In the calculation, the dominate frequency of the Ricker wavelet is 40 Hz, and the incidence angle ranges from 0 to 40 degrees. As can be seen in Figure 14, larger thickness of the gas hydrate-bearing reservoir corresponds to stronger reflection amplitudes from the BSR, Amplitude and phase of reflections vary with the thickness due to interference.