English (pdf)
Article in xml format
Article references
Automatic translation
Send this article by e-mail
Cited by SciELO 
Cited by Google 
Similars in
SciELO 
Similars in Google 
Permalink
Introduction
3D geomodeling is a crucial technique for understanding and defining the geometry of subsurface geology. The geological models represent the subsurface created by integrating geological and structural knowledge from multiple data sources. These models provide a quantitative basis for advanced, exhaustive studies (Brisson et al., 2023). Furthermore, they have played an essential role across multiple areas, including mineral exploration and the oil industry, urban planning, geological risk assessment, and geoscientific research (Cao et al, 2024). The quality and scale of the models are determined by the specific objectives, financial resources, available data, the complexity of the subsurface, and the interpolation algorithm used in the modeling (Ji et al, 2024).
A challenge in the construction and accuracy of geological models is the availability of data (Calcagno et al., 2008). In some cases, data are scarce, particularly, as in developing countries, where investment in subsurface exploration and data acquisition is insufficient. An example of this is Colombia, where is located the case of study. Although some information is available, it is limited. Despite these challenges, coherent geological models must still be developed. Traditionally, geological models are built using commercial software that provides a robust framework, such as Petrel (Gunnarsson, 2011), JewelSuite (Sram et al., 2015), or Leapfrog Geo (Braga et al., 2019). However, they impede research progress due to restricted access. Implementing open-source Python-based libraries helps democratize access to knowledge through their transparency and collaboration (Kavanagh, 2004). In Colombia, the energy transition is driving the need to develop geological models to identify areas suitable for the exploitation of natural resources (Thema and Roa-García, 2023). The use of accessible tools enables more researchers to develop, analyze, and evaluate geological models. This, in turn, fosters the development of new solutions to the challenges of increasing the share of renewable energies in the country's energy matrix.
The commonly used interpolation algorithms in 3D geomodeling are Inverse Distance Weighting (IDW), Linear, Spline, and Kriging (Yamamoto, 1998; Qu et al, 2021) . These algorithms have been implemented in open-source libraries such as GMSH (Liu et al., 2024; You et al, 2025) and GeoMeshPy (Dashti et al, 2023). However, although they have shown their usability in modeling complex geological structures and exploration applications, algorithms often do not accurately represent the true properties of the subsurface (Castro-Franco et al., 2017). Interpolation algorithms generally do not incorporate fundamental geological principles, which limits their ability to generate models that adequately reflect essential geological thinking. This geological sense has been incorporated into the potential-field interpolation method with a numerical approach (Almedallah et al., 2021; Fandel et al, 2021; Wu and Sun, 2021). This method employs field data and geological knowledge, using surface points and orientation data (Calcagno et al, 2008; de la Varga et al, 2019). The potential-field interpolation method has been implemented in an open-source Python-based library, allowing full manipulation of any algorithm and ensuring interoperability with different modeling algorithms. Furthermore, it uses a stochastic modeling approach, which allows for analyzing the uncertainty of geological units while considering the variability inherent to geological phenomena. It represents a high unexplored potential for cases where information is scarce and geological expertise is required (Pakyuz-Charrier et al, 2019).
We propose a framework that implements an open-source Python-based library to generate interoperable geologically meaningful models from scarce data. This framework strategically incorporates geophysical data from interpreted sources, such as seismic lines and borehole logs, to generate the inputs required to develop geological models. The framework was validated through three experiments: 1) comparison of the proposed method with state-of-the-art interpolation algorithms for constructing geological surfaces to evaluate the geological sense and accuracy; 2) simulation of scenarios with different scarcity rates to evaluate the robustness and performance of the algorithm; 3) geophysical modeling of gravimetric parameters to evaluate the model interoperability when integrating it with other modeling algorithms.
Proposed framework
The proposed framework is presented below (Figure 1). The potential-field interpolation method used, the 3D geomodeling process, and its reliability assessment through uncertainty analysis are described. Additionally, software validation is presented to demonstrate its accuracy, robustness, and interoperability.
Figure 1 Proposed framework with seismic and borehole input data to generate a 3D geological model, using the potential-field interpolation method. Furthermore, the validation process of the proposed framework is presented (accuracy, robustness, and interoperability).
Potential-field interpolation method
The potential-field interpolation method defines the interfaces I К (К = 1, 2, ...) between formations using geological surface data, guided by the orientation field inferred from the surface orientation data (Lajaunie et al., 1997). Both data types (surface points and orientation measurements) are co-kriged to interpolate a continuous 3D scalar potential field function T( x ) for any point x = (x, y, z) (Aug et al, 2005). The interface location defines the position of reference iso-values, i.e., the set of points x that satisfies T{pc) = I k aligned with the orientation field. The orientation data are represented as normal 3D unit vectors for each structural plane, determining the gradient of the potential field (Calcagno et al., 2008). The division of this space into distinct regions, based on the scalar field values, is referred to as a lithological block. The surface points can be assigned a formation or structural feature like faults (Wu and Sun, 2021 ). This approach allows for the representation of the geometry of the formations within a framework closely aligned with the geological sense.
3D geomodeling
The proposed framework implements GemPy (Figure 2), an open-source Python-based library for stochastic geological modeling and inversion, which uses the potential-field interpolation method (de la Varga et al., 2019).
The inputs for a GemPy model include (a) surface points with Cartesian coordinates positional (x, y, z) and (b) orientation measurements with cartesian coordinates positional (x, y, z), azimuth, dip, and polarity.
Uncertainty analysis
The open-source Python-based library Gempy implements a stochastic modeling approach. The stochastic approximation generates models with a probabilistic distribution of lithologies within each voxel after multiple iterations. This means that instead of assuming a deterministic distribution of lithologies, a stochastic approach incorporates variability and uncertainty (Bonakdari and Zeynoddin, 2022). The Monte Carlo simulation method for uncertainty propagation is employed to assess uncertainty within the geological model. This Bayesian approach propagates the initial uncertainty of the data in 3D geological models. Introduces geological variability using a normal distribution with a mean of 0 and a standard deviation of 30 (for this study), randomly perturbing the Z coordinate of surface points in each initial model, resulting in multiple statistically valid models (Pakyuz-Charrier et al, 2019; Silva-Cárdenas, 2024). To quantify model uncertainty, we compute the Shannon entropy (#) for each voxel in the model:
where Ƥi is the probability of lithological occurrence in each voxel (Brisson et al, 2023). Therefore, high entropy values indicate areas with greater geological uncertainty.
Validation of the framework
The framework validation was carried out through a series of experiments. We compared the proposed framework with different interpolation algorithms in 3D geomodeling, including kriging, linear, spline, IDW, and convergence interpolation, to identify the most consistent method to represent subsurface properties accurately. The framework’s robustness was then assessed by simulating multiple scenarios with varying scarcity rates ranging from 0 % to 90 %, with uncertainty estimated using a Monte Carlo simulation and the Shannon entropy model. Finally, the interoperability of the framework was tested by integrating the model with other modeling algorithms, such as geophysical modeling, using the open-source Python-based library SimPEG (Cockett et al., 2015). This integration demonstrates the ability of the open-source framework to produce results comparable to robust commercial software.
Case study
Study area
The proposed framework is tested for the Llanos basin in Vi chada, Colombia. The Llanos basin is a Cretaceous foreland basin created from bending subsidence caused by the tecto-nic load associated with the Cordillera Oriental (Bayona et al., 2007). Figure 3 shows that the Quaternary deposits cover the entire study area and host a sequence of sediments of the Cenozoic (Guayabo Formation, León Formation, and Carbonera Formation) deposited on a Precambrian crystalline basement (Parguaza Granite) (Bayona et al., 2007, 2008a, 2008b; Moreno-López and Escalona, 2015).
The sedimentary rocks in the study area gradually de-crease in thickness to the east, as supported by geological models based on aerogravimetric and aeromagnetometric data, which show a progressive decrease in the thickness of the sedimentary sequences as the distance from the foothills of the mountain range increases (Graterol, 2009).
Figure 3 Geological map and stratigraphy of the study area and survey dataset. The study area is characterized at the surface by Quaternary and Neogene units and at depth by the Cenozoic Guayabo, León, and Carbonera formations. The black lines and green points illustrate the spatial distribution of the input datasets.
Dataset
The data used in this study consists of eight 2D seismic reflection sections and four boreholes. The seismic data covers an area of approximately 9200 km2. In this case, the seismic lines were in time, so it was necessary to perform a time-to-depth conversion to precondition the data since GemPy works with Cartesian coordinates. The seismic and borehole interpretation horizons defi ne geological surfaces corresponding to the bases of the Quaternary deposits and the Guayabo Formation, León Formation, and Carbonera Formation. The formation tops in the boreholes are indicated as established markers from the interpretation of the petrophysical logs and the descriptions of the drill cuttings.
Figure 4 illustrates input data collection from interpreted seismic sections and borehole data. Cartesian coordinates of the surface points and the orientation measurements were extracted from the mapped horizons. For orientation measurements, the azimuth and dip values of the horizon were taken at some points.
Figure 4 Distribution of the interpreted seismic sections and the borehole (CPE5-SD). Examples of the input data extracted from the horizons and borehole markers are surface points and orientation measurements. The zoom area illustrates how the orientation is measured based on the inclination of the mapped horizon.
Results and Discussion
3D geomodeling
The potential-field interpolation method defines four interfaces and five formations. This stratigraphic sequence is shown in Figure 5, which presents the 3D geological model obtained by implementing the framework. We observe the tectonic uplift of the basement and thinning of sedimentary formations toward the southeast of the model. The 3D geological model is consistent with that described in the literature on the configuration of the subsurface in foreland basins, which present a tectonic uplift of the crystalline basement in response to an orogenic load (Catuneanu, 2004). Furthermore, this enables the dismissal of proposed hypothetical sedimentary models, such as the Vichada impact structure (Hernández et al., 2009, 2011; Torrado-Perez, 2019), which would create a crater basin not observed in the 3D model.
Figure 5 3D geological model of the study area is based on the potential-field interpolation method. The tectonic uplift of the Parguaza Granite and the thinning of sedimentary formations are exhibited in the southeast of the 3D geological model. The visualization has a vertical exaggeration of 60x.
The stratigraphic sequence comprises Cenozoic sedimentary units (Guayabo Formation, León Formation, Carbonera Formation) of the Llanos basin, overlying the Precambrian basement (Parguaza Granite). Furthermore, in studies such as Graterol (2009) , the uplift of the basement is inferred from aerogravimetric and aeromagnetic surveys conducted in the eastern Llanos Basin. The 3D geological model was generated with accuracy and geological meaning, as it aligns with previous studies in the area.
Uncertainty Analysis
An uncertainty analysis of the 3D geological model was performed, and the probability of occurrence of each geological unit was calculated using a Monte Carlo simulation to propagate uncertainty. Figure 6 shows the probability of lithology to block occurrences, indicating that variability is concentrated at the defined interfaces. The little variability in the input data implies that input uncertainty does not significantly affect the model'saccuracy. The Shannon entropy showed a mean entropy of 0.146 and a maximum entropy of 0.769 for the model.
Interpolation algorithms experiment results
The qualitative comparison of 3D surfaces is shown in Figure 7. The results of this experiment indicate that surfaces generated using traditional algorithms (kriging, spline, IDW) and the convergent interpolation employed in the commercial software Petrel exhibited overfitting to the data, resulting in irregularities without geological significance. Linear interpolation was the only traditional algorithm that did not show overfitting. However, it exhibited geometric patterns that were not of geological importance. In contrast, the potential-field interpolation method in our proposed framework accurately represented subsurface properties.
Figure 7 3D geological surfaces generated by different interpolation algorithms. The arrows highlight irregularity in the surfaces, resulting from overfitting and the emergence of geometric patterns without geological sense. The visualization has a vertical exaggeration of 10x.
To complement these observations, a quantitative smoothness analysis based on the RMS of the mean curvature (H) of each surface shows that the potential-field interpolation exhibits the highest curvature values among all methods (Table 1). A higher RMS(|H|) in the potential-field surface is consistent with a more realistic representation of morphological complexity of the geological context, while the lower values in traditional methods reveal a tendency to over-smooth the model and generate unrealistic planar patches in data scarce zones. This method modeled the 3D surfaces consistently with a geological sense without evidencing irregularities due to overfitting the interpolated data or geometric patterns.
Scarce data experiment results
The data reduction was performed by randomly selecting a fixed percentage of points from each layer. The selection followed a uniform probability distribution, ensuring that all points have the same likelihood of being removed. This approach guarantees that all layers retain a consistent proportion of data while keeping the removal process completely random. The Figure 8 presents the 3D visualization of the simulated scenarios, where it is observed that as the data scarcity rate increases, the surfaces lose resolution in their features. However, they still maintain their orientation and geological sense. The distortion observed in areas without data is attributed to the limited information available for a large modeling domain.
Figure 8 3D geological surfaces generated with different scarcity rates (sr). The model with 80 % and 90 % of sr exhibits distortions and low resolution of surface features. The visualization has a vertical exaggeration of 10, and red points represent the data used in each scenario.
A robust performance under scarcity conditions would be manifested in stable entropy, indicating that the framework maintains predictable and efficient behavior despite limitations.
Table 2 shows the entropy analysis as a function of the scarcity rate (sr), presenting the minimum, maximum, and mean entropy values for each simulated scenario. The results show that while the maximum entropy experiences a slight increase with the increase in the scarcity rate, the mean entropy remains notably stable. The observation that the maximum entropy slightly increases could indicate higher variability under data scarcity, but the stability of the mean confirms that these scenarios do not impact the overall performance of the framework. It is important to note that although our framework has strong performance in scenarios with scarce data with relatively homogeneous spatial distribution, areas with low data density tend to become overly smooth. This occurs because the potential-field algorithm generalizes the properties of regions with higher data density, thereby exerting greater control over the model geometry, as documented in previous studies (Aug et al., 2005; de la Varga et al, 2019; Pakyuz-Charrier et al, 2019). Consequently, these scarce data zones with heteroge-nous distribution may lead to geological representations that are unrealistic or biased. Therefore, under condi-tions of scarce data availability, it is essential to maintain the most homogeneous distribution possible.
Table 2 Entropy of the generated 3D geological models. Minimum, maximum, and mean entropy values as a function of the scarcity rate (sr) of each simulated scenario.
| Scarcity rate (%) | Entropy | ||
|---|---|---|---|
| Minimum | Maximum | Mean | |
| 0 | 0 | 0.769 | 0.146 |
| 40 | 0 | 1.095 | 0.167 |
| 80 | 0 | 1.098 | 0.152 |
| 90 | 0 | 1.098 | 0.141 |
The maximum entropy of the models is usually asso ciated with the contact surfaces between the formations, as shown in Figure 9. The metrics show that good performance was obtained even with a scarcity rate of 90 %, considering that the modeling domain is a large area of 9200 km2.
Interoperability experiment results
We evaluated the propagation of gravimetric parameters in geophysical modeling software. As shown in Figure 10, the results of gravimetric modeling are co parable and exhibit similar patterns of residual gravity anomalies.
Gravimetric modeling in both software revealed a negative residual gravity anomaly in the east of the study area, supporting the tectonic uplift of the basement in the basin as established by Graterol (2009) using aero-gravimetry data acquired directly in the study area.
Figure 10 Maps of residual gravity anomalies simulated by integrating the 3D geological model with density properties in A. Oasis Montaj and B. SimPEG software. The results are comparable and exhibit similar patterns of gravitational anomalies.
The discrepancies in the residual gravity anomaly va lues are caused by each model operator using a different approximation. Gravity model response in Oasis Montaj is based on the frequency-domain techniques published by Parker (1973) and Blakely (1995) . On the other hand, SimPEG is based on the numerical solution of gravity equations in the spatial domain, using numerical integration (Cockett et al, 2015; Heagy et al., 2017).
Conclusions
This study demonstrates the effectiveness of the proposed framework for 3D geomodeling using only open-source tools. Robustness enables it to perform well under extreme conditions, including a 90 % datascarci-ty rate. In addition to providing greater precision and geological realism in the representation of subsurface properties through the implemented potential-field in terpolation method, the framework's interoperability with other modeling algorithms enables its applicability to a wide range of analyses and applications. The results obtained with open-source Python-based libraries (GemPy and SimPEG) were comparable to those generated by commonly used commercial software that provides a robust modeling framework (Petrel and Oa sis Montaj), demonstrating that it is possible to achieve results of equivalent quality without the high costs asso ciated with accessing commercial software.
Our open-source framework promotes the democratization of knowledge, facilitating its access and applica tion both nationally and internationally. Our proposed open-source framework could serve as a key driver for exploring potential geological resources, thereby advancing the development of strategic sectors such as renewable energy.
Future work
In various industries, data comparable to that used as inputs in the proposed framework are also acquired, making the methodology potentially replicable across multiple contexts. The literature reports studies of 3D geomodeling with GemPy that provide a basis for mine ral exploration (Güdük et al., 2021; Sehsah et al., 2022; Liang et al., 2023), hydrocarbons (Stamm et al., 2019; Wu and Sun, 2021), and even aquifer characterization and modeling (Thomas et al., 2022; Haehnel et al., 2023). In this direction, future work aims to integrate additional information from shallow boreholes, which are commonly used in mining, to enrich the models and expand the framework's potential applications.
As future work on improving the framework, the flexibility offered by the Python-based implementation of the proposed framework enables the incorporation of gravimetric and magnetometric modeling equations directly into GemPy functions, thereby enabling direct simulation of potential-field responses. This development is feasible because GemPy is an open-source library that allows full access to and modification of its internal functionality. Finally, such integration would enable the framework to operate within a geophysical inversion scheme.
