Quantitative seismic pre-stack analysis of potential gas-hydrate resources in the Makran Accretionary Prism, offshore Iran

Marine and Petroleum Geology 48 (2013) 160e170

Contents lists available at ScienceDirect

Marine and Petroleum Geology journal homepage: www.elsevier.com/locate/marpetgeo

Quantitative seismic pre-stack analysis of potential gas-hydrate resources in the Makran Accretionary Prism, offshore Iran Ehsan Salehi a, Abdolrahim Javaherian a, b, *, Majid Ataee Pour c, Nasser Keshavarz Farajkhah d, Mojtaba Seddigh Arabani e a

Department of Petroleum Engineering, Amirkabir University of Technology, Tehran, Islamic Republic of Iran Institute of Geophysics, University of Tehran, Tehran, Islamic Republic of Iran Department of Mining and Metallurgical Engineering, Amirkabir University of Technology, Tehran, Islamic Republic of Iran d Exploration and Production Research Center, Research Institute of Petroleum Industry, Tehran, Islamic Republic of Iran e Department of Geophysics, Exploration Directorate, National Iranian Oil Company, Tehran, Islamic Republic of Iran b c

a r t i c l e i n f o

a b s t r a c t

Article history: Received 2 March 2013 Received in revised form 22 July 2013 Accepted 30 July 2013 Available online 11 August 2013

Gas hydrates are classified as a major unconventional resource and assumed to be a future potential energy source. High pressures, low temperature stability conditions of the hydrates restrict their occurrence to the permafrost and deep sea regions (from the continental slop to the abyssal). In these regions where the well information is most often absent or sparse, the amplitude-variation-with-angle (AVA) is a common method of hydrate/free gas assessment. In this method, various linear reflection coefficient equations are fitted to a series of data points to compute some pre-stack attributes. These attributes have the advantage of being an interface property rather than a layer property which makes them an easy to use characterization tool in lack of well data. In recent years, rock physics templates (nomograms) have been presented to estimate the hydrate/free gas saturations solely based on these interface properties. AVA inversion approach was applied to the seismic data in the Makran Accretionary Prism, offshore Iran, to estimate the bottom simulating reflector (BSR)-vicinity hydrate and gas saturations. This inversion process included four pre-stack attributes of the intercept, gradient, normal incidence (NI) and Poisson reflectivity (PR) and was limited to the BSR. Estimated BSR-vicinity saturations were ranged from 7%e8.5% for hydrate and 5%e6.5% for free gas at three representative locations. Accuracy of each attribute was analyzed by a synthetic modeling. In this modeling, the effective medium theory (EMT) was considered to calculate the physical properties of hydrate and gas-bearing sediments. Lithological properties and the acquisition geometry were set to the corresponding values of the field data. Results of the synthetic AVA inversion showed that if the hydrate saturation does not exceed 30%, the error of the estimation for pre-stack attributes would be in acceptable range. The expected estimation errors for actual attributes were derived from this synthetic modeling. Ó 2013 Elsevier Ltd. All rights reserved.

Keywords: Gas hydrate AVA inversion Pre-stack attribute Effective medium theory BSR Makran accretionary prism

1. Introduction Natural gas hydrates are crystalline solids in which the gas molecules (mainly methane) are trapped in a lattice of water molecules. The environmental and economical importance of these resources have been remarked by many authors, but their role as a potential future energy source attracts more attention (Collett, 2000, p. 123e136). As a result, the evaluation of gas hydrate

* Corresponding author. Department of Petroleum Engineering, Amirkabir University of Technology, Tehran, Islamic Republic of Iran. Tel.: þ98 21 64545131; fax: þ98 21 64543528. E-mail addresses: [email protected], [email protected] (A. Javaherian). 0264-8172/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.marpetgeo.2013.07.015

resources has grown significantly during the past two decades (Kvenvolden, 1993; Collett, 2008). Many resourceful techniques are developed to quantitatively assess these resources to estimate the amount of hydrates in marine sediments and the possible trapped free gas beneath them (Xia et al., 2000; Jin et al., 2002; Zillmer, 2006; Riedel and Shankar, 2012). Occurrence and stability of the hydrates require thermodynamic conditions present in the permafrost and deep sea regions (Kvenvolden and Barnard, 1983; Ecker and Lumley, 2001). In the latter case, well information is frequently absent or sparse (Pevzner et al., 2008; Shelander et al., 2010). In such a situation, seismic data is the main tool for the evaluation of these resources. Also a predrill estimation of hydrate saturation and distribution is sometimes required to locate the position of a well (Shelander et al.,

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

2012). Seismic processing velocity (stacking or migration) is a widely used quantity which many authors have tried to relate it to the physical properties (porosity, hydrate saturation etc.) by means of proper rock physic models (Ecker et al., 1997; Reister, 2003; Dewangan and Ramprasad, 2007; Schlesinger et al., 2012; Li et al., 2013). Besides the fact that the processing velocity is smooth and sometimes unrealistic, it does not consider the seismic amplitudes. The method of amplitude-variation-with-angle (AVA) eventuate the interface properties rather than the layer properties. Thus it can be linked to the physical properties only upon seismic data. AVA pre-stack attributes get more insight into the variations of elastic properties, especially across the bottom simulating reflector (BSR). For example, they might be used to delineate the potential hydrate reservoir (Fohrmann and Pecher, 2012) or to obtain information about Poisson’s ratio variation (Tinivella and Accaino, 2000) with no need to the well data. Methods of AVA inversion are mainly based on the linear approximation of Zoeppritz equations (1919) each of which accurate for a particular acoustic/elastic analysis. Different arrangement of the approximations may lead to different responses. Therefore, one should know about the possible sources of uncertainties. Recently, some authors have presented templates (nomograms) to estimate the hydrate and free gas saturations from AVA attributes (Cordon et al., 2006; Ojha et al., 2010). Nevertheless, computed attributes and estimated saturations are uncertain. Uncertainties are caused not only by the presence of noise in data (Cambois, 1998; Simm et al., 2000), but also, by inappropriate data processing and/or geological related impacts such as interfering and converted waves (Simmons and Backus, 1994, 1996; Bakke and Ursin, 1998). The accuracy of these attributes is also affected if the AVA assumptions are violated. To measure the accuracy of each attribute, a conceptual synthetic model, based on the effective medium theory (EMT) model was constructed (Helgerud et al., 1999). In this model, the hydrate crystals were considered as part of the load bearing matrix without any cementation effect which increase both P-wave and S-wave velocities. Elastic modulus, lithology and porosity values were taken from the reported values (Helgerud, 2001; Lee and Collett, 2001; Ojha and Sain, 2008). Different hydrate and gas saturations were considered which each pair form a different scenario. The acquisition geometry was also set similar to the corresponding field data. Using an elastic perturbational forward model, the seismic responses were derived and pre-stack attributes were calculated. The actual and estimated values for each attribute were graphically represented and the estimation error was measured by means of goodness-of-fit statistics. For the field data, in the next step, the calculation of AVA attributes was applied to three representative common reflection points (CRPs) and based on the EMT model the BSR-vicinity hydrate and free gas saturations were estimated. Upon the synthetic derived uncertainties, the amounts of expected errors for predicted saturations were noted. The presence of gas hydrate in the Makran Accretionary Prism has been investigated and proved by seismic pre- and post-stack analyses (Minshull et al., 1992; Grevemeyer et al., 2000; Kaul et al., 2000). In this study, BSR-vicinity hydrate/gas saturations in the Iranian sector of the Makran Accretionary Prism were estimated. For this purpose, several pre-stack angle-related attributes based on two different AVA approximations (of Zoeppritz equations) were calculated.

161

Figure 1. NortheSouth oriented seismic stack section in the Makran Accretionary Prism. The study of area is shown in the upper left part. The local uplift is caused by a trust fault which belongs to the continuum of an imbricate fault zone.

Figure 1 shows a NortheSouth oriented seismic line between CRPs 10,800 and 13,000. The BSR is conspicuous reflection in this section in about 2.4e2.6 s two-way time (TWT). This reflection is an indicator of bottom of the gas-hydrate stability zone (BGHSZ) which mimics the sea floor reflection with a reverse polarity. The local uplift is caused by a trust fault. The trust faults are dominant events at the accretionary wedge that form imbricate fault zone. There is no well data within the area, hence any hydrate analysis or characterization study is solely based on the seismic data. In such a situation, the AVA plays a major role in a qualitative and quantitative hydrate study (Muller et al., 2007). In this study, the AVA inversion was applied to the synthetic and field data. The estimation error of pre-stack attributes for synthetic data was assessed and the accuracy of different AVA approximations was determined via goodness-of-fit statistics. These results were also used to determine the inaccuracy of the estimated hydrate/gas saturations. The method of AVA inversion, the utilized goodness-of-fit statistics, the process of recovering true amplitude (for field data) and the method of saturation estimation are all represented in the next section.

2.1. AVA theory The AVA analysis also conventionally known as amplitudevariation-with-offset (AVO) has been widely utilized as an effective reservoir characterization tool (Ostrander, 1984; Fatti et al., 1994). This method is well described by the well-known Zoeppritz equations which completely determine the amplitudes of the reflected and transmitted plane waves at a single interface. In order to get a more intuitive understanding of the factors that control amplitude-variations-with angle/offset and also to simplify the computations, a series of linear approximations to the Zoeppritz equations have been developed. The Zoeppritz equations and their succeeding approximations establish a method to derive some useful pre-stack attributes via fitting the curves to the specific amplitudes. The Aki and Richards’ three-term approximation (1980, p. 153) provides a basis for the AVA inversion:


 
€ 1 2 1  4gsin2 q Dr=r þ 1=2cos2 q DvP =vP RzU 
 4gsin2 q DvS =vS ; =

2. Methodology The occurrence of gas hydrates in the Iranian sector of the Makran Accretionary Prism has been demonstrated through a detailed seismic attribute analysis (Hosseini Shoar et al., 2009).

(1)

where R is the reflection coefficient of P-wave, r is the average of density across an interface, vS is the average of S-wave velocity, vP is

162

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

the average of P-wave velocity, q is the average of the P-wave incidence and transmission angles and g ¼ v2S /v2P. D denotes the changes, assumed to be small, for the corresponding variables across the interface. Here the two following approximations to the Equation (1) and their corresponding pre-stack attributes are considered for an AVA inversion. Shuey (1985) proposed an alternative approximation of the Zoeppritz equations, as reported in Castagna and Smith (1994). This equation is totally equivalent to the Aki and Richards’s approximation (1980) but it is separated in such way that each term contributes in a different angular range hence giving a more tangible insight into the angle/amplitude relation:

RzI þ Gsin2 q þ Csin2 qtan2 q;

(2)

where,

I ¼ 1=2½DvP =vP þ Dr=r;

In the current study, a synthetic model is considered to investigate the effect of different hydrate and gas saturations on the prestack attributes of Equations (2) and (3). As it will be discussed later in Section 3.1, each hydrate/gas saturation pair makes a scenario. To quantitatively determine how much the estimated values are analogous to the actual ones, the goodness-of-fit statistics are computed for all scenarios. These statistics are root mean square error (RMSE) and R-square. RMSE is the average deviation of the predicted values from the response (actual) values. This statistic is considered as the standard error of the regression. In the case of synthetic result analysis, the individual difference between the actual and estimated attributes of each scenario is called the residual and RMSE serves to accumulate them into a single descriptive quantity as follows (Longley et al., 2005):

RMSE ¼

G ¼ 1=2½DvP =vP   4g½DvS =vS   2g½Dr=r; C ¼ 1=2½DvP =vP : The first term is the intercept which measures the reflection amplitudes at normal incidence, called acoustic impedance reflectivity; it is equivalent to the conventional definition of the seismic reflection coefficient (Badley, 1985, p. 7). The second and third terms are the gradient and the curvature and predominate at the intermediate and far angles, respectively. Gradient measures the variation rate of P-wave reflection with increasing of the incident angle. Gradient in conjugation with the intercept forms the basis of the AVO cross plot analysis (Castagna et al., 1998). Curvature is identical to the P-wave velocity reflectivity and is the major factor that controls the behavior of the reflection variation at large incident angles. Verm and Hilterman approximation (1995); consider the Poisson’s ratio as the major factor that controls the reflection behavior with angle variations; their approximation has the advantage of being based on the Poisson reflectivity. This approximation is originated from Shuey’s equation as follows:

i h i h RzNI 1  4gsin2 q þ PR sin2 q ;

2.2. Goodness-of-fit statistics

(3)

where NI is the normal incidence and is identical to the intercept. PR denotes the Poisson reflectivity and is defined as PR ¼ Ds/ (1  s)2, where s is the average Poisson’s ratio across the interface. The presence of outliers or non-Gaussian distributed errors in the amplitude profile may cause erroneous least square regression results. Therefore, the least square regression is considered to be non-robust, even only a few data samples are contaminated with the outliers (Walden, 1991). To cope with the impact of the outliers, some robust norms have been introduced into the seismic inverse problems (Bube and Langan, 1997; Guitton and Symes, 2003). This robust regression works by assigning a weight to each offset (or angle) within a CMP. This is achieved with an iteratively reweighted least square (IRLS) process (Andersen, 2008). For the initial iteration, all data points have an equal weight and therefore, the regression coefficients (pre-stack attributes) are computed with the ordinary least square method. For the next iterations, the points with greater errors in the previous iteration are assigned a lower weight and the coefficients are recomputed. This process is continued until the maximum number of iterations reaches or the regression coefficients vary within a pre-assigned interval.




1=n

X


b Y Y i i

2 0:5

;

(4)

^ i and Yi, in this case, are the actual and corresponding where, Y estimated attributes, respectively and n is the number of the data points which in this study is equal to the number of scenarios. It is usual to present this quantity in a normalized form, NRMSE. In this case, Equation (4) is divided by the range of actual values, (Ymax  Ymin). Another measure of goodness-of-fit, that is R-square, quantifies the amount of the data variation which is successfully accounted for by the estimation:

R-square ¼ 1 


X


b Y Y i i

2 .
X

 ðYi  Ymean Þ2 ;

(5)

where Ymean is the data mean. R-square most often varies between 0 and 1, with a value closer to the 1 indicating that the estimation is well and it explains the variation in the data.

2.3. True amplitude recovery A quantitative AVA inversion requires a true amplitude recovery. This goal was achieved by applying an amplitude recovery sequence to the raw amplitudes as follows: Geometrical spreading; in a layered medium, the signal’s amplitude decays due to the distribution of energy on the growing surface of the wave front. An equation by Ursin (1990) was used to compensate for offset-dependent geometrical spreading as a function of zero-offset travel time and the root-mean-square (rms) velocity. Attenuation; when a signal propagates into a medium, it losses part of its energy due to the energy mode conversion. This energy attenuation is a function of physical property of the medium and the frequency of the signal (Priest et al., 2006). Knopoff (1964) described the effect of this intrinsic attenuation on the signal’s amplitude after a traveled distance x:

AðxÞ ¼ A0 expðpx=Q lÞ

(6)

where A0 is the initial displacement amplitude, A(x) is the amplitude at the distance x, l is the dominant wavelength and Q is the quality factor. Q is often described in terms of the dissipation factor, Q1 (Lavergne, 1989, p. 18), which is the wave energy loss per cycle (O’Connell and Budanisky, 1978). To compensate for the attenuation effect of the BSR overlying layer, here Q ¼ 191 (measured by Sain and Singh, 2011) was considered for the P-wave.

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

Array directivity; array directivity filters the signal as a function of emergence angle on the surface. The filtering action varies with time and offset, so that it has an adverse impact on the AVA character of a reflector (Mazzotti and Ravagnan, 1995). The array effect correction, F, is accomplished using Sheriff and Geldart equation (1995, p. 247):

F ¼ sin½npðx=lÞsin q=fnsin½pðx=lÞsin qg;

(7)

where, n is the number of the receivers spaced within an interval x. q and l are incident angle and apparent wavelength, respectively. No source directivity correction is required if the source is considered as a point source (Shankar et al., 2005). Amplitude to reflection conversion; Warner approach (1990) is the most applied technique to convert the reflection amplitudes to the reflection coefficients. In the case of BSR, the reflection coefficient, RBSR, is retrieved as follows:

RBSR ¼ ABSR Am =A2p ;

(8)

where ABSR is the amplitude of the BSR and Ap and Am are the amplitudes of the sea floor primary reflection and its first order multiple, respectively. Ap and Am must be picked at the nearest offset where the water depth is relatively large and the sea floor is nearly flat. One assumption for Equation (8) is that the data are corrected for the spherical divergence and/or attenuation. If it is not the case, the right side of Equation (8) must be multiplied by 2 to account for the effect of the spherical divergence (Pevzner et al., 2008). Transmission loss; to correct the transmission losses due to the propagation across the sea floor, the amplitudes of all reflections underneath it are scaled by 1/(1  Rsf 2) (Fink and Spence, 1999), where Rsf 2 is the sea floor reflection coefficient. 2.4. Quantification of hydrate/gas saturations After the AVA inversion, two pre-stack attributes (for each equation) were calculated at each CRP that was utilized to estimate BSR-vicinity hydrate/gas saturations. The process of estimation included: (1) assuming a set of input physical properties such as lithology, porosity, contact number, (2) computing the elastic properties for a range of hydrate/gas saturations using EMT rock physics model, (3) calculating the desired AVA attributes (e.g. intercept and gradient) from the elastic parameters and finally (4) finding the nearest rock physics derived attributes to the AVA inverted attributes at each location. Corresponding saturations of the nearest pair were considered to be the BSR-vicinity hydrate/gas saturations for a specific CRP.

163

Table 1 The EMT modeling parameters including the elastic properties of different rock constituents (Helgerud, 2001). Modeling parameters

Values

Bulk modulus of quartz Shear modulus of quartz Density of quartz Bulk modulus of clay Shear modulus of clay Density of clay Bulk modulus of calcite Shear modulus of calcite Density of calcite Bulk modulus of pure hydrate Shear modulus of pure hydrate Density of pure hydrate Bulk modulus of methane Density of methane Bulk modulus of seawater Density of seawater Acceleration due to gravity Critical porosity Number of grain contacts

36.6 (GPa) 45 (GPa) 2.65 (g/cm3) 20.9 (GPa) 6.85 (GPa) 2.58 (g/cm3) 76.8 (GPa) 32 (GPa) 2.71 (g/cm3) 8.7 (GPa) 3.5 (GPa) 0.92 (g/cm3) 0.11 (GPa) 0.23 (g/cm3) 2.5 (GPa) 1.03 (g/cm3) 9.81 (m/s2)a 0.38a 9a

a

Lee and Collett (2001).

calculated using the HertzeMindlin contact theory (Mindlin, 1949). Below and above the critical porosity, these effective moduli were calculated via the upper and lower Hashin and Shetrikman bounds (1963). Finally, the saturated bulk modulus (Ksaturated) was calculated from the well-known Gassmann’s equation (1951) and saturated shear modulus (Gsaturated) was considered as equal to the dry frame shear modulus. The hydrate is then treated as an additional matrix component. When using HashineShetrikman bounds to calculate the dry moduli, this extra component stiffens the frame. It also reduces the porosity and hence decreases the weight of the pore fluid term on the bulk saturated moduli, when using Gassmann’s equation. Once the saturated bulk moduli are known, one can calculate the elastic velocities of rock. The required properties of the sediments as an input to the EMT modeling are listed in Table 1. The mineralogical composition was derived from a surface geological survey on the nearby area. After analyzing the specimens, it was surmised that the gas hydrate host sediments contained 70% quartz, 15% clay and 15% calcite. Based on this mineralogy and considering a hydrate saturation of 0%e90%, the P-wave velocities of sediments for a porosity of 40% and porosities of 40%  10% were calculated. The results are shown in Figure 2. To pick a nominated value for the porosity, these results

3. Examples 3.1. Synthetic model To build a set of proper synthetic seismograms, a series of plausible hydrate/free-gas models should be considered. Helgerud et al. (1999) introduced the EMT model to link the physical properties of a hydrate saturated rock to its elastic moduli and hence to the P-wave velocity, S-wave velocity and density. This theory has been used by many authors to model the AVA responses of the BGHSZ and gas or water saturated sediments (e.g. Matsumoto et al., 2011; Shelander et al., 2012). In this theory, the hydrate is considered as a load bearing part of a matrix with no cementing effect to reduce the initial porosity of the host sediments. At the critical porosity (according to Nur et al. (1998) this porosity varies from 0.36 to 0.43) the effective dry bulk modulus (KHM) and the shear modulus (GHM) were

Figure 2. Effect of the hydrate saturation on the P-wave velocity of a hydrate saturated sediment. Velocities were computed with the EMT model for three distinct porosities. Vc was set to 39% and hydrate saturations varied in a 0%e90% range. For this model, 70% quartz, 15% calcite and 15% clay were considered. Other modeling parameters are presented in Table 1.

164

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

were compared to the P-wave velocity and the hydrate saturation values of Ojha and Sain (2008). Comparisons showed that the middle curve in Figure 2 properly matched to their reported values (2100 m/s for 13% and 2160 m/s for 20% hydrate saturation). They also considered a porosity of 39% from Athy’s law (1930) for the BSR which was close to the nominated porosity of 40%. To build the corresponding synthetic seismograms, the simple, horizontally stratified model in Figure 3 was considered. As it is shown, a thick water column lays on three underlying layers with different pore contents. All elastic velocities and the density of sediments were calculated using the EMT. The hydrate saturation varies as 2%, 5%, 10%, 20%, 40% and 70% while the gas saturation varies as 0%, 3%, 6%, 10% and 25%. Each hydrate/gas saturation pair made a different scenario. Therefore, 30 different scenarios were considered for synthetic generation. Figure 3b shows an instance of 20% hydrate and 10% gas saturations. The depth of the BGHSZ was computed using a thermodynamic analysis. For this analysis, the geothermal gradient was taken 20  C/km. The sea floor temperature was set to 4  C and the gas hydrate was assumed to be entirely comprised methane. Using an elastic forward model algorithm, the above model was subjected to the computation of the interface’s reflection and the transmission coefficients. Here the applied forward algorithm was the wave-equation-based approach of Kennett (1985) while a zerophase Ricker wavelet with a 30 Hz dominant frequency had been selected for the source pulse. The events beyond the critical angles and the multiples were excluded for this modeling. Also, no interfering was taken into account for the calculated AVA responses. The acquisition geometry was set similar to the field data (Table 2). Although an AVA study of real data requires the amplitude correction process such as spreading correction, attenuation compensation or array directivity modification (Section 2.3), for a numerical modeling, it is supposed that all amplitudes have been recovered well and no further processing steps are required (Ruan et al., 2006). Figure 4a shows a synthetic CMP for the specific model of 20% hydrate saturation and 10% gas saturation. The BSR is a distinctive

Table 2 Seismic data acquisition parameters. Receiver interval Fold Total number of channels Number of streamer Streamer depth Source type/volume Shot interval Gun depth Recording length Sample rate Recording format

12.5 m 120 480 1 with length of 6100 m 6m Bolt air gun/2842 cu. inch 25 m 5m 7s 2 ms SEGD

reflector on this record. This is caused by an impedance contrast between the high velocity hydrate layer and the underlying low velocity gas layer. Preparing the data for the AVA analysis, the normal move out (NMO) correction and the NMO stretch mute were applied to the CMPs (Fig. 4b). Figure 4c shows the reflection variation of the BSR with the incident angle. These points were fitted to the linear Equations (2) and (3) to derive the pre-stack attributes. Also, the AVA responses of the bottom of the gas reflectivity (BGR) are plotted in this figure. This reflector nearly mirrors the shape of the BSR. In this study, the inversion is limited to the BSR reflectivity (incident angles of up to 25 ) to examine the effect of different hydrate and gas saturations on various pre-stack attributes. Figures 5e8 illustrate the actual and inverted attributes for all scenarios. The goodness-of-fit statistics are listed in Table 3. 3.2. Field data Four pre-stack attributes (intercept, gradient, NI and PR, respectively) were computed from three representative CRPs with a prominent BSR from the Makran Accretionary Prism. Locations of these CRPs are shown in the stacked section of Figure 1. The CRP numbers decrease landward with a 6.25 m distance. The acquisition parameters are listed in Table 2.

Figure 3. (a) A schematic stratified earth model used for the synthetic modeling. Depth of the bottom of the gas hydrate stability zone was determined by a thermodynamic modeling. All layers had the same lithology. Other modeling parameters are presented in Table 1. (b), (c) and (d) are the P-wave velocity, S-wave velocity and density for the 20% hydrate and 10% gas saturation scenario, respectively. All elastic parameters were computed using the EMT model.

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

165

The raw seismic data were processed in order to maintain accurate amplitudes. A preliminary velocity analysis by picking the rms velocity in the semblance window was done to get a compendious insight into the velocity structure. The geological complexity of this area and consequently the presence of the conflicting dips caused that a reflected event (e.g. BSR) not actually belong to a reflection point. To overcome the bad effects of dipping events with different stacking velocities, a Kirchhoff pre-stack time migration was applied to the data. Therefore, the data were transformed from the CMP domain to the CRP domain. Based on the represented sequence of Section 2.3, the true amplitudes were also recovered. All CRP traces are transformed from the offset to the angle domain using the following relation (Walden, 1991):

q ¼ sin1 ðxvP =tvrms Þ;

(9)

where x is the source-receiver offset, t is the total offset travel time, vP is the interval velocity of a particular layer and vrms is the rms velocity. Now the true reflectivity has been recovered, the CRPs have been transformed into the angle domain and four attributes of intercept, gradient, NI and PR have been computed for three representative CRPs (Table 4). In accordance with the synthetic data, AVA inversion of field data was also limited to incident angles of 25 (based on assumptions and limitations of different AVA approximations categorized in Li et al., 2007). The AVA inversion results for these CRPs will be discussed in the next section.

4. Results and discussion 4.1. Synthetic model

Figure 4. (a) The synthetic CMP gather based on the model of Figure 3. (b) The CMP of part (a), after NMO correction and NMO stretch mute. (c) Variation of the reflection coefficients with the angle of incidence. Negative and positive coefficient events are related to the hydrate/gas and the gas/brine reflections, respectively.

The synthetic CMP of the each scenario was inverted using Equations (2) and (3) to compute the discussed pre-stack attributes. The results for each attribute are shown in the Figures 5e8. These figures illustrate the actual and estimated attributes alongside their graphical comparison.

Figure 5. Intercept of the BSR calculated by the AVA inversion based on the Shuey’s approximation. (a) Actual intercept values for different scenarios computed by the EMT model. The color-coded surface is a linear interpolation among 30 distinctive scenario points. (b) Estimated intercept values using the AVA inversion, corresponding to the values in (a). (c) A graphical comparison of the actual and the estimated intercept values. Points of this figure were color coded based on the hydrate saturation. The dashed line is the best fit line. It is the line where the points will be on placed if the inversion exactly estimates the actual values. The solid line is the linear least square fit between the actual and estimated values and shows where the estimation deviates from the best fit line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

166

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

Figure 6. The gradient of the BSR calculated by AVA inversion based on the Shuey’s approximation. (a) Actual gradient values computed by the rock physics (EMT) model. (b) Estimated gradient values using AVA inversion, corresponding to the values in (a). (c) Graphical comparison of the actual and the estimated gradient values. Points enclosed by the circle are related to a 70% hydrate saturation which yields a positive gradient.

Figure 5 shows the intercept, calculated from Shuey’s approximation (Equation (2)). Figure 5a depicts that increasing the hydrate and free gas saturations causes an increase in the intercept value along the BSR (it becomes more negative). Appearance of hydrates in rocks increases while the gas existence drops the P-wave velocity, so that any increase in hydrate and/or gas leads to the enhancement of the acoustic impedance contrast at the BSR. From Figure 5a, it is obvious that the intercept is more sensitive to the gas saturation rather than the hydrate saturation, so that a little amount of gas results in a large negative intercept. It must be noticed that the reverse polarity of BSR (Fig. 1) is mainly due to this attribute and brighter amplitude of this reflector is attributed to the

higher amount of gas saturation. One indication of an increase in gas saturation along a BSR is the simultaneous increasing of the intercept and gradient (Yi et al., 2011). Figure 5b shows the estimated intercept values. From a qualitative point of view, the estimated values have an analogous trend compared to the actual ones (Fig. 5a), especially for lower amounts of the hydrate and gas saturations. A quantitative comparison of these two series is presented in Figure 5c. If the inversion process calculated the perfect estimation, all points would be along the best fit line (the dashed line). The solid line is a least-squares-fit between the actual and estimated points and shows a deviation from the perfect estimation. Points of this

Figure 7. NI attribute calculated by the AVA inversion based on the Verm and Hilterman’s approximation, across the BSR. (a) Actual NI values for different scenarios computed by the rock physics (EMT) model. (b) Estimated NI values using the AVA inversion, corresponding to the values in (a). (c) Graphical comparison of the actual and estimated attributes. Points of this figure were color coded based on the hydrate saturation. Points of this figure were color coded based on the hydrate saturation. The dashed line is the best fit line and the solid line is the linear least square fit between the actual and the estimated values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

167

Figure 8. PR attribute calculated by the AVA inversion based on the Verm and Hilterman’s approximation, across the BSR. (a) Actual PR values. (b) Estimated PR values using the AVA inversion, corresponding to the values in (a). (c) Graphical comparison of the actual and the estimated values. Points of this figure were color coded based on the hydrate saturation. The dashed line is the best fit line and the solid line is the linear least square fit between the actual and the estimated values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

figure were color coded based on the hydrate saturation. The RMSE and R-square for this attribute are 0.036 and 0.87, respectively. The measured goodness-of-fit statistics are listed in Table 3. A high value of R-square indicates that almost all of the data variation is accounted for by the estimation. The greatest deviation from the best fit line is occurred at a high saturation of hydrate. In this situation, the Aki and Richards’s assumption of a small property change across the interface is violated (Equation (1)). For example, considering the extreme scenario of 70% hydrate saturation and 25% gas saturation, the P-wave velocity varies from 3090 m/s to 1430 m/s across the BSR. If the rare case of 70% hydrate saturation is excluded, then the intercept estimation becomes less erroneous. The second attribute of the Equation (2) is the gradient. Figures 6a and 6b show the actual and inverted values of this attribute, respectively. This attribute predominates for middle angles. When the hydrate saturation is less than about 40%, the gradient values are negative (Fig. 6a). As the intercept values are also negative, these cause the class III AVO response based on the Rutherford and Williams (1989) classification. However, for hydrate saturations more than 50%, the gas hydrate crystals occupy more pore space and cause a significant S-wave velocity increase. In this case, the gradient values become positive. This leads to a decrease in reflection coefficients with an increase in the incident angle and results in the class IV AVO response (Castagna et al., 1998). For the hydrate saturation between these two ranges, the gradient is around zero and reflection coefficients are nearly steady within the middle angle range. Table 3 Goodness-of-fit statistics for the computed pre-stack attributes of synthetic modeling. Reflectivity equation

Pre-stack attribute

RMSE

NRMSE

R-square

Shuey approximation (1985)

Intercept Gradient NI PR

0.036 0.043 0.039 0.033

0.096 0.093 0.104 0.127

0.87 0.9 0.83 0.8

Verm and Hilterman approximation (1995)

A comparison of the actual and estimated gradient values shows that the both series have similar trends. The quantitative comparison of the gradient values in Figure 6c highlights that the estimated values match the real ones except for the case of a high hydrate or gas saturations. The depicted ellipse in this figure encloses the scenario of 70% hydrate saturation where the gradient is positive. The RMSE and R-square for this attribute were 0.043 and 0.9, respectively. The inversion results based on the Verm and Hilterman approximation (Equation (3)) are presented in Figures 7 and 8. Actual and estimated NI values are depicted in parts (a) and (b) of Figure 7. NI is equivalent to the intercept as mentioned before. In this regard, the results of Figure 7 are comparable to the ones in Figure 5. For this attribute, RMSE and R-square were 0.078 and 0.83, respectively. From Figure 8a, for constant gas saturation, PR gradually decreases with hydrate increasing. RMSE and R-square of PR estimation were 0.065 and 0.8, respectively. PR is a gas sensitive attribute (Verm and Hilterman, 1995) so that its variation is much larger when the gas is introduced to the pore space. Muller et al. (2007) showed that a relative increment in vS due to the hydrate saturation is larger than that in vP and the Poisson ratio decreases with hydrate increasing. Hence the largest PR corresponds to a case of the highest gas saturation while the hydrate saturation is of a lower amount. Similarly, Tinivella and Accaino (2000) indicated that Poisson’s ratio of hydrated sediment slightly decreases with increasing the hydrate saturation up to about 50%. After that, the Poisson’s ratio drops more rapidly. They related this behavior to the

Table 4 Pre-stack attributes values for three representative CRPs. Estimated amounts of hydrate and gas saturations are also listed. CRP number

Intercept

Gradient

NI

PR

Hydrate saturation

Gas saturation

12,315 12,355 12,482

0.1075 0.105 0.126

0.13 0.121 0.128

0.108 0.112 0.126

0.22 0.21 0.26

7%  2% 7%  2% 8.5%  2.4%

5%  1.3% 4.9%  1.3% 6.5%  1.5%

168

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

cementing effect of the hydrates at saturations higher than 50%. Therefore, a band-limited inversion of the PR section (in the absence of well data) yields a Poisson’s ratio section in which the hydrate distribution type can be identified. Based on the AVA inversion results represented in Figures 5e8 and statistics in Table 3, one can conclude that for this simple synthetic model, the inversion process overly computes the acceptable results. Higher saturations of hydrate and free gas cause greater errors. The possibility of greater errors in higher saturations is a violation of the assumption that the property change across a reflector is minor. Except the permafrost region, occurrence of sediments with hydrate saturations more than 40% is rare. If the higher hydrate saturation scenarios are excluded from the goodness-of-fit measurements, the amount of error will be decreased. For example, ignoring 40% and 70% hydrate saturations, improves the fit between the actual and estimated intercept values by 14%. Therefore, considering the common cases (hydrate saturations below 40%); the attribute estimation will be more precise. Table 3 also shows that using the intercept value to estimate normal incident reflectivity is more accurate than using NI. It is because that any probable error in vS/vP (in NI formula Equation (3)) does propagate into the estimation of NI and it causes a greater error of estimation for normal incident reflectivity.

4.2. Field data Scattered points in Figure 9 are the BSR reflectivities of three representative CRPs. Actually these points are super-gathers of five adjacent CRPs where the central CRP is marked in Figure 1. In Figure 9, two estimated BSR reflectivities, computed from

Figure 9. The BSR reflection variation with the angle of incidence in three different CRP locations. Circles are the actual BSR reflection values. Solid and dashed lines are the approximate reflection coefficients computed with Equations (2) and (3), respectively.

intercept-gradient and NIePR, are shown by solid and dashed lines, respectively. The computation of NI and PR requires an appropriate vS/vP estimation. Here a vS/vP value relevant to the hydrate and gas saturations is used instead of considering vS/vP ¼ 0.5. For this, estimated saturations from the intercept and gradient were related to the corresponding vS and vP using the EMT model. The results of AVA inversion for three CRPs are listed in Table 4. The estimation of hydrate/gas saturations was based on the process in Section 2.4, where porosity was derived from Athy’s law (1930) and lithology was set to 70% quartz, 15% clay and 15% calcite (as in the synthetic modeling). Here the BSR shows the class III AVA responses (with a small gradient). Although the BSR with a class IV AVA responses has been reported in the Makran Accretionary Prism (Ojha and Sain, 2007), the AVA responses in Figure 9 are the most common character of the BSR in this area. The gradient term in Equation (2) is a combination of the P-wave velocity reflectivity, S-wave velocity reflectivity and the reflectivity of density. Frequently, the reflectivity of density is smaller compared to the other two variables. This implies that when the contrasts of vP and vS across an interface have opposite signs, the reflectivity becomes larger with an increase in the angle. Considering the matrix support model for a hydrate reservoir, the vS contrast is negative (the same as the vP contrast) and smaller than that of the vP. In this situation, the gradient is still negative but its value is relatively small and the BSR reflectivity increases slightly with the incident angle. In the above statement, the lithology on both sides of the BSR is assumed to be identical. This is in contrast with a hydrate cementation model (Dvorkin and Nur, 1993) where a significant increase in S-wave velocity due to the hydrate presence causes a positive gradient and hence a class IV AVA responses. Another class IV possibility is the case of a high hydrate saturation (e.g. Fig. 6) which is rare. Based on the evidence, it is deduced that for this region the hydrate distribution has a matrix support form. Also, a relatively low saturation of hydrates causes a minor BSR reflectivity variation with the incident angle. The intercept values for CRPs 12,315 and 12,355 were almost equal and around 0.11. Moving to CRP 12482, this value slightly increases. Although the gas origin in the Makran Accretionary Prism is not well understood but one possibility of the higher intercept would be the lateral or vertical gas migration across the local high. Based on the synthetic results, the amount of error for the estimated intercept was minor and around 3.1%. Also, the gradient values of these CRPs did not vary much and ranged between 0.128 and 0.13. The error of the gradient estimation based on the synthetic results was 15% which was much larger than the intercept error. Using the computed values of the intercept and gradient, the amount of the BSR-vicinity hydrate and gas saturations were estimated (Table 4). These amounts were consistent with the observed BSR AVA responses. It must be noticed that the errors due to the gradient uncertainty were much larger than the intercept uncertainty. Also, the variation of hydrate saturation had a minor effect on these attributes (in accordance with Muller et al., 2007). For example, the hydrate variation from 5% to 10% altered the intercept by 8% and the gradient by 4%. This implied that the estimated hydrate saturation was more uncertain. The computed NI values for three CRPs were close to the intercept values listed in Table 4. For this NI range (based on the synthetic modeling results) the error was 4.2%. As discussed before, the greater error is correlated to the inappropriate vS/vP ratio in Equation (3). Also, the corresponding error for PR (z0.2) was 10%. Like gradient, PR also had a little dependency on the hydrate saturation so that varying the hydrate from 5% to 10% altered the PR value by only 1%. But this attribute had a lower estimation error and gave a physical/elastic insight to the hydrate resource.

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

5. Conclusions AVA intercept-gradient template is a gas hydrate characterization tool which can be utilized to investigate the presence of the free gas under the BSR and also used to directly estimate the hydrate and gas saturations. Other pre-stack attributes such as PR are frequently computed for hydrated sediments to gain a more physical insight into these resources. But the accuracy of these attributes is affected not only by non-attenuated noise in the data but also by violated assumptions in the AVA inversion. The results of the AVA inversion for the EMT-based synthetic modeling showed that if the hydrate saturation does not exceed 30%, then the inversion process can compute the pre-stack attributes quite well. Also, it was shown that the estimated intercept values had a smaller error than the NI values despite the fact that both of them had the same elastic meaning. Pre-stack attribute calculation for the three representative CRPs in the Makran Accretionary Prism indicated a BSR-vicinity hydrate saturation of 7%e8.5% and a gas saturation of 5%e6.5%. The averaged error of estimation were 2.1% and 1.4% of hydrate and gas saturations, respectively. The impact of hydrate saturation variations on the computed attributes is small and therefore its estimation has a higher uncertainty with respect to the gas saturation. The AVA pre-stack attributes establish an easy to use method of BSR characterization and hydrate and gas resource assessment.

Acknowledgments The authors wish to express their gratitude to the department of Geophysics of NIOC Exploration Directorate, Iran, for making available the seismic data and permission to publish this paper. We would like to thank to Research Institute of Petroleum Industry for helping with the data preparation.

References Aki, K., Richards, P.G., 1980. Quantitative Seismology, Theory and Methods. W. H. Freeman and Company. Andersen, R., 2008. Modern Methods for Robust Regression. Sage Publications Inc. Athy, L.F., 1930. Density, porosity and compaction of sedimentary rocks. AAPG Bulletin 14, 1e24. Badley, M.E., 1985. Practical Seismic Interpretation. International Human Resources Development Corporation. Bakke, N.E., Ursin, B., 1998. Thin-bed AVO effects. Geophysical Prospecting 46, 571e 587. Bube, K., Langan, R., 1997. Hybrid l1/l2 minimization with applications to tomography. Geophysics 62, 1183e1195. Cambois, G., 1998. AVO attributes and noise: pitfalls of crossplotting. In: 68th Annual International Meeting, SEG, Expanded Abstracts, pp. 244e247. Castagna, J.P., Smith, S.W., 1994. Comparison of AVO indicators: a modeling study. Geophysics 59, 1849e1855. Castagna, J.P., Swan, H.W., Foster, D.J., 1998. Framework for AVO gradient and intercept interpretation. Geophysics 63, 948e956. Collett, T.S., 2000. Natural gas hydrate as a potential energy resource. In: Max, M.D. (Ed.), Natural Gas Hydrate in Oceanic and Permafrost Environments. Kluwer Academic Publishers. Collett, T.S., 2008. Geology of marine gas hydrates and their global distribution. In: Offshore Technology Conference, OTC 19241. Cordon, I., Dvorkin, J., Mavko, G., 2006. Seismic reflections of gas hydrate from perturbational forward modeling. Geophysics 71, F165eF171. Dewangan, P., Ramprasad, T., 2007. Velocity and AVO analysis for the investigation of gas hydrate along a profile in the western continental margin of India. Marine Geophysics Research 28, 201e211. Dvorkin, J., Nur, A., 1993. Rock Physics for Characterization of Gas Hydrates: the Future of Energy Gases. U.S. Geological Survey Professional Paper, pp. 293e298. Ecker, C., Dvorkin, J., Nur, A., 1997. Estimating the Amount of Gas Hydrate and Free Gas from Surface Seismic. Stanford Exploration Project, Report 95, pp. 173e195. Ecker, C., Lumley, D., 2001. Seismic AVO Analysis of Methane Hydrate Structure. Stanford Exploration Project, Report 80, pp. 1e15. Fatti, J.L., Vail, P.J., Smith, G.C., Strauss, P.J., Levitt, P.R., 1994. Detection of gas in sandstone reservoirs using AVO analysis: a 3-D seismic case history using the geostack technique. Geophysics 59, 1362e1376.

169

Fink, C.R., Spence, G.D., 1999. Hydrate distribution off Vancouver Island from multifrequency single-channel seismic reflection data. Journal of Geophysical Research 104, 2909e2922. Fohrmann, M., Pecher, I.A., 2012. Analyzing sand-dominated channel systems for potential gas-hydrate-reservoirs using an AVO seismic inversion technique on the Southern Hikurangi Margin, New Zealand. Marine and Petroleum Geology 38, 19e34. Gassmann, F., 1951. Elasticity of porous media. Vieteljahrsschrift der Naturforschenden Gessellschaft 96, 1e23 (in German). Grevemeyer, I., Rosenberger, A., Villinger, H., 2000. Natural gas hydrates on the continental slope off Pakistan: constraints from seismic techniques. Geophysics Journal International 140, 295e310. Guitton, A., Symes, W., 2003. Robust inversion of seismic data using the Huber norms. Geophysics 68, 1310e1319. Hashin, Z., Shetrikman, S., 1963. A variational approach to the theory of elastic behavior of multiphase materials. Journal of the Mechanics and Physics of Solids 11, 127e140. Helgerud, M.B., 2001. Wave Speeds in Gas Hydrate and Sediments Containing Gas Hydrate: a Laboratory and Modeling Study (PhD thesis). Stanford University. Helgerud, M.B., Dvorkin, J., Nur, A., Sakai, A., Collett, T.S., 1999. Elastic wave velocity in marine sediments with gas hydrates: effective medium modeling. Geophysical Research Letters 26, 2021e2024. Hosseini Shoar, B., Sedigh Arabani, M., Javaherian, A., 2009. Application of seismic attributes for identification of gas hydrate bearing zone and free gas beneath it. In: First EAGE International Petroleum Conference and Exhibition, Iran, Shiraz, 4e6 May. Jin, Y.K., Lee, M.W., Collett, T.S., 2002. Relationship of gas hydrate concentration to porosity and reflection amplitude in a research well, Mackenzie Delta, Canada. Marine and Petroleum Geology 19, 407e415. Kaul, N., Rosenberger, A., Villinger, H., 2000. Comparison of measured and BSRderived heat flow values, Makran accretionary prism, Pakistan. Marine Geology 164, 37e51. Kennett, B.L.N., 1985. Seismic Wave Propagation in Stratified Media. Cambridge University Press. Knopoff, L., 1964. Q. Reviews of Geophysics 2, 625e660. Kvenvolden, K.A., 1993. A Primer on Gas Hydrate: the Future of Energy Gases. U. S. Geological Survey Professional Paper, pp. 279e291. Kvenvolden, K.A., Barnard, L.A., 1983. Hydrates of natural gas in continental margins. In: Watkins, J.S., Drake, C.L. (Eds.), Studies in Continental Margin Geology, vol. 34. Association of Petroleum Geology Memoir, pp. 631e640. Lavergne, M., 1989. Seismic Methods. Editions Technip. Lee, M.W., Collett, T.S., 2001. Elastic properties of gas hydrate-bearing sediments. Geophysics 66, 763e771. Li, Y., Downton, J., Xu, Y., 2007. Practical aspects of AVO modeling. The Leading Edge 26, 295e311. Li, L., Lei, X., Zhang, X., Sha, Z., 2013. Gas hydrate and associated free gas in the Dongsha Area of northern South China Sea. Marine and Petroleum Geology 39, 92e101. Longley, P.A., Goodchild, M.F., Maguire, D.J., Rhind, D.W., 2005. Geographic Information Systems and Science. John Wiley & Sons. Matsumoto, R., Ryu, B.J., Lee, S.R., Lin, S., Wu, S., Sain, K., Pecher, I., Riedel, M., 2011. Occurrence and exploration of gas hydrate in the marginal seas and continental margin of the Asia and Oceania region. Marine and Petroleum Geology 28, 1751e1767. Mazzotti, A., Ravagnan, G., 1995. Impact of processing on the amplitude versus offset response of a marine seismic data set. Geophysical Prospecting 43, 263e 281. Mindlin, R.D., 1949. Compliance of elastic bodies in contact. Journal of Applied Mechanics 16, 159e268. Minshull, T.A., White, R.S., Barton, P.J., Collier, J.S., 1992. Deformation at plate boundaries around the Gulf of Oman. Marine Geology 104, 265e277. Muller, C., Bonnemann, C., Neben, S., 2007. AVO study of a gas hydrate deposit, offshore Costa Rica. Geophysical Prospecting 55, 1e17. Nur, A., Mavko, G., Dvorkin, J., Galmudi, D., 1998. Critical porosity: the key to relating physical properties to porosity in rocks. The Leading Edge 17, 357e362. O’Connell, R.J., Budanisky, B., 1978. Measure of dissipation in viscoelastic media. Geophysics Research Letter 5, 5e8. Ojha, M., Sain, K., 2007. Seismic velocities and quantification of gas-hydrates from AVA modeling in the western continental margin of India. Marine Geophysical Research 28, 101e107. Ojha, M., Sain, K., 2008. Appraisal of gas-hydrate/free-gas from vP/vS ratio in the Makran accretionary prism. Marine and Petroleum Geology 25, 637e644. Ojha, M., Sain, K., Minshull, T.A., 2010. Assessment of gas-hydrate saturations in the Makran accretionary prism using the offset dependence of seismic amplitudes. Geophysics 75, C1eC6. Ostrander, W.J., 1984. Plane-wave reflection coefficients for gas sands at nonnormal angles of incidence. Geophysics 49, 1637e1648. Pevzner, R.L., Volkonskaya, A.L., Bouriak, S.V., Bocharova, A.A., Blinova, N.V., 2008. Prediction of properties of gas hydrate bearing sediments from the dynamic of reflected waves without borehole data. In: 70th EAGE Conference and Exhibition, Italy, Rome, 9e12 June. Priest, J.A., Best, A.I., Clayton, C.R.I., 2006. Attenuation of seismic waves in methane gas hydrate-bearing sand. Geophysics Journal International 164, 149e159. Reister, D.B., 2003. Using measured velocity to estimate gas hydrates concentration. Geophysics 68, 884e891.

170

E. Salehi et al. / Marine and Petroleum Geology 48 (2013) 160e170

Riedel, M., Shankar, U., 2012. Combining impedance inversion and seismic similarity for robust gas hydrate concentration assessments e a case study from the KrishnaeGodavari basin, East Coast of India. Marine and Petroleum Geology 36, 35e49. Ruan, A.G., Li, J.B., Chu, F.Y., Li, X.Y., 2006. AVO numerical simulation of gas hydrate reflectors beneath seafloor. Chinese Journal of Geophysics 49, 1665e1675. Rutherford, S.R., Williams, R.H., 1989. Amplitude-versus-offset variations in gas sands. Geophysics 54, 680e688. Sain, K., Singh, A.K., 2011. Seismic quality factors across a bottom simulating reflector in the Makran Accretionary Prism, Arabian Sea. Marine and Petroleum Geology 28, 1838e1843. Schlesinger, A., Cullen, J., Spence, G., Hyndman, R., Louden, K., Mosher, D., 2012. Seismic velocities on the Nova Scotian margin to estimate gas hydrate and free gas concentrations. Marine and Petroleum Geology 35, 105e115. Shankar, U., Sinha, B., Thakur, N.K., Khanna, R., 2005. Amplitude-versus-offset modeling of the bottom simulating reflection associated with submarine gas hydrates. Marine Geophysics Research 25, 29e35. Shelander, D., Dai, J., Bunge, G., 2010. Predicting saturation of gas hydrates using pre-stack seismic data, Gulf of Mexico. Marine Geophysical Research 31, 39e57. Shelander, D., Dai, J., Bunge, G., Singh, S., Eissa, M., Fisher, K., 2012. Estimating saturation of gas hydrates using conventional 3D seismic data, Gulf of Mexico Joint Industry Project Leg II. Marine and Petroleum Geology 34, 96e110. Sheriff, R.E., Geldart, L.P., 1995. Exploration Seismology. Cambridge University Press. Shuey, R.T., 1985. A simplification of Zoeppritz equations. Geophysics 50, 609e614. Simm, R., White, R., Uden, R., 2000. The anatomy of AVO crossplots. The Leading Edge 19, 150e155.

Simmons, J.L., Backus, M.M., 1994. AVO modeling and the locally converted shear wave. Geophysics 59, 1237e1248. Simmons, J.L., Backus, M.M., 1996. Waveform-based AVO inversion and AVO prediction-error. Geophysics 61, 1575e1588. Tinivella, U., Accaino, F., 2000. Compressional velocity structure and Poisson’s ratio in the marine sediments with gas hydrate and free gas by inversion of reflected and refracted seismic data (South Shetland Islands Antarctica). Marine Geology 164, 13e27. Ursin, B., 1990. Offset-dependent geometrical spreading in a layered medium. Geophysics 55, 492e496. Verm, R., Hilterman, F., 1995. Lithology color-coded seismic sections: the calibration of AVO crossplotting to rock properties. The Leading Edge 14, 847e853. Walden, A.T., 1991. Making AVO sections more robust. Geophysical Prospecting 39, 915e942. Warner, M., 1990. Absolute reflection coefficients from deep seismic reflections. Tectonophysics 173, 15e23. Xia, G., Sen, M.K., Stoffa, P.L., 2000. Mapping of elastic properties of gas hydrates in the Carolina trough by waveform inversion. Geophysics 65, 735e744. Yi, B.Y., Lee, G.H., Horozal, S., Yoo, D.G., Ryu, B.J., Kang, N.K., Lee, S.R., Kim, H.J., 2011. Qualitative assessment of gas hydrate and gas concentrations from the AVO characteristics of the BSR in the Ulleung Basin, East Sea (Japan Sea). Marine and Petroleum Geology 28, 1953e1966. Zillmer, M., 2006. A method for determining gas-hydrate or free-gas saturation of porous media from seismic measurements. Geophysics 71, N21eN32. Zoeppritz, K., 1919. Erdbebenwellen VIII B, on the Reflection and Penetration of Seismic Waves through Unstable Layers. Goettinger Nachr, pp. 66e84 (in German).