Geological Characterization of the 3D Seismic Record within the Gas Bearing Upper Miocene Sediments in the Northern Part of the Bjelovar Subdepression—Application of Amplitude Versus Offset Analysis and Artificial Neural Network

Abstract

As natural gas reserves are generally decreasing there is a need to successfully characterize potential research objects using geophysical data. Presented is a study of amplitude vs. offset, attribute and artificial neural network analysis on a research area of a small gas field with one well with commercial accumulations and two wells with only gas shows. The purpose of the research is to aid in future well planning and to distinguish the geophysical data in dry well areas with those from an economically viable well. The amplitude vs. offset analysis shows the lack of anomaly in the wells with only gas shows while the anomaly is present in the economically viable well. The artificial neural network analysis did not aid in the process of distinguishing the possible gas accumulation but it can point out the sedimentological and structural elements within the seismic volume.

1. Introduction

With a general decrease in the recoverable reserves of natural gas both onshore [1] and offshore [2], the need to locate new accumulations is more than evident. Exploration for hydrocarbon reservoirs relies on a very large number of input data for characterization of the subsurface. This work focuses on usage of attribute analysis on seismic volume, amplitude versus offset (AVO) analysis as well as searching for patterns within many extracted seismic attributes by self-organizing artificial neural networks (SOANNs). The seismic amplitude versus offset (AVO) analysis has become prominent as the direct hydrocarbon indicator in the last decade, aimed at characterizing the fluid content or the lithology of a possible reservoir and reducing the exploration drilling risk [3].

Seismic attribute analysis was performed on reprocessed and reinterpreted seismic volumes. Progress in characterization of hydrocarbon reservoirs is enabled by the analysis of seismic attributes. In saturated gas reservoirs the character of seismic signals changes. Results of the seismic attribute analysis indicate prospectives east and southeast from the discovered field. For the analysis of interception and gradient as primary AVO attributes, pre-stack data or common depth point (CDP) gathers data used by the method of summing traces from common depth points in the subsurface. Obtained crossplot of intercept and gradient can determine the appropriate AVO class that shows changes in amplitude relative to the offset.

SOANNs were used on the same dataset as the seismic attribute analysis which were extracted on the top of the targeted seismic horizon. Patterns indicate mostly sedimentological features that can be tied to the variability of the depositional environment.

2. Materials and Methods

2.1. Geological Settings of the Exploration Area

The Bjelovar Subdepression is located in the Croatian part of the Pannonian Basin (CPB), within the southwestern part of the Drava Depression (Figure 1). The total area of Drava Depression is approximately 12,000 km², of which 9100 km² is within the Republic of Croatia [4]. Geographically, it is surrounded by the Kalnik Mountain in the northwest, Bilogora Mountain in the northeast, Papuk Mountain, Ravna Gora Mountain and Psunj Mountain in the east and Moslavačka Gora Mountain in the south.

The total thickness of the Neogene and Quaternary deposits can be more than 7000 m in the middle part of the Drava Depression while in the Bjelovar Subdepression it can reach up to 3500 m thickness [5]. The analyzed gas field is located on the edge of the Bjelovar Subdepression (Figure 1) near the boundary with the Drava Depression (Croatian part of the Pannonian Basin System, CPBS). Three wells were drilled on the structure: Well-1, Well-2 and Well-3.

The formation of the Pannonian Basin System (PBS) began in Early Miocene when the Eurasian plate began to subduct under the Pannonian crustal fragment [8]. Numerous local transtensions as proto-areas of later depressions and subdepressions opened between the southern and northern boundaries of the PBS ([9,10]). The Croatian part of the PB (CPB) structural complex can be subdivided into three major megacycles [8,10] in relation to dominant tectonic events (Figure 2).

Major structural depressions, in the Croatian part of the Pannonian Basin are the Mura, Drava, Sava and Slavonija-Srijem depressions. The Bjelovar Subdepression is a part of the Drava Depression. Lithostratigraphy and regional geological relations are described in detail by several authors [8,11,12]. The focus of this research is located within Bilogora Fm. between well log markers Δ and D’ (Reservoir 1) and Reservoir 2 is located within Pepelana sandstones between well log markers Z and Δ (Figure 2).

2.2. Input Data

Seismic data reprocessing was done using up to date processing techniques and ended with Kirchhoff pre-stack time migration (PSTM, Figure 3a) flow. The input seismic data to the Kirchhoff pre-stack depth migration (PSDM) was the same set as for the time migration. The final pre-stack depth migration image resulted in better representation of the subsurface structures and stratigraphy, enabling greater interpretation confidence especially in the area of low signal to noise ratio.

PSDM shows enhancement in the image of the prospective area indicating anomalies present on the structural trap closure. Reliable and detailed analysis of thin deposits can only be performed on reprocessed and reinterpreted 3D volumes (Figure 3b). Location of the seismic section and structural maps of Reservoir 1 and Reservoir 2 are shown in Figure 3c,d.

2.3. Seismic Attribute and AVO Analysis

Seismic attributes are a compendium of all the measured, computed or observed information obtained from seismic data that can greatly contribute to better interpretation of the subsurface [14].

Any value derived from seismic data can be a seismic attribute using primary differentiated attributes, instantaneous amplitude, frequency and phase and all possible combinations are derived from them. By analyzing seismic attributes, detailed reservoir characterization and quantification can be obtained. The majority of attributes are derived from post-stack data, however there are a few pre-stack attributes. The first attributes developed were related to the 1D complex seismic trace and included: envelope amplitude, instantaneous phase, instantaneous frequency, and apparent polarity. Acoustic impedance obtained from seismic inversion can also be considered an attribute and was among the first developed [15].

Commonly used attributes in the exploration and detection process of possible hydrocarbons accumulations are variance, azimuth, dip, sweetness, spectral decomposition and AVO primary and secondary attributes.

The instantaneous amplitude, also known as reflection strength, is defined as the total energy of a seismic trace. High values of amplitudes are most often related to discordances, sudden lithological changes and boundaries caused by changes in the deposition environment. Instantaneous phase is a seismic attribute that gives equal weight to both strong and weak reflexes and is a good indicator of discontinuity, slope change and layer boundaries. Instantaneous frequency shows the time change of the current phase and is used in understanding the properties of reservoirs.

Sweetness is a seismic attribute used to highlight thick, clean reservoirs, along with hydrocarbons contained within [16]. Sweetness is calculated by dividing the instantaneous amplitude (amplitude envelope) by the square root of the instantaneous frequency. This attribute detects coarse-grained intervals. These intervals are highlights on the seismic section and are called sweet spots. They indicate the possible accumulation of gas or oil.

Amplitude versus offset (AVO) is primarily the variation in seismic reflection amplitude with change in distance between shot point and the receiver [3]. Simultaneous variations in amplitude are analyzed as a function of the incidence angle, also referred to as AVA (amplitude versus angle). By using linear regression or calculating the equivalent from seismic data, it is possible to calculate the line that describes how the amplitude changes with offset or angle using the intercept and the gradient as the primary attributes. Changes in amplitude can be displayed on their crossplot as single points.

Knott and Zoeppritz developed the theoretical work necessary for AVO theory [17,18], given the P-wave and S-wave velocities along with the densities of the two-bounding media. They developed equations for plane-wave reflection amplitudes as a function of incident angle.

Different authors later gave various simpler approximations of Zoeppritz’s equations. P. Richards and C. Frasier [19] expanded the terms for the reflection and transmission coefficients for a P-wave incident upon a solid–solid interface and simplified the result by assuming only small changes in elastic properties across the interface. This approximation was popularized in 1980 by K. Aki and P. Richards [20] and has since been commonly referred to as the Aki and Richards approximation.

Ostrander [21] was the first to introduce a practical application of the AVO effect, showing that a gas sand underlying a shale exhibited amplitude variation with offset.

Shuey [22] published a closed form approximation of the Zoeppritz equations. An approximation to the Shuey equation, which itself is an approximation of the Aki–Richards equation. It is essentially the same as the Shuey equation, with the simplification that the third term is ignored (Equations (1)–(3)),

R(Ө)R (0) + Gsin² Ө

(1)

where

$$\mathrm{R}\left(0\right)=\frac{1}{2}\left(\frac{\mathsf{\Delta }\mathrm{Vp}}{\mathrm{Vp}}+\frac{\mathsf{\Delta }\mathsf{\rho }}{\mathsf{\rho }}\right)$$

(2)

and

$$\mathrm{G}=\frac{1}{2}\frac{\mathsf{\Delta }\mathrm{Vp}}{\mathrm{Vp}}-2\frac{\mathrm{Vs}2}{\mathrm{Vp}2}(\frac{\mathsf{\Delta }\mathsf{\rho }}{\mathsf{\rho }}+2\frac{\mathsf{\Delta }\mathrm{Vs}}{\mathrm{Vs}})$$

(3)

where: R is the P-wave reflection amplitude, R(0) shows the zero-offset reflectivity, θ is the average of the incidence and transmission angles, Vp is the P-wave velocity, Vs is the S-wave velocity and ρ is the density of the reflector. ∆Vp is the P-wave velocity contrast across the interface, ∆Vs is the S-wave velocity contrast across the interface and ∆ρ is the density contrast across the interface.

Equation (1) is mostly used in reservoir exploration as one of the direct hydrocarbon indicators (DHI) and is known as Shuey’s two-term AVO equation. Data from the offset domain must be transformed into an angle domain to calculate the intercept and gradient using linear regression. The condition is that the equation is linear if the intercept is a function of the incidence angle. Amplitudes picked along a moveout-corrected event on a common-midpoint (CMP) gather plotted against sin²θ can be fitted to a straight line. The slope of the line gives the AVO gradient attribute and the ordinate at zero angle gives the AVO intercept attribute. This relationship is the one expressed by AVO crossplots, which usually show intercept R(0) (the zero-offset reflectivity) plotted against gradient G.

CDP gathers are transformed from the offset domain to the angle domain and analyzed to detect changes in amplitude relative to the offset because a sudden increase in amplitude may indicate gas saturation. Polarity changes indicating a particular AVO class are also observed. Intercept I (representing the segment on the y-axis showing the amplitude) and gradient G (showing the slope of the regression line) are calculated by the linear regression method. The diagram shows the amplitude versus angle of incidence sin²θ (Figure 4).

The generated attributes by the Shuey 2-terms method, intercept and gradient, are known as the primary AVO attributes of the linear approximation. From the primary attributes, secondary AVO attributes can also be calculated. Some secondary attributes are: product of intercept and gradient (indicates in a single volume the polarity and the angle dependence of the amplitude variations) and attribute Sign(intercept) × gradient (keeps the value of gradient but the result will be directly affected according to the combined values of intercept and gradient) which were used for the calculation in this case. Primary and secondary attributes were calculated from angle stacks (near angle stack; 0–10°, mid angle stack; 10–20°, far angle stack 20–30°). Angle stacks calculated from CMP gather with carefully selected angles for near, middle and far offsets.

AVO crossplots are a simple and elegant way of representing AVO data. A central concept that emerged from this work was the “fluid line,” or background trend, a hypothetical trend based on consideration of brine filled rock properties together with simplifications of the reflectivity equations.

If the intercept is plotted on the x-axis and the gradient on the y-axis, then for consolidated sand/shale rocks the top and base reflections form a trend from the upper left to the lower right quadrant of the crossplot that passes through the origin. When it was realized that data points for equivalent hydrocarbon-filled rocks plot to the left of this line, it became clear that normalizing the data against the fluid line might provide an optimum AVO indicator [24].

The classical works dedicated to AVO cross-plotting analysis [25,26] identified four types of gas sands based on their AVO characteristics. These can be summarized as follows (Figure 5):

Class I sands are characterized by higher acoustic impedances than those of the encasing shale. They are generally found in hard rocks at depth with lower porosities. They produce positive reflection coefficients, which decrease with angle.

Class II sands have almost the same impedance as the encasing shale. Hence, they have moderate positive or negative reflection coefficients, which results in strong negative coefficient at larger angle. In some circumstances this implies phase reversal between near and far angle.

Class III sands have lower impedances than the encasing shale. They are very common in the shallow subsurface. Class III reflection coefficients increase in absolute value with increasing angle of incidence.

Class IV (added by Castagna after Rutherford–Williams) sands may have identical normal incidence reflection coefficients as per Class III, but the magnitude of Class IV sand reflection coefficients decreases in absolute value with increasing angle of incidence.

2.4. Artificial Neural Networks

Artificial neural networks (ANNs) are tolls used widely in modern scientific purposes whilst their origin can be tracked to research [30] with numerous state of the art applications in geosciences [31,32,33]. Their employment in the last several decades in the task of reservoir engineering, enhanced oil recovery and subsurface characterization in general has seen large potential [34,35,36,37]. Fundamentally, there are two different types of ANN based on their learning principles. The first incorporates supervised learning when the ANN has a set of input data with desired outcomes on which they learn and form an algorithm. The second type are self-organizing artificial neural networks (SOANNs) with unsupervised learning principles where only an input data set is given, and the task of the network is to classify the data based on their similarity. One of the applications of SOANN is a cluster analysis using Kohonen training [38] as lithology prediction from various geophysical data [39,40,41]. Earlier research in similar stratigraphic intervals and geographical extent were limited to supervised ANNs for the purpose of lithology [14,42] and hydrocarbon saturation prediction [43]. Cluster analysis using SOANN for the purpose of this research was performed using StatSoft Statistica.

3. Results and Discussion

The sweetness seismic attribute is a very useful tool for proper description of the depositional environment, reservoir quality and lithofacies discrimination. It is very useful in fluvial systems for the identification of isolated sand bodies because they generate stronger and broader reflections than the surrounding shales. It becomes less useful in environments where sands and shales are highly interbedded or where contrasts in acoustic impedance between sands and shales are low. Sweetness improves the imaging of sand intervals or bodies and identifies oil and gas prone places called “sweet spots” [44].

On the RMS amplitude maps on the volume of sweetness attribute along the top of horizon of the Reservoir 1 (Figure 6), a small anomaly is visible in the area of the positive Well-1, while the other two negative wells are de facto in the area of lower amplitudes. The area southeast of Well-1 shows the location of highly increased amplitude which may possibly indicate gas saturation. The RMS amplitude map shown in the same way on the volume of the sweetness attribute along the top of Reservoir 2 (Figure 6) shows a positive Well-1 surrounded by a highly visible anomaly while the other two wells are in the low amplitude zone. The perspective area to the southeast is in the zone of higher amplitudes. It is possible that there is gas saturation but further consideration and other analyses are needed.

AVO analysis on the exploration area done on available pre-stack depth migration seismic volume. AVO attributes analyzed via crossplot intercept/gradient. Normally, amplitude decreases with offset. Increasing amplitude with distance from the source (offset) can be used to predict the gas saturation areas. The AVO anomaly is present at all three well locations and is also present to the south of Well-1. The anomaly can be seen in all four classes which may indicate possible gas saturation. A limitation of AVO analysis using only P-energy is its failure to yield a unique solution, so AVO results are prone to misinterpretation. One common misinterpretation is the failure to distinguish a gas-filled reservoir from a reservoir having only partial gas saturation (e.g., “fizz water “) [45]. It is nearly impossible to distinguish “fizz water sand units” from “commercial gas sand units by just using P-wave impedance or P-P reflectivity data. This is because “fizz water sand units” and “commercial gas sand units” have similar P-wave impedances (or the difference is so small that it is unresolvable seismically) [46].

The S-wave for all three wells was not recorded. In that case, sometimes, we can try using very far (beyond 30°) offset data. At far offset, seismic reflections are heavily influenced by not only the P-wave impedances, but also the S-wave impedances, the P-S impedance ratio, and density contrasts. As a result, very far offset seismic data could contain information on the amount of gas saturation, giving potential to distinguish “fizz water sand units” from “commercial gas sand units”. This of course depends on the quality of very far offset data and many other challenging factors [46].

On the seismic profile of the volume of the intercept (primary AVO attribute calculated from Angle stacks) through the Well-1 within the volume of interest (220 ms -XL888-XL948) the volume Sign(intercept) × gradient is overlapped where between the top of horizons of Reservoir 1 and Reservoir 2, inside the window clearly visible increase in amplitude in green, yellow and red, which indicates the possible presence of gas (Figure 7a). The Intercept/gradient crossplot shows a deviation from the background trend or AVO anomalies in all 4 classes (Figure 7b) On Well-1 along the top of horizons of Reservoirs 1 and 2. Furthermore, extracted AVO anomalies at one point shown in crossplot intercept/gradient as a pink triangle for Reservoir 1 which shows that the point falls within the third quadrant and shows AVO class III. The blue triangle shows Reservoir 2 also within the AVO class III.

The Figure 7c shows a map of RMS amplitudes on the volume Sign(intercept) × gradient between the top of horizons of reservoirs 1 and 2 with a calculated window of 10 ms above Reservoir 1 and below Reservoir 2. Furthermore, increasing amplitude of Well-1 is visible, confirming the presence of gas.

The same procedure was also performed for Well-2 and Well-3, the drilling of which did not confirm the presence of gas as in positive Well-1, as well as for potential area southeast of Well-1.

A window with a volume of interest of 180 ms (XL 916–950) was taken for Well-2 (Figure 8a). For Well-3 (Figure 9a), a window of 150 ms (XL 880–915), and for the anomaly area southeast of Well-1 (Figure 10a), a window with a size of 200 ms (XL 848–910). Within the volume of interest, an increase in amplitudes in green, yellow and red is shown.

Furthermore, the other two wells and the anomaly area on the southeast show the AVO anomaly at one point, the pink triangle for Reservoir 1 and the blue for Reservoir 2 located overall within the third quadrant, AVO class III (Figure 8b, Figure 9b and Figure 10b).

RMS amplitude maps at the Sign(intercept) × gradient volume between the tops of the horizons of both reservoirs it is evident that an anomaly exists, but around the two negative wells, not at their locations (Figure 8c and Figure 9c). To the southeast, in the perspective area, the anomaly is clearly visible (Figure 10c).

The AVO anomaly is also visible on two negative wells 2 and 3 in all 4 classes, but where class III which represents the classic bright spot is very visible. A limitation of AVO analysis using only P-energy is its failure to yield a unique solution, so AVO results are prone to misinterpretation. One common misinterpretation is the failure to distinguish a gas-filled reservoir from a reservoir having only partial gas saturation (“fizz water”). The S-wave was not recorded. Very small amounts of gas at two wells have the same Vp/Vs ratio as with 90% gas saturation, thus the AVO anomaly does not differ between 5 and 90% of saturation gas.

When there are no records of the S-wave, sometimes very far offset seismic data could give the potential to distinguish “fizz”-gas from “commercial gas sand units”. Far angle stacks are calculated with angles of 20–30 degrees. However, we also calculated a version with slightly larger angles of 24–36 degrees. The figure shows the extracted RMS amplitude on the far angle stack (with larger angles) volume along the top of Reservoir 1 (Figure 11a), and the figure shows the top of Reservoir 2 (Figure 11b) in the same way. Amplitude increases are present on both reservoirs on Well-1 as well as southeast in the potential area. Negative wells are located near anomalies. We cannot determine with accuracy that these anomalies do not show only very small amounts of gas and more detailed analyses should be made.

SOANN analysis was performed using TIBCO STATISTICA software on the top of Reservoir 1 horizon through the available 3D volume in order to recognize the patterns that could implicate the variation of lithology, gas accumulation anomalies or structural elements. Almost 170,000 data points were extracted from the surface which relates to Reservoir 1 top with associate seismic attribute values of envelope, instantaneous frequency, sweetness, RMS amplitudes, chaos and cosine of phase. Two outcomes with topological cluster maps of 2 × 2 and 3 × 3 dimensions were made to visualize the features on top of the Reservoir 1 horizon. Limit for the learning process was set for 1000 iterations using Kohonen training. Learning process results present relatively small amounts of errors in the classification of seismic data within the categories (

Table 1. Self-organizing neural network learning process result.

Neural Network Architecture	Training Error (%)	Test Error (%)
2 × 2 SOANN	5.2	8.4
3 × 3 SOANN	3.4	4.8

).

SOANN with 2 × 2 neural network architecture/topological map resulted in data being divided into 4 classes and 3 × 3 architecture. Data points with associated classes were imported into petrel where there were upscaled within a simple model representing the top of the Reservoir 1 horizon. All three wells fell in the same class category which pointed out that this kind of approach can hardly distinguish gas presence within the area of the research (Figure 12). However, structural features such as faults (linear features of category distribution) and sedimentological features can be observed which presents a possible path for further research for SAONN approaches within 3D seismic volumes.

4. Conclusions

AVO analysis on the exploration area done on available pre-stack depth migration seismic volume. AVO attributes analyzed via crossplot intercept/gradient. Normally, amplitude decreases with offset. Increasing amplitude with distance from the source (offset) can be used to predict the gas saturation areas. The AVO anomaly is present at all three well locations and is also present southeast of Well-1. The data were extracted at one point for each reservoir and each well, and are located within the third quadrant on the crossplot intercept/gradient and show AVO class III.

AVO class III shows an increase in amplitude with distance and shows a classic bright spot. This also means that there is a large difference in acoustic impedances between the reflectors. The reason may be sudden lithological changes as well as very small amounts of gas.

Only a small percent of the gas will cause a significant drop in the effective fluid modulus, and consequently a significant drop in the bulk modulus in the rock. Therefore, P-wave seismic can usually only determine the presence of gas, but not saturation.

The S-wave was not recorded for any wells. Therefore, we cannot confirm the amount of saturation with certainty. However, other methods such as the analysis of the seismic attribute sweetness and the extracted RMS amplitude on the far angle stack can help us. Two negative wells are located outside or on the edge of anomalies. Further studies will confirm whether commercial gas saturation exists within the larger anomaly.

Distinguishing “fizz water” from “commercial gas” saturations with surface seismic and borehole data remains a challenging problem. The solution and help can be in obtaining and using quality S-wave data, scaling and attenuation in data processing and analysis of other seismic attributes.

Although SOANN did not yield the desired result of distinguishing the locations of possible gas accumulations, the structural and sedimentological features which can be observed on top of Reservoir 1 horizon present a potential for future SOANN application.