Estimation of Gas Hydrate Saturation Using Constrained Sparse Spike Inversion: Case Study from the Northern South China Sea

Bottom-simulating reflectors (BSRs) were observed beneath the seafloor in the northern continental margin of the South China Sea (SCS). Acoustic impedance profile was derived by Constrained Sparse Spike Inversion (CSSI) method to provide information on rock properties and to estimate gas hydrate or free gas saturations in the sediments where BSRs are present. In general, gas hydrate-bearing sediments have positive impedance anomalies and free gas-bearing sediments have negative impedance anomalies. Based on well log data and Archie’s equation, gas hydrate saturation can be estimated. But in regions where well log data is not available, a quantitative estimate of gas hydrate or free gas saturation is inferred by fitting the theoretical acoustic impedance to sediment impedance obtained by CSSI. Our study suggests that gas hydrate saturation in the Taixinan Basin is about 10 - 20% of the pore space, with the highest value of 50%, and free gas saturation below BSR is about 2 - 3% of the pore space, that can rise to 8 - 10% at a topographic high. The free gas is non-continuous and has low content in the southeastern slope of the Dongsha Islands. Moreover, BSR in the northern continental margin of the SCS is related to the presence of free gas. BSR is strong where free gas occurs.


Bottom-simulating reflectors (BSRs) were observed beneath the seafloor in the northern continental margin of the South China Sea (SCS). Acoustic impedance profile was derived by Constrained Sparse Spike Inversion (CSSI) method to provide information on rock properties and to estimate gas hydrate or free gas saturations in the sediments where BSRs are present.
In general, gas hydrate-bearing sediments have positive impedance anomalies and free gas-bearing sediments have negative impedance anomalies. Based on well log data and Archie's equation, gas hydrate saturation can be estimated. But in regions where well log data is not available, a quantitative estimate of gas hydrate or free gas saturation is inferred by fitting the theoretical acoustic impedance to sediment impedance obtained by CSSI. Our study suggests that gas hydrate saturation in the Taixinan Basin is about 10 -20% of the pore space, with the highest value of 50%, and free gas saturation below BSR is about 2 -3% of the pore space, that can rise to 8 -10% at a topographic high. The free gas is non-continuous and has low content in the southeastern slope of the Dongsha Islands. Moreover, BSR in the northern continental margin of the SCS is related to the presence of free gas. BSR is strong where free gas occurs.
(Key words: Gas hydrate, Free gas, Saturation, Acoustic impedance, South China Sea)

INTRODUCTION
The hydrate-bearing sediment has an anomalously high compressional velocity, and gas hydrate saturation can be estimated when a reference velocity curve is known (Tinivella 1999). On the other hand, a low compressional velocity zone may indicate the presence of free gas in the pore space. The strong acoustic impedance contrast between sediment containing gas hydrate and sediments with free gas is evident on marine seismic reflection profiles as bottom simulating reflectors (BSRs, e.g., Shipley et al. 1979;Hyndman and Spence 1992). The elastic properties of pure gas hydrate are well known, but little is known about the elastic properties of hydrate-bearing sediments. The Poisson's ratio and the Lamé parameter term λρ are estimated from the P-impedance and S-impedance. McMechan (2002, 2004) proposed that the hydrated sediments have high elastic impedance, high P-impedance, high S-impedance, high λρ , and slightly lower Poisson's ratio compared to those of the surrounding non-hydratebearing sediment. The free gas-bearing sediment has low elastic impedance, low P-impedance, non-anomalous background S-impedance, low λρ and low Poisson's ratio. The anomalously high acoustic impedance is related to gas hydrate saturation.
The South China Sea (SCS) lies in the juncture of Eurasia, Pacific and Indo-Australia plates. The tectonics and geological evolution have been investigated for its potentially rich hydrocarbon reserves and unique marginal sea basin setting (Lee and Lawver 1995;Lüdmann and Wong 1999;Hayes 1980, 1983). Recent studies have focused on gas hydrate and deep-water hydrocarbon exploration. Milkov (2000) showed that structures such as active fault systems, mud volcanoes and other tectonic pathways control gas migration and gas hydrate accumulation. The tectonic regime in the northern SCS is mainly controlled by faults trending in NE and NEE directions and BSRs are observed to the south of Taiwan on reflection seismic records and at the Xisha Trough (Chi et al. 1998(Chi et al. , 2006McDonnell et al. 2000;Wu et al. 2005;Liu et al. 2006). The structures of the SCS continental margin are characterized by diapirs, troughs, seamounts, seabed valleys, fan channels and slope slumps (Wang et al. 2005;Yang et al. 2006).
Acoustic impedance profiles can be obtained by Constrained Sparse Spike Inversion (CSSI), which provides a quantitative approach to determine the presence of gas hydrate and free gas and to estimate their saturations. The gas hydrate saturations in the Taixinan Basin and in the Pearl River Mouth basin west of the Dongsha Islands are inferred by using CSSI on two seismic profiles 0102 and 97301, respectively (Fig. 1).

DATA DESCRIPTION
The two multi-channel seismic (MCS) reflection profiles we used for this study (Fig. 1) were shot by the Guangzhou Marine Geological Survey. Seismic line 0102 was shot in July, 2002 by R/V "Fen Dou Wu Hao" using a 2400-m long streamer with 192 channels (trace interval 12.5 m) and a coherent airgun source with a total volume of 8 × 20 in 3 shooting at every 25 m. The sampling interval is 1 ms. Streamer depth is 6 m and source depth is 2 m. Seismic line 97301 was shot in April, 2001 by R/V "Tan Bao Hao" using a 3000-m long streamer with 240 channels (trace interval 12.5 m) and the shot interval being 50 m. The sample interval is 2 ms and the source depth is 8 m. The data were processed by the Institute of Processing, Guangzhou Marine Geological Survey. The processing sequences mainly include pre-filtering, true amplitude recovery, noise attenuation, deconvolution, velocity analysis (first), residual statics, velocity analysis (second), multiple attenuation, DMO velocity analysis, DMO stack and noise attenuation.

Constrained Sparse Spike Inversion
Acoustic impedance provides rock property information and has been used to discern rock types and as a direct hydrocarbon indicator (Latimer et al. 2000). CSSI can be used to derive acoustic impedance variation based on seismic data and low-frequency impedance information obtained from well logs (Lindseth 1979;Jason Geosciences Workbench 2003). If there is no well log data, stacking velocities can be used to derive acoustic impedance, though the reliability of the results depends on the quality of seismic data. In this study, VelMod package (Jason Geosciences Workbench 2003) was used to convert stacking velocities into interval velocities, impedance can also be obtained. Pseudo well logs were created by sonic or P-wave velocity. Inversion of seismic data for acoustic impedance is performed using CSSI package of Jason Geosciences Workbench (2003). The input data include time-migrated seismic data, the seismic wavelet, interpolated control horizons and interpolated low-frequency impedance trend. Then InverTrace package was used to derive band-pass P-impedance from the time-migrated seismic data, wavelet, pseudo wells and time gate information. The lowfrequency information of P-impedance is obtained from the low-frequency trends of the well log impedance, which can be interpolated and extrapolated from the control horizon to the whole study area. The band-passed P-impedance can be merged with low frequency impedance to obtain the acoustic impedance. CSSI finds acoustic impedance values within prescribed constraints by minimizing the objective function (Jason Geosciences Workbench 2003) : at time i , where r i is the reflection coefficient, d i is the recorded seismic acoustic trace, s i is the synthetic trace, p and q are empirically determined exponents, and λ is the data mismatch weighting factor. The factors p and q can be set from the user interface. The p factor is 0.9 by default. For p in the range between 0 and 1, the solution will be sparser; for p > 1, the solution is more band-limited. The q factor is 2 by default, as seismic noises are generally assumed to be Gaussian in nature (Jason Geosciences Workbench 2003). The factor λ is used to control the balancing of the misfit norms. A low λ emphasizes the reflection misfit and will result in acoustic impedance with a few sharp contrasts, little detail and high residuals. A high λ emphasizes the mismatch term between seismic data and synthetic trace data and will result in detailed acoustic impedance and low residuals. In CSSI, λ factor is selected to strike an optimum balance between generating acoustic impedance traces. During generating acoustic impedance data, the λ factor had better demand high signal-to-noise, high correlation misfit, low seismic misfit and low reflectivity misfit. In seismic profile 0102 inversion, we set λ = 24, p = 1, q = 2. In seismic profile 97301 inversion, we set λ = 18, p = 1, q = 2. P-wave velocity, porosity, resistivity, density and impedance information can be obtained from Site 1148 of ODP Leg 184 (Fig. 2). The acoustic impedance profiles can then be obtained respectively (Figs. 3, 4).

Water-Saturated Porosity
If we assume that the sediment without gas hydrate and free gas is under uniform   compaction, the relationship between acoustic impedance (I) and water-saturated porosity ( Φ) will be a smooth baseline curve along which Φ decreases and I increases with increasing depth (Lu and McMechan 2002). Φ will decrease in the presence of gas hydrate in the sediment while I will increase. However, if there is free gas in the sediment, both Φ and I will decrease. In a gas hydrate-bearing zone, the deviation of water-saturated Φ− I from the baseline curve indicates gas hydrate saturation, which is the same as the deviation of resistivity from R 0 baseline in Fig. 2b.
To obtain the relationship between water-filled porosity and acoustic impedance, we use well log data at Site 1148 of ODP Leg 184 in the northern SCS. The gas hydrate-bearing zone can be determined by the anomalies of resistivity, acoustic impedance and compressional velocity. We can obtain Φ− I curve by constrained least-squares fitting of the porosity and acoustic impedance from site 1148. Φ = 0.299I 3 + 1.752I 2 -33.745I + 131.159 .
(2) Φ is porosity of neutron log in %, I is acoustic impedance in 10 6 kg m -2 s -1 . The water-filled porosity profile of seismic profile 97301, which passes Site 1148, can be obtained by equation 2 (Fig. 4b). Since there is no well data available near seismic profile 0102, we used Hamilton's empirical formula (Hamilton 1976) to calculate the porosity of this profile, assuming that the background porosity of the whole seismic profile is the same (which might not be true with terrigenous sediment in some areas): Φ = 0.72 -0.8165z + 0.361z 2 . (3) Φ is porosity, z is depth below seafloor in km.

GAS HYDRATES OR FREE GAS SATURATION
Two different approaches are adapted to determine gas hydrate and free gas quantities in seismic profiles 0102 and 97301, respectively. Gas hydrate saturation of seismic line 97301, where log data are available, can be obtained by Archie's equation (1942). For seismic profile 0102, since there is no well log data available, we use acoustic impedance contrast between water-saturated and gas hydrate-bearing or free gas-bearing sediments to derive their saturations.

Estimation of Saturation from Impedance
The compressional wave velocity equation of Domenico (1977) was used to calculate water-saturated acoustic impedance (see explanation of symbols in Table 1). and It can be applied in the case of (1) full water saturation; (2) water and gas hydrates in the pore space; and (3) free gas and water in the pore space (Tinivella 1999). Gas hydrate-bearing sediment can be determined from its anomalous acoustic impedance value. If the local geological setting is simple, the discrepancies between the inverted P-impedance by CSSI Profile and the theoretical P-impedance for water-saturated sediment are interpreted to be due to the presence of gas hydrate (positive anomalies) or free gas (negative anomalies). By comparing sediment impedance obtained by CSSI and the theoretical acoustic impedance calculated by equations 4 and 5 for water-saturated sediment, we can determine gas hydrate and free gas zones. Saturations of gas hydrate and free gas are estimated by fitting the theoretical acoustic impedance data (equation 4) to the P-impedance data obtained by CSSI. Only if we give a certain gas hydrate or free gas saturation value, and the calculated impedance by equations 4 and 5 is close to the result from CSSI, we can say that the sediment is gas hydrate or free gas saturated. If they are not accordant, we can modify the saturation value until they are accordant with each other. The shear moduli of dry sediment without hydrate can be calculated using modified Hashin-Shtrikman-Hertz-Mindlin theory (Dvorkin and Nur 1996). The sediment rock consists of a mixed mineralogy, the bulk and shear moduli of the rock can be determined using the Hill's average formula (Dvorkin and Nur 1996). The coupling factor k in equation 5 describes the degree of coupling between pore fluid and solid frame. Numerically, it ranges from one to infinity. In this study, we assumed that gas hydrate was part of solid, so it reduced porosity of the sediments and affected the bulk and shear muduli. The components and physical properties are given in Table 2. Figure 3b shows the acoustic impedance distribution of seismic line 0102 obtained by CSSI. We applied equations 3, 4, 5 and used the physical property values in Table 2 to obtain saturations of gas hydrate (Fig. 3c) and free gas (Fig. 3d). Figure 3c shows that gas hydrate distribution is non-uniform and blocky, and the corresponding saturation is about 10 -30% of the pore space. High values of 50% of the pore space appear around the central part of the profile. The free gas distribution is layered and discontinuous. Free gas saturation is about 2 -3% of the pore space, and can be as high as 8 -10% of the pore space near the central position of the profile (local topographic high).

Estimation of Saturation from Resistivity
Like ice, gas hydrates are electrical insulators (Pearson et al. 1983). Gas hydrate-bearing sediments have high electrical resistivities in comparison with water-saturated formation. If we assume that resistivity anomalies are caused by the presence of gas hydrate in the sediments and the pores are filled with water and gas hydrate then gas hydrate saturation S h is Table 2. Components and physical properties of framework. Table 1. List of parameters in eq. (5) and material properties parameters.
given by Archie's equation (1942): where S w is water saturation, R 0 is the resistivity of the sediment saturated with water, R t is the measured formation resistivity, and n is empirical constant (Paull et al. 1996;Collett and Ladd 2000). For gas hydrate-bearing sediments, n = 1.9386 (Pearson et al. 1983). If we assume that the background resistivity R 0 of the whole seismic line is the same (which can be determined at Site 1148 of ODP Leg 184 by a least-squares linear fit with depth) i.e., where z is the depth below seafloor in meters, then gas hydrate saturation above BSR can be obtained from equation 6 and equation 7. The relationship between gas hydrate saturation S h , water-filled porosityφ , the resistivities of water-saturated formation ( R 0 ), and the pore water ( R w ) is given by: where R 0 , m, a, R w and n are all determined empirically for a particular environment; hence, the only unknown being φ .
If we use m = 2.56, a = 1.05 (Collett and Ladd 2000) and n = 1.9386 (Pearson et al. 1983), we can obtain the gas hydrate saturation section estimated from equations 2 and 8 (Fig. 4d). Distribution of gas hydrate along seismic line 97301 is layered and continuous, which is quite different from that of seismic line 0102 in Taixinan Basin. Gas hydrate saturation is about 10 -20% of the pore space; in general, it increases to about 30% of the pore space near BSR.

The Nature of BSRs
Analyses of seismic and well log data in ODP Leg 164 indicate that free gas exists beneath the BSRs (Ecker et al. 1998(Ecker et al. , 2000Yuan et al. 1996;McMechan 2002, 2004), yet some doubt remains that free gas is present in the underlying sediments of GHSZ (Dillon and Paull 1983;Hyndman and Spence 1992).
Free gas bearing sediments have low P-impedance (Lu and McMechan 2004). We inferred no free gas exists in the sediment along seismic line 97301 since P-impedance is high (Fig. 4b). Moreover, the BSR is very weak in the seismic profile. The structure of seismic line 0102 is complex and the strata below BSR in CDP 11700-8900 exhibits discontinuity (Fig. 5). We suggest that a mud diapir lies between CDP 10300-8900. There is a rock density inversion and fault occurrence, and possible gas accumulations in the deep subsurface as suggested from the free gas saturation profile (Fig. 3d), all this being favorable for the formation of mud diapir (Milkov 2000). On profile 0102, the BSR is strong where free gas is present. This supports the idea that BSR is related with the presence of free gas.

Accuracy and Limitations
The accuracy of saturation estimation is greatly affected by selected parameter values, especially porosity values. If the water-filled porosity used is higher than that of the actual porosity, gas hydrate saturation obtained from the given water-filled porosity will be lower than the actual values. While if the water-filled porosity used is lower than that of the actual porosity, gas hydrate saturation obtained from the given water-filled porosity will be higher than the actual values. Sea. Note the complex mud diapiric structure, active faults and BSR interpreted on this profile. See Fig. 1 for profile location.
Gas hydrate saturation in seismic line 97301 is based on Archie's equation and an indirect empirical estimation of water-filled porosity from acoustic impedance. The values of R w and R 0 in equation 8 were determined by ODP Leg 184 Site 1148 in the SCS, whereas n, m, and a were assumed to be the same as those for the Blake Ridge (Pearson et al. 1983), those assumptions could be the most significant source of error.
For the region, a lack of well log data, and the quality of the seismic data may greatly influence inversion accuracy. The saturation in seismic line 0102 is based on the discrepancy between acoustic impedance obtained from CSSI and the assumed water-filled acoustic impedance. The acoustic impedance obtained is constrained by stacking velocity in the absence of well logs. The porosity values used are from the empirical function of Hamilton (1976), which is not applicable everywhere. Porosity was used to estimate saturations and the more accordant with formation porosity, the higher the result accuracy. Moreover, components and physical properties (Table 2) in the northern SCS are not applicable everywhere, although the method is applicable to any other data.

CONCLUSIONS
Acoustic impedance profiling derived from CSSI allows us to identify the presence of gas hydrate where positive impedance anomalies occur and the presence of free gas where negative impedance anomalies are observed. The impedance anomalies permitted us to estimate gas hydrate and free gas saturation based on well log data.
Gas hydrate distribution at Site 1148 of ODP Leg 184 appears to be layered and continuous, and the saturation is about 10 -20% of the pore space. While gas hydrate distribution in the Taixinan Basin is patchy and non-continuous, and saturation is about 10 -30% of the pore space, with the highest value of 50%. Free gas is present in this area, with saturation of 2 -3% of the pore space. Below GHSZ, BSR is stronger in areas where free gas exits than where it does not.