Delineation of gas hydrate reservoirs in the Ulleung Basin using unsupervised multi-attribute clustering without well log data

Accepted Manuscript Delineation of gas hydrate reservoirs in the Ulleung Basin using unsupervised multiattribute clustering without well log data Jaewook Lee, Joongmoo Byun, Bona Kim, Dong-Geun Yoo PII:

S1875-5100(17)30310-4

DOI:

10.1016/j.jngse.2017.08.007

Reference:

JNGSE 2257

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 27 April 2017 Revised Date:

7 August 2017

Accepted Date: 8 August 2017

Please cite this article as: Lee, J., Byun, J., Kim, B., Yoo, D.-G., Delineation of gas hydrate reservoirs in the Ulleung Basin using unsupervised multi-attribute clustering without well log data, Journal of Natural Gas Science & Engineering (2017), doi: 10.1016/j.jngse.2017.08.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Delineation of gas hydrate reservoirs in the Ulleung Basin using unsupervised multi-attribute clustering without well log data Jaewook Leea, Joongmoo Byuna, Bona Kima and Dong-Geun Yoob a RISE Lab., Hanyang University, Seoul, Korea

(KIGAM), Daejeon, Korea

Abstract

RI PT

b Petroleum and Marine Research Division, Korea Institute of Geoscience and Mineral Resources

SC

To target a well location for gas hydrate reservoirs, many previous studies have used only some discontinuous indicators such as the bottom simulating reflector (BSR) and dimming or acoustic

M AN U

impedance (AI). However, when drilling the first well, there are no well log and core analysis data for AI inversion and geological analysis. Moreover, these methods may have some risk of finding gas hydrate zones and overestimating the reservoir distribution. First of all, we inverted AI using seismic data and root-mean-square (RMS) velocity for AI inversion excepting well logs. In this inversion, the low-frequency impedance variations were estimated from RMS velocity by the product of the interval

TE D

velocity calculated by the Dix equation and the bulk density from this interval velocity. Then, this low frequency information was integrated with the seismic frequency information obtained from a reflectivity series by sparse-spike impedance inversion. To prevent overestimation due to the use of AI

EP

alone, we focused on another rock property, shear impedance (SI), which indicates the ‘rigidity’ of rocks, as gas hydrate consolidate sediments and increases the value of SI noticeably compared to that

AC C

in surrounding non-reservoirs. To estimate this property, we used this inverted AI and partial stack seismic data using the two-term Fatti’s equation. This amplitude variation with an offset (AVO) equation excludes the density term, which is too sensitive to noise, and has only the two terms of AI and SI. As a result, we applied K-mean clustering, which is the method of unsupervised machine learning, to delineate a more accurate and quantitative distribution of hydrated reservoirs from areas with higher impedances and higher values of two additional attributes (RMS amplitude and instantaneous frequency) compared to surrounding formations. In conclusion, we verified that this workflow is useful for identifying the distribution of potential reservoirs and aiding in well-site

ACCEPTED MANUSCRIPT location determination. Keywords: Quantitative seismic interpretation, multi-attribute analysis, Ulleung Basin, gas hydrate,

1. Introduction

RI PT

unsupervised clustering, machine learning

With a focus on deep-sea small fields and unconventional resources in complex reservoirs, the exploration and development of reservoirs have become increasingly difficult and risky. Although

SC

these reservoirs are easily detectable due to improvements in seismic exploration, it is more difficult

M AN U

to estimate their hydrocarbon reserves and develop them because of the higher geological or geophysical complexity than in previous reservoirs found from direct hydrocarbon indicators (DHI) such as bright spots (Sheriff and Geldart, 1995; Brown, 2010). To overcome this problem, the techniques of quantitative seismic interpretation (QSI), based on rock physics theory, have been developed for both better interpretations to detect prospective resources and accurate estimation to

TE D

determine whether those resources are commercial (Hilterman, 2001; Avseth et al., 2006; Mavko et al., 2009). Using techniques such as impedance inversion and seismic attributes, we can delineate the distribution of reservoirs and quantify reservoir properties away from the well site.

EP

Many studies of gas hydrate reservoirs have determined the existence of hydrates by detecting distinct events, such as the bottom-simulating reflector (BSR). However, these indicators are not

AC C

sufficient to determine the reservoir distribution and several seismic attributes have been analyzed to overcome the lack of information, such as the root-mean-square (RMS) amplitude and instantaneous frequency (Lee et al., 2009; Lee et al., 2016). Nevertheless, this traditional method has questionable accuracy (Tsuji et al., 2009) and the candidate zones need to be examined to determine the acoustic and shear impedance of the target well location for each zone. Then, after drilling an exploratory well, the inversion of acoustic impedance (AI) and shear impedance (SI) was estimated using sonic and density logs (Inamori et al., 2006). Furthermore, reservoir parameters, such as the porosity and concentration of gas hydrate, were estimated based on the relationship between the rapid increase in

ACCEPTED MANUSCRIPT velocity of the hydrated sediments and the concentration of hydrates (Lu et al., 2002; Chatterjee et al., 2016). However, many fields lack well log data because of the immense expense of drilling and the significantly high risk of field development. In these cases, applicable QSI techniques for determining

RI PT

sweet spots may be limited and difficult. The main goal of this research was to delineate the gas-hydrate reservoir distribution and to identify sweet spots by adding a step that estimates these impedances without the use of well logs and clustering the results of the seismic attribute analysis, based on the unsupervised machine learning,

SC

objectively and quantitatively. In this study, these methods were applied to the case of unconventional gas hydrate reservoirs in the Ulleung Basin, Korea. After a geological analysis of the target basin in

M AN U

the study area for selecting the effective attributes, we first identified candidate zones by detecting the distribution of hydrocarbon reservoirs and structural discontinuities and examined their feasibility and stability of development using multi-attribute analysis. In addition, we determined the sweet spot, as a well location, by delineating the approximate continuity and thickness of the reservoirs in each

TE D

candidate zone by a three-step inversion to estimate the impedance information and unsupervised classification using K-mean clustering from four key attributes. The estimated impedances in this study were verified by comparison with other results including velocity information such as existing

EP

vertical seismic profile (VSP) data, two well logs (UBGH2-2-2, UBGH1-09) located within the 3D seismic range, and inverted impedance models from a typical inversion using these logs from previous

AC C

studies (Tak et al., 2013; Jeong et al., 2014).

2. Geological setting

The study area was the hydrate-bearing region in the Ulleung Basin, which is located in the southwestern part of the East Sea. The basin is a bowl-shaped back-arc basin that deepens northeastward, and the seafloor is about 2000–2300 m beneath the sea surface and almost flat. The basin formed during the Late Oligocene to Early Miocene from crustal extension related to southward movement of the Japan Islands (Tamaki et al., 1992; Kim et al., 2007). Faulting and folding were

ACCEPTED MANUSCRIPT generated by a change of the stress regime from tension to compression at the latest Neogene (Lee et al., 2001), and upward movement of gas-rich fluids and formation of hydrate-bearing layers also occurred (Ryu et al., 2013). In this basin, many siliciclastic mass-transported deposits (MTD) were accumulated since the late Miocene and distal turbidites and hemipelagic sediments were deposited

RI PT

after the Pleistocene (Lee and Suk., 1998). Gas hydrates and associated gas are dominant in the turbidite and hemipelagic sediments, and these sediments show seismic indicators such as the BSR, strong reflections, and chimneys (Yi et al., 2011; Bahk et al., 2013).

SC

A 3D seismic survey was conducted in the central-northern part of the basin by the Ulleung Basin Gas Hydrate Drilling Expedition (UBGH) projects in 2007 and 2010. In addition, well data including

M AN U

sonic and density logs and VSP data from five wells were obtained in the survey area (Fig. 1). As shown by these data, gas hydrate-bearing reservoirs have obvious characteristics, such as higher Pwave velocity and shear modulus. They can be categorized into two types of gas hydrates, pore-filling and chimney-type. The pore-filling gas hydrates almost always occur in turbidite sands, whereas the

TE D

chimney-type hydrates occur mostly in fracture-dominated zones. Earlier studies investigated the quantitative estimation of shear impedance and hydrate saturation using well log data and the supervised neural network method. However, this workflow requires multiple well log data, and the

EP

cost of obtaining these data can be enormous. Therefore, an effective method without a reliance on well data could be meaningful for determining well location.

AC C

In this case, research was focused on pore-filling gas hydrates with low structural discontinuity and on indicators such as BSR and strong reflections in 3D seismic data. However, these indicator events are discontinuous, and it was difficult to identify reservoirs using only these indicators. Therefore, the aims of this research were the identification of reservoirs and the determination of sweet spots using multi-seismic attribute analysis, including impedance inversion.

ACCEPTED MANUSCRIPT 3. Method 3.1 Multi-attribute analysis Pore-filling gas hydrates can generate several seismic indicators, such as BSR and strong reflections (Fig. 2). However, it is not easy to fully delineate gas hydrate and free gas reservoirs because of their

RI PT

discontinuities and ambiguities. Therefore, 3D seismic multi-attribute analysis was conducted to identify potential gas hydrate-bearing zones for determining a sweet spot. Gas hydrate-bearing layers have peculiar elastic characteristics, such as sharply higher P- and S-wave velocities caused by

SC

consolidation of unconsolidated sediments by gas hydrate (Stoll and Bryan, 1979). Moreover, free gas-bearing layers beneath gas hydrate layers can be detected by a sharp decrease in both of these

M AN U

velocities as well as the attenuation of high-frequency components. Furthermore, from the viewpoint of geological structures, lateral continuity and vertical thickness should be considered along with the avoidance of faults, chimneys, and fractures. Seismic attributes can reveal these features more successfully than original seismic data, so 3D seismic multi-attribute analysis is highly efficient in the identification of candidate zones for gas hydrates. Thus, five geologically and physically meaningful

TE D

attributes were chosen: two amplitude attributes (RMS amplitude and instantaneous frequency) for detecting hydrates and three geometric attributes (variance, curvature, and dip illumination) for

EP

avoiding discontinuities.

In addition, it is important to identify structural discontinuities such as faults, chimneys, and

AC C

fractured zones because gas leakage from these discontinuities can reduce the stability of gas production. Therefore, in this case, three typical attributes, variance, structural dip, and curvature, were applied to measure these geologic features. As the opposite of coherence, variance is a type of 3D attribute that emphasizes lateral changes of AI by calculating the trace-to-trace variability in a specified sample interval, whereas RMS amplitude indicates the vertical variation of AI in a given interval. Therefore, faults, fluvial channels, chimneys, and salt domes could be highlighted. Furthermore, structural dip helps to improve the calculated result of variance and reduce the uncertainty of interpretation. Structural dip is an important geological feature, and many attributes such as consistent dip, dip deviation, local structural dip, and dip illumination are related to this

ACCEPTED MANUSCRIPT feature. For this reason, dip illumination (Aqrawi et al., 2011) was applied as an effective way to emphasize the edge of dipping reflectors, such as those of chimney and column structures, in the candidate zones.

RI PT

3.2 Three-sequential impedance inversion The purpose of impedance inversion is to build layer-based reservoir models for better interpretation and to estimate the probabilities of hydrocarbon contents. Well logs are crucial to the building of an impedance model for three reasons: deterministic wavelet extraction, low frequency modeling, and

SC

result verification.

M AN U

However, in the case of no available well log data, other available data such as 3D seismic and RMS velocity data may solve the problem of no well log data. In this study, the wavelet for inversion was extracted by a statistical method from near-stack seismic data instead of by a typical deterministic method from well logs. Although the seismic data have a lower resolution than well logs, the statistical wavelet estimation has advantages compared to a deterministic wavelet extraction because

TE D

deterministic phase mismatches use well log data from a short logging sequence and inaccurate phase extrapolation away from well locations. Furthermore, we can avoid AVO effects by using near-stack data, rather than full-stack data, which are generated from the summation of near-stack and far-stack

EP

data. The minimum-phase wavelet was extracted from the specific time interval (2950–3150 ms) of the near-stack data acquired from the region including Zones 1 and 3 by a statistical method using an

AC C

autocorrelation of each trace.

Secondly, most inversion techniques reveal an absolute impedance model based on information from well log data. These techniques obtain broadband frequency components from well log data and solve an absolute impedance model, but other techniques without well log data usually create only a relative impedance model. A relative impedance model shows high and low relative impedances corresponding to different lithologies or pore fluids, but these are not comparable to absolute impedances because of the lack of a low-frequency initial model from well log data (Cooke and Cant, 2010). Given this case, RMS velocity and some rock physics equations were used to solve this

ACCEPTED MANUSCRIPT problem rather than well log data. To overcome the lack of information from well logs, we integrated three algorithms for AI and SI inversion: the band-limited algorithm for establishing the pseudo-AI model; the sparse-spike

RI PT

algorithm for improving the inverted AI model; and simultaneous AVO inversion for estimating SI from the partial stack data and inverted AI. After the inversion, the result for the inverted impedance models without well log data was verified by measuring against the acoustic impedance log at

UBGH2-2-2 and the shear impedance log at UBGH1-09, which is located in the center of the 3D

3.3 Unsupervised clustering

M AN U

between the inverted model and the impedance log.

SC

seismic volume. The accuracy and suitability could be evaluated from the very strong correlation

Since the birth of the neural network, there have generally been two categories of many neural network methods: supervised methods and unsupervised methods. In conventional cases using well logs, supervised methods provide a train estimation model that estimates reasonable output data from

TE D

given input data. Assuming that the same model is applicable to similar inputs to calculate outputs, the model is established by a decrease in the error between the trained data and the original data. These methods can predict facies computed from the well logs and estimate properties such as gas hydrate

EP

saturation based on already-known values from the log data. According to previous studies, estimation of gas hydrate saturation is possible based on three well log data. Unlike conventional cases with well

AC C

logs, we could not estimate exact properties such as gas hydrate saturation because there was no exact value for subsurface layers. However, classification by unsupervised methods is meaningful because it can efficiently subdivide our input data, such as two attributes and two inverted impedances, into several groups.

After inverting AI and estimating SI, we had to delineate the quantitative reservoir distribution for determination of optimal well location by the significant characteristics of higher RMS amplitude, instantaneous frequency, and AI and SI that could characterize gas hydrate layers. To integrate the four types of data and classify some groups from them, we applied K-mean clustering (Lloyd, 1982),

ACCEPTED MANUSCRIPT which is a kind of unsupervised neural network for data classification, to extract reservoir bodies. From among the groups, we can select a specific group including characteristic properties, such as a higher AI and SI. Then, this group could be regarded as the potential reservoir distribution and this could be verified through the BSRs in the seismic data and the higher AI/SI layers in the inverted

RI PT

model. In this study, we identified four attributes based on K-mean clustering, which is a popular

SC

method for cluster analysis.

4. Identification of candidate zones

M AN U

4.1. Amplitude attributes

Among several amplitude attributes, RMS amplitude is a popular statistical quantity for measuring the variation in a specified time interval. From the seismic traces, RMS amplitude indicates the magnitude of impedance variation. Therefore, a higher RMS amplitude means a higher variation of pore fluid contents and lithologies. Gas hydrate causes a steep increase in both P- and S-wave velocity,

TE D

so a high RMS amplitude implies a contrast between hydrated sediments and surrounding lithologies/infill pore fluids. The time interval of the measurement was determined by the interval in which BSRs mainly occur (2950–3150 ms). The blue area in the RMS amplitude map in Fig. 3

EP

indicates values of RMS amplitude higher than those of the surrounding area and shows the BSR and strong reflections. The yellow area in the figure indicates a lower RMS amplitude and shows chimney

AC C

and column structures. By focusing on the BSR and the strong reflection by pore-filling gas hydrates, five regions with a higher RMS amplitude including the indicators (BSR and strong reflections) were selected as candidate zones (Fig. 4). As stated previously, the occurrence of gas hydrate results in a rapid increase of velocity and impedance because it fills the pore space and consolidates the grains. As the result of the RMS amplitude calculation suggests, the regions of higher RMS amplitude effectively show BSR and strong reflection. However, there are many other factors aside from pore fluid contents affecting the amplitude, and an additional attribute also should be analyzed. One additional effective attribute usually used to detect gas sands is instantaneous frequency.

ACCEPTED MANUSCRIPT Instantaneous frequency is a kind of instantaneous attribute from complex seismic traces derived by a Hilbert transform. Unlike the frequency of the wavelet, this attribute indicates the mean amplitude of the wavelet and provides information about thickness and lithology parameters. Therefore, it is commonly used to delineate thin bed tuning or abnormal attenuation due to hydrocarbons in pores

RI PT

(Thakur and Rajput, 2011). In the case of gas hydrate reservoirs, free gas layers below the BSR usually attenuate high frequency components more than hydrated layers. Thus, this facilitated visualization of the boundaries between upper hydrated layers and lower free gas layers for the five

SC

candidate zones. Using these attributes, candidate zones were ranked by evaluating the degree of lateral and vertical variation, to estimate the distribution of hydrated layers with higher frequency

M AN U

values above free gas layers with lower frequency values indirectly. In this evaluation, Zone 3 had the best continuity of high instantaneous frequency zones. Following Fig. 5, we assumed the red zones with a higher value above BSR to be hydrate-bearing layers and the blue zones with a lower value below BSR to be gas-bearing layers. Interpreting the multiple sections of the candidate zones with this attribute, the continuity and thickness of the red zones were largest in Zone 3. The other candidate

TE D

zones have much thinner and more discontinuous red zones than Zones 1 and 3. After an analysis of the two amplitude attributes, geometric attributes such as variance, curvature, and dip illumination

EP

were applied to measure the structural continuity for the candidate zones. 4.2 Geometric attributes

AC C

As shown in Fig. 6, structural discontinuities such as faults and chimneys are emphasized by these two attributes. Moreover, for better interpretation, curvature was calculated to delineate and highlight fault surfaces. Variance usually allows detection of the position of faults, but it is limited in describing the edge of the footwall and the hanging wall (fault plane). These edges can be visualized using K1 curvature, which is the largest positive curvature, and K2 curvature, which is perpendicular to K1 curvature (Roberts, 2001; Klein et al., 2008). In this case, the formations are almost flat, but the discontinuities at fault tips and chimneys have higher values of curvature than the surrounding areas. As a result, the combination of these attributes allows a better interpretation of structural

ACCEPTED MANUSCRIPT discontinuities (Fig. 7). Following Fig. 7, the structural discontinuities in the five candidate zones were ranked, and Zones 1 and 3 were shown to have fewer discontinuities, such as fault zones and chimneys, than the other

RI PT

candidates. Zones 1 and 3 could hold plausible sweet spots for well location (Table 1). Given this attribute analysis, impedance inversion could aid the delineation of the extent and thickness of gas-

SC

hydrate and free-gas reservoirs and identify sweet spots of Zones 1 and 3.

5. Delineation of reservoir distribution

M AN U

5.1. Establishing the pseudo-AI model

To solve the problem of the absence of low frequency contents in seismic data during impedance inversion, low frequency information was estimated from RMS velocity rather than well logs. Similar to the advantage of statistical wavelet extraction, an interval velocity calculated from RMS velocity

TE D

can be sufficient to provide the information and overcome the short and narrow coverage of the sonic log despite the lower resolution of that log. Subsequently, an initial low frequency model was built based on this information and the low frequency component of seismic data from a recursive

EP

algorithm (Lavergne, 1975). Band-limited inversion based on this algorithm can be used to estimate AI information. This information was used as a pseudo-AI log to obtain the inverted AI model.

AC C

Following the flowchart for establishing the pseudo-AI information (Fig. 8), the Dix equation was first used to convert RMS velocity to the interval velocity of each layer. Then, the AI information near 0 Hz was calculated from the product of the interval velocity and the bulk density estimated from this velocity (Mavko et al., 2009). In this case, the range of estimated P-wave velocity in the time interval (2900–3200 ms) was 1500.0–1676.6 m/s, so Gardner’s equation (Gardner et al., 1974), which is valid for the velocity range of 1500–6100 m/s, could be applied. As a result, the AI information with low frequency near 0 Hz was calculated.

ρ ( g / cc) =1.74Vp (km / s)0.25 (2-2)

ACCEPTED MANUSCRIPT Band-limited inversion can be used to integrate low frequency information from RMS velocity and mid-frequency information from seismic traces based on a recursive algorithm. The inverted AI information is roughly equivalent to AI logs for well sites seismic trace locations (Becquey et al., 1979). This information can yield lithology or pore-fluid variations without well control. In this case,

RI PT

the addition of the low frequency contents from seismic data can improve the inverted result (Fig. 9), allowing a better understanding than can be obtained from the results of conventional relative impedance inversion, the changes of lithology and pore fluid, and the visualization of the actual AI

SC

variation. Thus, the inverted AI model had structural details from seismic data and could serve a role

5.2. Improving the inverted AI model

M AN U

similar to that of AI logs.

Although the inverted AI model from band-limited inversion shows geologic features like the original seismic data, geologically meaningless events, such as random noise and side lobes of the wavelet, are prone to generate non-existent reflectivity series and cause erroneous interpretations.

TE D

Therefore, a sparse-spike inversion can help reduce these undesirable events. This inversion method yields a model having the smallest non-zero reflectivity series with the major structural features. This algorithm honors the assumption of layered geology (Li, 2002) and is suitable for this field.

EP

Categorizing large events and small events of a series from actual reflectivity, this inversion assumes that only the large events of the actual reflectivity series are geologically meaningful and should be

AC C

used to invert the model because the magnitudes of the nonsensical events in the series are usually small (Russell, 1988). Supposing that there are few considerable changes of lithology, the Linear Programming Method (Oldenburg et al., 1983), minimizing the L1-norm of the reflectivity coefficients of the model as a constraint, was chosen to decrease the effect of outliers such as the meaningless events. Consequently, the inverted model has a more simple and sparse variation of AI than usual band-limited inversion. The objective function (2-3) is shown below:

J = w1 (T −W ∗ r ) + w2 ∗ ∑ ri , (2-3) where T, W, and r indicate the original seismic data, the wavelet used in this inversion, and the

ACCEPTED MANUSCRIPT inverted reflectivity coefficients, respectively, and w1 and w2 are the weighting factor of each term. The first term is the difference between the inverted results and the input seismic data, and the second term is a constraint that controls the number of reflectivity coefficients. Optimization of the inversion parameters was computed from the value having the lowest RMS error between synthetic and original

RI PT

seismic data. After the parameter determination, the resulting model was compared to the previous model from band-limited inversion.

Compared with the previous inverted result, the layering was more obvious and the interpretation

SC

was improved by the suppression of meaningless AI variations (Fig. 10). Moreover, the inverted model (Fig. 11) was verified by computing the correlation values with an AI log at well as UBGH2-2-

M AN U

2 (first blind well), and with another inverted model from a model-based inversion using the well log and VSP data from a previous study (Tak et al., 2013). The correlation coefficients of these were considerable (UBGH2-2-2 AI log: 0.85; model-based inversion: 0.92). Consequently, this inverted result could be useful for interpreting AI variation in the 3D seismic area and estimating SI from these

TE D

AI and partial stack data.

Using this inverted model, a sweet spot could be identified using the lateral continuity and vertical thickness of the hydrate-bearing sands of Zones 1 and 3. Measuring these values of Zone 3 against

EP

those of Zone 1 in the sections of instantaneous frequency and the AI model (Fig. 12), Zone 3 seems more suitable for the location of the first well because of the comparatively thicker and more

AC C

continuous reservoirs that are presumed to have a higher instantaneous frequency and AI than surrounding lithologies.

5.3. Estimating SI from inverted AI and partial stack data After obtaining the AI model from the two inversion steps, estimating SI from these AI and partial stack data helps prevent the overestimation of reservoir distribution for the optimization of well location in Zone 3. Estimating the distribution of reservoirs for exact well location with only a higher AI has some risk because there are many factors other than pore fluids affecting the AI. Therefore, better prediction of reservoirs could be obtained from the areas with both a higher AI and higher SI

ACCEPTED MANUSCRIPT compared to surrounding lithologies. Shear impedance indicates the rigidity of rock, and hydrates consolidate sediments and increase their rigidity, as well as their AI (Fig. 13). The estimation of SI from AI and partial stack data (near-stack and far-stack) is based on AVO theory, and Fatti’s Equation (Fatti et al., 1994) was chosen from among many other equations. Fatti’s equation is a revised form of

RI PT

the equation of Aki and Richards (Aki and Richards, 1980) in terms of AI, SI, and density, and clearly shows the effects of the three properties on the three terms. We used a two-term equation including only AI and SI, rather than the three-term equation, because density is quite sensitive to random noise

SC

and the term of this property may decrease the accuracy and robustness of the inverted result. This equation is available at incident angles of less than 30 degrees, and it is sufficient for the application

M AN U

of two partial stack data (near-stack: 0–15o; far-stack: 16–30o). The two-term equation (2-4) is shown below:

R(θ ) ≈

∆Z p 2Z p

(1+ tan 2 θ ) − 8(

∆Z p 2Z p

)2

∆Z s 2 sin θ (2-4) 2Z s

TE D

where R(θ) is a reflectivity function of the incidence angle, Zp is AI, and Zs is SI. We recognized two terms in this equation: the first term is practically driven by AI, and the second is driven by SI. Using this, we determined the angle to the average angle of the partial stack data (near-stack: 7.5 o; far-stack:

EP

22.5 o). As a result, we finally generated an SI cube from two partial stack data and the inverted AI as in a conventional simultaneous inversion using the AI log.

AC C

For verification, we compared this SI to the estimated SI log at well UBGH1-09 (second blind well), calculated by the three-phase Biot-type Equation (TPBE) model that could quantify the gas hydrates from the two seismic impedances, AI, and SI in unconsolidated sediments (Jeong et al., 2014). Two results had a noticeably high correlation coefficient (0.82). From Fig. 13, we could easily see that higher AI layers did not always have SI values higher than the surrounding layers. This means that using only AI may produce a fallacious reservoir distribution and that layers having both a higher AI and higher SI are more likely to be promising reservoirs. Afterward, integrating the characteristics of the two attributes and two impedances, we tried to estimate a more accurate quantitative reservoir

ACCEPTED MANUSCRIPT distribution for determining well location. 5.4. Estimating reservoir distribution by an unsupervised clustering K-mean clustering is a type of unsupervised training algorithm for dividing input data into K sets of

RI PT

output by minimizing the sum of distance functions from each point of data to each center of K groups. The specified value of K should be less than or equal to the number of input data. In the process of dividing the data objects into K groups, the similarity of the objects in each group is

SC

increased. Following a classification by K-mean clustering (Fig. 14), there were four classes: Class 1 (green), Class 2 (blue), Class 3 (pink), and unclassified residuals (yellow). Among the four classes,

M AN U

Class 3 highlights the layers with four higher values of our selected seismic features around the BSR, and we could estimate this class as having promising reservoirs. Then, we isolated the bodies showing Class 3 from the 3D seismic data, and these extracted bodies were used to provide a better interpretation and target a well location, considering the lateral continuity and vertical thickness of these bodies. Furthermore, before sighting a well location, it will also be helpful to consider the

TE D

geological information such as sedimentation patterns, channel-levee systems, grain-size variation and

EP

gas-migration pathways with the result.

6. Conclusions

Locating gas hydrate reservoirs by geophysical indicators such as BSR and enhanced reflections is

AC C

possible; however, these indicators are mostly discontinuous and have limited value for determining the approximate extent of reservoirs. Thus, a multi-attribute analysis was conducted to identify candidate zones for sweet spot determination. The selection of the attributes was focused on geological and geophysical meanings reflecting the field geology and elastic characteristics of the reservoirs rather than on a mathematical meaning. In this study, amplitude attributes such as RMS amplitude and instantaneous frequency were used to reveal elastic features of the reservoirs because hydrate-bearing layers have a higher RMS amplitude due to the variation of acoustic impedance and the significant contrast of frequency due to the attenuation of high frequency components. In addition,

ACCEPTED MANUSCRIPT structural attributes, such as variance, structural dip, and curvature, were applied to emphasize discontinuous structures, such as faults, chimneys, and fractured zones, because these structures could reduce production stability. Next, two sequential impedance inversion techniques (band-limited inversion and sparse-spike inversion) without well logs were conducted to estimate the lateral

RI PT

continuity and vertical thickness of potential reservoirs in each candidate zone. The wavelet for the inversion was extracted statistically, and RMS velocity information was used to establish an initial low frequency model. Moreover, the extent and thickness of reservoirs was estimated and a sweet spot

SC

was successfully determined through the use of these two sequential impedance inversion methods without well log data. In addition, by estimating SI from the AI model and partial stack data based on

M AN U

AVO theory, we were able to isolate the reservoirs from the seismic section using K-mean clustering and target the optimal well location. In conclusion, the optimized workflow for the case without well logs in the target field will be used to identify sweet spots as well as the distribution of reservoirs where only seismic data exist. Consequently, the methodologies developed in this research may be

productivity.

TE D

beneficial for reducing production risk, establishing a more optimized drilling plan, and enhancing

Acknowledgement

AC C

EP

This work was supported by the Human Resources Program in Energy Technology of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) granted financial resource from the Ministry of Trade, Industry & Energy, Republic of Korea (NO. 20164010201120). This research was also supported by the Basic Research Project (GP2012-026) of the Korea Institute of Geoscience and Mineral Resources (KIGAM) funded by the Ministry of Science, ICT and Future Planning of Korea. Moreover, we wish to thank Hampson-Russell from CGG and Petrel from Schlumberger for providing an academic license.

References

Aki, K., Richards, P. G., 1980, Quantitative seismology: Theory and methods, W. H. Freeman and Co. Aqrawi, A. A., Boe, T. H., Barros, S., 2011, Detecting salt domes using a dip guided 3D Sobel seismic attribute, 81st Annual International Meeting, SEG, Expanded Abstracts, 1014–1018. Avseth, P., Mukerji, T., Mavko, G., 2006, Quantitative seismic interpretation, Cambridge University Press. Bahk, J.J., Kim, D.H., Chun, J.H., Son, B.K., Kim, J.H., Ryu, B.J., Torres, M., Riedel, M., Schultheiss, P., 2013. Gas hydrate occurrences and their relation to host sediment properties: results

ACCEPTED MANUSCRIPT from second Ulleung Basin gas hydrate drilling expedition, East Sea, Korea. Marine and Petroleum Geology, 47, 21–29. Becquey, M., Lavergne, M., Willm, C., l979, Acoustic impedance logs computed from seismic traces, Geophysics, 44, 1485–1501.

RI PT

Brown, A. R., 2010, Dim Spots in Seismic Images as Hydrocarbon Indicators, AAPG Search and Discovery Article, 40514. Chatterjee, R., Singha, D. K., Ojha, M., Sen, M. K., Sain, K., 2016, Porosity estimation from prestack seismic data in gas-hydrate bearing sediments, Krishna-Godavari basin, India. Journal of Natural Gas Science and Engineering, 33, 562-572.

SC

Cooke, D., Cant, J., 2010, Model-based Seismic Inversion: Comparing deterministic and probabilistic approaches, Canadian Society of Exploration Geophysicist Recorder, 35, 28-39.

M AN U

Fatti, J., Smith, G., Vail, P., Strauss, P., Levitt, P., 1994, Detection of gas in sandstone reservoirs using AVO analysis: a 3D Seismic Case History Using the Geostack Technique, Geophysics, 59, 1362-1376. Gardner, G.H.F., Gardner, L.W., Gregory, A.R., 1974, Formation velocity and density – the diagnostic basics for stratigraphic traps, Geophysics, 39, 770-780. Hilterman, F., 2001, Seismic Amplitude Interpretation: 2001 Distinguished Instructor Short Course, Society of Exploration Geophysicists.

TE D

Inamori, T., Saeki, T., 2006, Delineation of gas hydrate-bearing sediments by multi seismic attributes using 3D seismic survey in the eastern Nankai Trough area, 8th SEGJ International Symposium, 1–6. Jeong, T., Byun, J., Choi, H., Yoo, D.G., 2014, Estimation of gas hydrate saturation in the Ulleung basin using seismic attributes and neural network, Journal of Applied Geophysics 106, 37-49.

EP

Kim, G.Y., Yoo, D.G., Lee, H.Y., Lee, Y.J., Kim, D.C., 2007, The relationship between silica diagenesis and physical properties in the East/Japan Sea: ODP Legs 127/128. Journal of Asian Earth Science, 30, 448-456. Klein, P., Richard, L., James, H., 2008, 3D curvature attributes: a new approach for seismic interpretation, First Break, 26, 105-111.

AC C

Lavergne, M., 1975, Pseudo-diagraphies de Vitesse en offshore profond, Geophysical Prospecting, 23, 695–711. Lee, G. H., Suk. B. C., 1998, Latest Neogene-Quaternary seismic stratigraphy of the Ulleung Basin, East Sea (Sea of Japan), Marine Geology, 146, 205–24. Lee, G. H., Kim, H. J., Han, S. J., Kim. D. C., 2001, Seismic stratigraphy of the Ulleung Basin in the East Sea (Japan Sea) back-arc basin, Marine and Petroleum Geology, 18, 615–634. Lee, J., Byun, J., Kim, B., Yoo, D., 2016, Sweet spot identification for gas hydrate reservoirs in boreholefree field, 7th SEG/EAGE International Geosciences Student Conference Lee, M. W., Collett, T. S., Inks, T. L., 2009, Seismic-attribute analysis for gas-hydrate and free-gas prospects on the North Slope of Alaska, Natural gas hydrates— Energy resource potential and associated geologic hazards: AAPG Memoir 89, 541 – 554.

ACCEPTED MANUSCRIPT Li, Q., 2002, Sparse Spike Inversion and the Resolution Limit, CSEG Geophysics 2002. Lloyd., S. P., 1982, Least squares quantization in PCM, IEEE Transactions on nformation Theory 28(2), 129–137. Lu, S., McMechan, G. A., 2002, Estimation of gas hydrate and free gas saturation, concentration, and distribution from seismic data, Geophysics, 67(2), 582-593.

RI PT

Mavko, G., Mukerji, T., Dvorkin, J., 2009, The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media, 2nd ed., Cambridge University Press. Oldenburg, D.W., Scheuer, T., Levy, S., 1983, Recovery of the acoustic impedance from reflection seismograms, Geophysics, 1318–1337.

SC

Roberts, A., 2001, Curvature attributes and their application to 3D interpreted horizons, First Break, 19, 85-99.

M AN U

Russell, B. H., 1988, Introduction to seismic inversion methods: SEG Course Notes Series No. 2, Society of Exploration Geophysicists. Ryu, B.J., Collett, T.S., Riedel,M., Kim, G.Y., Chun, J.H., Bahk, J.J., Lee, J.Y., Kim, J.H., Yoo, D.G., 2013, Scientific results of the second gas hydrate drilling expedition in the Ulleung Basin (UBGH2), Marine and Petroleum Geology, 47, 1–20. Sheriff, R. E., Geldart, L. P., 1995, Exploration Seismology, 2nd ed., 415–418 Stoll, R.D. Bryan, G.M., 1979, Physical properties of sediments containing gas hydrates, Journal of Geophysical Research, 84, 1629–1634.

TE D

Tak, H., Byun, J., Seol, S.J., Yoo, D.G., 2013, Zero-offset vertical seismic profiling survey and estimation of gas hydrate concentration from borehole data from the Ulleung Basin, Korea, Marine and Petroleum Geology, 47, 204-213.

EP

Tamaki, K., Suyehiro, K., Allan, J., Ingle Jr., J.C., Pisciotto, K.A., 1992. Tectonic synthesis and implications of Japan Sea ODP drilling. In: Tamaki, K., Suyehiro, K., Allen, J., McWilliams, M., Proc. O DP, Sci. Results, 127/128, Ocean Drilling Program, College Station, TX, 1333–1348. Thakur, N. K., Rajput, S., 2011, Exploration of Gas Hydrates, Geophysical Techniques, Springer.

AC C

Tsuji, Y., Fujii, T., Hayashi, M., Kitamura, R., Nakamizu, M., Ohbi, K., Saeki, T., Yamamoto, K., Namikawa, T., Inamori, T., Oikawa, N., Shimizu, S., Kawasaki, M., Nagakubo, S., Matsushima, J., Ochiai, K. and Okui, T., 2009, Methane-hydrate occurrence and distribution in the eastern Nankai trough, Japan: Findings of the Tokai-oki to Kumano-nada methane-hydrate drilling program, Natural Gas Hydrates Energy Resource and Associated Geologic Hazards, 228-246. 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 AVO characteristics of the BSR in the Ulleung Basin, East Sea (Japan Sea), Marine and Petroleum Geology, 28, 1953-1966.

ACCEPTED MANUSCRIPT

Table 1. Multi-attribute analysis results for the candidate zones.

Instantaneous Frequency

2

SC

Zone 1

Geometric Attributes (Discontinuity) 1 (The least discontinuities) 5 (The most discontinuities) 2

5 (Most Discontinuous and thinnest)

3

4

4

TE D

1 (Most continuous and thickest)

EP

Zone 5

Higher than surrounding areas

AC C

Zone 4

3

M AN U

Zone 2 Zone 3

RI PT

RMS Amplitude

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 1. (a) Map of the Ulleung Basin in the East Sea, Korea (modified from Ryu et al., 2013). The yellow box represents the study area in the central p art of the basin, shown in (b). (b) A more detailed map of the study area showing the extent of the 3D seismic survey and the locations of the five wells where logs were acquired.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 2. 2D seismic section including well UBGH2-6 (N–S). This region was selected to identify gas hydrate-bearing sand layers within the turbidite facies. The seismic section shows the bottom simulating reflector (BSR) and strong reflections (modified from Ryu et al., 2013).

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 3. Map of root-mean-square (RMS) amplitude in the time interval of 2950–3150 ms. The blue area indicates values of RMS amplitude higher than those of the surrounding area and shows BSR and strong reflections. The yellow area indicates lower values of RMS amplitude and shows chimney and column structures.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 4. Candidate zones that have higher values of root-mean-square (RMS) amplitude and include bottom simulating reflector (BSR) and strong reflections.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 5. Amplitude, instantaneous frequency, overlapped section, and evaluation result of continuity and thickness of estimated gas hydrate-bearing layers (red) for each candidate zone.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 6. Overlapped section of variance and structural dip (dip illumination) showing discontinuous features in the candidate zones.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 7. Overlapped section obtained by color addition of red, green, and blue (RGB) components for variance (red), K1 curvature (green), and K2 curvature (blue).

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 8. Flowchart of the process used to establish the pseudo-acoustic impedance model by three sequential inversion methods to obtain acoustic impedance (AI) and estimate shear impedance (SI).

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 9. Inverted impedance section at the crossline 242 including the location of well UBGH 2-2-2 from the band-limited inversion.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig 10. Improved inverted impedance section at the crossline 242 including the location of well UBGH 2-2-2 from the sparse-spike inversion.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Fig 11. (a) Inverted section from the sparse-spike inversion without well log data and (b) inverted section from conventional model-based inversion along with the result of a comparison with the acoustic impedance log of UBGH 2-2-2. The inverted model from this methodology was verified by high correlation values with the acoustic impedance log (0.85) and the inverted model from model-based inversion (0.92).

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

EP

Fig 12. (a) Sections of original amplitude, instantaneous frequency, and acoustic impedance at Zone 1, (b) those at Zone 3, and (c) the comparison result of the multi-attribute analysis and impedance inversion between Zone 1 and Zone 3. Hydrate-bearing layers can be estimated

.

AC C

as having a higher positive amplitude over the bottom simulating reflector (BSR), higher instantaneous frequency and acoustic impedance (AI)

The layers showing these features in Zone 3 are mostly thicker than those in Zone 1. Zone 3 is more suitable for the first well site to obtain well log data than any other candidate zone.

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

Fig 13. Sections of the (a) acoustic impedance (AI) and (b) shear impedance (SI) model at the inline 245 from the simultaneous inversio n by comparison with the AI log of well UBGH1-09 and the inverted result from the sparse-spike inversion. In these models, hydrate-bea ring layers can be estimated as the layers with higher AI and SI.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

Fig 14. (a) Result of the classification by K-mean clustering that indicated four cl asses: Class 1 (green), Class 2 (blue), Class 3 (pink), and unclassified residuals (y ellow). (b) Class 3 with seismic data including bottom simulating reflector (BSR)

RI PT

and enhanced reflection, and (c) Class 3 with an inverted acoustic impedance (AI)

model that corresponds to a higher AI than surrounding lithologies. Class 3 coul

AC C

EP

TE D

M AN U

SC

d be used to estimate hydrated layers.

2

ACCEPTED MANUSCRIPT Highlights 1. To complement the preexisted interpretation results using BSRs, we performed multi-attribute analysis for delineation of gas hydrate reservoirs in the Ulleung Basin. 2. To solve the absence of well log data, we conducted three-step sequential impedance inversion to estimate AI and SI.

RI PT

3. To identify the reservoir distribution objectively and quantitatively, we applied the method of unsupervised clustering to the attributes.

AC C

EP

TE D

M AN U

SC

4. The proposed workflow helps to identify the reservoir distribution and aid in first well-site location determination.