Introduction

Gas hydrate is regarded as a promising new energy natural resources in the twenty-first century by scientists. It has the advantages of high energy, extensive distribution, large scale, shallow burial and so on (Kvenvolden 1993; Ojha et al. 2010). The reserves of gas hydrate may be 2 times that of the total of gas, oil and coal we known. Many countries have carried out a lot of exploration and experiment study of gas hydrate for decades, such as America, Japan, Canada, etc. Researches show gas hydrate in land area is found only in the lower part of the permafrost area, in contrast the amount of gas hydrate in sea area is very large (Collett et al. 2009; Sloan 2003). Due to the limitation of temperature and pressure, the formation and distribution of gas hydrate in marine sediments are diversified and complex (Cook et al. 2010). Hence, accurately exploring the location of gas hydrate is challenging, which has been widely concerned by researchers.

There are various exploration methods of gas hydrate, for example, geophysical techniques, geochemistry techniques (Edwards 1997, 2005; Latychev and Edwarrds 2003), geomorphic survey (Kvenvolden and Kastner 1990), underwater imaging survey (Milkov et al. 2003) and so on. At present, seismic identification techniques are still the most effective and practical exploration method in general use. Bottom simulating reflectors (BSR) is an important symbol to identify gas hydrate in seismic profiles, which are approximately parallel to the seafloor, opposite to the polarity of seabed reflector and are usually accompanied by blank zones (Shipley et al. 1979). However, sometimes BSR may not be present where the hydrate occurs. Therefore, in some cases, it is hard to identify gas hydrate just through the identification of BSR. When gas hydrate filled in some pore space of seawater-saturated sediments, the compressional and shear wave velocity will increase (Stoll and Bryan 1979). However, other minerals (such as limestone, opal, pyrite and so on) can also cause an increase in velocity. Liu et al. (2009) proposed to identify gas hydrates through the ratio between compressional wave velocity increment and the shear wave velocity increment (△Vp/△Vs). They discovered the △Vp/△Vs of gas hydrate decreases with increasing saturation, while the △Vp/△Vs of other minerals increase with saturation increasing. This approach has worked in some hydrate areas, but accurate velocity data are needed.

In this paper, we experiment some physical properties of a suite of different minerals and gas hydrate bearing sediments, finally proposing a new approach to identify marine gas hydrate based on the properties of compressional wave velocity and density. The physical properties of sediments can be predicted by many rock physics models (Lee et al. 1996; Ecker et al. 1998; Carcione and Tinivella 2000; Gei and Carcione 2003; Jakobsen et al. 2000). In our study, we obtain compressional wave velocity and density via the effective medium theory, which is widely applied to simulate gas hydrate bearing sediments (Helgerud et al. 1999; Lee and Collett 2013; Ecker et al. 2000; Chand et al. 2004). In addition, we use the data of borehole measurements in the Dongsha area, of the South China Sea, and the Hydrate Ridge, along the Oregon continental margin, to test our new approach.

Theory

We use the effective medium theory (see Appendix 1) to model the unconsolidated marine sediments. Assuming a part of the pore space of seawater-saturated sediments is filled respectively by gas hydrate, clay, sand and other authigenic minerals, which are regarded as a part of solid frame sediments here. Moreover, we will further discuss other gas hydrate morphology in greater detail. The related parameters are listed in the Table 1. The matrix minerals constituent of sediments are set as 5% quarts, 35% calcite and 60% clay (Ecker et al. 2000; Mavko and Mukerji 1996; Wang et al. 2013; Liu et al. 2009). We finally analyze the physical properties of different minerals-filled sediments.

Table 1 Calculation parameters

The compressional wave velocity of different minerals-filled sediments varies with filling minerals’ saturation at different primary porosity (50%, 60%, 70% and 80%) modeled by effective medium theory, as shown in Fig. 1. The saturation refers to the volume fraction of filling minerals in the pore space of sediments. The velocity of all minerals-filled sediments increases with saturation increasing. The variation tendencies are consistent, making hard to identify gas hydrate from other minerals. When the saturation is 0%, it indicates the sediment is 100% seawater-saturated sediment.

Fig. 1
figure 1

Compressional wave velocity of different minerals-filled sediments as a function of saturation. The primary porosity of sediments are a 50%, b 60%, c 70% and d 80%, respectively

Figure 2 shows the density of different minerals-filled sediments changes with the saturation increasing at different primary porosity (50%, 60%, 70% and 80%). The density of gas hydrate-filled sediment decreases with the saturation increasing. However, for other minerals-filled sediments, the density variation tendencies increase. The reason for this phenomenon is that gas hydrate filled part of pore space replacing seawater. Gas hydrate density is slightly lower than that of seawater, whereas the density of other minerals is exceeding that of seawater. Nonetheless, although there are differences in density, we cannot solely distinguish gas hydrates from sediments by this property, which is poorly reliable in some degree.

Fig. 2
figure 2

Density of different minerals-filled sediments as a function of saturation. The primary porosity of sediments are a 50%, b 60%, c 70% and d 80%, respectively

From the above analysis, neither compressional wave velocity nor density can be used directly to discriminate gas hydrate. Therefore, we propose an attribute (Z) in order to distinguish gas hydrate from other minerals in sediments. The Z should have the characteristic of opposite variation tendencies with saturation increasing for gas hydrate and other minerals filled sediments. The most effective of combinations of properties is \(Z = V{p^{0.5}} /\rho^{5}\), as illustrated in Fig. 3. The \(Vp^{0.5} /\rho^{5}\) increases with saturation for gas hydrate, but shows a decreasing-trend in other minerals. When the primary porosity is certain, there is a threshold value at the saturation of 0 that indicates the \(Vp^{0.5} /\rho^{5}\) of 100% seawater-saturated sediment. The \(Vp^{0.5} /\rho^{5}\) caused by filling hydrate in pore space is larger than the threshold value, nevertheless, the \(Vp^{0.5} /\rho^{5}\) caused by filling other minerals is smaller than the threshold value. Moreover, the threshold value is related to the primary porosity, and the higher the porosity, the higher the threshold value and the greater the difference between gas hydrate and other minerals. As a consequence, the \(Vp^{0.5} /\rho^{5}\) can be utilized to distinguish gas hydrate in marine sediments.

Fig. 3
figure 3

\(Z = Vp^{0.5} /\rho^{5}\) of different minerals-filled sediments as a function of saturation. The primary porosity of sediments are a 50%, b 60%, c 70% and d 80%, respectively

Application

To validate our method of identifying gas hydrate in marine sediments, we choose the Dongsha area of the South China Sea and the Hydrate Ridge along the Oregon continental margin as our study areas. Both sites have been discovered to contain gas hydrates as confirmed via logging-while-drilling (LWD) and wireline log (Sha et al. 2015; Tréhu et al. 2003). In 2007 and 2013, the Guangzhou Marine Geologic Survey carried out two gas hydrate explorations in South China Sea, respectively GMGS1 and GMGS2, and obtained hydrate samples (Zhang et al. 2015; Wang et al. 2006). In addition, a series of seismic data and LWD data were also obtained during the Ocean Drilling Program (ODP) Leg 204 in Hydrate Ridge (Tréhu and Bangs 2001; Tréhu et al. 2002, 2003). Logging data of GMGS2 Sites 05 and 08, and Leg 204 Sites 1244E and 1245E are applied in our study (Fig. 4). Non-hydrate zones are used to fit compressional velocity (\(V_{p0}^{{}}\)) and density (\(\rho_{0}\)) of 100% seawater-saturated sediments by least square method, then the threshold value \(Z_{0}\) (\(Z_{0} = V_{p0}^{0.5} /\rho_{0}^{5}\)) can be derived.

Fig. 4
figure 4

Location of GMGS2 and Leg 204 expedition area

GMGS2 site 05

At Site 05, the logging data shows that there are velocity and density anomalies between 193 and 206.5 m below mudline (BML), where gas hydrate has been discovered (Sha et al. 2015). In Fig. 5a, b and c, the velocity, density and resistivity well logs of site 05 are shown, and the area highlighted in yellow represents gas hydrate zone. The non-hydrate zones are used to fit the background trends of velocity and density with red lines, which are denoted as \(V_{p0}^{{}}\) and \(\rho_{{0}}\). We compute the results of \(\Delta Z_{{}} = Vp^{0.5} /\rho^{5} - Vp_{0}^{0.5} /\rho_{0}^{5}\), when \(\Delta Z < 0\), we set it to − 0.007 in order to better display in the figure. If the calculated value \(Vp^{0.5} /\rho^{5}\) is greater than the threshold value, then we can judge there is gas hydrate, otherwise, no gas hydrate is in here. From Fig. 5d, we can easily find most of points (\(\Delta Z > 0\)) are concentrated in the hydrate zone, and very few points are in other places. The nonconforming points account for 0.68% of the total points. Therefore, the result of our method at site 05 agrees with the well logging result basically.

Fig. 5
figure 5

At Site 05: a compressional wave velocity (blue line) of logging data and a linear background trend (red line); b density (blue line) of logging data and a linear background trend (red line); c resistivity of logging data; d the result of \(\Delta Z_{{}} = V_{p}^{0.5} /\rho^{5} - V_{p0}^{0.5} /\rho_{0}^{5}\) (blue dots). The dots (\(\Delta Z < 0\)) are set as − 0.007 for clarity. Gas hydrate bearing zone is marked by yellow area

GMGS2 site 08

The gas hydrate accumulates in two layers at Site 08, which are marked by yellow in Fig. 6. The upper layer is at 9–22 m BML with low saturation. And the lower layer is at 65–98 m BML with obviously high velocity, high resistivity and low density, respectively. In addition, there is a carbonate layer within the depth interval at 58–62 m BML (Sha et al. 2015). We calculate the background trends on the basis of the non-hydrate zones and get the results of \(\Delta Z\). If \(\Delta Z < 0\), we set it to -0.19 for displaying clearly. The data points mainly occur at gas hydrate zones, especially, in the lower layer the results coincide well with the gas hydrate zone. However, the results of the upper layer are not so good as the lower layer. We infer that it is due to the low saturation of hydrate at the upper layer. The unmatched points account for 3.43% of the total points. In summary, the test results conform to the rule at Site 08.

Fig. 6
figure 6

At Site 08: a compressional wave velocity (blue line) of logging data and a linear background trend (red line); b density (blue line) of logging data and a linear background trend (red line); c resistivity of logging data; d the result of \(\Delta Z_{{}} = V_{p}^{0.5} /\rho^{5} - V_{p0}^{0.5} /\rho_{0}^{5}\) (blue dots). The dots (\(\Delta Z < 0\)) are set as − 0.19 for clarity. Gas hydrate bearing zones are marked by yellow area

Leg 204 site 1244E

According to previous reports, the presence of gas hydrate is located in 40–130 m BML at Site 1244E by direct sampling, geochemical core analyses, downhole logging data and so on (Tréhu et al., 2003). In Fig. 7, the gas hydrate zone is marked in yellow and we use least square method to obtain background trends shown in red lines. Similarly, the points (\(\Delta Z < 0\)) are set as − 0.007 for the sake of clarity. The test results of \(\Delta Z\) indicate that the majority of points (\(\Delta Z > 0\)) exist at hydrate interval. However, there is some points (\(\Delta Z > 0\)) between 130–135 m BML, we infer the hydrate position in data may not be accurate or there may be errors of fitting. The nonconforming points account for 1.63% of the total points. In general, the results of test fit relatively well with actual condition at Site 1244E.

Fig. 7
figure 7

At Site 1244E: a compressional wave velocity (blue line) of logging data and a linear background trend (red line); b density (blue line) of logging data and a linear background trend (red line); c resistivity of logging data; d the result of \(\Delta Z_{{}} = V_{p}^{0.5} /\rho^{5} - V_{p0}^{0.5} /\rho_{0}^{5}\) (blue dots). The dots (\(\Delta Z < 0\)) are set as − 0.007 for clarity. Gas hydrate bearing zone is marked by yellow area

Leg 204 site 1245E

At Site 1245E, gas hydrate occurs from ~ 50 to 131 m BML with high resistivity and velocity (Tréhu et al. 2003). We draw the velocity, density and resistivity well logs in Fig. 8a, b and c, and the hydrate location is shown in yellow. Likewise, the background trends of velocity and density are fitted. We calculate the value of \(\Delta Z\) and set the value of points (\(\Delta Z < 0\)) as − 0.01 just for displaying. From Fig. 8d, in hydrate zone there are most points (\(\Delta Z > 0\)) and a few discrete points (\(\Delta Z > 0\)) appear in non-hydrate zone. The nonconforming points account for 1.42% of the total points. The results can basically indicate the gas hydrate zone, which is consistent with well log measurement.

Fig. 8
figure 8

At Site 1245E: a compressional wave velocity (blue line) of logging data and a linear background trend (red line); b density (blue line) of logging data and a linear background trend (red line); c resistivity of logging data; d the result of \(\Delta Z_{{}} = V_{p}^{0.5} /\rho^{5} - V_{p0}^{0.5} /\rho_{0}^{5}\) (blue dots). The dots (\(\Delta Z < 0\)) are set as − 0.01 for clarity. Gas hydrate bearing zone is marked by yellow area

Discussion

In this new method, the results of application are generally in agreement with actual hydrate position at GMGS2 Site 05, GMGS2 Site 08, Leg 204 Site 1244E and Leg 204 Site 1245E. Nonetheless, we also find some issues. At Site 08, the hydrate zone is not fully identified by using \(Vp^{0.5} /\rho^{5}\) in the upper layer, whose saturation is relatively low. Hence, the velocity and density anomalies are correspondingly not obvious, which may create the observed errors. Moreover, the existence of background line fitting errors of velocity and density logs might even be possible.

The values of a and b

The new attribute is based on the different physical properties caused by the difference in the internal structure of hydrate and other minerals. It is worth mentioning that the a and b of \(Vp^{{\text{a}}} /\rho^{b}\) are not fixed values. Due to the small velocity difference and big density difference between hydrate and other minerals, we can reduce the effect of the velocity difference and increase the effect of the density difference in order to get a characteristic attribute. We suggest the values range can be: 0 < a < 1 and b > 1. When a = 0.5 and b = 5, the attribute can generally satisfy most of high porosity unconsolidated conditions. For example, we find different primary porosity have different effects on \(Vp^{0.5} /\rho^{5}\) from Fig. 3. When the primary porosity is 60%-80%, \(Vp^{0.5} /\rho^{5}\) presents a better ability to distinguish hydrate compared to low porosity. We can make a concrete analysis of the a and b according to porosity of each specific study area before identification process. In our tests, the attribute \(Vp^{0.5} /\rho^{5}\) has been proved to be adaptive for hydrate identification in unconsolidated marine sediments.

The effects of lithology on the identification

The lithology factor is crucial to be considered in the identification of gas hydrates, playing an important role in target reservoir formation and distribution. In order to research the effects of lithology, we compute how the \(Z = Vp^{0.5} /\rho^{5}\) of different minerals-filled sediments varied with saturation in different lithology (host sediments): (1) pure clay; (2) 50% clay and 50% quartz; (3) pure quartz (Fig. 9). We set the primary porosity of sediments is 70%. From Fig. 9, we can discover the lithology just have a slight effect on the identification. The main effects are on the threshold values of \(Vp^{0.5} /\rho^{5}\), which are 0.1666, 0.1629 and 0.1599, respectively.

Fig. 9
figure 9

\(Z = Vp^{0.5} /\rho^{5}\) of different minerals-filled sediments as a function of saturation in different lithology: (1) pure clay; (2) 50% clay and 50% quartz; (3) pure quartz. The primary porosity of sediments is 70%

The effects of morphology on the identification

In nature, there are two basic distribution morphologies of gas hydrate observed: pore-filling and fracture-filling. The pore-filling gas hydrate occurs in pores and can be regarded as isotropic and homogeneous medium, whose elastic properties are calculated by effective medium theory in our study. However, the fracture-filling gas hydrate generally accumulates in fracture as nodules layers or veins and can be regarded as anisotropic medium. Therefore, we use the two-laminated media model (see Appendix 2) that is an anisotropic rock physics model proposed by Lee and Collett (2009) to simulate the fracture-filling gas hydrate sediments. We set the incident angle as 0° and 90°, which represent horizontal and vertical cracks to a vertical borehole, respectively (Lee and Collett 2009). They can generally represent two extreme conditions. Moreover, gas hydrate can act as a part of pore fluid and solid frame component, which are also two scenarios for pore-filling gas hydrate model. Then, we compute how the \(Z = Vp^{0.5} /\rho^{5}\) of different morphologies of gas hydrate bearing sediments varied with saturation (Fig. 10). For different morphologies, the values of Z have the same trend with saturation. The gas hydrate can be identified from other minerals bearing sediments.

Fig. 10
figure 10

\(Z = Vp^{0.5} /\rho^{5}\) of different morphologies of hydrate bearing sediments as a function of saturation. The primary porosity of sediments are a 50%, b 60% and c 70%, respectively

In theory, this approach can also calculate gas hydrate saturation. The values of \(V_{p}^{0.5} /\rho^{5}\) increase with increasing saturation. The mineral components of study areas are required in order to compute the theoretical values of \(V_{p}^{0.5} /\rho^{5}\) varied with saturation. Then the saturation can be obtained by comparing the actual values with theoretical values of \(V_{p}^{0.5} /\rho^{5}\). Unfortunately, we presently don’t have the accurate mineral components of study area. Moreover, the saturation accuracy depends on velocity and density data in actual application.

Conclusion

In order to accurately identify gas hydrate in unconsolidated marine sediments, we research a series of attributes of different minerals-filled sediments and finally propose a new attribute \(V_{p}^{0.5} /\rho^{5}\) based on the effective medium theory. The \(V_{p}^{0.5} /\rho^{5}\) of gas hydrate-filled sediments increases with saturation, while the \(V_{p}^{0.5} /\rho^{5}\) of other minerals-filled sediments decrease with saturation. Thus, we use this difference to identify gas hydrate by comparing the computed \(V_{p}^{0.5} /\rho^{5}\) with threshold values. Well logging data at GMGS2 Site 05, 08 and ODP Leg 204 Site 1244E, 1245E are used to verify our method. The test results are generally agree with the actual conditions. However, few unmatched points are caused by unobvious velocity and density anomalies in low hydrate saturation position or background line fitting errors.

This proposed approach provides a new perspective to identify gas hydrate in unconsolidated marine sediments, which aims to improve identification accuracy. If velocity and density data can be derived by inversion in seismic data, the gas hydrate would be accurately detected by previous methods combined with this new approach.