Seismic Data Processing and the Characterization of a Gas Hydrate Bearing Zone Offshore of Southwestern Taiwan

1 CAS Key Laboratory of Marginal Sea Geology, South China Sea Institute of Oceanology, Chinese Academy of Sciences, Guangzhou, China 2 Guangzhou Institute of Geochemistry, Chinese Academy of Sciences, Guangzhou, China 3 Guangzhou Center for Gas Hydrate Research, Chinese Academy of Sciences, Guangzhou, China 4 Graduate School of Chinese Academy of Sciences, Beijing, China 5 Sinopec International Petroleum Exploration and Production Corporation, Beijing, China * Corresponding author address: Dr. Hui Deng, CAS Key Laboratory of Marginal Sea Geology, South China Sea Institute of Oceanology, Chinese Academy of Sciences, Guangzhou, China; E-mail: nhsdenghui@163.com Various seismic attributes of gas hydrate bearing sediments were analyzed in the accretionary prism offshore of southwestern Taiwan utilizing seismic imaging, velocity analysis, AVO analysis, and AVO inversion of large offset seismic data. A bottom-simulating reflector (BSR) is clearly observed on the seismic section with a reversed polarity compared to that of the seabed reflection. Instantaneous amplitude sectioning clearly shows lateral variations of the BSR. The zero-phase waveform of the BSR is distinct and the weak reflectors above the BSR can be observed on the instantaneous phase section. AVO analysis shows the absolute value of the negative BSR amplitude increasing with offset. A low P-wave interval velocity layer was found below the strong BSR by detailed velocity analysis. Both the P (normal incident P-wave reflection coefficient) and G (AVO gradient) values are highly negative for the strong BSR on the P and G sections, and they lie in the third quadrant of the P and G cross-plot section. The P+G (reflectivity of Poisson’s ratio) value is also negative on the P+G section and the P-G (normal incident S-wave reflection coefficient) value is approximately zero on the P-G section along the same strong BSR. All the seismic characters described above suggest that a gas hydrate layer exists together Terr. Atmos. Ocean. Sci., Vol. 17, No. 4, December 2006 782 with a free gas layer below it along strong BSRs in the area offshore southwestern Taiwan. (

Various seismic attributes of gas hydrate bearing sediments were analyzed in the accretionary prism offshore of southwestern Taiwan utilizing seismic imaging, velocity analysis, AVO analysis, and AVO inversion of large offset seismic data.A bottom-simulating reflector (BSR) is clearly observed on the seismic section with a reversed polarity compared to that of the seabed reflection.Instantaneous amplitude sectioning clearly shows lateral variations of the BSR.The zero-phase waveform of the BSR is distinct and the weak reflectors above the BSR can be observed on the instantaneous phase section.AVO analysis shows the absolute value of the negative BSR amplitude increasing with offset.A low P-wave interval velocity layer was found below the strong BSR by detailed velocity analysis.Both the P (normal incident P-wave reflection coefficient) and G (AVO gradient) values are highly negative for the strong BSR on the P and G sections, and they lie in the third quadrant of the P and G cross-plot section.The P+G (reflectivity of Poisson's ratio) value is also negative on the P+G section and the P-G (normal incident S-wave reflection coefficient) value is approximately zero on the P-G section along the same strong BSR.All the seismic characters described above suggest that a gas hydrate layer exists together

INTRODUCTION
Rifting and sea-floor spreading have led to the formation of the Cenozoic sedimentary basins and petroleum accumulation in the northern margin of the South China Sea.The area offshore of southwestern Taiwan is the arc-continent collision zone between the Southern China continental margin and the Taiwan-Luzon island arc.Abundant sediments were transported from Taiwan and were convolved into an accretionary prism by the subduction-collision processes.The accretionary prism is bounded by the Manila Trench to the west and by the North Luzon Trough to the east.Complex imbricated structure zones are formed within the accretionary prism (Huang and Yuan 2002).
Comprehensive geological and geophysical surveys have been carried out for academic research offshore of southwestern Taiwan (Huang and Yin 1990;Huang et al. 1992;Reed et al. 1992;Liu et al. 1997;Huang et al. 2001).The seismic proxy of gas hydrate, i.e., the bottomsimulating reflector (BSR), has been recognized and classified in association with the tectonosedimentation geology.The accretionary prism offshore of southwestern Taiwan has been suggested to be favorable and have bright prospects for the formation of gas hydrate (Chi et al. 1998;Shyu et al. 1998;McDonnell et al. 2000;Chi et al. 2006;Liu et al. 2006).
In 2001, seismic reflection data were acquired in the area offshore of southwestern Taiwan (line L, Fig. 1) using a long multi-channel streamer.This survey was supported by the Chinese Basic Research Priorities Program.BSR is clearly observed in the accretionary prism on the migrated seismic profile (Fig. 2).Song et al. (2004) and Deng et al. (2005) found a BSR in this dataset after conventional seismic data processing and analyzing.The seismic data of line L have a maximum source-receiver offset of 3237.5 m, and present 30-fold with high signal-tonoise ratio.This makes it possible to conduct extensive velocity analyses, seismic imaging and AVO (Amplitude-Versus-Offset) processing.In this paper, we describe the data acquisition first, and show different processing and analyses performed for various seismic attributes of BSR, and finally, the characters of the gas hydrate bearing zone are presented.

SEISMIC DATA AND ACQUISITION
The seismic data were collected in 2001 by the seismic vessel Tanbao of the Guangzhou Marine Geological Survey.The seismic line runs across the accretionary prism offshore of southwestern Taiwan (line L, Fig. 1).The source is composed of two arrays of Sleeve guns with a total capacity of 3000 in 3 at pressure of 2000 psi.The shot interval is 50 m.The number of channels is 240 with interval of 12.5 m.The near offset is 250 m and the maximum sourcereceiver distance is 3237.5 m. 10 second of data is recorded with a sampling rate of 2 ms.The air guns and the streamer are deployed at 7 and 9 m in depths, respectively.The seafloor relief varies sharply and the water depth is from 900 to 3500 m across line L (Figs. 1, 2).AB section is relatively shallow and smooth.Along this section, the seismic data have a high signal-to-noise ratio and the frequency bandwidth is from 10 to 70 Hz.The BSR is evident and very clear (Fig. 3).The seismic data in the AB section thus are suitable for highresolution seismic imaging, accurate velocity analyses and AVO analyses.

SEISMIC IMAGING AND VELOCITY ANALYSES
High resolution seismic imaging processing was conducted for the seismic character study of gas hydrate.The keys to the seismic processing are seismic amplitude preservation and detailed velocity analyses.The NMO (normal move-out) corrected CDP (common depth point) gathers and velocity models were also obtained for the AVO inversion.
The main processing sequences include: geometry definition, trace editing, spherical divergence correction, velocity analysis, NMO, stacking, and migration.The stacked and migrated profiles are shown in Figs. 4 and 5, respectively.Complex trace analysis was used to the migrated data to produce the instantaneous amplitude and instantaneous phase profiles (Figs. 6, 7).Seismic data were processed using the Focus processing system.
The depths of the BSR vary from 0.4 to 0.5 s below the sea floor.As the sea floor multiple appears 2 s below the sea floor, the BSR is not interfered by multiples.Exponential amplitude compensation was applied in the seismic amplitude preservation processing.Many processing procedures, including filtering, DMO (dip move-out), sea floor multiple suppression, deconvolution and amplitude balance, were not used in this study in order to avoid artificial influence on the waveform and amplitude of the BSR.
To obtain accurate seismic velocity information, seismic velocity analysis was carried out using semblance and velocity panel techniques (Sheriff and Geldart 1995).The RMS (rootmean-square) velocity was picked at every 20th CDP and AGC (automatic gain control) was applied before velocity analysis so as to help identifying and picking velocity values of the   weak reflectors (Fig. 8).The RMS velocities were converted into the P-wave interval velocities using Dix's equation (Dix 1955) (Fig. 9).
BSR depth varies from 1800 to 2100 m over the section AB and it is far smaller than the maximum source-receiver distance (3237.5 m).With high signal-to-noise ratio, the seabed and BSR of the section AB are continuous and strong reflectors (Fig. 8); therefore, the RMS velocity errors of these reflectors are small.However, the reflectors beneath the BSR are discontinuous and weak, and their RMS velocity errors are relative large.
The picking error of the RMS velocity and the thickness of the reflectors are the main factors that affect the P-wave interval velocity accuracy.The picking error of the RMS velocity (the minimum RMS velocity increment picked on the computer screen during interactive velocity analysis) is about ± 3 m s -1 .The thickness of the reflectors below the BSR varies from 50 to 100 m.According to Hajnal and Sereda (1981), the maximum P-wave interval velocity error below the BSR is about ± 50 m s -1 , which is far smaller than the P-wave interval velocity variations across the BSR (Fig. 9).

AVO ANALYSES
AVO analysis has been widely used to delineate physical properties around a reflection boundary and to analyze whether free gas is trapped beneath the hydrated layer (Hyndman and Spence 1992;Katzman et al. 1994;Andreassen et al. 1995;Ecker et al. 1998;Carcione and Tinivella 2000;Tinivella and Accaino 2000;Diaconescu et al. 2001).Generally, accurate AVO analysis of a BSR requires correction for directivity of both the air guns array and the hydrophones array (after spherical divergence correction).No directivity correction of the air guns array was applied to the data, because the source we used (two Sleeve guns distributed cross-line) can be treated as a point source.Only directivity correction of the hydrophones array was made.
The hydrophones array is linear and the elements are spaced at equal intervals along the streamer.The hydrophones array attenuation (F) was described by Sheriff and Geldart (1995): where n is the number of hydrophones (in our case, 14) separated by ∆x = 0.9615 m; λ is the wavelength; and θ is the incident angle. (2 where Vp is the P-wave interval velocity of the medium above the interface; V rms is the RMS velocity; and t(x) is the picked time at x offset.Directivity correction to reflection amplitude is given by the following formula: where A i is the uncorrected amplitude, and A j is the corrected amplitude.The values of corrected amplitude for the BSR and the seabed are displayed in Figs.10a  and b, respectively, using equation (3).
In Fig. 10a, the maximum offset (3237.5 m) corresponds to an incident angle of about 43 degrees by equation (2), and the offset of 2400 m corresponds to an incident angle of 30 degrees, which is the critical angle for the AVO inversion in this paper.

AVO INVERSION
The variation of the P-wave amplitude with offset can be used as a direct hydrocarbon indicator (Ostrander 1984).The Zoeppritz's equations (Zoeppritz 1919) describe theoretically the reflection amplitude changes with the incident angle.BSR is generally interpreted as the base of the gas hydrate stability zone (BGHS) (Shipley et al. 1979), which separates the sediments containing gas hydrate above and the sediments containing free gas below.The elastic property contrasts between the gas hydrate bearing and the free gas bearing sediments lead to BSR amplitude variations with incident angles.
The Zoeppritz equations were approximated ( θ 30 ≤ o ) by Shuey (1985) as: R( ) where R( ) θ is the reflection coefficient; θ is the incident angle; P is the P-wave reflection coefficient at normal incidence; and G is the gradient.
The AVO inversion was performed based on equation (4).The key processing sequences are: data input, pre-conditioning, ray tracing, AVO quality control, angle gathers generation, and AVO attributes inversion.The input data were NMO-corrected CDP gathers and seismic velocity model data.The offsets were converted into the incident angles by ray tracing, and the CDP gathers beyond 30 degree of the incident angle were excluded before AVO inversion.The quality control assured the reliability of AVO inversion.Inversion was done by the Probe system, and the resulting P, G, P + G, and P -G attributes were obtained (Figs. 11,12).Here P -G is the S-wave reflection coefficient at normal incidence, and P + G is the Poisson's ratio reflectivity (Verm and Hilterman 1995).

BSR
BSR appears to be a zero-phase reflection in shot gathers and CDP gathers, and its polarity is reversed compared to the seabed reflection.It was found that the BSR shows increasingly negative amplitude with offset in a CDP gather (Fig. 3b).
Better images of the sub-bottom shallow structure and the BSR were obtained by seismic image processing.BSR is recognized clearly in both the stacked and migrated sections.It is 0.38 to 0.58 s below seabed, and sub-parallels seabed.The depth of the BSR increases as the water   depth increases.The BSR cuts across dipping sedimentary reflections between CDP4600 and CDP5160.The weak reflection or blank reflection zone occurs between the seabed and the BSR, and the thickness of the weak reflection zone increases from west to east.The amplitude and continuity of the BSR vary laterally.The amplitude is stronger between CDP2240 and CDP3000 (Figs. 4, 5).
The amplitude of the BSR appears more clearly in the instantaneous amplitude section (Fig. 6).The weak reflectors between the seabed and the BSR can be observed in the instantaneous phase section (Fig. 7).
For the strong BSR between CDP2240 -3000 (see Figs. 4, 5, and 6), the amplitude varies nonlinearly with offset (Fig. 10a).The BSR amplitude increases very little or even decreases at near offset, but increases rapidly at middle and far offsets, especially beyond 2400 m (or 30 degree of the incident angle).In contrast, amplitude variations with offset are not distinct for weak BSR (CDP3000 -5160).

P-wave Interval Velocity
For strong BSR (CDP2240 -3000), the P-wave interval velocities (Vp) change from high to low across the BSR (Fig. 9).Vp increases gradually from about 1700 to 2000 m s -1 downward from the seabed, then sharply decreases to between 1500 and 1700 m s -1 below the BSR and forms a low velocity zone.The thickness of the low velocity zone is between 40 and 80 m.Vp then increases rapidly below the low velocity zone (Fig. 9).For the sections where weak BSRs are observed (CDP3000 -5160), no low velocity zone below the BSR was identified.

AVO Attributes
P-attribute section represents the P-wave zero-offset reflection coefficient, which reveals the variation of P-wave impedance.The BSR is negative in the P profile (Fig. 11a), indicating that the impedance below the BSR is smaller than that above the BSR.The G-attribute section represents the P-wave zero-offset reflection coefficient variations with incident angles (Fig. 11b).The G-attribute values along BSR in section AB are also negative, implying that the absolute value of the BSRs amplitude increases with incident angles.
Poisson's ratio reflectivities along the BSR between CDPs 2600 -3000 are also highly negative as the P + G attribute section reveals (Fig. 12a), implying that the values of the Poisson's ratio below the BSR is far smaller than that above the BSR.The S-wave reflection coefficient is close to zero in the P -G attribute section between CDPs 2600 -3000 (Fig. 12b), suggesting that the S-wave interval velocity contrast is very small across the BSR in this section.

DISCUSSION
Within the section AB where strong BSRs are observed between CDPs 2600 -3000, the absolute value of the BSR amplitude increases with offset (Fig. 10).If the amplitude variation is considered to be caused mainly by the rock physical property variations across the BSR, three types of BSR amplitude variation may be produced according to Hyndman et al. (2001): type I, only hydrated layer above the BSR; type II, both hydrate layer above the BSR and free gas layer below the BSR; and type III, only free gas layer below the BSR.However, in the section AB, the AVO anomaly of the BSR belongs to the type III (Rutherfort and Williams 1989) and the BSR's value lies in the third quadrant of the crossplot of the P and the G sections (Castagna and Swan 1997;Castagna et al. 1998) (Fig. 13c), thus the sediments should contain a hydrate layer above the BSR and free gas below the BSR.
In the section AB, the P-wave interval velocity is decreased across the BSR (Fig. 9), the Poisson's ratio of the sediments decreases sharply below the BSR as revealed in the P + G attribute section (Fig. 12a), and the S-wave reflection coefficient nearly equals zero in the P -G attribute section (Fig. 12b).This implies that the S-wave interval velocity contrast is small across the BSR.The seismic characters established support the presence of a gas hydrate layer above the BSR and a free gas layer below the BSR (Domenico 1976;Minshull and White 1989;Carcione and Tinivella 2000;Tinivella and Accaino 2000;Lu 2003).
As per integrated analyses above, the gas hydrated layer above the BSR and the free gas layer below the BSR may most probably occur in the place where strong BSRs are observed (e.g., CDP2240 -3000).Where weak BSRs are presented in section AB (CDP3000 -5160), the P-wave interval velocity reversal is not clearly observed, and the trend of amplitude variation with offset at BSR is not clear.In this case, only the hydrated layer may exist above the BSR.
According to Deng et al. (2005), the depth of the BSR calculated from the seismic data of line L is consistent with that derived based on the temperature-pressure curve of gas hydrate formation.Many gas hydrate surveys, such as ODP (ocean drilling program) Leg 164 and ODP Leg 204 (Wood and Rupple 2000;Trèhu et al. 2003), have proved the presence of the gas hydrate above BSR.On the migrated profile of line L (Fig. 2), the waveform, polarity, amplitude variation of the BSR, and the weak reflections above the BSR are clearly visible, which is comparable to those of BSRs in other convergent margins of the world, such as in the Peru margin and in the northern Cascadia margin, where gas hydrate samples have been collected (von Huene and Pecher 1999;Hyndman et al. 2001;Pecher et al. 2001).This further suggests that gas hydrate may be present in the accretionary prism offshore of southwestern Taiwan.

CONCLUSIONS
Seismic imaging techniques, detailed velocity analyses, and AVO inversion have been applied to a large offset multichannel seismic dataset that runs over the accretionary prism offshore of southwestern Taiwan.Various seismic attributes derived from these analyses indicate that both gas hydrate and free gas layers may be present in places where strong BSRs are observed, and only gas hydrate bearing zones may occur above weak BSRs in the shallow subbottom sediments of the accretionary prism of line L. Comprehensive analyses of types of seismic attributes in the study strengthen the reliability of identifying gas hydrate bearing zones in this area.

Fig. 1 .
Fig. 1.Location of the survey line in the area offshore southwestern Taiwan.

Fig. 2 .
Fig. 2. The migrated seismic profile of line L. The AB box shows the section where various analyses have been performed.

Fig. 7 .
Fig. 7.The instantaneous phase profile of section AB.

Fig. 9 .
Fig. 9. P-wave interval velocity structures at various CDP locations along section AB.

Fig. 10 .
Fig. 10.Amplitude vs offset display of the BSR (a) and the seabed (b) for section AB (CDPs 2445 -2905).The dots are corrected amplitude of BSR and seabed.The dashed lines are illustrative, which are used to delineate the BSR and seabed amplitude variations with offsets.

Fig. 11 .
Fig. 11.The P-attribute section (a) and the G-attribute section (b) for part of the section AB.

Fig. 12 .
Fig. 12.The P + G attribute section (a) and the P -G attribute section (b) for part of the section AB.

Fig. 13 .
Fig. 13.The P-attribute section (a), the G-attribute section (b), and their crossplot (c) for part of the section AB.The data of the crossplot are sampled from the rectangle area in the P and G sections above.