Delineation of shallow channel geometry and infill lithology using Spectral decomposition 1 and seismic attributes : A case study from the North Sea Basin , Netherlands 2

9 In the Southern North Sea, 3D seismic data had been widely acquired to explore for hydrocarbons, but 10 interpretations of these datasets until now focus mainly on the deep exploration targets of the petroleum 11 companies. Less attention is given to shallow sediments. But these sediments often contain channels that can 12 serve as potential reservoir units. Thus the mapping and identification of these shallow channels and defining 13 their infill lithology is important. In this study, seismic spectral decomposition technique has been used to 14 delineate shallow thin channel geometry in a 3D seismic data acquired in the Dutch sector of the North Sea. The 15 concurrent interpretation of curvature and coherence cubes with seismic facies analysis based on reflection 16 terminations and geometry, amplitude and continuity enables the discrimination between shale versus sand filled 17 channels. The results of the spectral decomposition show two distinct low sinuosity channel features in NNE18 SSW direction but becomes diffuse towards the North. The strong negative curvature anomaly along the 19 channels’s axes observed in the most negative curvature attribute implies that the sediments within the channels 20 have undergone more compaction. These strong negative curvature anomalies are interpreted to be due to 21 differential compaction of shale filled channels. 22 23


Introduction
In the Southern North Sea, an extensive 3D seismic data had been acquired to explore for oil and gas in the Upper-Jurassic and Lower Cretaceous by the petroleum companies.But for most 3D seismic datasets, interpretations are focused mainly on the deep exploration targets of the petroleum companies.Less attention is given to shallow (younger) sediments.But these shallow (younger) sediments often contain channels that could serve as potential reservoir units.The mapping and identification of these shallow channels and defining their infill lithology in a fluvio-deltaic system is thus important in exploration and production.This is because these channels could serve as tools for shallow geohazard analysis (Selvage et al., 2012;Bouanga et al., 2014).For instance, in a gas filled shallow stratigraphic trap, the gas can be a hazard and a risk when drilling a borehole (Schroot and Schuttenhelm, 2003;Stuart and Huuse, 2012;Qayyum et al., 2013).Additionally, the occurrence of a shallow gas filled stratigraphic trap can be an indication for the presence of deeper hydrocarbon reserves, and thus an exploration tool.Besides, some of the shallow gas fields are even big enough to be considered as commercial gas fields (Stuart and Huuse, 2012).
Channels are generally visually close to or below seismic resolution (Caldwell et al., 1997;Tetyukhina et al., 2010), so thin to their surrounding geometry that their subtleties are nearly invisible in traditional seismic data.Thus, delineating thin reservoir sands from conventional seismic data had always been a challenge.Recent innovations such as coherence technology (Marfurt et al., 1998) and other edge sensitive attributes (Luo et al., 1996) are common methods employed in mapping boundaries of these geological subtle targets (channels).
Although coherence images and edge sensitive attributes reveal channels edges, a key limitation in these techniques is that they cannot delineate the channel's thickness (Chopra and Murfurt, 2006).
Spectral decomposition is a recent seismic interpretation technology that reveals otherwise hidden geological information and thus is being used extensively as an excellent tool for mapping channels (Partyka et al., 1999).In spectral decomposition, reflection from a thin bed has a peculiar expression in the frequency domain that gives an indication of the temporal bed thickness.It is a powerful seismic imaging and mapping tool that provides the interpreter useful quantitative information for determining bed thickness (Partyka et al., 1999), visualization of stratigraphy (Marfurt and Kirlin, 2001) and detection of hydrocarbon (Castagna et al., 2003;Sinha et al., 2005) to a level that was previously impossible.Spectral decomposition is also an effective tool in enhancing geohazard analysis as it is sensitive to wavelet, reflectivity, tuning and attenuation changes (Selvage et al., 2012).In spectral decomposition, the seismic data is converted from the time domain to the frequency domain and decomposed into frequency components.Studying the individual frequency components and comparing their responses provides significant insight into the subsurface geology.The time-frequency mapping process is a non-unique process; as a result, there are several methods for carrying out time-frequency analysis of non-stationary signals.Popularly used spectral decomposition methods include Fast Fourier Transform (FFT), Continuous Wavelet Transform (CWT), S-transform (ST), and Matching Pursuit decomposition (MPD).It is important to note that each method has its strengths and weaknesses (Chakaborty and Okaya, 1995;Leppart et al., 2010).
Several published works have discussed thoroughly various aspects of the post-stack time migrated 3D seismic dataset provided by dGB Earth Sciences in the F3 block.Some of these studies carried out in the study area are on delineation of geological features using spectral decomposition (Erlangga et al. 2013) spectral analysis (Honorio et al. 2014), porosity prediction from seismic inversion (Mojeddifar et al., 2015) etc.
Although these studies and others provide very rich literature on various aspects of the dataset, studies on prediction of the infill lithology of fluvio-deltaic channels particularly in the shallow (younger) sediments distinguishing channels filled with sand from those filled with shale are lacking.In this study, we have used attribute-assisted interpretation workflow to study the shallow channel geometry and infill lithology.We have applied FFT and CWT methods of spectral decomposition and seismic attributes (coherence and curvature) to a 3D seismic data set acquired in the upper Cenozoic fluvio-deltaic system in the block F3 in the North Sea basin to delineate shallow thin channel geometry and distinguish between intrachannel shale versus sand lithologies.
The intrachannel lithologies predicted using the seismic attributes were validated using well logs available in the area in conjunction with other lines of evidence.

Geological Setting
The study area (F3 block) is located in the Dutch sector of the North Sea.Much of the entire North Sea region in the Cenozoic era was a thermally subsiding epicontinental basin that was confined by land masses (Sorensen et al., 1997).Sedimentation rates during the Neogene outpaced the subsidence rate, resulting in rapid deposition and shallowing of the North Sea basin.An extensive fluvio-deltaic system (Eridanos delta) prevailed in the basin during the late Cenozoic Period (Ziegler, 1990;Overeem et al., 2001), draining the Fennoscandian High and the Baltic Shield.According to Overeem et al., (2001), the Eridanos drainage developed due to the Neogene uplift of the Fennoscandian Shield and accelerated subsidence of the North Sea Basin that occurred at the same time.The drainage system (Eridanos delta) started when the Scandinavian Shield was uplifted during the Oligocene (Rohrman et al., 1995).Nonlin.Processes Geophys.Discuss., https://doi.org/10.5194/npg-2017-62Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 27 October 2017 c Author(s) 2017.CC BY 4.0 License.Sales (1992) reported that the uplift rate increased during the late Miocene and also in the early Pliocene (Ghazi, 1996).Due to the late Miocene uplift, high sediment influx filled the northern offshore regions of the Dutch Sector.The increasing sediment load resulted in a differential load throughout the region.As a result, the buried Permian Zechstein salt started moving in the region forming several localized unconformities within the Pliocene interval that are underlain by salt domes.
The Cenozoic succession consists of two main packages, separated by the Mid-Miocene Unconformity (Steeghs et al., 2000) (Fig. 1).The lower package consists predominantly of relative fine-grained gradational Paleogene sediments (Steeghs et al., 2000), while the package above the unconformity is largely a progradational deltaic sequence that are made up of coarse Neogene sediments.The package above the unconformity can be subdivided into three sequences (Units 1, 2, and 3) corresponding to three phases of delta evolution (Fig. 1).Generally, in this package, conspicuous large-scale sigmoidal bedding pattern, downlap, toplap, onlap and truncation structures are observed.The base of Unit 2 is the zone of interest for this study.
Unit 2 is the delta foreset with a coarsening upward sequence (Tetyukhina et al., 2010) and the age of this unit is estimated to be Early Pliocene.

Spectral Decomposition Methods
The thin-bed tuning effect is the reason for the application of spectral decomposition method on a seismic data.The thin-bed tuning effect occurs when reflections from top and down layers have a constructive interference.In this instance, the peak amplitude response will occur at ¼ wavelength of the dominant period and layer thickness less than this value will not be detected in the seismic section (Saadatinejad et al., 2011).Laughlin et al. (2002) illustrated the relationship between tuning thickness and frequency using a wedge synthetic seismic model (Fig. 3).Fig. 3A shows a thickness increase from 0 to 30m at the left side, while the right side indicates amplitude tuning in three different frequencies.Fig. 3B shows the basis of the spectral decomposition technique, high frequency (36Hz, green colour) delineates the thinner part of the paleo channel and the low frequency (15Hz, red colour) shows the thicker part (Chopra and Marfurt, 2007).

Fast Fourier Transform (FFT)
The Fourier transform ) ( F  of a time-domain seismogram f(t) is expressed mathematically as: where t is time.Although a non-stationary signal when converted into the frequency domain via the Fourier transform method gives the overall frequency behaviour of the signal; such a transformation is not adequate for analysing seismic data (a non-stationary signal), whose frequency content is not constant but varies with time.
By taking short segments of the signal which are considered stationary parts (i.e windowing the signal) and then performing the Fourier transform for each segment, provides the frequency content of the signal at that time period (Chakraborty and Okaya, 1995;Zabihi and Siahkoohi, 2006).When this time window is shifted appropriately, it is possible to extract the frequency content of the signal and thus produce a 2-D representation of frequencies versus time.This 2-D representation is commonly known as short-time Fourier transform (STFT).The implementation of FFT is based on Short Window Discrete/Fast Fourier Transform.
The STFT is given by the inner product of the signal

 
where the window function  is centered at time t = τ, with τ being the translation parameter, and   is the complex conjugate of  .

Continuous Wavelet Transform (CWT)
The continuous wavelet transform (CWT) introduced by Morlet et al., (1982) is another method used to analyze the time-frequency content of a signal.Unlike the STFT where the window function has a fixed length, the CWT uses a variable window length.If the length of the interval on which the window function is non zero increases, the time resolution decreases, and the frequency resolution increases.On the other hand, when the length of the interval decreases, the time resolution increases and the frequency resolution decreases.The foregoing means that by increasing the frequency resolution, the time resolution will decrease and vice versa (Mallet, 1999).
The wavelet transform consists of wavelets which are functions defined as


, that have zero mean, which is localized in both time and frequency (Sinha et al., 2005).Each wavelet basis is generated by dilating and translating a two parameter function known as the mother wavelet, ) (t  .Given a wavelet basis, we can represent all functions in the basis by translations and scalings of the mother wavelet, ,  and  are the dilation or scale and translation parameters.In Eq. 3, as the value of  increases, the wavelet is compressed, its spectrum dilates and the peak frequency shifts to a higher value. Conversely, as the wavelet is scaled such that it dilates, the value of  decreases, its spectrum is compressed and the peak frequency shifts to a lower value (Chopra and Marfurt, 2015).The CWT is defined as the inner product of the family of wavelets

Materials and Method
The dataset used in this study is a post-stack time migrated open source 3D seismic dataset that was acquired in the F3 block covering an area of approximately 16 × 23 km 2 to explore for oil and gas in the Upper Jurrassic and Lower Cretaceous made publicly available by dGB Earth Sciences through Opendtect share seismic data repository (Aminzadeh, and de Groot, 2006).The F3 is a block is located in the North-eastern part of the Dutch sector of the North Sea.The 3D seismic data is made up of 650 inlines and 950 crosslines with a line spacing of 25m in both inline and crossline direction.The sampling rate is 4ms with a total data length of about 1.8s.
Figure 2 shows a vertical seismic section (in line 250) with gamma ray logs overlayed on the section in the respective well locations (F02-1, F06-1, F03-2 and F03-4).The study area is the horizon (in red) at the base of Unit 2 located between 800ms and 1100ms (Fig. 2).Because the original F3 dataset is noisy, only Dip-steered Median Filtered dataset was used as input data in this study.Since this study is aimed at delineating channels and to distinguish between sand-fill from shale-fill channel system in the study area, time slices of the 3D seismic volume was carried out between 1000 ms and 1055 ms, and the geological features in each time slice was analysed.In time slice 1007 ms (Fig. 4a), a channel like feature having a NNE-SSW pattern was observed, and extended towards the southern part in time slices 1028, 1037 and 1055 ms respectively (Fig. 4c, d and d).At time slice 1055ms, only a small part of the feature was observed and become indiscernible in time slice 1064 ms (Fig. 4d).From the above, it implies that a channel system existed between 1007 and 1055 ms which shows maximum prominence at 1028 and 1036 ms.
Based on this, the horizon shown in red was picked for analysis (Fig. 5A and B).
We also analysed the amplitude spectrum of the data (Fig. 5C).We define three dominant frequencies from the amplitude spectrum for RGB multi-colour display.The three frequencies were chosen such that they represent the low (28 Hz), middle (42 Hz) and high (60 Hz) frequency of the seismic bandwidth around the horizon.The three frequencies were then output as a single RGB blended full colour image.This is important because mixing outputs of different frequencies enables us to analyse results that depict different geological features related to different geometrical scales simultaneously i.e higher frequencies reveal features of more detailed character, whereas lower frequencies those which are more coarse.In using the CWT technique for spectral decomposition, the choice of the wavelet is important as it affects the output result (Castagna and Sun, 2006;Chopra and Marfurt, 2015).As a result, we tested all three mother wavelets i.e Morlet, Gaussian and Mexican-Hat wavelets to determine which wavelet would give better resolution.Figure 6 shows comparison of the resolution of the channel features of the three mother wavelets at 42 Hz frequency.Although the resolution of the Mexicat-Hat and Gaussian mother wavelets appears to be similar (Fig. 6), the Gaussian wavelet resolution of the channel features is slightly superior to that of the Mexican-Hat wavelet.Hence the Gaussian wavelet was chosen as the ideal mother wavelet for this data.In carrying out the spectral decomposition, we applied the low, middle and high frequencies on both the FFT and CWT (Fig. 7).In the spectral images shown in Fig. 7, we observed that some parts of the channel are resolved better at low frequency, some at mid frequency and others at higher frequency.However, at 42 Hz frequency, a clearer image of the channel is observed.This frequency at which the channel geometry is clearly pronounced is the tuning frequency.Since different parts of the channel features are resolved at different frequencies, it was considered that the channel geometry can be obtained in its complete form when the different frequencies are stacked together.Thus, a stacked frequency volume was obtained by summing up the 28Hz frequency (low frequency), 42Hz frequency (mid frequency) and 60 Hz frequency (high frequency) volumes (Fig. 8A).We used RGB colour-blending technique to display the multiple spectral components in a single 'full colour' image.Figure 8B displays RGB colour blending of the frequencies 28 Hz (red), 42 Hz (green) and 60 Hz (blue).In addition to FFT and CWT, we also generated coherence attribute image.The coherence attribute also known as similarity is an effective tool in determining lateral changes in the waveform and enables the mapping of lateral changes that might occur in stratigraphy (Mai et al., 2009).
Sandstones and claystones are the main infill rock types of paleo channels (Cao et al., 2015).Each of these lithotypes has different responses to compaction, hence the lateral changes that accompany compaction of a rock volume as it lithifies in a channel depends on the type of lithology (Torrado et al., 2014;Chopra and Marfurt, 2012).Because shales compact faster than sand (Torrado et al., 2014), a channel filled with sand and suspended in a shaly matrix may look like a mound or 'structural' high while a channel that is filled with shale 311 in a sandier interfluve may appear as structural low (Chopra and Marfurt, 2007a;Chopra and Marfurt, 2007b;312 Chopra and Marfurt, 2012).comparison is to enable verify potential differences between the methods.In Fig. 7, two significant distinct channel features are observed, particularly by the frequency of 42 Hz and are shown using red circles.The general direction of the channels observed is NNE-SSW and exhibit low sinuosity.Though the channels shape and low sinuosity are revealed by both algorithms, the CWT algorithm enhances the channel features better, the FFT results were considered rather poor.We attribute this behaviour to the length of the time window used for the FFT decomposition process.This is because the length of the time window used for the FFT decomposition is of paramount importance and the output of the FFT decomposition is always dependent on this characteristic (Kwietniak et al., 2016).
In using different frequencies, features of different scales are separated.High frequencies delineate smaller features, whereas low frequencies delineate more coarse structures.This effect is visible by comparing Figs.7 a,     b and c.Fig. 7e, in its left southern part, there is a more pronounced sign of the channel continuity (indicated with red circle) that was barely visible with lower frequencies especially with the FFT algorithm.We believe that the channel's width and thickness, are smaller in this part than those which are visible for lower frequencies.Additionally, the channel on the right (Fig. 7f) appears to have a branch towards the southern part (indicated with a red circle).This feature was not observed in the FFT algorithm and was barely visible even in the CWT at lower frequencies, but become pronounced at higher frequencies.This suggests that the main channel is relatively thicker than its branch.
Mixing outputs of different frequencies enables us to analyse results that depict different geological features related to different geometrical scales simultaneously i.e higher frequencies reveal features of more detailed character, whereas lower frequencies reveal those which are more coarse.Fig. 8 shows a RGB colour blended full colour image using 28 Hz (in red), 42 Hz (in green) and 60Hz (in blue).The RGB colour blending also  indicates differential compaction and corresponds to the incised channels in the coherence image, consistent with the observation of Chopra and Marfurt, (2012) and Torrado et al., (2014), who reported that differential compaction can also be seen in seismic amplitude section.We have used the coherence attribute to enhance channel edges (Fig. 10A).Fig. 10 B and C show the most-positive and the most-negative curvature volumes computed from the picked horizon.The most-positive curvature delineates the likely levees and flanks of the channels, while the most-negative curvature indicates the channel axis (thalwegs) Chopra and Marfurt, (2012).
Figure 10C shows a strong negative curvature anomaly along the channel axis (blue), implying that sediments within the channels have undergone more compaction.In the coherence image shown in Fig. 10A, the channels become progressively thinner towards the North, eventually fall below tuning and become indistinct.Corendering (i.e simultaneously displaying) coherence and the most negative curvature (Fig. 11A) shows how the curvature and coherence attributes can be used in a complimentary manner, in that in the lower part of the Figure where the channels are relatively thick, the coherence edges and channel axises mimic each other.That is, the curvature anomalies correlate to the channel geometry delineated in the coherence image and one can trace the channel geometry further on the most-negative curvature image even though the coherence channel edge anomalies become diffuse or disappear towards the North.Figure 11B shows co-rendered most-positive (red) and most-negative curvature (blue).Fig. 11 shows a strong negative curvature anomaly along its axis (blue).This implies that sediments within the channel have undergone more compaction (Chopra and Marfurt, 2012;Torrado et al., (2014).These strong negative curvature anomalies are interpreted to be due to differential compaction of shale filled channel.Levees and channel edges are delineated as mounds or ridges and give rise to strong positive curvature anomalies (red).In order to validate our prediction of the infill lithology, we overlay the gamma ray logs of the respective wells on the seismic section (Fig. 2) with the horizon of interest shown in red.This horizon is at the bottom of Unit 2.
According to Fig. 2, the horizon of interest is at an interval rich in sediments with gamma ray values generally greater than 70 API (Luthi, 2001).This is consistent with the observation of Tetyukhina et al., (2010) who reported that shaly sediments with higher impedance values were deposited mainly at the bottomset of the clinoform system in the F3 block of the Southern North Sea basin.

Conclusions
The delineation of subtle geological features such as channels has always been a challenge.This is because channels are subseismic, so thin to their surrounding geometry that their subtleties are nearly invisible in traditional seismic data.But these geobodies are important because they serve as hydrocarbon traps.In an attempt to delineate shallow thin sand reservoirs in the F3 block, we have carried out spectral decomposition analysis using the Fast Fourier Transform (FFT) and Continuous Wavelet Transform (CWT) on a 3D seismic data acquired in the F3 block.We assessed the relative performance of the FFT and CWT on the data.We used a red-green-blue (RGB) colour-blending technique to display the composite full colour image to enhance better resolution of the channels features.In order to determine the infill lithology of these channels, we have also used the coherence and curvature attributes in a complementary manner.While the coherence attribute was used to enhance channel edges, the curvature attribute was used to discriminate between intrachannel shale versus sand lithologies based on differential compaction of the channel relative to its edges.
The results show two distinct almost linear channels trending in the NNE-SSW direction.The CWT algorithm remarkably delineated the channel geometries in a much better way than the FFT algorithm.While the most positive curvature defines the likely levees and overbank deposits and as well as the flanks of the channel, the most negative curvature defines the channel axis.The curvature anomalies also correlate to the channels geometry obtained from the coherence attribute.The strong negative curvature along the axis of the channel is due to differential compaction of a channel filled likely with shale.

Fig. 3 .
Fig. 3. Relation between tuning thickness and frequency (Lui, 2006) and (B) results of spectral decomposition at 36Hz, 15Hz, maps and channel thickness detected with variable frequency (Chapra and Marfurt, 2007).Low -frequency slices indicate high thickness (green lines) and high frequency slices indicate low thickness (red lines).

Where.
is the complex conjugate of  and ) , (   W F is the time-scale map (scalogram).At each scale (i.e., for each value of ) the kernel wavelet is scaled by a factor 1/  and translated by  to produce the wavelet Nonlin.Processes Geophys.Discuss., https://doi.org/10.5194/npg-2017-62Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 27 October 2017 c Author(s) 2017.CC BY 4.0 License.Useful wavelets commonly used in wavelet transform are Morlet, Gaussian and Mexicat-Hat.

Fig. 4 .
Fig. 4. Time slices (a) 1007ms (b) 1028 ms (c) 1037 ms (d) 1055 ms (e) 1064ms (red oval shape indicate the channel like feature.This feature is not delineated in time slice 1064ms Fig. 5. (A) Display of inline 250 and the mapped seismic horizon shown in red (B) colour blended display of the mapped horizon shown in A (C) Amplitude Spectrum of inline 250 showing the seismic bandwidth Nonlin.Processes Geophys.Discuss., https://doi.org/10.5194/npg-2017-62Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 27 October 2017 c Author(s) 2017.CC BY 4.0 License.

Fig. 6 .
Fig. 6.Comparison between the three mother wavelets at 42 Hz frequency.The Gaussian wavelet gives better resolution of the channel features (indicated with red arrows) than the Morlet and the Mexican Hat wavelets Nonlin.Processes Geophys.Discuss., https://doi.org/10.5194/npg-2017-62Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 27 October 2017 c Author(s) 2017.CC BY 4.0 License. 313

Fig. 9 .
Fig. 9. Display showing incised channels on a coherence slice and its seismic amplitude signature.The seismic signature of the incised channel is seen as a sag at the position of the yellow arrows.We interpret the sag over the channels to indicate that they contain more shale than the surrounding matrix.

Fig. 9
Fig.9shows a coherence attribute image showing the channel features through a seismic section perpendicular to the thalweg of the channels.The sag observed in the vertical seismic section (yellow arrows) Fig. 10.(A) Coherence image (B) Most positive curvature and (C) Most negative curvature, the yellow arrows delineate channel edges seen in the Coherence.Note that these channels can be followed further on the most negative curvature image.Notice the strong most-negative curvature anomaly along the channel axis (blue).We interpret the most negative curvature anomaly to be due to differential compaction over shale-filled channels.

Fig. 11 .
Fig. 11.(A) Co-rendered coherence and most negative curvature.(B) Co-rendered mostpositive and mostnegative curvature volumes where moderate curvature values are rendered transparent.Sediments within the channel have undergone more compaction and give rise to a strong negative curvature anomaly along its axis (blue).Levees and channel edges appear as ridges and give rise to strong positive curvature anomalies (red).