Pearson et al. (1983)applied the three-phase time-average equation (Eq. 1) to hydrate bearing rocks and qualitatively explained the acoustic properties of hydrate bearing strata in the consolidation medium: 1 V p = ϕ ( 1 − S ) V w + ϕ S V h + ( 1 − ϕ ) V m . (1)

In Formula 1, V p is the P-wave velocity of hydrate bearing formation. V h is the P-wave velocity of pure hydrate. V w is the P-wave velocity of the fluid; V m is the P-wave velocity of the matrix. Φ is porosity; S is the saturation of the hydrate.

Similar to the three-phase time-average equation given by Pearson (1983), the three-phase equation Wood and Stoffa Pand Shipley (1994) of hydrate bearing sediments can be defined as 1 ρ V p 2 = ϕ ( 1 − S ) ρ w V w 2 + ϕ S ρ h V h 2 + ( 1 − ϕ ) ρ m V m 2 . (2)

In Formula 2, ρ h is the density of pure hydrate, ρ w is the density of the fluid, ρ m is the density of the matrix, and the volume density of the formation ρ = ( 1 − ϕ ) ρ m + ( 1 − S ) ϕ ρ w + S ϕ ρ h . (3)

In 1996, M.W. Lee et al. used the weight method of the three-phase time-average equation and the three-phase Wood equation used by Nobes et al. (1986)to estimate the velocity of deep-sea hydrate bearing sediments, and the three-phase weight equation was written as 1 V p = W ϕ ( 1 − S ) n V p 1 + 1 − W ϕ ( 1 − S ) n V p 2 . (4)

In Formula 3, V p 1 is the P-wave velocity calculated by Wood’s equation; V p 2 is the P-wave velocity calculated by the time-average equation. W is the weighting factor; n is the constant that simulates the variation of diagenesis with hydrate saturation.

The higher rock velocities estimated by time-average equations require the artificially lower matrix velocities. Based on the matrix velocity of 4.37 km/s and the assumption that hydrate is not present, the relationship between formation velocity and porosity is obtained by using the three-phase time-average equation, three-phase Wood equation, and three-phase weight equation respectively, as shown in Figure 1. The curve obtained by the weighted equation is located between the time-average equation and the Wood equation.

When the hydrate saturation is 0, when W > 1, the weighting equation is inclined to Wood’s equation, and when W < 1, the weighting equation is inclined to the time-average equation. Since (1-s)≤1, the weighting equation rapidly approaches the time-average equation as n increases (as shown in Figure 2). Therefore, when applying the three-phase weighting equation, a flexible method is provided by using the weighting factor and the exponential term, which is more suitable for consolidation (the time-average equation is more suitable) or suspension conditions (Wood’s equation is more suitable). The blue dots in Figure 3 are the actual data of a well in the study area. By changing the size of N, different data in different depth sections can be fitted.

The consolidation factor W, which controls the degree of consolidation, is determined by using the formation data without hydrate and fitting with the forward curve. By comparing the data from Well A, it can be seen that the fitting effect is best when W = 1.5.

Finally, using the W and N obtained by analysis, the P-wave velocity in the study area can be forward modeled, as shown in Figure 4. Among them, the blue curve is the measured curve, and the red curve is the forecast curve. The overall trend is generally consistent, but the prediction error of the free gas layer is large, so further stratified modeling may be needed.

As shown above, the most popular hydrate rock physical models are based on the six models proposed by Zhang (2012) (Figure 5). The six models refer to the contact cementation model, particle-enfolded model, skeleton particle support model, pore filling model, nodules fracture filling model, and fracture filling model respectively. In the contact cementation model, rock grains are regarded as freely accumulated spheres, and the hydrates bond the spheres together in contact. In the model of particle-enfolded, the rock particles are regarded as freely accumulated spheres, and the hydrate grows around the particles, acting as cement. In the skeleton particle support model, the hydrate is regarded as the bearing particles on the framework. In the pore filling model, hydrates are treated as particles or fluids to fill pores. In the nodular fissure filling mode, hydrate is deposited or filled in the fissure as nodules, which belongs to uneven distribution. The doping mode is to treat the hydrate as evenly distributed in the rock matrix, similar to the ice layer of the tundra. The filling mode of nodule fracture is to deposit hydrate as a nodule or fill it in the fracture, which belongs to a kind of heterogeneous distribution. However, the fracture filling model is not common in this work area, so the follow-up content has not carried out research on this model.

According to the analysis of mineral composition and porosity data of several Wells in the study area, the mineral composition in the study area is relatively single, mainly concentrated in clay (illite), quartz, and calcite. As shown in the multi-well histogram analysis (Figure 6), the porosity ranges from 30% to 55%, illite from 10% to 45%, quartz from 10% to 40%, and calcite from 0% to 15%.

Based on the mineral analysis results, we constructed measurement plates of five classical rock physical models and compared them with the actual data (Figure 7).

According to the rock physics template analysis, the hydrate enrichment form in this study area is the most consistent with the trend of the third type (skeleton particle support type) and the fourth type (pore filling type), which confirms the hydrate enrichment form in this area from the side. However, the measured values of longitudinal and S-wave velocity in the study area are slightly lower than the template curve, which may be caused by the incorrect elastic parameters of mineral components.

According to the previous study, most of the existing petrophysical models can be divided into two categories: one category regards hydrate as pore fluid, and the other category regards hydrate as skeleton mineral. Neither of the two models can well simulate the hydrate enrichment in the study area.

Through inductive analysis, we think that the two kinds of modeling ideas have some disadvantages. When hydrate is completely used as pore fluid, the modeling process uses Gassmann’s theory for fluid replacement, which results in the S-wave velocity of such models not changing with the hydrate content (as shown in Model 4). However, through statistical analysis, the shear modulus of hydrate is 2.57 MPa. It is precise because the shear modulus of hydrate is non-zero that its content must affect the S-wave velocity of the mixture. For the other types of models (for example, Model 3 and Model 5), which treat hydrate as skeleton minerals, the relationship between hydrate and pore fluid is ignored, resulting in an equal proportion of S-wave velocity with the increase of hydrate content. According to the template analysis, the correlation between hydrate saturation and velocity in this study area is between the above two categories. The S-wave velocity increases with the increase of hydrate content, but the increase rate is lower than that of the template of Model three and Model 4. Therefore, the existing hydrate models cannot meet the needs of hydrate modeling in this study area, and how to choose an appropriate equivalent medium theory to simulate hydrate is the core and key of modeling in this study area.

As shown in Figure 8, the two-phase modeling idea constructed in this study is between the existing pore-filled model and the skeleton supported model, which means that part of hydrate exists in the skeleton form, while the other part is enriched in the form of pore fluid. Based on the above ideas, the two-phase rock physical model we constructed is shown in Figure 9. As shown in the virtual frame in the figure, we divided the hydrate into two parts: the skeleton that plays a supporting role and the fluid in the pore. The hydrate, which acts as the skeleton, is combined with calcite, quartz, and illite to form the skeleton of the model, and they are fused by VRH averaging theory to obtain the dry rock skeleton. The hydrates enriched in the pores, together with water and gas, are added to the dry rock skeleton as pore fluid. Based on the above analysis, since the shear modulus of hydrate is non-zero, we use the BK solid substitution theory to simulate hydrate in pores. BK theory has been widely used in the simulation of rock physical model of unconventional shale in recent years, which is used to simulate the kerogen with non-zero shear modulus. In this study, we apply this idea to hydrate simulation. Gas and water are replaced by conventional Gassmann’s theory, and the equivalent hydrate mixture is finally obtained.

Based on the rock physical model of the two-phase hydrate constructed above, we carried out S-wave velocity prediction on Well E in the study area to test the feasibility of the model. The results are shown in Figure 10. The curve in the figure on the left is the S-wave velocity, the curve on the right is the P-wave velocity, the red curve represents the measured curve, and the blue curve represents the predicted curve. The error between the predicted velocity curve and the measured velocity curve is small, which proves that the rock physical model of double phase hydrate is available.

The main petrophysical model formulas used in this model are as follows:

Gassmann's equation is based on an isotropic assumption, which can not meet the needs of anisotropic shale petrophysical modeling. Brown and Korringa derived the Gassmann equation in anisotropic form. The Brown–Korringa formula can be used to describe the theoretical relationship of the equivalent modulus of an anisotropic rock skeleton when it is saturated with fluid. The specific formula is as follows: S i j k l d r y − S i j k l s a t = ( S i j α α d r y − S i j α α 0 ) ( S k l α α d r y − S k l α α 0 ) ( S α α β β d r y − S i α α β β 0 ) + ( β f l − β 0 ) ϕ , (16) where S i j k l d r y is the equivalent elastic flexibility tensor of rock skeleton, S i j k l s a t is the equivalent elastic flexibility tensor of saturated fluid rock, S i j α α 0 is the equivalent elastic flexibility tensor of constituent minerals, β f l and β 0 compressibility of pore fluid and mineral, ϕ is the porosity.

Based on the constructed model, we have carried out shear wave velocity prediction for well E in the work area, as shown in Figure 10. The left side is the shear wave velocity, the right side is the longitudinal wave velocity, the red curve is the measured curve, and the blue curve is the predicted curve. The error between the predicted results and the measured results is small, which proves the feasibility of the model.

Then we applied the previously proposed two-phase hydrate petrophysical model to other Wells in the study area, as shown in Figure 11. The predicted S-wave velocities are shown in the figure on the left, with the red curve representing the measured S-wave velocities and the blue curve representing the predicted S-wave velocities. The results show that the fitting effect of the model is good when the target well contains only a hydrate layer. However, when the target well contains not only the hydrate layer but also the free gas layer, the fitting effect of the model is still good for the hydrate layer, but poor for the free gas layer.

According to the analysis, the petrophysical model of double phase hydrate mainly focuses on the modeling of hydrate and does not consider the enrichment form of free gas too much. As shown in Figure 12, the longitudinal and S-wave velocity relationship between the hydrate layer and free gas layer was calculated. Red represents the hydrate layer and blue represents the free gas layer. As can be seen from the figure, the S-wave velocities at different depths in the hydrate layer are different, while the S-wave velocities at different depths in the free gas layer have almost no difference. It is verified from the side that hydrate saturation has a direct effect on the S-wave velocity due to its non-zero shear modulus, while free gas saturation has little effect on the S-wave velocity. Therefore, for the free gas layer, we need to build a targeted rock physical model of the free gas layer to improve the prediction accuracy of S-wave velocity.

Since the physical properties of the free gas layer and hydrate layer are quite different, it is difficult to simulate the characteristics of both by using the same petrophysical model. Therefore, we construct the idea of a segmented modeling method of hydrate + free gas. For the hydrate layer, we also use the rock physical model of double phase hydrate proposed above, while for the free gas layer, we build a new modeling idea. Since the free gas layer does not contain hydrate, we do not need to consider hydrate simulation. The specific idea is as follows: first, the main skeleton minerals such as calcite, quartz, and illite were fused with the VRH average theory to obtain a solid mixture. Then, the empty pores are added to the solid mixture background based on DEM theory. Finally, the Gassmann fluid substitution theory is used to add water and free gas to the dry rock skeleton to obtain the final equivalent medium Figure 13.

Based on the idea of piecewise modeling, we used different petrophysical models to predict the S-wave velocity for hydrate and free gas respectively, and then combined the prediction results of different intervals. The prediction results are shown in Figure 14. In the figure, the predicted S-wave curve based on the rock physical model of two-phase hydrate is shown on the left side, and the predicted S-wave curve based on the modeling method for free gas is shown in the middle. The red represents the measured curve, and the blue represents the predicted curve. It can be found that the prediction result of the free gas layer S-wave has been significantly improved by using the segmented modeling method, and the predicted curve is basically consistent with the measured curve.