Skoči na glavni sadržaj

Izvorni znanstveni članak

https://doi.org/10.17794/rgn.2025.5.13

SEED-GUIDED, FAULT-AWARE HORIZON TRACKING VIA GREEDY PROPAGATION FOR GEOLOGICAL MODELLING

Ana Brcković orcid id orcid.org/0009-0004-1915-3266 ; Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb, Pierottijeva 6, 10000 Zagreb, Croatia *
Marko Cvetković orcid id orcid.org/0000-0002-4555-6083 ; Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb, Pierottijeva 6, 10000 Zagreb, Croatia
Josipa Kapuralić orcid id orcid.org/0000-0002-2613-2859 ; Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb, Pierottijeva 6, 10000 Zagreb, Croatia
Ivan Medved ; Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb, Pierottijeva 6, 10000 Zagreb, Croatia

* Dopisni autor.


Puni tekst: engleski pdf 4.334 Kb

str. 167-178

preuzimanja: 155

citiraj

Preuzmi JATS datoteku


Sažetak

The Gola Gas field in the northwestern Drava Basin lies within the structurally complex petroleum system. Its sandstone reservoirs exhibit spatial heterogeneity due to deltaic fan sedimentation, necessitating robust interpretation tools. We present a fully Python-based seismic horizon extraction and 3D geological modelling workflow tailored to this challenge. The method leverages the greedy propagation method to track horizons with optimal reflector continuity, overcoming limitations of manual interpretation and traditional auto-pickers in faulted or noisy zones. An initial approximation is provided via seed-guided interpolation. Around this surface, a seismic sub-volume is extracted and flattened to align reflections. Within this volume, the greedy algorithm identifies a locally optimal path, maximizing continuity and amplitude strength. Our implementation extends the cost-based framework with slope guidance, local correlation, and robust outlier filtering to ensure geologic realism and repeatability. Surfaces are converted to 3D meshes using interpolation and smoothing, and the results are integrated into an interactive 3D geological model of the Gola field. For the first time in the research area, the model was built entirely in Python, enabling automated horizon extraction, geobody visualization, and seismic interpretation within a flexible, open-source environment.

Ključne riječi

seismic horizon; geological model; greedy propagation; dynamic programming; Gola field

Hrčak ID:

337429

URI

https://hrcak.srce.hr/337429

Datum izdavanja:

21.10.2025.

Podaci na drugim jezicima: hrvatski

Posjeta: 500 *




1. Introduction

The Gola Gas field, located in the northwestern Drava Basin, lies within the "Deep Drava" petroleum system, a structurally complex zone shaped by syn-rift tectonics from the Early to Middle Miocene (Lučić et al., 2001; Pavelić, 2001; Royden, 1988) (see Figure 1). Miocene sediments, unconformably overlying Mesozoic basement, include bioclastic limestones, lithic sandstones, and fine-grained conglomerates deposited in marine and deltaic environments. In the NW part of Gola, sandy conglomerate sequences and bioclastic limestones dominate, with evidence of high-density turbidity currents depositing coarse detritus on delta slopes (Lowe, 1982). Gas reservoirs occur mainly in Lower and Middle Miocene clastic units at depths of 1800–3950 m (Saftić et al., 2003), with lateral and vertical heterogeneity reflecting variable depositional settings—from continental to shallow marine. These units exhibit significant lithological variability and reservoir quality, with bioclastic carbonates forming key gas-bearing intervals. Although some conglomerates and sandstones lack fossils, their sedimentological features suggest subaqueous deposition in a proximal prodelta to shallow marine setting (Tadej, 2011). The importance of the selected research area lies in the fact that the deeper subsurface is known for its gas bearing calcitic reservoirs in Early to Middle Miocene and has therefore been explored (Tadej, 2011). However, it has been established that the shallow, Upper Miocene reservoir sandstones of the Gola field were deposited within the sedimentary channels and lobes of prograding delta fans and are classified as stratigraphic traps (Brcković et al., 2024). The importance of delineating sandstone reservoirs lies in the heterogeneity of their spatial distribution and complex sedimentation in deltaic fan environment. The reservoir properties, such as porosity and permeability, have been determined using the core measurements (Tadej, 2011). However, due to intertwined layering of sandstones and shales, the reservoir connectivity is still to be determined and described. The reservoir sandstones have been established in well data (Brcković et al., 2024; Tadej, 2011). However, limitations of traditional seismic interpretation methods provides the space for the application of machine learning methods assisted interpretation. This is why the aim of the article is to present a reproducible Python-based framework combining local greedy propagation method for seismic interpretation, leading to the geological modelling of the research area. Horizons themselves are not isolated artifacts but are essential inputs to stratigraphic modelling and the basis for understanding reservoir architecture. To illustrate the broader geological utility of the extracted horizons, independently derived geobodies were integrated with stratigraphic surfaces to construct a more comprehensive subsurface model.

image1.png

Figure 1. Geographical position of the main research area (Gola field), broader area with seismic data containing the main research area, wells position (well top data) on a map of central Europe with Panonnian Basin System borders delineated.

Seismic horizon interpretation is a basic step in geological modelling framework, which is essential for characterizing structural and stratigraphic frameworks in the subsurface. Traditionally, horizons are identified manually by interpreters who follow coherent reflectors in seismic sections. While this approach can provide geologically valid results, it is often subjective, labor-intensive, and limited in resolution and consistency, especially in large 3D seismic datasets or in structurally complex areas with faults, noise, or stratigraphic discontinuities (Yan & Wu, 2021). Recent methodological innovations have aimed to automate horizon picking to overcome these limitations. Among them, slope-based methods use local dip estimates to propagate horizon picks (Wu & Fomel, 2018). While robust in relatively continuous settings, they struggle with abrupt terminations of horizons such as those formed by faulting. To address this, Wu & Fomel (2018) proposed multigrid correlations to enhance continuity, and Bugge et al. (2019) introduced non-local trace matching to better capture lateral reflection relationships in discontinuous settings. A major leap in horizon extraction has come with the use of deep learning. For instance, Bi et al. (2021) developed a volume-to-volume convolutional neural network that estimates a relative geologic time (RGT) volume from seismic data. This RGT volume can simultaneously capture horizons and faults by treating them as isosurfaces and discontinuities, respectively. The approach leverages synthetic training data and structural similarity-based loss functions to generalize across various seismic settings (Bi et al., 2021). Others similarly used U-Net architectures for direct horizon extraction from seismic amplitudes, showing that neural networks can handle complex geometries and discontinuities when trained appropriately (Tschannen et al., 2020) . Despite the efficiency of deep learning, it remains computationally expensive, data-hungry, and often difficult to interpret or control. As a complementary alternative, dynamic programming (DP) offers an efficient, deterministic method for extracting optimal seismic paths and surfaces. Yan & Wu (2021) proposed a DP-based refinement method that starts from a coarse initial horizon and searches for the optimal amplitude-consistent surface within a defined window. This approach does not rely on continuous slopes alone, making it particularly adept at following reflectors through structurally disturbed zones such as faults or unconformities. Moreover, it integrates user guidance through sparse control points or seed values, which enhances interpretability and flexibility (Yan & Wu, 2021). In contrast to deep networks, dynamic programming is not data-driven but instead exploits local correlation structure and global optimization to ensure geological realism. It is computationally lightweight and interpretable, making it a valuable tool for semi-automatic workflows or for use alongside machine learning where ground truth is limited. By combining these approaches—traditional slope estimation, non-local matching, and dynamic programming—modern seismic interpretation frameworks can achieve unprecedented precision and reliability in horizon tracking. This enables more robust geobody extraction, reservoir modelling, and subsurface analysis, especially in complex geologies where manual or purely slope-based methods fail to trace seismic horizons adequately.

Recent research in 3D geological modelling reveals a strong trend toward combining open-source tools, machine learning, and structured geological logic. Many studies focus on hierarchical or stochastic modelling strategies, often relying on geological rules or learning from data (De La Varga et al., 2019; Feng et al., 2024; Marquetto et al., 2024). A wide range of studies emphasize hybrid methodologies that combine traditional interpretation with computational models, often using decision trees, support vector machines, or unsupervised clustering to improve reservoir prediction and structural modelling accuracy (Chen et al., 2024; Otmane et al., 2025). While Shahbazkia et al. (2025) further demonstrated the robustness of hybrid approaches like Quick Invariant Signature and Dynamic Time Warping (DTW) for handling noisy seismic data. These papers show a broader shift toward modular, flexible modelling workflows. Our approach in this article is simple, modular, deterministic and fully Python-based. It combines interpretable seed-guided surface extraction with structured hexahedral meshing. Our emphasis was on bridging semi-automated interpretation with geobody extraction and volumetric mesh creation, and doing so entirely with accessible, open-source libraries.

2. Materials and methods

The methodology implemented for seismic horizon extraction implements the greedy propagation method (GP) to refine or generate geologically consistent surfaces that follow seismic reflectors with optimal continuity. This approach was motivated by the limitations of both manual interpretation, which is subjective and time consuming, and traditional automated methods, which often struggle in the presence of structural discontinuities like faults or noisy reflections. Traditional methods include but are not limited to the methods such as Maximum Amplitude Picking, Zero-Crossing or Peak/Through Picking, Tracking Using Instantaneous Attributes and Dip-Steered Auto-picking (Gradient-Based), which have been tested throughout our research as well (see Table 1). The proposed methodology applies a greedy algorithm in line with principles of greedy optimization strategies (Lu et al., 2016; Nie et al., 2023).

For seismic horizon picking in Gola field, several seismic attributes were utilized to enhance reflector continuity and improve the reliability of the greedy approach. The core attribute used during tracking was the seismic amplitude, specifically focusing on peak or trough continuity, depending on reflector polarity. However, amplitude alone can be ambiguous in areas of structural complexity, so this was supplemented with the instantaneous phase and local similarity (correlation) attributes to guide the algorithm through low signal-to-noise regions. Amplitude values were extracted directly from the seismic cube and used to identify local maxima or minima as target reflection events. To improve lateral continuity, local similarity (computed via a windowed normalized cross-correlation) was applied in the inline and crossline directions. This measure acted as a confidence mask during DP pathfinding, preferring paths with consistent waveform shapes. Slope estimates were optionally calculated using gradient filters, allowing the algorithm to enforce realistic dip constraints, especially important near faults. All attribute calculations and tracking were implemented in Python, allowing full reproducibility and seamless integration with the automated interpretation framework. These attributes collectively ensured robust and geologically valid horizon extraction even in structurally deformed or noisy seismic sections. The process begins with an initial approximation of the horizon, typically interpolated from sparse seed points or derived through a separate coarse tracking mechanism. This preliminary surface does not need to be precise; it merely guides the subsequent refinement. Around this surface, a sub-volume of seismic amplitudes is extracted, which is flattened relative to the initial horizon to align local reflections. Within this sub-volume, the dynamic programming algorithm searches for a path that maximizes the global consistency of reflection events, typically focusing on peaks, troughs, or zero crossings, while accounting for local continuity and amplitude strength. The DP algorithm ensures a globally optimal solution by avoiding local minima that might mislead simpler tracking methods. A crucial strength of this approach lies in its ability to follow reflectors through complex structures, including faults and noisy regions, by treating horizon tracking as a global optimization problem. This is supported by the previous work by Yan & Wu (2021), who demonstrated how dynamic programming can overcome the shortcomings of slope-based and correlation-based horizon tracking. Their formulation interprets the horizon as an optimal path through a 2D or 3D cost function, where the cost is inversely related to desirable seismic features, such as high correlation or strong amplitude alignment. Our implementation extends this idea by also incorporating local slope guidance, where available, and by using efficient neighbourhood correlation schemes to provide local predictive depth shifts. The method has been enhanced to preserve determinism, ensuring repeatable results, and to filter out physically implausible outliers in the resulting surface. Additionally, by transforming the extracted horizon into a point cloud and applying robust outlier filtering and smoothing (via interpolation and Gaussian filtering), we generate a final continuous surface that is both geologically plausible and visually consistent. Complementing this, recent advances in deep learning-based approach-es, as presented by Bi et al. (2021), offer a data-driven alternative using Relative Geologic Time (RGT) estimation. These methods predict RGT volumes from seismic data using encoder-decoder neural networks, allowing simultaneous interpretation of horizons and faults. While such methods are powerful, especially when large training datasets are available, our dynamic programming method provides a transparent, interpretable, and user-guided alternative that does not rely on black-box models or extensive training data—making it particularly suitable for exploratory interpretation or projects with limited labelled datasets. Together, these methodologies establish a robust, reproducible, and interactive horizon extraction workflow that balances automation with interpretability and adapts well to varying geologic scenarios and data qualities.

The starting point of this workflow is the function that extracts horizon (see Appendix A.1), which implements a robust seed-based greedy algorithm. It propagates from one or more user-defined seed points, each defined by an inline, crossline, and depth value, and traverses neighboring traces to determine the best-matching horizon location. The match is based on correlation of trace segments within a depth window to ensure that the local waveform shape is preserved. Importantly, the method does not match each trace to its immediate neighbor. Instead, it compares all segments to the reference pattern taken from the seed trace, promoting geological consistency across a larger area. Local slope estimates help steer the search direction to follow natural reflector trends, while a depth penalty discourages large vertical jumps, preventing the algorithm from locking onto strong but unrelated reflectors. These combined strategies ensure that the extracted horizon honors both signal similarity and geological structure. After the horizon is extracted, it is typically noisy or contains small-scale deviations caused by trace-level inconsistencies or interference. The next step uses the filtering function to identify and remove extreme outlier values in the horizon surface. This step uses a z-score analysis to quantify how much each depth value deviates from the statistical trend of the rest of the surface. Traces that deviate beyond a user-defined threshold are masked, removing the influence of extreme values without erasing genuine geological variation. This smoothed horizon can now be safely used for surface reconstruction. To move from a scattered set of horizon points to a continuous surface, the function was defined (see Appendix A.2) that converts the 2D DataArray (with inline, crossline, and depth) into a point cloud. It can apply an optional amplitude threshold filter to reject regions with low signal strength, reducing noise propagation into the surface. The resulting dataset contains structured coordinates of interpreted horizon points, which form the basis for gridding and smoothing. The core gridding and smoothing are performed with smoothing function (see Appendix A.3) that transforms irregular horizon point clouds into continuous surfaces. This function supports flexible interpolation methods and applies rigorous quality control, including filtering by z-score and global depth deviation. Users can opt to fill missing gaps using the nearest-neighbor interpolation and apply a Gaussian smoothing filter to gently remove noise while preserving regional structure. This creates clean, geologically plausible surfaces that are ready for visualization, modelling, or comparison with known horizons.

The seismic interpretation workflow developed here revolves around the automation and enhancement of one of the most fundamental tasks in subsurface characterization, which is horizon mapping. Through numerous modifications (see Table 1), it was determined which was the best approach for horizon extraction in our research area. This workflow is structured to bring repeatability, robustness, and geological plausibility into the interpretation process by combining dynamic programming, local dip guidance, amplitude-based correlation, and surface smoothing techniques. The process ultimately produces high-confidence horizon surfaces that reflect subsurface structure with minimal manual intervention and are suitable for subsequent geological modelling. To create geological model and interactive visualization workflow in Python, a combination of specialized scientific, geospatial, and visualization packages was used (see Appendix B). Data was uploaded and processed by using mostly numpy (Harris et al., 2020), xarray (Hoyer & Hamman, 2017), pandas (The pandas development team, 2024) and scipy (Virtanen et al., 2020). Rasterio was used for reading and interpreting the DMR surface from GeoTIFF format (Gillies et al., 2013). To build a 3D geological model, a structured grid was created by layering known surfaces from top to bottom: starting with the digital elevation surface (DMR), followed by a sequence of interpreted horizon surfaces, and finally a user-defined bottom surface. A regular grid was constructed horizontally (X and Y directions), and for each layer, the corresponding depth (Z value) was calculated using surface data. Between each pair of adjacent surfaces, the model volume was filled to represent realistic geological intervals. This layered grid structure reflects the true geometry of the subsurface and can serve as a foundation for simulations, property modelling, or further geological analysis. For the interactive geological model, pyvista was utilized (Sullivan & Kaszynski, 2019).

Table 1. Results of modification for automatic horizon picking functions, with increasing version number meaning better geological plausibility

image2.png

Together, these functions form a modular and adaptable workflow for seismic horizon interpretation. The workflow emphasizes clarity and determinism, with each step explicitly grounded in seismic amplitude patterns and geological structure. It is equally suited for isolated surface extraction or large-scale horizon mapping. Compared to manual picking, which is time-consuming, user-biased, and often inconsistent if multiple people are working on a bigger area, this approach provides rapid, reproducible, and explainable results in Gola field area and even in the wider are of research (see Figure 1). This methodology has been validated on the wider area of our Gola field, making it a powerful foundation for modern geological modelling.

In order to reliably showcase the spatial distribution of extracted horizon surfaces, a geological model of the research area has been built.

3. Results

3.1. Greedy propagation guided horizon extraction

The proposed methodology of local greedy propagation initially based on a dynamic programming approach, has been proven to reliably and semi-automatically extract seismic horizons throughout the seismic volume. The goal was to extract seismic horizon surfaces in order to build an interactive 3D geological model of the research area using free and open-source Python tools.

The development of the horizon extraction workflow followed a stepwise refinement from simple amplitude-based picking to more geologically realistic and robust methods (see Table 1). The earliest versions (v0–v1) focused on trace-wise strategies such as selection of maximum amplitude or inclusion of local dip information to follow structurally plausible trends. While these methods offered fast results, they were structurally inaccurate and unreliable. Since they were based on the amplitude values and with no restraint, the horizon line was visibly connecting similar amplitude strengths on different locations and depths which lead to unexpected peaks in the line (see Figure 1). Subsequent versions (v2–v3) utilized Dynamic Time Warping (DTW) to improve reflector alignment across traces and better handle lateral reflector variability, while versions v4–v5 incorporated smoothing and attribute-based masking to enhance surface continuity and suppress noise from potentially chaotic or faulted zones. The improvement in horizon tracking seemed more stable but the depth was translated and thus the horizon detection was unreliable (see Figure 1). Subsequent versions (v6–v8) introduced dynamic programming (DP) as a global optimization strategy, in which the horizon was treated as a path across traces that minimizes a cost function combining amplitude and smoothness. DP proved to be a powerful tool for enforcing lateral continuity and global structural coherence. It greatly reduced the erratic behaviors seen in previous versions and the selected horizon is tracked well (see Figure 1). However, a critical limitation persisted: DP has no notion of geological relevance. It finds the best path, based on signal, but not necessarily the desired reflector. Seismic data often contains multiple reflectors of similar amplitude, and not all of them are stratigraphically or reservoir relevant. Without a way to guide DP toward a specific reflection of interest, the method could select unintended features, especially in layered or noisy intervals. This challenge was resolved by shifting to a seed-guided approach (v9–v11), which means that one or more known horizon points can be defined, either from well data or from visual inspection of strong reflectors. The final method (v11) implemented greedy seed-based propagation. From each seed point, the method expands outward by choosing the best-matching depth position in adjacent traces, guided by a reference pattern and informed by local dip. This design not only respects structural trends but ensures that only the desired reflector is followed, solving the ambiguity that DP struggled with. This ensured the possibility of extracting horizons for different well tops from well data, in total six of them (see Figure 2).

The decision to adopt greedy local propagation over global optimization was grounded in both methodological and practical considerations. While DP ensures global continuity, it is computationally expensive, especially in 3D volumes or large 2D sections, and it lacks the flexibility to adapt to faults, reflector terminations, or seed constraints. In contrast, greedy propagation is computationally lightweight and directionally flexible. Just by comparing the time needed for execution of the code for one (same) horizon extraction, it is visible how much easier the GP approach is for the end user (in our case around ten times faster). By expanding from seeds using local rules and updating only what is needed, it achieves high performance and scalability, especially when working in interactively defined zones. It also adapts naturally to discontinuities, when a fault or reflector break is encountered, correlation falls below a set threshold, and the propagation halts. This makes the method fault-aware without requiring explicit fault detection. This is a key strength, as many workflows either ignore faults or require a separate, often complex, fault-picking step. This is especially interesting for our research, as there are only few known and detected faults, while discontinuities in reflections are usually due to the stratigraphical and structural complexity of the area.

image3.png

Figure 2. Results of modification to automatic horizon picking functions with some of the modification (versions) from Table 1 being shown. Higher version numbers show better geological reliability in an automatic horizon extraction (more observable on the magnified panels in green).

image4.png

Figure 3. Automatically picked horizons that represent different well tops spatially propagated over the seismic volume

To ensure the reliability of the extracted horizons, a thorough visual validation process was carried out using an interactive 2D seismic viewer (see Figure 3, Figure 4). This custom-built tool allowed for a detailed inspection of both inline and crossline slices, overlaid with the interpreted horizon surfaces (see Appendix A.4). Each horizon derived from well tops was tracked and systematically verified across the full seismic volume. By interactively scrolling through inlines and crosslines, interpreters were able to confirm that the extracted surfaces consistently followed the expected reflector patterns and aligned with geological features observed in the well log data. This visual quality control step was crucial in identifying and correcting any inconsistencies, and it reinforced the accuracy and geological realism of the interpretation workflow.

Using the same sparse set of seed points (between one and five seeds per desired horizon) derived from the test dataset, the proposed methodology successfully scaled to a broader seismic volume that encompassed the original area (see Figure 1, Figure 4). Despite being calibrated on a smaller subset (Gola field, 7 x 7 x 7 km seismic volume), the greedy seed-based horizon tracking approach maintained geological consistency and structural coherence when applied to the larger dataset (broader seismic covered area containing Gola field, 25 x 25 x 10 km seismic volume) (see Figure 1). This outcome highlights the method's robustness and generalizability. It demonstrates that a carefully selected seed, whether from well data or visual interpretation, can guide horizon extraction effectively across extended spatial extents, even when the seeds have local distribution. Such transferability is particularly valuable for regional interpretation, where limited ground truth exists but broader coverage is desired.

image5.png

Figure 4. Horizon Z picking results – tracking was done based on dynamic programming (DP) with the initial seeds from well data (well top position). The results were visually validated and tracking is within the seismic resolution as well as within the expected manual interpretation of the experts.

image6.png

Figure 5. Horizon Z picked in a wider area that includes our research area, Gola field

3.2. Interactive geological model (Python-based)

One of the key outcomes of this work is the construction of an interactive 3D geological model of the research area (Gola field) (see Figure 5). This model integrates multiple data sources and is built on the extracted seismic horizon surfaces from this study, publicly available Digital Model of Relief (DMR) for topographic context, and geobodies previously extracted through semi-automated analysis (unpublished manuscript). These constructive elements (horizons, geobodies) were processed, interpreted, and visualized using fully open-source Python tools. Unlike the usual modelling that relied on proprietary software or merely visualization of the externally prepared models, this project demonstrates a complete open-source workflow from uninterpreted 3D seismic data to geological insight. The model provides an intuitive, high-resolution view of the subsurface architecture while the Python-based workflow offers established workflow for new data generation and insertion (see Appendix B). This is particularly valuable given that the area remains actively investigated, with limited previous modelling attempts. By offering both transparency and interactivity, the model enhances accessibility for geoscientists and stakeholders, while also giving a foundation for future updates as new data become available.

image7.png

Figure 6. Geological model of the Gola field retrieved with a full python framework, containg DMR (uploaded), horizon surfaces (automatically extracted with dynamic programming) and 14 geobodies presenting potential reservoir sandstones with different petrophysical properties (extracted with unsupervised learning methods (unpublished manuscript).

4. Discussion

The developed methodology is centered around a semi-automated, seed-guided horizon tracking function that applies a greedy expansion logic (Lu et al., 2016). It mimics the cognitive process of human interpreters by matching local waveform patterns and following anticipated reflector dips. Specifically, the function extracts a waveform pattern from a seed location and performs trace-by-trace pattern matching using correlation within a defined depth window. A vertical deviation penalty ensures geological continuity, while optional slope steering biases the vertical search window to favor structurally consistent reflectors. Faults are not explicitly detected as some methodology suggestes (Nayak & Nayak, 2025), yet the method exhibits natural fault-tolerance.

Discontinuities are respected due to the local nature of the correlation window, a strict minimum correlation threshold, and the use of a visited mask to prevent backtracking or looping. This reflects how human interpreters hesitate or revise their interpretations near fault zones.

Originally motivated by the dynamic programming (DP) principle, the implementation evolved into a greedy expansion approach for both memory efficiency and geophysical realism. This concept has also been emphasized in seismic applications of dynamic programming by Yan & Wu (2021) and in general algorithmic optimization strategies like those described by Goldner et al. (2013) and Nie et al. (2023). The greedy approach focuses only on reflectors of interest, with no need to explore the entire volume. This stepwise evolution of the method is also reflected in the functional versions developed during this study. Early versions relied on basic amplitude picking and slope guidance. The final version implemented a greedy seed-guided expansion with slope and pattern constraints that offer the best balance between interpretive control, geological consistency, and computational efficiency. The visual validation step ensured that each surface consistently honored the original seed and followed the targeted log marker or reflector over the full seismic volume. Notably, the method generalized well, even when seeds were placed in a smaller test dataset, they were successfully used to extract geologically viable horizons in a much larger surrounding volume.

Strengths of the workflow lie in its reproducibility, interpretability, and full automation potential. With only a single seed, surfaces can be propagated through large volumes, as can be seen in our results. First, we adapted the methodology for Gola field and the results were geologically plausible, especially given the seismic resolution that is around 30 m. Similar integration of automated horizon tracking into regional stratigraphic frameworks has been explored in recent works as well (Chang et al., 2024; Feng et al., 2024; Marquetto et al., 2024). The method is especially attractive to both academia and industry due to its ability to run interactively in modern Python environments. The integration with geobody extraction (unpublished manuscript), 3D surface reconstruction, and topographic data further allows complete subsurface modelling, all within a free and open framework. However, some challenges could be acquainted. They are primarily related to parameter tuning (e.g. depth window, correlation threshold, dip weighting), which can impact stability in structurally complex or noisy areas, but are easily modified. In regions with multiple competing reflectors or low signal-to-noise ratios, the method may still require expert supervision to adjust parameters or validate seed placement. Ambiguity in such areas is inherent to seismic interpretation and not unique to this method.

Ultimately, an interactive geological model was built, constructed entirely using open-source tools, from surface horizons to previously extracted geobodies. This model represents the first fully open-source 3D geological interpretation of the Gola field area, offering valuable insights for research, education, and ongoing industry exploration.

5. Conclusions

This study presents a fully open-source workflow for semi-automatic seismic interpretation and 3D geological modelling, applied for the first time to the Gola gas field. At its core is a seed-guided horizon extraction method that mimics human interpretation logic by applying local pattern matching, vertical continuity constraints, and optional dip steering. The approach balances geological realism with computational efficiency through a greedy expansion algorithm, as a robust, scalable tool. The final implementation successfully tracks geologically consistent surfaces through structurally complex zones without requiring full-volume analysis.

Our results show that even a single seed can initiate reliable surface propagation across large seismic volumes, enabling reproducible, geologically plausible interpretation. When combined with geobody extraction, topography, and interactive visualization, the method forms a complete modelling workframe. Although some parameter tuning is still necessary in complex geological settings, this limitation is minor thanks to the method’s transparency and adaptability. The resulting interactive 3D model marks the first fully open-source structural representation of the Gola field and serves as a valuable resource for research, education, and exploration.

Funding

This research received no external funding.

Author’s contribution

Ana Brcković (Assistant, MSc): conceptualization, data curation, formal analysis, investigation methodology, software, visualization and writing – original draft. Marko Cvetković (Associate professor, PhD): conceptualization, data curation, methodology, resources, supervision, validation and writing – review & editing. Josipa Kapuralić (Assistant professor, PhD): methodology, resources, validation and writing – review & editing. Ivan Medved (Assistant professor, PhD): methodology, resources, software and writing – review & editing.

All authors have read and agreed to the published version of the manuscript.

References

 

Bi, Z., Wu, X., Geng, Z., & Li, H. 2021Deep Relative Geologic Time: A Deep Learning Method for Simultaneously Interpreting 3-D Seismic Horizons and Faults. Journal of Geophysical Research: Solid Earth. 126(9)https://doi.org/10.1029/2021JB021882

 

Brcković, A., Orešković, J., Cvetković, M., & Marić-Đureković, Ž. 2024Enhancing the Understanding of Subsurface Relations: Machine Learning Approaches for Well Data Analysis in the Drava Basin, Pannonian Super Basin. Applied Sciences (Switzerland), 14(14).https://doi.org/10.3390/app14146039

 

Bugge, A. J., Lie, J. E., Kjelsrud Evensen, A., Faleide, J. I., & Clark, S. 2019Automatic extraction of dislocated horizons from 3D 2 seismic data using non-local trace matching 3 4 Right-running head: Seismic horizon correlation. http://library.seg.org/

 

Chang, P. Y., Puntu, J. M., Lin, D. J., Amania, H. H., Chen, W. S., & Lin, A. T. shun. 2024Application of machine learning and resistivity measurements for 3D apparent geological modeling in the Yilan plain, Taiwan, at the SW Tip of the Okinawa trough. Geoscience Letters. 11(1)https://doi.org/10.1186/s40562-024-00339-5

 

Chen, Z., Wang, R., Xu, B., & Zhu, J. 2024Research on Oil and Gas-Bearing Zone Prediction and Identification Based on the SVD–K-Means Algorithm—A Case Study of the WZ6-1 Oil-Bearing Structure in the Beibu Gulf Basin, South China Sea. Energies. 17(22)https://doi.org/10.3390/en17225771

 

De La Varga, M., Schaaf, A., & Wellmann, F. 2019GemPy 1.0: Open-source stochastic geological modeling and inversion. Geoscientific Model Development. 12(1):1–32. https://doi.org/10.5194/gmd-12-1-2019

 

Feng, Y., Wen, G., Shang, J., Wen, S., & Wu, B. 2024Research on 3D geological modeling based on boosting integration strategy. Ore Geology Reviews. 171:https://doi.org/10.1016/j.oregeorev.2024.106157

 

Goldner, E., Silva, P. M. ´ A., & Gattass, M. 20132D Horizon Tracking Using Dynamic Programming.

 

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C.,… Oliphant, T. E. 2020Array programming with NumPy. In Nature. 5857825:357–362. Nature Research. https://doi.org/10.1038/s41586-020-2649-2

 

Hoyer, S., & Hamman, J. 2017xarray: N-D labeled Arrays and Datasets in Python. Journal of Open Research Software. 5(1):10https://doi.org/10.5334/jors.148

 

Lowe, D. R. 1982Sediment gravity flows: II: Depositional models with special reference to the deposits of high-density turbidity currents. Journal of Sedimentary Petrology. 52:279–297

 

Lu, G., Xiong, Y., Ding, C., & Wang, Y. 2016An optimal schedule for urban road network repair based on the greedy algorithm. PLoS ONE. 11(10)https://doi.org/10.1371/journal.pone.0164780

 

Lučić, D., Saftić, B., Krizmanić, K., Prelogović, E., Britvić, V., Mesić, I., & Tadej, J. 2001The Neogene evolution and hydrocarbon potential of the Pannonian Basin in Croatia. Marine and Petroleum Geology. 18(1):133–147. https://doi.org/10.1016/S0264-8172(00)00038-6

 

Marquetto, L., Jüstel, A., Troian, G. C., Reginato, P. A. R., & Simões, J. C. 2024Developing a 3D hydrostratigraphical model of the emerged part of the Pelotas Basin along the northern coast of Rio Grande do Sul state, Brazil. Environmental Earth Sciences. 83(10)https://doi.org/10.1007/s12665-024-11609-y

 

Nie, J., Graizer, V., & Seber, D. 2023A greedy algorithm for wavelet-based time domain response spectrum matching. Nuclear Engineering and Design. 410:https://doi.org/10.1016/j.nucengdes.2023.112384

 

Nayak, B., & Nayak, P. 2025Optimal fault detection from seismic data using intelligent techniques: A comprehensive review of methods. Journal of Groundwater Science and Engineering. 13(2):193–208. https://doi.org/10.26599/JGSE.2025.9280049

 

Otmane, M., Imtiaz, S., Jaluta, A. M., & Aborig, A. 2025Boosting Reservoir Prediction Accuracy: A Hybrid Methodology Combining Traditional Reservoir Simulation and Modern Machine Learning Approaches. Energies. 18(3)https://doi.org/10.3390/en18030657

 

Pavelić, D. 2001Tectonostratigraphic model for the North Croatian and North Bosnian sector of the Miocene Pannonian Basin System. Basin Research. 13(3):359–376. https://doi.org/10.1046/j.0950-091x.2001.00155.x

 

Royden, L. 1988Late Cenozoic tectonics of the Pannonian Basin system. In L. H. Royden & F. Horvath (Eds.), , editor. The Pannonian Basin: A study in basin evolution. AAPG Memoir 45. AAPG.;

 

Saftić, B., Velić, J., Sztanó, O., Juhász, G., & Ivković, Ž. 2003Tertiary subsurface facies, source rocks and hydrocarbon reservoirs in the SW part of the Pannonian Basin (Northern Croatia and south-western Hungary). Geologia Croatica. 56(1):101–122. https://doi.org/10.4154/232

 

Gillies, S. and others. 2013Rasterio: geospatial raster I/O for Python programmers. https://github.com/rasterio/rasterio

 

Shahbazkia, H. R., Khosravani, H. R., Pulatov, A., Hajimani, E., & Kiazadeh, M. 2025A Comprehensive Comparative Study of Quick Invariant Signature (QIS),Dynamic Time Warping (DTW), and Hybrid QIS + DTW for Time Series Analysis. Mathematics. 13(6)https://doi.org/10.3390/math13060999

 

Sullivan, C., & Kaszynski, A. 2019PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK). Journal of Open Source Software. 4(37):1450https://doi.org/10.21105/joss.01450

 

Tadej, J. 2011Evolution of the Early and Middle Miocene sedimentary environments in the north-western part of the Drava depression based on the well analysis data (Razvoj ranomiocenskih i srednjomiocenskih taložnih okoliša sjeverozapadnog dijela Dravske depresije na temelju podataka iz dubokih bušotina). Faculty of Mining, Geology and Petroleum Engineering.

 

The pandas development team. 2024pandas-dev/pandas: Pandas. Zenodo.;

 

Tschannen, V., Delescluse, M., Ettrich, N., & Keuper, J. 2020Extracting horizon surfaces from 3D seismic data using deep learning. Geophysics. 85(3)https://doi.org/10.1190/geo2019

 

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., … Vázquez-Baeza, Y. 2020SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods. 17(3):261–272. https://doi.org/10.1038/s41592-019-0686-2

 

Wu, X., & Fomel, S. 2018Least-squares horizons with local slopes and multi-grid correlations. http://library.seg.org/

 

Yan, S., & Wu, X. 2021Seismic horizon extraction with dynamic programming. Geophysics. 86(2):51–62. https://doi.org/10.1190/geo2020-0039.1

Acknowledgements

The authors would like to thank INA-Industrija nafte d.d. for providing data. All the analyses, interpretation and modelling have been done using Python.


This display is generated from NISO JATS XML with jats-html.xsl. The XSLT engine is libxslt.