International Conference on Gas Hydrates (ICGH) (6th : 2008)

SEISMIC MODELING OF HETEROGENEITY SCALES OF GAS HYDRATE RESERVOIRS Huang, Jun-Wei; Bellefleur, Gilles; Milkereit, Bernd Jul 31, 2008

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


59278-5666.pdf [ 714.07kB ]
JSON: 59278-1.0040992.json
JSON-LD: 59278-1.0040992-ld.json
RDF/XML (Pretty): 59278-1.0040992-rdf.xml
RDF/JSON: 59278-1.0040992-rdf.json
Turtle: 59278-1.0040992-turtle.txt
N-Triples: 59278-1.0040992-rdf-ntriples.txt
Original Record: 59278-1.0040992-source.json
Full Text

Full Text

Proceedings of the 6th International Conference on Gas Hydrates (ICGH 2008), Vancouver, British Columbia, CANADA, July 6-10, 2008.  Seismic modeling of heterogeneity scales of gas hydrate reservoirs Jun-Wei Huang∗ Department of Physics University of Toronto 60 St. George Street, Toronto, ON, M5S 1A7 CANADA Gilles Bellefleur Geological Survey of Canada 615 Booth Street, Ottawa, ON, K1A 0E9 CANADA Bernd Milkereit Department of Physics University of Toronto 60 St. George Street, Toronto, ON, M5S 1A7 CANADA ABSTRACT The presence of gas hydrates in permafrost regions has been confirmed by core samples recovered from the Mallik gas hydrate research wells located within Mackenzie Delta in the Northwest Territories of Canada. Strong vertical variations of compressional and shear velocities and weak surface seismic expressions of gas hydrates indicate that lithological heterogeneities control the lateral distribution of gas hydrates. Seismic scattering studies predict that typical horizontal scales and strong velocity contrasts due to gas hydrate concentration will generate strong forward scattering, leaving only weak energy to be captured by surface receivers. In order to understand the distribution of gas hydrates and the scattering effects on seismic waves, heterogeneous petrophysical reservoir models were constructed based on the P-wave and S-wave velocity logs. Random models with pre-determined heterogeneity scales can also be used to simulate permafrost interval as well as sediments without hydrates. Using the established relationship between hydrate concentration and P-wave velocity, we found that gas hydrate volume content can be determined by correlation length and Hurst number. Using the Hurst number obtained from Mallik 2L-38, and the correlation length estimated from acoustic impedance inversion, gas hydrate volume fraction in Mallik area was estimated to be 17%, approximately 7x108 m3 free gas stored in a hydrate bearing interval with 250,000 m2 lateral extension and 100 m depth. Simulations of seismic wave propagation in randomly heterogeneous models demonstrate energy loss due to scattering. With the available modeling algorithm, the impact of heterogeneity scales on seismic scattering and optimum acquisition geometries will be investigated in future studies. Keywords: gas hydrate, Monte Carlo simulation, random media, seismic scattering, permafrost, heterogeneity, petrophysics ∗  Jun-Wei Huang: Phone: +1 416 978 5177 Fax +1 416 978 7606 E-mail:  INTRODUCTION Natural gas hydrates, a type of inclusion compound or clathrate, are composed of gas molecules trapped within a cage of water molecules. The presence of gas hydrates in permafrost regions and in deep ocean along continental slopes have been confirmed by core samples recovered from boreholes. However, neither the geological nor the seismic model of the gas hydrate reservoirs is fully understood. Significantly elevated P- and S-wave velocities from borehole logs indicate the cementation of sediment due to the formation of gas hydrates. However cross-well survey [1] and zero-offset Vertical Seismic Profiling (VSP) data [2] indicate strong attenuation in the gas hydrate bearing zone. Guerin et al. [3] explained the acoustic wave dissipation from a rock physics perspective. The high velocity and high attenuation properties of gas-hydrate-bearing sediments can also be attributed to scattering and elastic wave mode conversions [4]. Regardless of the chemical composition, geologic age and tectonic history of the earth media, power spectra of all available data such as sonic logs uniformly decay as 1/f0.5~1.5, where f is the spatial frequency [5]. This type of data was referred to as 1/f-noise and is characterized by their fractal nature or scale-invariance [6]. Using windowed Fourier Transform the power spectrum of P-wave velocity logs from Mallik 2L-38 present decay with different decay rates corresponding to different lithologies (Figure 1). The depth window is 100 m with 10 m step for each Fourier Transform. Sediments without hydrates (600 m – 900 m deep) show a faster decay rate than hydrate interval and permafrost interval. Von Karman style autocorrelation algorithm is a widely used approach throughout earth science to model earth phenomena with fractal nature due to its self-invariant and band-limited properties. From observations of a von Karman style autocorrelation function and a bimodal probability density function in the P-wave velocity logs from MITI Nankai Trough Post Survey well #1, Kamei et al. [7] constructed a random heterogeneous acoustic model and successfully explained the frequency dependence of the Bottom Simulating Reflector (BSR) as well as the scattering attenuation. The same stochastic features were  observed in the borehole logs from Mallik 2L-38 and Mallik 5L-38. To include energy dissipation due to elastic wave mode conversion, we modified Kamei’s algorithm [7] and constructed 2D and 3D elastic heterogeneous random models. Here, we provide an algorithm on the model construction and show how such models can provide first order gas hydrate volume estimates and explain the high attenuation at the seismic frequency range. Results from this study demonstrate that within the stability zone, the lithological heterogeneities (e.g., lithology, porosity, fluid saturation, pore pressure etc.) partly control the distribution of gas hydrates.  Figure 1. The space and spatial frequency decomposition of Mallik 2L-38 P-wave velocity logs with a 100 m depth window and 10 m step. The color code shows the logarithm of the power spectrum.  MODEL CONSTRUCTIONS A multi-dimentional stationary random field is characterized by an autocorrelation function (ACF) and a probability density function (PDF). A spectral-based approach was adopted to simulate the randomly heterogeneous model [8]. In our statistical analysis, two borehole logs from Mallik 2L-38 and Mallik 5L-38 are available. Mallik 2L38 covered from the permafrost down to the bottom of gas hydrate stability zone, with 2 meter depth sampling interval, while Mallik 5L-38 only covered the gas hydrate bearing zone with 0.15 meter sampling interval. Due to the wide coverage, Mallik 2L-38 was used for the model construction. However Mallik 5L-38 is considered more statistically representative of the local area than Mallik 2L-38. For the first step of model construction, a linear trend was identified and then subtracted from the processed well logs of Vp and Vs in the sediment strata (Figure 2). Velocity logs were scaled by  estimated standard deviations. Secondly, a bimodal probability density function was obtained by fitting the histogram of well logs from the gas hydrate strata: 890 m – 1144 m (Figure 3). The bimodal probability function was considered as a linear combination of two weighted Gaussian distribution functions [7]: ⎛ (V − μ1 )2 ⎞ (1 − w1 ) ⎛ (V − μ2 ) 2 ⎞ (1) w1 p(V ) = + exp − exp −  of Vp and Vs was assumed. In our case, the PSDF is  where w1 is the weight factor, V is the random parameter related to the velocity, μ and σ are mean velocity and standard deviation, respectively. Least square fitting provide σ1 values. To ensure a zero mean, weight factor w1 can be decided with the two mean values: w1=μ2/(μ2-μ1), and to have a unit variance, the second variance σ2 is derived by: σ 22 = ⎡⎣ μ 2σ 12 − ( μ2 − μ1 )( μ2 μ1 + 1) ⎤⎦ μ1 . The mean values  k = k x2 ax2 + k y2 a y2 + k z2 az2  2πσ 1  ⎜ ⎝  2σ 12  ⎟ ⎠  2πσ 2  ⎜ ⎝  2σ 22  ⎟ ⎠  2 D : S d (k ) =  4πν ax a y  (1 + k 2 )  ν +1  , k = k x2 ax2 + k y2 a 2y and  3  3D : S d (k ) =  ( 4π ) 2 ax a y az 3 2 ν +2  (1 + k )  3⎞ ⎛ Γ ⎜ν + ⎟ 2⎠ ⎝ ⋅ , Γ (ν )  (3)  where kx,ky,kz are the wave number component, and ax,ay,az are the characteristic scales in 3-D. More general expressions can be found in Lord [11].  and standard deviations are listed in Table 1. The third step was to fit the autocorrelation of the Vp logs and Vs logs by a von Karman style autocorrelation function [9]: C ( r ( x) ) =  Gv ( r ( x) ) rν Kν ( r ) , = ν −1 Gv (0) 2 Γ (ν ) 2  2  ⎛ x ⎞ ⎛ y ⎞ ⎛ z ⎞ r ( x) = ⎜ ⎟ + ⎜ ⎟ + ⎜ ⎟ ⎜ ⎟ ⎝ ax ⎠ ⎝ a y ⎠ ⎝ az ⎠  (2)  2  where x is a vector in the multi-dimentional random field, x=[x,y] in 2D and x=[x,y,z] in 3D, and Kν is the modified Bessel function of the second kind, Γ(ν) is the Gamma function, and ν is the Hurst number describing the roughness of the medium. The vertical characteristic scale az and the Hurst number ν for both P-wave and S-wave velocity logs were obtained by a least square curve fitting (Figure 4). The same statistical analysis was applied to Mallik 5L-38 logs, and the estimated vertical correlation length and Hurst number is listed in Table 1 for comparison. Due to the original data processing including median filters and averaging [10], the Hurst number estimated from Mallik 2L-38 may be overestimated. The desired power spectrum density function (PSDF) of the medium, Sd(K), was obtained from ACF equation (2) by a Fourier transform. The intermediate Gaussian random field g(x) was then generated from the PSDF which is followed by a Cumulative Distribution Function (CDF) mapping to construct the non-Gaussian bimodal random field: Vp and Vs [8]. A strong positive correlation  Figure 2: Two linear trends were identified in the sediment with and without gas hydrates (Mallik 2L-38). Constant values were added to the shallow permafrost interval and beneath the gas hydrate interval due to the absence of data. Density measurements above 370 m are unreliable due to the casing problem.  Figure 3: Least square curve fitting of velocity logs histograms, in which μ and σ are the mean values and standard deviation, respectively. “Vp-linear-mean/StvP” means P-wave velocity logs subtracted from the linear trend and the mean value, then divided by estimated standard deviations. The scaling was used to ensure zero mean and unit variance.  As an example, we simulated a 2D elastic random medium and selected two synthetic well logs with  zero mean shown in Figure 5 to compare with the field well logs. We can see that the random heterogeneity simulation did not reproduce the same logs but generated pseudo-logs exhibiting the same stochastic features as the field logs.  hydrate-bearing sediment intervals are identified from Mallik 2L-38 logs. Again, the random models exhibited the same stochastic features as the well logs. In the example, the lateral correlation length for permafrost and sediment were assumed to be 200 m, while the gas hydrate interval has anisotropic correlation length: 100 m along X direction and 50 m along Y direction. Table 1. List of the best fitting parameters for zero mean logs. The standard deviations, and vertical characteristic scales are scaled up by estimated standard deviations.  Best Fitting of Mallik 2L-38 P-wave velocity logs from gas hydrate interval μ1(m/s) μ2(m/s) σ1(m/s) σ2(m/s) az(m) ν  Figure 4: Least square curve fitting of the autocorrelations of scaled velocity logs. The vertical correlation length in meters and the Hurst number were given by the fitting results. The strong random vibrations at large lags indicate a weak correlation.  -400 500 128.4 355.3 Best Fitting of Mallik 2L-38 S-wave from gas hydrate interval μ1(m/s) μ2(m/s) σ1(m/s) σ2(m/s) -150 200 46.0 196.2 Best Fitting of Mallik 5L-38 P-wave from gas hydrate interval  7.9 0.59 velocity logs  μ1(m/s) -156  az(m) 19.3  μ2(m/s) 553  σ1(m/s) 345.7  σ2(m/s) 350.3  az(m) ν 4.3 0.94 velocity logs ν 0.28  Best Fitting of Mallik 5L-38 S-wave velocity logs from gas hydrate interval μ1(m/s)  μ2(m/s)  σ1(m/s)  σ2(m/s)  az(m)  ν  -101  358  223.6  219.1  21.3  0.20  Figure 5. Zero mean field well logs were compared with the simulated well logs. The linear trend and the mean values are subtracted. Instead of reproducing exactly the same logs, the simulated logs only represent the same stochastic features as the field logs.  The same statistical analysis was applied to the permafrost interval and the sediment interval without gas hydrates. As expected, both intervals showed Gaussian distributions. The Cumulative Distribution Functions (CDF) of models and logs are shown in Figure 6. The mismatch between the model and the density logs is as a result of the casing problem in the borehole which caused the density measurements in the shallow part to be unreliable (Figure 2). Multi-dimensional (2D and 3D), multi-variable (acoustic and elastic) models are thus available for gas hydrate volume estimation and further seismic scattering analysis. Figure 7 shows an example of 3D P-wave velocity petrophysical model with 2km length and width, 1.1km depth. Permafrost, sediment, and gas-  Figure 6. Cumulative Distribution Functions of models and logs for (from the top to the bottom) permafrost, sediment, and gas-hydrate-bearing sediment. Note the bimodal distribution of gas hydrate interval has a different CDF from Gaussian distribution field. The mismatch for density was caused by the unreliable logs in the shallow part.  approximately 7x108 m3 under 273.15 temperature and 101.325 kPa pressure.  Figure 7. A 3D P-wave velocity model based on the Pwave logs from Mallik 2L-38. Three intervals: permafrost (PM), sediments (SD) and gas-hydrate-bearing sediments (GH) were simulated individually with different correlation length. Only for the gas hydrate interval, the correlation length along X direction (100 m) is different from along Y direction (50 m).  K  Figure 8. Empirical relationship between hydrate concentration and P-wave velocity established by Lee [13] for Mallik2L-38 case. In this rock physics model, the volume clay content was assumed to be 30%.  GAS HYDRATE VOLUME ESTIMATION In Figure 8, the preliminary relationship between the gas hydrate concentration and the P-wave velocity was established based on the assumption of constant porosity (33%) and clay volume content (30%) [12]. If the heterogeneous distribution of gas hydrates is governed by the fractal nature of lithology, the volume content of gas hydrate is then determined by the correlation length and Hurst number. Figure 9 shows the volume fraction of gas hydrate as a function of correlation length for a specific Hurst number. We used a Hurst number of 0.58, estimated from Mallik 2L-38 P-wave velocity logs. The gas hydrate volume content increased from 16% to 18%, as correlation length increased from 50 m to 1000 m. For constant correlation length (ax=ay=500), an increasing Hurst number results to an increase of gas hydrate content, with lower and upper limits of 16.5% to 22.5% (Figure 9). Due to the limited number of boreholes in Mallik area, additional information is needed to obtain a volume estimate of gas hydrate in Mallik area, such as lateral correlation length and gas hydrate random distribution region. Acoustic impedance inversion of industry’s 3D surface seismic data provided the lateral extension of the gas hydrate interval in Zone C (identified in [13]) about 0.25 km2, and the lateral correlation length about 500 m [13]. If we assume that gas hydrates identified in zone C have a 100 m vertical distribution range, free gas stored in gas hydrate format is  Figure 9. The variance of gas hydrate volume fraction with horizontal correlation length (top) and with the Hurst number (bottom). The color code shows the volume fraction of gas hydrate in a 3D randomly heterogeneous model.  SEISMIC SCATTERING Wave attenuation is an important physical property of hydrate-bearing sediments that is rarely taken into account in site characterization  with seismic data. Estimates of Q-factor (inverse of attenuation) from seismic data usually neglect effects from scattering and assume that attenuation is entirely due to intrinsic rock properties [14]. Here we show that part of the strong attenuation can be explained by scattering. A simple Finite Difference numerical experiment was carried out on a 3D model consisting of permafrost, sediments without gas hydrates with infinite long correlation lengths (layered), and heterogeneous gas-hydratebearing interval (Figure 2). Beneath the gas hydrate zone, two homogeneous layers were added to generate reference reflections. A finite difference modeling code was used to solve 3D elastic wave equations [15]. The P-wave velocity in gas hydrate zone varies from 1909 m/s to 5128 m/s and its contrast to the linear background is 28% in a root mean square sense. The correlation lengths for the hydrate interval were 100 m in X direction and 50 m in Y direction. The offset was 160 m. In our numerical experiment, wave-number (k) associated with the source and correlation length (a) of the model make ka≈10. Thus using Wu’s definition [16] the scattering regime covered both large angle (scattering attenuation) and small angle scattering regime. The energy loss due to scattering was demonstrated in Figure 10. Even if the back scattered energy (up-going wave) was compensated, the scattering within the heterogeneous gas hydrate reservoir contributed up to 10dB energy loss. Seismic scattering and the effects on conventional seismic data processing will be further analyzed using the available random models.  Figure 10. Separated P and S wave components of near offset synthetic VSP data demonstrated up to 10 dB energy loss due to scattering in the heterogeneous gas hydrate reservoir with the compensation of up-going energy.  CONCLUSIONS Randomly heterogeneous elastic 2D and 3D multivariable models of permafrost, sediments with and without gas hydrates were constructed based on the spectral approach, which honored the statistical features represented by logs. Assuming that the gas hydrate heterogeneous distribution is determined by the fractal nature of earth media, we can estimate the gas hydrate volume and thus free gas in the gas hydrate zone. It was found that the increasing Hurst number will increase the gas hydrate volume more significantly than the increasing correlation length. However, care must be taken when estimating the Hurst number from smoothed logs and additional information is a must to provide gas hydrate volume estimate. Energy loss due to scattering was demonstrated from a simple near offset numerical experiment. It implied that the scattering due to certain scale heterogeneity cannot be ruled out in the study of the attenuation observed in the field data. Further analysis will be conducted in this direction. In order to simulate scattering accurately, numerical experiments in 3D models must be used [17]. REFERENCES [1] Pratt, R.G., Hou, F., Bauer, K., and Weber, M.H. : Waveform tomography images of velocity and inelastic attenuation from theMallik 2002 crosshole seismic surveys; in Scientific Results from theMallik 2002 Gas Hydrate Production Research Well Program, Mackenzie Delta, Northwest Territories, Canada, (ed.) S.R. Dallimore and T.S. Collett; Geological Survey of Canada, Bulletin 585. [2] Milkereit, B., Adam, E., Li, Z., Qian, W., Bohlen, T., Banerjee, D., and Schmitt, D.R.: Multi-offset vertical seismic profiling: an experiment to assess petrophysical-scale parameters at the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate production research well; in Scientific Results from the Mallik 2002 Gas Hydrate Production Research Well Program, Mackenzie Delta, Northwest Territories, Canada, (ed.) S.R. Dallimore and T.S. Collett; Geological Survey of Canada, Bulletin 585, 2005. [3] Guerin, G., and Goldberg, G., Modeling of acoustic wave dissipation in gas hydrate-bearing sediments. Geochemistry Geophysics Geosystems, Vol.6, Number 7,2005. [4] Huang, J.W., and Milkereit, B., Seismic imaging of gas hydrate distribution – a case study,  Extended abstract at 2006 CSPG-CSEG-CWLS convention, 2006. [5] Holliger, K., and Goff, J., A Generic Modle for the 1/f-Nature of Seismic Velocity Fluctuations, in Heterogeneity in the Crust and Upper Mantle: Nature, Scaling, and Seismic Properties, edited by John A. Goff, and Klaus Holliger, Kluwer Academic/Plenum Publisher, New York, 2003. [6] Mandelbrot, B., The Fractal Geometry of Nature, Freeman, New York, 1983. [7] Kamei, R., Hato Masami, and Matsuoka, T., 2005, Random heterogeneous model with bimodal velocity distribution for Methan Hydrate exploration, exploration Geophysics, 36, 41-49. [8] Yamazaki, A., and Shinozuka, Digital generation of non-Gaussian stochastic fields, J. Eng. Mech., ASCE, 114, 1183-1197, 1988. [9] Goff, J.A., and Jordan, T.H., Stochastic modeling of seafloor morphology: Inversion of sea beam data for second-order statistics: Journal of Geophysical Research, 93, 13589-13608, 1988. [10]. Collett, T.S., Lewis, R.E., Dallimore, S.R., Lee, M.W., Mroz, T.H., and Uchida, T.,: Detailed evaluation of gas hydrate reservoir properties using JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well downhole well-log displays, in Scientific Results from JAPEX/JNOC/GSC Mallik 2L-38 Gas Hydrate Research Well, Mackenzie Delta, Northwest Territories, Canada, (ed.) S.R. Dallimore, T. Uchida, and T.S. Collett; Geological Survey of Canada, Bulletin 544, p. 295-311, 1999. [11] Lord, R.D., 1954, The use of the Hankel Transform in Statistics I. General Theory and Examples, Biometrika, Vol. 41, No. 1/2, pp.44-55. [12] Lee, M.W., and Collett, T.S., Comparison of Elastic Velocity Models for Gas-Hydrate-Bearing Sediments, in Natural Gas Hydrates Occurrence, Distribtuion, and Detection, edited by Charles K. Paull, and William P. Dillon, American Geophysical Union, Washington, D.C., 2001. [13] Bellefleur, G., Riedel, M., and Brent, T., Seismic characterization and continuity analysis of gas-hydrate horizons near Mallik research wells, Mackenzie Delta, Canada, The Leading Edge, May, 2006. [14] Bellefleur, G., Riedel, M., Brent, T., Wright, F., and Dallimore, S.R., Implication of seismic attenuation for gas hydrate resource characterization, Mallik, Mackenzie Delta, Canada. Journal of Geophysical Research, Vol. 112, 2007.  [15] Bohlen, T., Parallel 3-D viscoelastic finite difference seismic modeling, Computers & Geosciences, 28, 887-899, 2002. [16] Wu, R.S., Seismic wave scattering, invited article for the "Encyclopedia of Geophysics", ed. By D. James, Van Nostrand Reinhold and Comp., 1166-1187, 1989. [17] Bohlen T., and Milkereit, B, 2001, Parallel 3D Viscoelastic Finite-Difference Seismic Modeling, Extended abstract at EAGE 63rd Conference & Technical Exhibition — Amsterdam, The Netherlands.  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items