Abstract
The concept of a transition zone, known as Wolf ramp, was incorporated into the seismic interpretation of a 3D seismic survey situated within the Baltic Basin (Northern Poland). Within the survey area, there exists one formation, the Pasłęk Formation, (Lower Silurian—Llandovery), that exhibits a linear change of velocity. This characteristic—linear change of velocity—causes a reflection coefficient (i.e., seismic amplitude) produced at such a boundary to be frequency dependent. The Pasłęk Formation was considered to be a potential shale gas reservoir and it was necessary to determine its structural position and thickness. The formation is challenging for robust seismic interpretation on the migrated seismic section—it does not manifest a stable reflection coefficient, and the amplitude contrast associated with the borders of the formation is low. There is no impedance contrast that would produce a reflection of high amplitude at the top or base of the formation which excludes determination of the formation thickness, hence the estimation of reservoir volume. Within a 3D dataset, there exists only one well with complete logs that were used for the analysis. The Pasłęk Formation is a flat-lying layer that continues itself far beyond the 3D survey. It is present in wells in the vicinity of the study area. These wells lay within other 3D or 2D datasets, but the quality of the seismic is poor, and similar seismic analysis is not possible. Nevertheless, these wells were incorporated in the research to reason about the possible link between the existence of transition zone and mineral content. The method used for recognition of transition zone is spectral decomposition and spectral analyses. The integrated studies enabled us to find a link between the Wolf ramp and mudstone-claystone interval of the Silurian age and give a new example of a transition zone which is present in shale plays. The transition zone concept might be applied for shale plays identification and analysis.
1. Introduction
The concept of a transition zone, a Wolf ramp, is relatively old—the original work of Alfred Wolf describing a ramp comes from the second volume of Geophysics [1]. The transition zone is defined as a layer of a given thickness that separates two half-spaces of different and constant velocity values: (upper layer) and (lower layer). The relation between velocities is not essential. The transition zone is characterized by a linear change of velocity with depth, hence the name. In this model, the density is neglected [1]. The most important characteristic of this layer is that such a sequence produces a reflection coefficient that is a function of the frequency of the elastic wave for a zero-offset seismic trace.
The paradigms were extended to a zone of linearly changing velocity and density within 60-80. Works of Gupta [2] give solid analytical solutions for variations of density and velocity for elastic waves in liquid media. The idea of the transition zone is applicable in many geological settings: Justice and Zuba [3] used it for permafrost analysis. In the paper, they presented 1D modeling with the use of convolution to find a seismic signature link to permafrost analysis. These models show that the transition zone is frequency dependent and that permafrost resembles such a characteristic.
Recently, the concept was rediscovered and described with a modest approach by Liner and Bodmann [4] where authors use the transition zone model for interpretation of sea bed. This paper contains a repeated analytical solution of Wolf. In the paper, the authors present concept of applying spectral decomposition for indicating reflectivity dispersion present for normal-incidence reflection of a linear velocity transition zone. Many different factors may be the reason for linearly changing velocity with depth; one of them is gas saturation. Gómez and Ravazzoli [5] studied this relation to carbon dioxide transition layers. The work gives the application of the theory for CO2 storage. They started with Wolf solution and added factors related to fluid properties—bulk density, bulk modulus. The numerical solution lays in accordance with the analytical solution, and the effect of carbon dioxide saturation influences the model. This technique enables more precise monitoring of carbon dioxide sequestration sites.
The model of a transition zone and its attenuating properties was the subject of a research project [6] for indicating shale-bearing deposits within the Polish shale gas concessions. The Pasłęk Formation was one of several prospecting horizons for potential shale gas exploration. However, the formation was challenging to map with the use of the migrated seismic section. The problem was due to its elastic properties and the reflectivity series—the top and bed of the formation cannot be easily identified on the seismic section. Hence, there was a need to incorporate other methods for seismic interpretation of this interval. It was proposed that the interval under analysis can be considered to be a transition zone, and that frequency analysis might enable its interpretation. The main focus of the research was to apply different spectral decomposition algorithms so that the interpretation of the transition zone would be possible. In this article, we focus mainly on the applicability of the transition zone for seismic interpretation of shale reservoir. We will show a new example of the transition zone and evaluate spectral decomposition algorithms and their performance.
2. Materials
2.1. Theory of a Transition Zone
The transition zone is a layer that separates two intervals of a given velocity, and its velocity characteristic is described by the velocity of the upper layer, lower layer, and distance between them. The simplified model of a transition zone is presented in Figure 1.
In Figure 1, the layers are indicated by the velocity of the first layer, , the velocity of the second layer, , and the transition zone of the defined thickness . The density changes are neglected for the initial definition. The relation between and is negligible, i.e., can be lower or higher than . By the definition of only three parameters, the three-layer geological model is created. The critical property of the transition layer is the fact that the reflection coefficient produced in the seismic reflection process is frequency dependent. The reflection coefficient of the transition zone layer is defined as [1]
In Equation (1), a reflection coefficient for absolute amplitude is opposite to the displacement reflection coefficient [5]. The precise derivation of the reflection coefficient is presented in many papers, starting with the original paper by Wolf [1], but also more current [4, 7]. In Equation (1), is a velocity ratio of the upper and lower part of the transition zone, frequency is in hertz, transition zone thickness, in meters, and and depend on frequency by
Parameter is defined according to the upper-layer velocity .
Such defined reflection coefficient (Equation (1)) is a function of frequency, which implies that in respect of the convolutional model of seismic trace, the same is true for seismic amplitude value. For example, for a transition zone of thickness 50 meters, the maximal spectral energy will concentrate for frequencies below 20 Hz (Figure 2).
The function presented in Figure 2 depicts the behaviour of the absolute seismic amplitude for a case when is twice bigger than . For other cases, the exact solution will be different, but two of properties of the function will be kept: (1) the existence of periodicity, i.e., notch width and peak position, here of 30 Hz, and (2) the decaying value of spectral amplitude. The exact notch position is specific for the chosen initial parameters. With the different thickness and velocity ratio, the function will keep the periodicity, but the notch width will vary [5].
For such a defined transition zone, it is, therefore, crucial to analyze low frequency range of the seismic data. For this reason, spectral analysis is an accurate tool for transition zone interpretation: given the relative thickness of a layer (geological information from the well), one can model the exact behaviour of the amplitude versus frequency response by incorporating Equation (1) and analysis of the frequency maxima (Figure 2). Spectral decomposition volume related to such frequency maxima should reveal the layer under analysis and enable its thickness estimation.
2.2. Motivation for Application of Wolf Ramp Concept
The convolutional model of a seismic trace says that the observed seismic amplitude is relative to reflection coefficient. Assuming a stationary, zero-phase wavelet, the change in seismic amplitude is proportional to the reflection coefficient. The reflection coefficient that is a step function will produce distinguishable and distinct peak (Figures 3(a) and 3(b)), even in the presence of a tuning effect (Figures 3(c) and 3(d)).
Every deviation from a step function that is associated with velocity and/or density gradients will result in a modified seismic reflection. Such examples are shown in Figures 3(e)–3(n), and the resulting reflection coefficient might be extended in time (Figures 3(e)–3(g)) or, if the thickness of a transition zone is lower, the resulting reflection will be rotated in phase and scaled in terms of amplitude values. The exact shape of reflection will be governed by the model and interference of the wavelet. The classical Wolf ramp is depicted in Figures 3(f) and 3(g). It can be noticed that reflections produced by Wolf ramp manifest lower amplitudes, even though the impedance change is the same as for other models. Additionally, the reflection produced by a Wolf ramp is extended in time. The character of the reflection changes when the transition zone has larger thickness (compare Figures 3(f) and 3(g)), so in terms of layer of nonuniform thickness across the study area, seismic signature will vary in relation to the thickness interval. Nonetheless, the application of frequency analysis, to which the transition zone is very sensitive, might give additional information that would result in more accurate seismic interpretation.
2.3. Modeling for Transition Zone Recognition
Before we proceeded with real data example, we performed a modeling step. To do so, we used full-waveform modeling and a synthetic geological model of a Wolf ramp. For a model, we used 0-phase Ricker wavelet of specific dominant frequencies. The model is created to simulate the seismic response from a layer of linear velocity change. The model consists of the transition zone of a thickness of 50 meters; the / ratio is the same as in analytical solution (Figure 2) in order to compare the results and verify if the transition zone would be visible on the migrated seismic section. In Figure 4, the resulting traces are presented that were created with different dominant frequencies. The trend is visible—with the increasing frequency, the reflection from a transition zone decreases. The amplitude scale for all panels is consistent; for each frequency, ten zero-offset traces are presented. For frequencies around 30 Hz, it is almost impossible to depict any seismic reflection, even in a complete synthetic scenario, with no noise added. This might suggest that even the transition zone of significant thickness would not be visible on the migrated seismic section.
In comparison, in the numerical solution (Figure 2 vs. Figure 4), it is difficult to indicate the exact amplitude periodicity. Nonetheless, according to the modeling performed by Gomez and Ravazzoli [5], the response from the transition zone might as well reveal only a decaying character of amplitude for higher frequencies. In this respect, the results obtained by full-waveform modeling lay in accordance with the theoretical solution.
2.4. Geological Setting
A potential transition zone was indicated within the Silurian sediments—the Pasłęk Formation that lays within the Polish part of Baltic Basin, Poland (Figure 5). The survey area is situated in the Pomeranian Voivodeship at the border of the Wejherowo and Puck counties. The seismic survey was designed within the concession of Polish Oil and Gas Company, and the area that is covered by the 3D seismic survey is 12.65 km2. The survey area is situated in the western part of the Peri-Baltic Syneclise, within the range of the southeastern slope of the Łeba elevation.
The lithostratigraphic profile of the area is represented by the Precambrian crystalline basement, deposits of the Eocambrian, Cambrian, Ordovician, Silurian, Zechstein, Triassic, Jurassic, Cretaceous, and Cenozoic deposits. The clastic series creates two main complexes: one is of the Caledonian orogeny and encompasses deposits of the Cambrian to Silurian, and the other is associated with the Laramian phase and includes deposits from Permian to Cretaceous. These two complexes are divided by the Variscidian gap that ranges from the Devonian to Carboniferous. The lithostratigraphic profile of the Paleozoic sequence is presented in Figure 6, and a more detailed description of the region can be found in recently published papers [10, 11]. Within this sequence, there were two possible shale gas reservoirs—the Sasino claystone formation of Upper Ordovician and the Pasłęk Formation of Lower Silurian (Aeronian-Telychian).
The Pasłęk Formation consists of dark-grey and grey claystones that are laminated with the greenish and black claystones [12]. The confirmed thickness of the formation might reach up to 60 meters in the Baltic Basin. Most recently, the deposits are described as rhythmic alternations of black, laminated mudstones and greenish, bioturbated mudstones [11]. The formation is considered to be noncalcerous and showing variable to the high degree of bioturbation [11]. The top of the formation is gradational, and this lithological characteristic causes the acoustic impedance of the layer to be gradational hence producing very weak seismic reflection. The specific geological position altogether with the sedimentary history probably explains the source of linear changes in velocity. For the methods and the seismic interpretation, the formation can be just treated as a superficial layer that exhibits linear changes of velocity with depth which can be stated after analysis of well logs (Figure 7). The goal is to verify whether the Pasłęk Formation can be treated as a transition zone and if yes to interpret it, applying a frequency analysis within the 3D seismic survey in order to determine its thickness away from the wells.
(a)
(b)
(c)
2.5. Potential Shale Gas Prospects
The Pasłęk Formation can be found in in the marine and terrestrial region of the Peri-Baltic Syneclise. The complete cores of the formation are available in three wells. The formation is very rich in organic debris, and a wide variety of graptolites was found in the deposits. The Lower Silurian complex can be classified as a potential shale gas reservoir. While drilling of well Darżlubie-IG1 (1973), the Lower Silurian deposits, represented by the Pasłęk Formation, the hydrocarbon content in the drilling fluid increased from 2.28% to 9.43% [13]. More recent well data from other wells estimated the reservoir parameters with the total organic content (TOC) of the Llandovery profile to be at the averaged level between 1 and 6% TOC [14]. Within the study area, the volumetric kerogen content (VKER) of the Pasłęk Formation manifests stable level of about 4-5% [10]. The hydrogen index (HI) varies greatly across the East European Craton, and it is very difficult to indicate the exact type of kerogen by analytical methods due to the high maturity of the sediments [14]. Nonetheless, the Pasłęk Formation manifests slightly lower maturity than the neighboring Jantar and Sasino Formations that both manifest higher TOC. The kerogen type, most confidently, can be determined as type II kerogen of algae-marine origin with high generation potential [14]. In terms of burial and thermal history, Lower Silurian strata reach temperatures of 120° and the maturity modelling has shown that an increase of temperature and burial was continuous from Late Silurian to the Middle Carboniferous [15]. The detailed rock physics modelling of elastic properties [16] classified the formation under analysis to shales/shales with organic matter and hydrocarbons.
3. Methods
Reflection coefficients produced by the transition zone exhibit frequency-dependent attenuation (Figure 2). To interpret the interval, we chose methods of spectral decomposition as the primary tool for analyzing the spatial continuity of the formation. Spectral decomposition enables the definition of the amplitude component associated with the given frequency range. For the transition zone, the most significant result should be associated with the lower frequency range (low frequency anomaly). Firstly, the instantaneous frequency was computed, and then spectral decomposition followed.
For spectral decomposition a seismic volume with relative amplitude preservation (RAP volume) should be preferably used. Such processing would not modulate both amplitude relation and frequency relation [17, 18]. More invasive processing sequences may introduce some frequency range [10] that will affect the response for a transition zone. For these reasons for spectral decomposition, we used a seismic volume without amplitude modulations and gain.
For analysis, three algorithms of spectral decomposition are used: fast Fourier transform (FFT), continuous wavelet transform (CWT), and complete ensemble empirical mode decomposition (CEEMD). The results of spectral decomposition are consistent but differ in details. These differences are highlighted in the next section.
FFT algorithm is based on a linear transform and performs Fourier analysis in a given time window. This algorithm performs well in situations, where the approximate time span of the event can be defined beforehand. For this reason, it is crucial to use information from a borehole and verify the potential thickness of the interval in question. For our research, we have information from 3 wells. Additionally, the layer under analysis is not much involved tectonically—we can assume that it is a flat-lying formation that does not manifest any structural deformations.
CWT algorithm works with the library of wavelets that are modified Gaussian wavelets—called Morlet wavelets. The decomposition process is based on the scaled and shifted Morlet wavelets (unlike the FFT that uses sine functions). CWT algorithm can reach a higher temporal resolution and does not require any information on the scale of the geological feature. By applying CWT, the results are reliable for small features (insignificant thicknesses) and more massive objects (significant thicknesses) which are especially helpful in the case of no geological information from a borehole.
CEEMD algorithm is based on Huang decomposition [19] and uses the sifting process: the local minima and maxima of a seismic trace are found and then based on the defined points, and the upper and lower envelope is computed. Next, the mean of the envelope is computed, and this is called the first intrinsic mode functions (IMFs). The first IMF is subtracted from the signal, and the whole process repeats. In order to secure mode mixing phenomena, a defined portion of Gaussian noise is added to the signal before the process starts. The decomposition continuous until the remaining signal is monotonic. The advantage of the decomposition is that it does not need a predefined set of the decomposition kernels, and the IMFs are purely designed to fit the specific seismic trace.
The parameters used for the specific decomposition are presented in Table 1. The parameters were tested, and the presented set is optimal. Tests were run on the seismic traces from the vicinity of wells and adjusted so that to reach in this control, positions required resolution verified by the synthetic seismograms and real traces. For the FFT algorithm, we tested window lengths between 15 and 50 ms. For the CWT algorithm, other types of wavelet were applied: Ricker and Gaussian; for the CEEMD realization number, it was increased up to 50 that substantially elongated computation time. The results were comparable, and we decreased the number of realizations to the point when the results exhibited the same rate of details, and the computation time was acceptable, resulting in 20 being optimal.
4. Results and Discussion
4.1. Results of Spectral Analysis
In Figure 8, the normalized amplitude spectra computed by FFT and CWT decomposition are presented. These spectra were computed for a single trace that lays in a location of well L-1 for the time position corresponding to the middle point in a transition zone.
(a)
(b)
Spectra shown in Figure 8 are similar: they both indicate approximately the same dominant frequency, that is, 20 Hz, and for the FFT method, it is higher, around 30 Hz. The higher amplitudes for lower frequency content agree with the synthetic example (Figure 2). For the CWT results, higher spectral dynamics are visible, and two frequency peaks are visible (20 and 60 Hz). The second peak has twice lower amplitude. Similar characteristics are visible for the FFT method; nevertheless, the rate of changes is not as prominent; the amplitude drop for the second notch is at the level of 20%. However, the periodicity of the spectral peaks is the same (notch width is app. 40 Hz).
Such a behaviour of the frequency spectra is typical for Wolf ramp that was shown for the modelled data [1, 4, 5] In real data example, the shape of the frequency spectra is affected by many factors, e.g., attenuation, interference, and the existence of the multiples; hence, the function presented in Figure 8 is not smooth. Nevertheless, the lower than averaged dominant frequency (dominant frequency of the survey is 35 Hz) might suggest that the interval under analysis can be understood as a transition zone.
Instantaneous frequency (IF) manifests the rate of changes in the phase of the signal, and it is linked to the maximal spectral amplitude [20]. The results of IF are shown in Figure 9(b), in comparison with the migrated seismic section (Figure 9(a)).
(a)
(b)
In Figure 9(a), the top of the Pasłęk Formation is marked by a light blue line. The reflection that corresponds to the top of the formation diminishes from left to right, losing its strength, which complicates the interpretation and excludes the application of the automatic horizon picking. In Figure 9(b), the interval under analysis is similarly indicated, by the frequency drop that corresponds to the transition zone is more continuous. The formation manifests itself by lower values of instantaneous frequency, at the level between 15 and 25 Hz. It also can be seen that the top of the transition zone resembles very low frequencies. The presented low frequency anomaly has a continuous character and can be visible within all span of a 3D seismic survey. The anomaly has symmetric character around the time of 1755 ms, and its position correlates with the interval, where the velocity is approximated by a linear function. The bottom of the anomaly is clearly defined by a significant increase in the frequency value (indicated by black arrows), and this change indicates the top of the Ordovician (O3). Another interesting effect that can be seen is associated with decreasing thickness of the top of the Ordovician—Prabuty formation, and it is marked by a blue arrow in Figure 9(b). The Prabuty formation undergoes tuning effect that can manifest itself by a sudden change in frequency value. Such effect was modelled and explained by Zeng [21], and we present a real data example of how changes in thickness can be interpreted on the migrated seismic profile (see Figure 9(b), blue arrow in the left-hand side).
A similar comparison is presented in Figure 10. After computing decomposition with the use of the FFT method, it is possible to extract dominant frequency—the frequency value for which the amplitude value reaches its maximum (Figure 10(b)).
(a)
(b)
In the case of spectral decomposition performed with the use of the FFT algorithm, results for a transition zone are consistent with the previously presented IF. However, the low frequency anomaly has higher values of frequency. Also, there exists an increase in frequency (app. 45 Hz) around the OrV seismic horizon; this also occurs in IF display. Nevertheless, the top of the Ordovician similarly manifests itself by a rapid frequency increase (50 Hz). The low frequency anomaly that corresponds to the transition zone has different morphology, and the FFT decomposition reveals more details into it. These differences are indicated by black arrows in Figure 10(b). Spatial comparison of the two frequency attributes proves the continuity of the anomaly within a 3D survey.
The CEEMD algorithm, in comparison with the previously presented, requires significantly more computation time [22]. The result of dominant frequency after the CEEMD for the same inline is presented in Figure 11(b).
(a)
(b)
The results of the CEEMD algorithm are consistent with the IF values, and the low frequency anomaly is visible. The anomaly is more continuous than the result of the FFT decomposition but reveals similar features in terms of its spatial thickness. The frequency values, however, are closer to those indicated by instantaneous frequency (compare Figures 9(b) and 11(b)).
The last decomposition algorithm that was incorporated into the analysis was based on the CWT method. Its results are presented in Figure 12.
(a)
(b)
The results of the CWT algorithm, although consistent with the other decomposition results, reveal slightly different behaviour. The frequency values are similar to FFT and IF, but two morphological objects can be visible. In Figure 12(b), there are regions that can mislead the interpretation (marked with black arrows—Figure 12(b), 1st arrow from left indicating low anomaly overlapping with the top of Ordovician and 2nd arrow indicating abrupt change in anomaly geometry) that exceed the range of the Pasłęk Formation and overlap with the top of Ordovician—the Jantar formation. For this reason, CWT results might be somehow misguiding the automatic interpretation of the base of the Pasłęk Formation. Nonetheless, the top of the formation is more coherent and consistent in comparison with other algorithms revealing the simillar morphology for the top of the formation.
With the results of the spectral analysis, it was possible to estimate the thickness of the Pasłęk Formation (Figure 13). It was constructed by autopicking the top and bed of the low frequency anomaly based on decomposition results. The resulting map is taken as an average value from the thickness estimated by FFT, IF, and CEEMD algorithms. The CWT results were excluded from the computation since in many places, the anomaly was smeared out for more than 150 ms, the temporal distance much greater than a documented thickness of the Pasłęk Formation. The map in Figure 13 shows temporal thickness reaching a maximum of 60 ms and diminishing to almost zero (white region). The distribution of the low frequency anomaly exhibits the specific pattern and is not random, e.g., there exist areas of higher and lower thicknesses of the low frequency anomaly. There is a coincidence with the thickness of the Pasłęk Formation and its structural position. The higher value of temporal thickness is localised in lower structural position, and more elevated areas are associated with lower values of thickness. Such a trend links the existence of lower frequency anomaly to the geological setting (i.e., the structural position of the deposits). Hence, the presence of Wolf ramp manifested by the frequency anomaly can be used as a guideline for seismic interpretation.
4.2. Discussion
We associate the low frequency anomaly with the Pasłęk Formation, which exhibits a linear change of velocity with depth. In Figure 14, the sonic and density logs together with lithological models are presented. The well L-1 lays within the 3D data; two others are situated outside the seismic data coverage. The top of the Pasłęk Formation is marked dashed black line which approximates the linear character of velocity changes for the interval (as shown in Figure 7(b), density value is constant). It can be observed that the formation corresponds to the increased value of clay mineral (VIL) content (green area in Figure 14). Moreover, the change in clay mineral content also shows a specific relation—it increases with depth through the Pasłęk Formation (Figure 14, red dashed line).
(a)
(b)
(c)
In Figure 15, we present the clay mineral content (fraction) for the three wells. The depth of the sample is indicated by colour. For the clay mineral content, a distinctive pattern can be observed—the deeper part of the interval corresponds to the higher clay mineral content, that is also associated with lower velocities. The velocity decreasing with depth makes the Pasłęk Formation the inverted velocity layer. The relationship between the clay mineral content and velocity resembles the linear behaviour. For this reason, we speculate that the velocity decrease observed in the Pasłęk Formation is associated with the increase in clay mineral content.
(a)
(b)
(c)
5. Conclusions
Given the above analysis, we can classify the Pasłęk Formation as an example of a transition zone. The interval is characterized by a low frequency anomaly that was indicated by instantaneous frequency and by spectral decomposition. For the expected thickness of interest, it was crucial to analyze low frequency content of the seismic data as the existence of transition zone manifests itself in the range of lower frequencies. For transition zone interpretation, we propose an application of various decomposition algorithms, so as to compare the results and incorporate advantages of several decomposition methods.
The decrease in velocity for the Pasłęk Formation is associated with the increased fraction of clay mineral content, and the existence of a transition zone is governed by lithofacies changes within the Pasłęk Formation.
With the application of spectral decomposition, it was possible to map the top and bed of the formation that otherwise were difficult to be interpreted. With these results, the map showing the thickness distribution of the interval was created. The thickness of the formation aligns in a specific pattern and is not random, which gives reason to believe that the distribution of the formation was controlled by a geological factor; here, the factor is most likely linked to the structural position.
With the presented methodology and application of transition zone concept, it was possible to interpret the potential shale play reservoir. This enabled us to indicate the thickness of the potential shale play formation that is crucial for estimating the reservoir volume.
Data Availability
Research data are not shared. The data owner is Polish Oil and Gas Company that generously shared the data with the authors for scientific and didactic purposes.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This research has been supported by AGH University of Science and Technology in Kraków (grant numbers 11.11.140.645 and 16.16.140.315). The article is the result of research conducted in connection with a project: seismic tests and their application in the detection of shale gas zones. Selection of optimal parameters for acquisition and processing in order to reproduce the structure and distribution of petrophysical and geomechanical parameters of prospective rocks was part of the program Blue Gas—Polish Shale Gas; (BG1/GASLUPSEJSM/13). We would like to thank CGG for their provision of seismic interpretation software through the University Software Grant Program. We would like to thank Polish Oil and Gas Company for sharing the data for scientific and didactic purposes and for permission for publication of the results. Authors would like to express appreciation for two anonymous reviewers for valuable comments and insight. We thank our colleague, Gabriel Ząbek, for his help in the preparation of the chosen figures.