Two-Phase Modeling Technology and Subsection Modeling Method of Natural Gas Hydrate: A Case Study in the Shenhu Sea Area

It is found that natural gas hydrate is not only a pore-filling material but also exists in the reservoir in the form of rock skeleton particles. Therefore, the traditional petrophysical simulation method cannot well describe the physical properties of natural gas hydrate reservoir. At the same time, the physical properties of the hydrate layer and its associated free gas layer are quite different, so it is difficult to fit the physical properties of the two media using traditional modeling methods. The two-phase modeling technology used in this paper is the equivalent medium modeling technology based on BK solid substitution theory and Gassmann fluid substitution theory, which simulates hydrate particles in rock skeleton and hydrate filling in pores, respectively. The forward simulation results show that the two-phase simulation technology of natural gas hydrate can well fit the P-wave and S-wave velocity information of the medium. At the same time, the equivalent medium model of the free gas reservoir is established by using only Gassmann fluid substitution theory. The practical application shows that the subsection modeling method can well solve the problem of the too large difference between the two sets of reservoir physical properties and make the calibration results of forward modeling synthetic records more accurate.


INTRODUCTION
Natural gas hydrate is a new type of clean energy that has been paid attention to since the 20th century. It is widely distributed in the continental shelf, deep-water basin, and permafrost area and has the characteristics of extremely rich resources. In the 1990s (Paull and Matsumoto, 2000;Wood and Ruppel, 2000;Boswell R et al., 2012;Yamamoto K, 2014), China gradually carried out the investigation of natural gas hydrate and completed the leap from exploration to experimental exploitation in more than 20 years (Zhang H Q et al., 2007;Wang X J et al., 2014;Yang et al., 2015;Yang et al., 2017a). Especially in the South China Sea, a large number of natural gas hydrate investigations and two rounds of experimental exploitation have been carried out (Li et al., 2018;Sha et al., 2019;Zhang et al., 2017;Zhang et al., 2018;Zhang et al., 2020a;Zhang et al., 2020b;Su et al., 2020;Yang et al., 2017b;Yang et al., 2020). At present, the investigation method of natural gas hydrate is still dominated by the seismic method, but the seismic exploration method can only describe the possible range and spatial distribution of ore body and cannot quantitatively characterize the distribution characteristics of hydrate saturation in the ore body. It is necessary to establish the mathematical expression relationship between P-wave, S-wave velocity, density, and other physical parameters obtained by seismic data inversion and saturation through the method of petrophysical modeling. In this way, the 1dimensional gas hydrate physical parameters obtained from logging data are extended to 2-dimensional and 3-dimensional space to complete the prediction of gas hydrate saturation (Worthington, 2008;Huang et al., 2012;Pan et al., 2019;Dong et al., 2020;Liu et al., 2020;Betlem et al., 2021;Yang et al., 2021).
With the development of natural gas hydrate, many petrophysical models suitable for natural gas hydrate have been developed.
Natural gas hydrate in sediments shows the properties of both fluid and solid skeleton, which have been found in laboratory research (Buffett and Zatsepina, 2000). With the efforts of many scholars, a variety of petrophysical theories suitable for simulating hydrate deposition have been proposed, including the weighted empirical formula (Lee et al., 1996), equivalent medium theory (Zhang and Toksoz, 2012), three-phase Biot type theory (Lee 2008), modified Biot Gassmann theory (Lee 2002), etc. The threephase Biot type equation proposed by Carcione and Tinivella (2000) assumes that the formation is composed of three phases of sediment, hydrate, and pore fluid to calculate the elastic wave velocity in the hydrate stability zone. Lee (2008) introduces parameters to describe pore filling and contact behavior on the basis of the three-phase Biot type; The equivalent medium theory and modified Biot Gassmann theory are essentially porefilling models for the simulation of natural gas hydrate bearing sediments . Zhang Yuwen et al. (2004) based on Biot's two-phase medium theory and aiming at the three hydrate deposition models proposed by Ecker, studied the variation law of the velocity and attenuation of fast P-wave, slow P-wave, and S-wave with frequency in hydrate bearing strata with and without dissipation. According to the research results of Dai et al. (2008); Zhang et al. (2011) summarized the hydrate petrophysical models based on equivalent media into five categories, and quantitatively analyzed the relationship between reservoir wave velocity and Poisson's ratio with hydrate saturation. Gao et al. (2012) used the improved Biot Gassmann (bgtl) model proposed by Lee to estimate the saturation of gas hydrate in well A by using P-wave velocity for the unconsolidated deep-water sedimentary strata with high porosity and silty clay in Shenhu sea area of the South China Sea. Pan et al. (2014) compared the saturation predicted by the effective medium model, the improved boit Gassmann model, and the simplified three-phase equation under the same occurrence state. It is found that the hydrate saturation predicted by the effective medium model and the improved Biot Gassmann model is more reasonable than that predicted by the simple three-phase equation (stpe).
Based on BK solid replacement theory and Gassmann fluid replacement theory, this paper establishes the equivalent medium two-phase modeling technology to characterize the two phase states of fluid filling and rock skeleton of natural gas hydrate respectively. This technology is applied to the natural gas hydrate work area in the Shenhu sea area of the South China Sea. At the same time, the subsection modeling method is used to fit the natural gas hydrate layer and free gas layer. The test shows that this method and technology has good applicability to the petrophysical modeling of natural gas hydrate in the Shenhu sea area, can well characterize the physical properties of natural gas hydrate and free gas reservoir in this work area, and can provide strong support for the quantitative study of natural gas hydrate saturation.

TRADITIONAL ROCK PHYSICAL MODELING METHODS
There are many types of hydrate petrophysical models. According to the previous research results and the actual situation of the study area, this study attempts to implement the existing model and analyzes the applicability of the existing model according to the actual data. 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:

Empirical Formula Class (Three-Phase Weight Equation)
( 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 ( 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 In 1996, M.W. Lee et al. used the weight method of the threephase 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 In Formula 3, V p1 is the P-wave velocity calculated by Wood's equation; V p2 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-  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 timeaverage 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.

Equivalent Medium Model (Five Classical Hydrate Models)
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 FIGURE 3 | Selection of formation consolidation factor W for three-phase weighting equation.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 884375 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.

TWO-PHASE MODELING TECHNOLOGY AND SEGMENTED MODELING METHOD 3.1 Petrophysical Modeling of Dual-Phase Hydrate
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:

Voigt-Reuss-Hill Boundary
Voigt-Reuss-Hill boundary is composed of Voigt upper bound, Reuss lower bound, and hill average. Voigt boundary is the upper FIGURE 7 | Rock physics template for five classic hydrate models.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 884375 limit of the VRH boundary, also known as equal strain average. Its results describe the average stress-strain relationship when each phase of the mixture is assumed to have equal strain. The specific formula is as follows: In the formula, M V and M i are the elastic modulus of the mixture and phase i, respectively, ϕ i is the volume fraction of phase i.
Reuss average is the lower bound of the VRH average, also known as the equal stress average, which assumes that each phase has the same stress. The specific formula is as follows: The meaning of each parameter in the formula is the same as that of the Voigt formula.
Hill pointed out that the arithmetic weighting of Voigt and Reuss can be used to predict the equivalent modulus of rock. The specific formula is as follows: It can be seen that the Voigt upper limit is the arithmetic weight of the elastic parameters of each phase, and the Reuss average is the arithmetic average of the reciprocal of the elastic parameters of each phase. Hill average is to average the equivalent results of the two. This averaging is based on the assumption of stress averaging or strain averaging. Therefore, when applying Voigt-Reuss-Hill upper and lower limits, it is necessary to assume that each component of the mixture is identical and the rock is linear and elastic.

Gassmann Equation
Gassmann's theory describes the propagation of seismic waves in porous media with saturated fluid, which belongs to the lowfrequency model. The equation constructs the relationship between seismic wave velocity, pore fluid, and mineral skeleton. Gassmann formula is widely used to calculate the change of elastic modulus caused by the change of fluid in pores (i.e. fluid substitution). The relationship between the elastic modulus of rock skeleton and the elastic modulus between porosity and matrix is as follows: In the formula, K dry , K m , K ϕ They are the equivalent elastic modulus of dry rock skeleton, matrix, and rock with pore space. ϕ is the porosity. Under low-frequency conditions, the relationship between the rock equivalent elastic modulus of fully saturated fluid and porosity is as follows: ( 9 ) FIGURE 9 | Flow chart of two-phase hydrate modeling. In the formula, K sat is the equivalent elastic modulus of rock in saturated fluid, K , ϕ can be expressed as In the formula, K fl is the elastic modulus of the fluid. Simultaneously from the abovementioned formula, eliminating K , which is the basic form of the Gassmann equation. Based on the basic equation, Xu and White (1995) changed the abovementioned formula to Gassmann equation assumes that the change of pore fluid does not affect the equivalent shear modulus of rock (μ sat ) , that is, the equivalent shear modulus of saturated fluid rock is equal to the shear modulus of rock skeleton (μ dry ), as shown in the following formula: Using the abovementioned formula, the equivalent elastic modulus of saturated fluid rock can be calculated, and when combined with the density equation, the P-wave velocity and the S-wave velocity of rock can be obtained.
The assumptions that Gassmann's theory needs to meet are as follows: 1) The rock skeleton is homogeneous and isotropic 2) All pores in the rock are interconnected, and there is no chemical reaction between pore fluid and skeleton; that is, fluid and skeleton are two independent systems 3) The model needs to meet the low-frequency assumption; that is, the pore fluid (without viscosity) can flow fully at the half wavelength of seismic wave, and the pore pressure is always in equilibrium 4) The fluid is completely saturated

Brown-Korringa Replacement Model
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: where S dry ijkl is the equivalent elastic flexibility tensor of rock skeleton, S sat ijkl is the equivalent elastic flexibility tensor of saturated fluid rock, S 0 ijαα is the equivalent elastic flexibility tensor of constituent minerals, β fl 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.

Segmented Modeling Method
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 nonzero 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.

OPTIMIZATION OF CALIBRATING DOWNHOLE SEISMIC BASED ON MODELING
Through analysis, it is found that the P-wave velocity of some Wells in the study area is discontinuous. In order to solve this problem, conventional commercial software is used to complete the missing velocity value through the interpolation method, but the prediction accuracy of velocity value based on the interpolation method cannot be guaranteed. Therefore, we  tried to use the constructed rock physical model to improve the accuracy of velocity prediction and finally improve the effect of calibrating downhole seismic. As shown in Figure 15, the logging curve in the left figure shows the absence of P-wave velocity at the boundary between the hydrate and free gas layers. The missing velocity is simply interpolated by commercial software, and the forward trace is calculated based on the interpolation velocity. However, there are some differences between the forward trace and seismic data, which may be caused by inaccurate velocity prediction. On the right of Figure 16 are the P-wave velocities calculated based on the model proposed in this study. Blue represents the predicted velocity and red represents the measured velocity. The synthetic records are obtained by using the velocities calculated by the rock physical model. It can be seen from the figure that in the P-wave velocity missing section, the forward modeling record has been significantly improved, and the reflection interface between the hydrate layer and the free gas layer has been clearly displayed. On the whole, the S-wave velocity prediction based on the rock physical model can improve the incomplete forward logging profile caused by the lack of velocity curve to some extent, and finally improve the precision of calibrating downhole seismic.

CONCLUSION
Through the application of actual data, two conclusions can be drawn: 1) The two-phase modeling technology based on BK solid substitution theory and Gassmann fluid substitution theory can take into account both the existing forms of hydrate as rock skeleton particle and pore-filling material. The equivalent medium model obtained by this technology can be used to fit more accurate longitudinal and S-wave velocities. 2) The physical properties of the gas hydrate layer and its associated free gas layer are very different, so the same modeling idea cannot be used. The segmenting modeling method can fully consider the physical properties of the two layers, which can greatly improve the effect of forward synthesis record calibration.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials; further inquiries can be directed to the corresponding author.