Acoustic Velocity Log Numerical Simulation and Saturation Estimation of Gas Hydrate Reservoir in Shenhu Area, South China Sea

Gas hydrate model and free gas model are established, and two-phase theory (TPT) for numerical simulation of elastic wave velocity is adopted to investigate the unconsolidated deep-water sedimentary strata in Shenhu area, South China Sea. The relationships between compression wave (P wave) velocity and gas hydrate saturation, free gas saturation, and sediment porosity at site SH2 are studied, respectively, and gas hydrate saturation of research area is estimated by gas hydrate model. In depth of 50 to 245 m below seafloor (mbsf), as sediment porosity decreases, P wave velocity increases gradually; as gas hydrate saturation increases, P wave velocity increases gradually; as free gas saturation increases, P wave velocity decreases. This rule is almost consistent with the previous research result. In depth of 195 to 220 mbsf, the actual measurement of P wave velocity increases significantly relative to the P wave velocity of saturated water modeling, and this layer is determined to be rich in gas hydrate. The average value of gas hydrate saturation estimated from the TPT model is 23.2%, and the maximum saturation is 31.5%, which is basically in accordance with simplified three-phase equation (STPE), effective medium theory (EMT), resistivity log (Rt), and chloride anomaly method.


Introduction
Gas hydrate mainly exists in the seafloor and polar permafrost [1], and it owns the cage structure of solid crystal, which is formed by water molecules and natural gas (usually dominated by methane). The formation of gas hydrate needs a low-temperature and high-pressure environment, and the concentration of methane must exceed its solubility in the pore water. So gas hydrate is commonly distributed in the water depth greater than 300 m in the continental slope belt [2,3]. The submarine gas hydrate reserves mainly depend on the distribution of gas hydrate area, the thickness of gas hydrate stability zone, the porosity of sedimentary layer, the saturation of gas hydrate, and so on. However, the accurate estimation of gas hydrate reserves is very difficult due to the lack of research on the determination of gas hydrate distribution and gas hydrate saturation. As gas hydrate is rich in methane, it is associated with a series of scientific issues [4], including the global carbon cycle [5], global temperature changes [6], the sea-level rise [7], and future energy supply [8]. Therefore, researches on determining gas hydrate distribution and estimating gas hydrate saturation have become the focus of the scientists all over the world. In May 2007, the gas hydrate samples and various log data of gas hydrate zone were firstly obtained in Shenhu area, South China Sea, which made a significant breakthrough in exploration of gas hydrate in China. Meanwhile, it provided a great convenience for investigating the properties of the gas hydrate reservoir [9].
Geophysical logging is an important tool for evaluating gas hydrate saturation, and valuable information can be obtained by studying the resistivity and acoustic velocity log data. According to the log data in Shenhu area, South China Sea, Chinese scholars have utilized some theoretical models and empirical formulas to estimate gas hydrate saturation [10][11][12][13][14], and the results of these researches have greatly promoted the process of studying gas hydrate saturation by using log data in China. However, previous studies mainly 2 The Scientific World Journal focused on the estimation of the gas hydrate saturation, and the discrepancy of the estimation of gas hydrate saturation was large because of the application of various methods. Besides, previous studies did not systematically study the relationship between any two of the gas hydrate saturation, elastic wave velocity, and sediment porosity, which led to incomplete understanding of the log data.
Foreign scholars have carried out researches on the evaluation of the marine gas hydrate saturation earlier. A variety of theoretical models or experimental models have been proposed to estimate gas hydrate saturation, such as Wyllie et al. [15] time average equation with the seismic velocity [16][17][18], the effective medium theory [19][20][21][22][23], Biot-Gassmann theory model [24][25][26][27], compression wave (P wave) velocity of thermal-elastic theory [28], the three-phase equation (TPE) [29][30][31][32][33][34], and velocity model theory based on the two-phase theory (TPT) model [35]. Moreover, Tinivella et al. [36] made a research to compare the TPT model with the TPE model for evaluating gas hydrate saturations in marine sediments, and the comparison showed that the two theoretical approaches were in very good agreement. Based on this, the TPT model has been applied to verify the TPE model and estimate the gas hydrate and free gas saturations in several different areas [37][38][39].
In this work, based on the log data at site SH2 in Shenhu Area, we first establish gas hydrate model and free gas model by applying elastic wave velocity numerical model of the TPT method, then study the dependence of the P wave velocity on gas hydrate saturation, free gas saturation, and sediment porosity, and finally choose the gas hydrate model to estimate gas hydrate saturation at site SH2.

Numerical Simulation of Elastic Wave Velocity by the TPT Model
The TPT model [41,42], which supposes that the rock solid part is composed of rock matrix and gas hydrate and that the rock pore fluid is composed of free gas and water, can be used to study the elastic wave velocity model about the elastic characteristics of marine sand-shale reservoirs. Based on this theory and using the reported in Tinivella [35], the velocity relation between P wave velocity ( p ) and shear wave (S wave) velocity ( s ) is as follows: Where eff is the effective porosity, is the average rigidity of the skeleton, m is the average density, f is the density of the fluid phase, is the proportionality coefficient, is the coupling factor, and b , f and m are the compressibility of the solid phase, the fluid phase, and the matrix, respectively.

Geological Setting
Shenhu area is considered as one of the occurrences of gas hydrate, which is located in the middle of the northern slope of the South China Sea, between the Xisha Trough and the Dongsha Islands and Baiyun Sag of Zhu II depression of the Pearl River Mouth Basin ( Figure 1). Shenhu area experienced a geological evolution process similar to the northern margin of the South China Sea and eventually formed the regional sedimentary sequences in which marine sediments were the dominant composition [43]. Taking the Cenozoic sedimentary strata in the Shenhu area, for example, it is about 1000-7000 m, and the organic matter content is 0.46%-1.9% [44][45][46], which can provide the material base for gas hydrate. In recent years, some studies have confirmed that mud volcanoes, seafloor slips, mud diapers and other special structural units beneficial to the formation of gas hydrate are widely being developed in Shenhu area [47], and the bottom simulating reflections (BSRs) have also been identified using various geophysical methods in the northern South China Sea by Guangzhou Marine Geological Survey, China Geological Survey [48]. During April-June 2007, eight sites were drilled in Shenhu area (see Figure 2), among which, sites of SH2, SH3, and SH7 in water depth of 1105 to 1423 m were determined to contain gas hydrate in recovered core samples. The thickness of gas hydrate stability zone was about 10 to 25 m [9,49], and the sediment lithology in and above the zone was silt and silty clay respectively, according to the core data.

Data Set Used in the Study
As previous researchers have done much work on site SH2 [10-14, 40, 50], many valuable references can be used to obtain better results and verify the reliability of our research, so we select site SH2 as the research object. The water depth of site SH2 is 1232 m, and the maximum drilling depth is 245 mbsf. Figure 2 shows the conventional logs of site SH2. The measurement depth range is from 50 to 245 mbsf, and the measurement projects include caliper, density, natural gamma ray, acoustic, and resistivity logging. The occurrence of gas hydrate reservoir in this site appears to be "response characteristics of two high and two low" in the log curve, namely high resistivity, high natural gamma ray and low density, low acoustic time, especially for resistivity and acoustic logs. Besides, when the layer of well diameter changes, caliper curve can be used as the effective parameter to identify gas hydrate reservoir because the abnormality of other logs has nothing to do with the well condition.
Based on the previous research methods for the thickness of gas hydrate stability zone [51][52][53] and combined with the  analysis of conventional log data, the gas hydrate stability zone of site SH2 is determined to be at the depth of 195 to 220 mbsf [40].

Methodology
Seafloor sediments containing gas hydrate are generally composed of rock grain, gas hydrate, water, and natural gas. In order to research the characteristics of gas hydrate reservoir, gas hydrate model and free gas model have been established in this section, and based on these two models, the numerical simulation method and the TPT are used to study the dependence of the elastic wave velocity on sediment porosity, gas hydrate saturation and free gas saturation.

Establishment of Gas Hydrate
Model. The gas hydrate model assumes that the sediments are composed of rock grain, gas hydrate and water, and gas hydrate, is in the pore space, which is regarded as a part of the rock matrix. Supposing that s , h , w , and g represent the volume percentage of rock grain, gas hydrate, water, and free gas in the sediments respectively, the gas hydrate model can be expressed as where is sediment porosity.

4
The Scientific World Journal Gas hydrate saturation ( h ) and water saturation ( w ) can be written as The volume percentages of rock grain in solid phase ( s ) and gas hydrate in the solid phase ( h ) can be written, respectively, as

The Parameter Determination for Numerical Simulation of Gas Hydrate Model Based on the TPT.
In order to apply the TPT to gas hydrate model, some parameters in (1) should be known. The parameters can be determined by Tinivella and Schon's derivation formula [35,54].
(1) Effective porosity ( eff ) can be written as (2) Average density of sediments ( m ), density of the solid phase ( b ), and density of the fluid phase ( f ) can be written as where s is density of the rock grain, h is gas hydrate density and w is water density.
(3) Assume that the solid compressibility lies between the Voigt and Reuss averages [54]. b can be written as where s is rock grain compressibility and h is gas hydrate compressibility.
(4) Compressibility of the fluid phase ( f ) is where w is water compressibility.
(5) Compressibility of the matrix ( m ) indicates the compressibility of sediments without water, and it can be calculated by the following equation: where p is pore compressibility. The algorithm [35] to calculate p is: where 0 is the sediment porosity at the sea bottom and d is differential pressure.
(7) Shear modulus ( ) indicates average rigidity of the skeleton, and it can be calculated by the following equation: where h is gas hydrate rigidity and sm is the shear modulus of solid matrix with gas hydrate [55], which can be calculated by the following equation: where smKT is Kuster and Toksoz's shear modulus [56] and sm0 is the shear modulus of solid matrix without gas hydrate.

Implementation Steps of Numerical Simulation of Gas
Hydrate Model Based on the TPT. Consider the following.
(1) Given the , or calculate it by the empirical formula.

Establishment of Free Gas
Model. Free gas model assumes that the sediments are composed of rock grain, water, and free gas, and it can be expressed as The Scientific World Journal 5 Water saturation ( w ) and free gas saturation ( g ) can be expressed as (1) eff can be written as (2) m , b , and f can be written as where g is free gas density.
(4) Assume that the fluid compressibility lies between the Voigt and Reuss averages [54]. f can be written as where g is the compressibility of free gas. (5) m is calculated by the same equation as that used in gas hydrate model. (6) is calculated by the same equation as that used in gas hydrate model. (7) can be written as where is rock grain rigidity. (8) Variation range of is 1-∞.

Implementation Steps of Numerical Simulation of Free
Gas Model Based on the TPT. Consider the following.
(1) Given the , or calculate it by the empirical formula.

Estimation of Sediment Porosity.
Sediment porosity is a key parameter for the estimation of gas hydrate saturation for both gas hydrate and free gas models. Therefore, appropriate log data should be selected to estimate the sediment porosity at site SH2. The log data that can be applied to determine the sediment porosity include density, acoustic, resistivity, and neutron logging.
When it comes to the determination of sediment porosity by acoustic log, the data need to be corrected by the regional core data because seafloor sediments are always loose silt and silty clay. However, the core data is always insufficient, and it is difficult to determine the compaction correction coefficient, so acoustic log is not available. When it comes to the determination of sediment porosity by resistivity log, the Archie formula should be used to calculate porosity, and the Archie constants and formation water resistivity need to be known for the Archie formula [57]. However, the previous two parameters are generally determined by some empirical equations, and the estimation error of the sediment porosity is significant. As for the determination of sediment porosity by neutron log, it usually cannot be realized because of the lack of the neutron log data.
Compared with the resistivity, the acoustic and the density logs are less affected in gas hydrate reservoirs, and can generally reflect the situation of sediment porosity, so we select the density log data to estimate the sediment porosity in this study. The intensity of scattering gamma ray can be measured by the density log, which reflects electron density of the strata and volume density of rock ( b ). The estimation of by density log data can be expressed as [58] where ma is matrix density and f is fluid density. Considering the effect of shaly sediments, (29) can be written as [58] = where sh is the volume content of the shale, SH is the content index of the shale, sh is the density of the shale, GR is the value of natural gamma log in the research interval, GR min is the value of natural gamma log in pure sandstone interval, GR max the is value of natural gamma log in pure mud interval, and GCUR is the Hilchie index, which is 3.7 in the Tertiary of North America and 2 in old stratum [59]. Equations (29) and (30) are used to calculate the sediment porosity at site SH2, and processing parameters can be set as follows: ma = 2.65 g/cm 3 , f = 1.04 g/cm 3 , sh = 2.70 g/cm 3 [12,35]. The porosity calculated by (29) and (30) is close to each other (Figure 3), which varies in the range of 30% to 55%, and the average value is 45%. The result indicates that sediments at site SH2 are of high porosity.

Results
Gas hydrate model of the TPT is used to forward stimulate P wave velocity of sediment formation in different gas hydrate saturation conditions in depth of 50 to 245 mbsf at site SH2. Table 1 shows the values of the main parameters used to evaluate the velocity. When S h = 0, we can get the P wave velocity of water saturated sediments forward stimulated by gas hydrate model based on the TPT. From Figure 4, the tendency between actual curve of P wave velocity log and P wave velocity curve of the saturated water condition is almost consistent in the gas hydrate interval (above 195 mbsf), so the model and its parameters are rational for numerical simulation in this study. The difference between actual P wave velocity of the log and P wave velocity of saturated water condition reflects the value of gas hydrate or free gas saturation, which can be used to qualitatively identify the gas hydrate reservoir. The specific response characteristics are as follows: the possibility of containing gas hydrate is dominant when the actual P wave velocity of the log is higher than P wave velocity of saturated water condition; the possibility of containing free gas is dominant when actual P wave velocity of the log is lower than the P wave velocity of saturated water condition [35]. In depth of 195 to 220 mbsf, actual P wave velocity of the log is significantly higher than the P wave velocity of saturated water condition, so this interval is the gas hydrate stability zone. In depth of 220 to 245 mbsf, the actual P wave velocity of the log has an increase relative to the P wave velocity of saturated water condition. However, Figure 2 does not indicate the increase of resistivity. Without the coring analysis data in this interval, whether the abnormality is caused by gas hydrate or not can not be ascertained, and should be researched in further study. When S h gradually increases, P wave velocity made by forward stimulation also increases; when S h > 15%, P wave velocity increases significantly; when S h = 30%, P wave velocity curve is located at the right of the actual P wave velocity logging curve. The above results indicate that the basic range of gas hydrate saturation is 0-30% in depth of 50 to 245 mbsf at site SH2.
In order to study the dependence of P wave velocity on sediment porosity, gas hydrate saturation, it is assumed that values of gas hydrate saturation increase from 0 to 1 in the interval of 0.1. Using gas hydrate model of the TPT to model the corresponding P wave velocity of the previous gas hydrate saturations, the relation surface of previous three properties can be formed as Figure 5 shows SH2, the porosity presents a decreasing trend except for the abnormality caused by borehole conditions in some intervals, and P wave velocity of forward stimulation (S h = 0) slowly increases from 1743 to 1795 m/s. But with the increase of gas hydrate saturation, the increase rate of P wave velocity is obviously accelerated, and P wave velocity (burial depth is 51 mbsf) increases from 1743 to 3961 m/s. From the above analysis, the general rule between P wave velocity of forward stimulation and sediment porosity, gas hydrate saturation at site SH2 is the smaller the sediment porosity, the greater 8 The Scientific World Journal  the P wave velocity; the higher the gas hydrate saturation, the greater the P wave velocity. This result is basically in accordance with the research result made by Tinivella [35]. Similarly, it is assumed that values of free gas saturation increase from 0 to 1 in the interval of 0.1 in order to study the relation between P wave velocity and free gas saturation. Using free gas model of the TPT to model the corresponding P wave velocity of the previous free gas saturations, the relation surface of P wave velocity, sediment porosity, and free gas saturation can be formed as Figure 6 shows. In depth of 50 to 245 mbsf at site SH2, with the increase of free gas saturation, P wave velocity (burial depth is 51 mbsf) decreases from 1773 to 597 m/s (The velocity 597 m/s is obtained supposing 100% free gas saturation), and the decrease rate is significant. Considering the depth effect, the decrease rate of P wave velocity slows down with the increase of burial depth. From the above analysis, the general rule between P wave velocity of forward stimulation and free gas saturation at site SH2 is the higher the free gas saturation, the lower the P wave velocity. This result is also basically in accordance with the research result made by Tinivella [35]. The estimation of gas hydrate saturation for gas hydrate reservoir evaluation has an important significance. In order to estimate gas hydrate saturation in sediments, it is necessary to associate the P wave velocity of the log with the P wave velocity of gas hydrate model based on the TPT. Given an initial gas hydrate saturation, based on gas hydrate model of the TPT, the difference between P wave velocity of forward simulation and the actual P wave velocity of the log in this saturation can be acquired. If the difference is in the range of allowable error, the saturation can be treated as the actual saturation; if the difference does not satisfy the error requirement, the value of gas hydrate saturation should be modified until meeting the error precision.
Using gas hydrate model of the TPT to inverse gas hydrate saturation at site SH2, the values of the main parameters are listed in Table 1, and the inversion result is shown in Figure 7. In the interval of 50 to 90 mbsf at site SH2, the range of gas hydrate saturation is 0-17.5%, and the average value is 4.8%. As the shallow sediments are influenced by variation of borehole conditions, the estimation error of gas hydrate saturation in this interval is significant, which should be noticed during the analysis. In the interval of 90 to 195 mbsf, the range of gas hydrate saturation is 0-18.9%, and the average value is 7%. In the interval of 195 to 220 mbsf, the range of gas hydrate saturation is 7-31.5%, and the average value is 23.2%. With the increase of burial depth, gas hydrate saturation gradually increases and finally reaches the peak value of 31.5% in 208 mbsf. Then the gas hydrate saturation decreases slowly The Scientific World Journal with the increase of burial depth, and the range of gas hydrate saturation is 0-25.8% in the interval of 220 to 245 mbsf, with an average value of 15.5%.

Discussion
It is very important to determine the porosity for the evaluation of gas hydrate saturation. In order to analyze the accuracy of porosity estimation by density log data, we use resistivity log data and combine the Archie formula [57] to make a comparison between the estimation results, and the comparison results are shown in Figure 8. The porosity estimated by resistivity log data generally changes in the range of 30 to 50%, and the average value is 43% [12,64]. In the interval of 50 to 195 mbsf, the curves of density porosity and resistivity porosity are approximately coincident, while the former fluctuates due to the borehole effect. In the interval of 195 to 220 mbsf, the sediment contains gas hydrate, and the curve of resistivity porosity decreases significantly compared with the curve of density porosity due to the significant increase of resistivity (Figure 2), so the porosity calculated by resistivity log data needs to be corrected to exclude the influence of the increase of skeleton components. In the interval of 220 to 245 mbsf, the curves of density porosity and resistivity porosity are approximately coincident again. These results indicate that using density log data to estimate the porosity in the gas hydrate stability zone at site SH2 is relatively more reliable. According to previous studies, the range of core porosity by laboratory analysis in this interval was 40-55%, which was almost coincident with the estimation of porosity by density log data in this study [13]. As Figure 8 shows, the core porosity distribution corresponds with the curve of density porosity of the well, and this result proves that the porosity estimated by density log data can meet the requirement for evaluating gas hydrate saturation at site SH2.
In order to verify the accuracy of gas hydrate saturation estimated by the TPT, we compare gas hydrate saturation in this study with that estimated by Wang et al. [50] in the occurrence of gas hydrate (195 to 220 mbsf) at site SH2 ( Figure 9). The curve of gas hydrate saturation made by the TPT first increases and then decreases as the burial depth increases, which reaches the peak value of 31.5% in the 208 mbsfs and then decreases gradually. The peak value of gas hydrate saturation by the TPT is slightly smaller than the resistivity log (Rt) method (40.5%), the simplified three-phase equation (STPE) method (41%), and the effective medium theory (EMT) method (38.5%). However, the curve trend of gas hydrate saturation estimated by the TPT and other three methods is basically consistent, and the difference only lies in the amplitude of the curve, which indicates that using the TPT method to estimate gas hydrate saturation at site SH2 is available. The average value of gas hydrate saturations calculated by chloride anomaly method is 25%, and the peak value is 45% [9]. The peak value of gas hydrate saturation estimated by the TPT is relatively lower than that estimated by chloride anomaly method, but the average values estimated by the two methods are basically same, and most distribution dots of gas hydrate saturations obtained by chloride anomaly method correspond with the curve of gas hydrate saturation estimated by the TPT method ( Figure 9). This also indicates that using the TPT method to estimate gas hydrate saturation at site SH2 is available.

Conclusions
In summary, the relationships between P wave velocity and gas hydrate saturation, free gas saturation, and sediment porosity at site SH2 are studie, respectively, by virtue of elastic wave velocity of numerical stimulation based on the TPT, and gas hydrate model and free gas model are established to estimate the sediment porosity in order to determine the gas hydrate saturation of the research area. Some conclusions can be drawn as follows.
(1) Using the difference between P wave velocity of saturated water condition and actual P wave velocity of the log, whether the sediment contains gas hydrate or not can be identified quickly. In the interval of 195 to 220 mbsf at site SH2, the actual P wave velocity of the log increases significantly relative to P wave velocity of forward stimulation in saturated water condition, so this interval is determined to contain gas hydrate.
(2) By virtue of elastic wave velocity of numerical stimulation based on the TPT, combined with log data, the dependence of P wave velocity on gas hydrate saturation, free gas saturation, and sediment porosity at site SH2 can be analyzed, respectively. In the interval of 50 to 245 mbsf, as sediment porosity decreases, P wave velocity gradually increases; as gas hydrate saturation increases, P wave velocity increases; as free gas saturation increases, P wave velocity gradually decreases.
(3) The log data can be used to calculate gas hydrate saturation of the whole well, and the availability is better than the coring data. The average value of gas hydrate saturation estimated by the TPT is 23.2%, and the peak value is 31.5%, which is basically in accordance with the values estimated by the STPE model, the EMT model, the Rt model and chloride anomaly method.