Analytical Modelling of Geological Structures with Digital Terrain Surfaces, Three-Variable Functions, Level Surfaces, Graphics Processing Units (GPUs) and Compute Unified Device Architecture (CUDA)

 

by J. Mateo-Lázaro^1*, J. Castillo-Mateo^2,3, J.Á. Sánchez-Navarro¹, A. García-Gil⁴, R. Martínez-Cebolla⁵

¹ Department of Earth Sciences, University of Zaragoza (Spain)

² Department of Statistical Methods, University of Zaragoza (Spain)

³ Department of Statistical Science, Duke University. USA

⁴ Geological Survey of Spain (IGME-CSIC), Madrid

⁵ Urban and Territorial Planning. Government of Aragon

*  Corresponding author: jesmateo@unizar.es

Abstract

Here we present a tried-and-tested application that allows for the construction of geological structure models based on piecewise three-variable functions and their associated level surfaces. Once the master equations have been adjusted, their intersection with specific surfaces can be obtained. This is therefore an analytical approach. Surfaces that intersect analytical geological structures can be analytical, also defined by equations, or non-analytic, such as those defined by a mesh or data grid, such as the digitally defined terrain surface. All of this can be done very quickly using Graphics Processing Units (GPUs), specifically with the Compute Unified Device Architecture (CUDA). This proposal has applications in multiple fields such as geological research, mining, geotechnics, the environment, and others interested in building three-dimensional geological models.

Keywords

3D-Geological Models, GPUs, CUDA, Analytical Three Variable Equations, SHEE software

Note:

Papers published in this special issue of the European Geologist journal have undergone a thorough peer-review process but have not been copy-edited. Authors bear full responsibility for the linguistic accuracy of their contributions.

This work is licensed under a Creative Commons Attribution 4.0 International License.

---

1. Introduction

In the field of work that concerns us, there are currently applications such as Leapflog, GOCAD or GemPy with multiple examples of applications in the literature (Brooks et al, 2011; Hoffman and Gelman, 2011; Schneider et al., 2011; De la Varga et al., 2019)

The interpretation of folds has been addressed since ancient times (Smoluchowski, 1909). The first computer-based geological maps appeared in the early 1970s (Rhind, 1971). The use of three-dimensional (3D) visualization for the reconstruction of geologic structures has become an area of ​​research and development of new methodologies and the visualization of 3D geologic maps began with (Caumon et al., 2013). They proposed a fast and simple tool for dividing a geologic model into one or more arbitrary planes in real time using the OpenGL stencil buffer.

Hypsometric hill-shading is a common way to represent a DEM in a GIS interface. Each hypsometric range can be easily interpreted as a succession of level surfaces belonging to a three-variable “Φ” function that shows the position of the strata. Three variable nonlinear functions can be introduced with components (coefficients, terms, degrees) that allow them to be adjusted to different characteristics of the structural geometry, such as the position of a point in the lithological column. Mathematical operators such as translations, rotations, and scaling can also be applied. The rotations around the axes are equivalent to traditional parameters of Structural Geology, such as strike, dip, and plunge.

Other geometric functions can also be added, of which there are many examples, such as parabolic cylinders, quadrics in general, Gaussian, trigonometric functions, and many others. Another interesting case is that of implicit functions, which do not depend on themselves and are solved using indirect methods such as Newton-Raphson, which is highly effective, bisection and the false-regular, which are suitable for different cases.

The main innovations of the SHEE application’s geological structures module are the first two, but there are other minor ones: (1) For any point in the domain, no numerical operations (calculations or data storage) are performed unless required by the application. For example, in visualization, the entire volume of the domain is not computed, only the sections that need to be visualized, which is a clear advantage in terms of numerical efficiency. (2) The surfaces to be calculated can be defined analytically, for example, plane equations for sections, but also data grids such as an irregular terrain surface. A new feature for the latter case is that no additional numerical operations are required beyond reading the irregular grid data and applying the equation of the three-variable function. (3) The CUDA application is a relatively recent development. On the one hand, it is a widely used platform, and on the other, it is likely the first application that allows the representation of geological structures, both as a preliminary step before adjustment and in real time, given its excellent computational efficiency.

2. Background

For this work, SHEE software interface was used. This software was developed at the Department of Earth Sciences at the University of Zaragoza (Mateo-Lázaro et al., 2014, 2016a, 2016b). Figure 1 shows a screenshot of the program interface. This application is programmed in Visual Basic and uses OpenGL as its main graphics component, specifically for digital terrain models by generating textures with hypsometric colors.

To visualize geological structures, a C++ library supported by the CUDA platform was compiled. This library recalculates the new values ​​of the OpenGL textures. To solve some types of geometric structures, iterative mathematical methods were used for root resolution.

---

---

3. Application

Hypsometric shading is a common way to represent a DEM in a GIS interface. Each color in the table represents a range of altitudes on the terrain surface. Each hypsometric color can be interpreted as representing a succession of horizontal layers, and instead of viewing the “z” altimetry of the DEM, one can view a three-variable “Φ” function (Equation 1) showing the position of strata (Figure 2).

---

---

If the altimetric variable “z” is multiplied by a scalar “e”, a new function “Φ” is obtained (Equation 2) by symbolizing geological characteristic, resulting in this case in a greater thickness for each layer (Figure 3a).

The equation for a monoclinal series of inclined strata (Figure 3b) would be: Ø(x,y,z)=a•x+b•y+c•z, where x, y and z are the positions of an arbitrary point in space, and a, b and c are the components of a vector perpendicular to a plane. This equation represents an infinite number of inclined and parallel planes (Equation 3).

---

---

3.1. Geometrical Functions

The next step is to define “Φ” as a more complex function. For the case of a parabolic cylinder structure (Figure 4), it is necessary to introduce the equation of the parabola in the XZ plane (Equation 4) with a focal length of “P“ with A=±1/4P².

---

---

Other three-variable functions, as well as combinations of multiple functions, can also be used, to fit and represent a three-dimensional geological structure with sufficient fidelity. Figure 5 shows Sine and Exponential functions (Equation 5, 6).

---

---

3.2. Vergence

In addition to the aforementioned parameters, the vergence angle “v” can be introduced, thus obtaining new expressions corresponding to Equations. 7a, 7b in the case of the XZ plane. Figure 6 shows the structure of Figure 5a by introducing Equations. 7.

---

---

3.3. Rotation

Rotations about the axes are equivalent to the traditional parameters of geological structures. Strike is the angle of rotation about the Z axis, and dip is the rotation about the X or Y axis, i.e., angles rZ, rX, and rY, respectively. Figure 7 shows the structure of Figure 5b by introducing Equations 8.

---

---

3.4. Combining Functions and Non-Parallel Layers

Other structures, such as interference folding, can be obtained by combining functions. In Figures 8a and 8b two sine functions, one along the X axis and other along the Y axis, are combined.

Other example is the layout of nonparallel layers as shown in Figure 8c where the exponential equation is used with the parameter “A” (Equation 6) depending on the location and other new parameters (a, b, and c) of the Equation 9.

---

---

3.5. General Equation

Equation 10 is a general equation, which represents all cases above discussed where “φ” corresponds to a combination of functions such as those seen.

3.6. Implicit Functions

An implicit function (e.g., Eq. 11) is characterized by its self-dependence, which means it cannot be solved using direct methods. In the case of Fig. 6, a displacement proportional to “ɸ” can be introduced along the X axis along the tangent plane relative to the inflection point of the folds. Fig. 9 shows this effect using a system of imbricated thrusts.

---

---

3.7. Division into Subdomains

The study domain can be divided into subdomains bounded by conditionals. There are several types of conditionals that can be implemented. In Figure 10, two subdomains are represented: the lower one presents two formations (brown and blue hues) and the upper one is composed of orange, green, and red layers.

Two ways to operate and interrelate subdomains also arise. For the oldest subdomain, the corresponding three-dimensional function can be selected (Equation 12a) or the latest function can be added to the oldest function (Equation 12b), so that the lower subdomain is affected by two episodes of deformation. 

---

---

4. Conclusions and Outlook

This work demonstrates the suitability of combining techniques to investigate the geometry of subsurface geological structures based on data such as their intersection with the terrain. The combination of geological structure modelling with digital terrain surfaces, analytical three-variable functions, and GPU/CUDA-based computation represents an innovative approach. This integration enables applications in a wide range of fields.

In addition to the geometric approach, the perspective could be expanded with equations of physical processes. An example would be sedimentation basin models, whose parameters represent physical processes such as the spatially distributed sedimentation rate, proportional to the depth of the water body and the distance from the basin edge, or subsidence rates, among many others. The visualization methods presented could also be applied in other fields of Earth Sciences. Heat sources, radioactive sources, etc., can be defined using geometric functions (for example, spheres) whose intersection with the ground surface can be visualized to generate intensity maps on that surface.

Piecewise or domain geometric definition techniques are not exclusive to Geology; they are also used in other disciplines. For example, in civil engineering, the definition of the axis of the layout of a highway is not defined with a single analytical equation but rather consists of a progression of sections or pieces of equations, or alignments, circles and straight lines, between which there is a smooth transition of the curvature radius by means Clothoid Functions.

---

References

1. Brooks, S., Gelman, A., Jones, G., Meng, X.-L., 2011. Handbook of Markov Chain Monte Carlo. CRC Press, New York, http://dx.doi.org/10.1201/b10905.
2. Caumon, G.; Gray, G.; Antoine, C.; Titeux, M.O. 2013. Three-Dimensional Implicit Stratigraphic Model Building From Remote Sensing Data on Tetrahedral Meshes: Theory and Application to a Regional Model of La Popa Basin, NE Mexico. IEEE Transactions on Geoscience and Remote Sensing, 51(3). 1613-1621.
3. De la Varga, M., Schaaf, A., Wellmann, F., 2019. GemPy 1.0: open-source stochastic geological modeling and inversion. Geosci. Model. Dev. 12 (1), 1–32. http://dx. doi.org/10.5194/gmd-12-1-2019.
4. Hoffman, M.D., Gelman, A., 2011. The No-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. http://dx.doi.org/10.48550/ARXIV.1111. 4246. 11 K.
5. Mateo-Lázaro, J., Sánchez-Navarro, J.A., García-Gil, A., Edo-Romero, V. (2014). 3D-geological structures with digital elevation models using GPU programming. Computers & Geosciences, 70, 147-153. DOI: http://10.1016/j.cageo.2014.05.014
6. Mateo-Lázaro, J., Sánchez-Navarro, J.A., García-Gil, A., Edo-Romero, V., Castillo-Mateo, J., 2016a. Modelling and layout of drainage-levee devices in river sections. Eng. Geol. 214, 11–19. https://doi.org/10.1016/j.enggeo.2016.09.011
7. Mateo-Lázaro, J, Sánchez-Navarro, JA, García-Gil, A, Edo-Romero, V., 2016b. Flood Frequency Analysis (FFA) in Spanish catchments. J. Hydrol. 538, 598–608. https://doi.org/10.1016/j.jhydrol.2016.04.058
8. Rhind D.W., 1971. The production of a multi-colour geological map by automated means: Nachr. aus den Karten und Vermessungeswesen, Heft 52, 47-52.
9. Schneider, A., Gerke, H.H., Maurer, T. (2011). 3D initial sediment distribution and quantification of mass balances of an artificially-created hydrological catchment based on DEMs from aerial photographs using GOCAD. Physics and Chemistry of the Earth, Parts A/B/C, Volume 36, Issues 1–4, Pages 87-100.
10. Smoluchowski, M.S., 1909. Uber ein gewisses Stabilitats problem der Elastizitatslehre und dessen Beziehung zur Entstehung von Faltengebirgen. Akad. Wissensch. Krakau, Math. Kl., 3-20.

---

This article has been published in European Geologist journal 60 – 5th IPGC Special Edition 1