We offer a versatile
workflow to convert geological models built
with the Paradigm^{™}
GOCAD^{©} (Geological Object
Computer Aided Design) software into the
open-source VTU (Visualization Toolkit unstructured grid) format for
usage in numerical simulation models. Tackling relevant
scientific questions or engineering tasks often involves
multidisciplinary approaches. Conversion workflows are needed as
a way of communication between the diverse tools of the various
disciplines. Our approach offers an open-source, platform-independent, robust, and comprehensible method that is potentially
useful for a multitude of environmental studies. With two
application examples in the Thuringian Syncline, we show how
a heterogeneous geological GOCAD model including multiple layers and
faults can be used for numerical groundwater flow modeling, in our case employing the OpenGeoSys open-source numerical toolbox for groundwater flow simulations. The
presented workflow offers the chance to incorporate increasingly
detailed data, utilizing the growing availability of computational power
to simulate numerical models.

To gain understanding of phenomena in our environment, it is common practice
to develop models

Overview map of the Thuringian Syncline region: border of the GOCAD
geological model (blue line), elevation and rivers (light blue) in the
simulation model of the Unstrut catchment until the Oldisleben gauge (border
as black line), position of cross section A–B (dashed line); background
catchment map from the BFN (2015), digital elevation data from the SRTM (Shuttle
Radar Topography Mission)

In the following, we present a workflow to convert data between two types of
models, i.e., from a

GOCAD is a sophisticated tool to model complex geological structures; see for
example

Workflow of GO2OGS for GOCAD to OGS mesh conversion; Setup A shows
a model of different sedimentary layers; Setup B shows a model of the
Thuringian Syncline (legend given in Fig.

Legend for geological units of Setup B; for abbreviations, see
Table

GOCAD was deployed within the Influins (INtegrierte FLuiddynamik IN
Sedimentbecken) project to create a complex, 3-D, geological structure model of
the Thuringian Syncline including several layers and multiple faults (see
Figs.

We selected OpenGeoSys as a numerical modeling tool, while the presented
workflow to utilize the hydrogeological information from GOCAD (GO2OGS) is
not necessarily limited to the usage of this specific numerical modeling
tool. OpenGeoSys

Geological models can be built in several data formats; in this case a GOCAD stratigraphic grid format (SGrid) was generated for further usage. OpenGeoSys makes use of the VTK framework to read the mesh files in the VTU (Visualization Toolkit unstructured grid) format.

As a link between geological and visualization sciences,
ParaViewGeo

Numerical simulation studies on the hydrogeology in the Thuringian Syncline
have been conducted by

As far as the authors are aware, similar groundwater flow simulation codes,
e.g., GMS

The conversion workflow to transform a geological structure mesh (SGrid) into
a usable simulation mesh (VTU) is sketched in
Fig.

In order to foster and encourage the interested scientific audience, the
equivalent source code for the individual algorithms can be found in an
open-access repository

For support, contact thomas.fischer@ufz.de, or info@opengeosys.org

, and is available to the community viaAfter a brief description of the underlying mathematical equations, we introduce meshing terms and element quality measures. Based on this, we present an algorithm for the conversion of the GOCAD SGrid ➀ to an open data format mesh ➁. By evaluating the element quality of ➁, we propose a reconstruction algorithm to obtain a corrected mesh ➂. Finally, in a third section, we utilize ➃ to delineate a mesh that is used in an application of a regional-scale groundwater flow modeling ➃ that makes use of the heterogeneous information initially obtained by the geological modeling in ➀.

Partial differential equations are used as models for natural processes. For
instance, a linear elliptic partial differential equation of second order
defined on a domain

In the following, we want to apply the finite element method

We used volume elements for the partition of a polyhedral approximation of
the problem domain

The FEM performs arithmetic operations based on the mesh cells. Due to the
finite representation of real numbers as floating point numbers in
a computer, the arithmetic operations are not exact, and the operations'
(intermediate) results are rounded. For example, two nodes can be so close to
each other that during the computation of their distance, cancellation will
occur; i.e., the difference of the operands of the subtraction contains
a smaller number of accurate (significant) digits than the original operands.
Additionally, the element quality depends on the physical process that should
be simulated. For example, in transport process simulations, the Peclet
number has to be acknowledged

In general, mesh quality influences both, the accuracy and efficiency of
a FEM. Therefore, the quality of the cells has to satisfy specific quality
criteria.

OpenGeoSys can read and write the VTU file format of the well-documented
Visualization Toolkit

A GOCAD SGrid data set is comprised of several files. A central file manages all important metadata. This file (extension.sg) may already include the necessary data to describe the geological model and its properties, or may also refer to external data containers. In the following we will describe the methodology to read the GOCAD SGrid data ➀ and write to an open data format ➁.

Algorithm 1 contains the main steps performed to convert
the GOCAD SGrid data to an open data format. The algorithm was implemented
within the OpenGeoSys project; i.e., it uses its data structures to read and
convert the data. The algorithm as well as the definition and implementation
of data structures used within the algorithm can be found in GO2OGS:

Intermediate results of Algorithm 1:

Executing steps 1–3 of the algorithm results in a structured grid (see
Fig.

We employed Algorithm 1 on Setup A, i.e., the small-scale block domains.
After converting the hydrogeological models, the resulting meshes of
➁ could immediately be used as an input for numerical flow
simulations with OpenGeoSys (see

Schematic representation of implementation of split nodes into a
structured mesh

Due to the fact that Setup B shows a more complex structure, the application
of Algorithm 1 did not yield a mesh that is already suitable for FE
simulations. Therefore, in the following, we will only attend to further
processing work on Setup B. Specifically, the domain of Setup B is comprised
of 13, non-continuous stratigraphic units. Furthermore, the model domain
contains 54 vertical faults. Converging stratigraphic units and faults that
are realized through the integration of split nodes (compare
Fig.

To evaluate the overall quality of the mesh of Setup B, we calculated the
edge ratio criterion (compare Eq.

To increase the mesh quality of the original GOCAD SGrid mesh, we will use
the method of reconstructing the mesh, converting ➁ to
➂ (see Fig.

For reconstruction, one could employ TetGen

Mesh elements at non-continuous geological units, vertical cross
section A–B (see Fig.

Histogram of aspect ratio classes for reading GOCAD SGrid mesh data (gray) and reconstructed mesh data (green); sample elements for selected aspect ratios.

For the hexahedra resampling-based reconstruction, we use the implementation provided by Algorithm 2.

In the first step, employing GO2OGS:

In the second step, cells intersecting fault surfaces are marked using
GO2OGS:

Exemplary comparison of meshes before

In the last step, in GO2OGS:
, cells marked as
invalid in the first step are removed. The resulting unstructured grid is
stored in the open data format VTU ➂. We called this grid
“VTU

The selection of the resolution for the reconstruction influences several
aspects, including (1) the element quality, (2) the number of cells needed
for the domain approximation, (3) the level of detail of the geological
information, and (4) the approximation quality for the investigated processes
of the FEM. Issue (1) is relevant to ensure that the generated elements are
not corrupted (i.e.,

Cross section A–B for comparison of a coarse and a fine
reconstruction of geological units to the original GOCAD model, vertical
exaggeration

In the preceding sections, we introduced a workflow to convert a heterogeneous GOCAD mesh, as a result of geological modeling, into an input file for numerical simulations (here: OpenGeoSys) for the application of the FEM method. In the following, we want to show how we used the converted mesh of Setup B in a catchment-scale modeling, primarily as a proof of concept for the conversion workflow. As a secondary result, we expect a first insight into the significance of flow regimes of the different hydrogeological units and faults. This should provide a basis for further investigations within the framework of the Influins project and future work in the study area.

We want to simulate steady-state, confined groundwater flow in a catchment of
the Thuringian Syncline. The chosen study area, the Unstrut catchment until
the Oldisleben gauge (compare Fig.

For a solution of Eq. (

For the model at hand, i.e., Setup B, we employed Algorithm 2 on ➁
and were able to produce a reconstructed mesh ➂ considering all
issues raised in Sect.

To define the border of the numerical model, we had to cut the geological
model of the Thuringian Syncline ➂ with the catchment boundary of
the Unstrut River. Therefore, we established the following workflow: (a) the
catchment boundary was derived from the SRTM 90

With specified boundaries, we could parametrize
Eq. (

Isotropic permeability values

We assumed common constant values for dynamic viscosity

As the Basement (“G”) layer from the GOCAD model does not significantly
contribute to the flow regime (J. Geletneky, Thüringer Landesanstalt für
Umwelt und Geologie, TLUG, personal communication, 2015), we neglected this material group in the setup. The respective
elements were removed following a similar approach as described in
Sect.

Hydraulic head distribution in the simulation model, vertical
exaggeration

Assuming that surface and subsurface catchments are identical in this case,
and that the Basement (G) layer is practically impermeable, we can set
no-flow boundary conditions after Eq. (

A homogeneous groundwater recharge was assumed on the top of the model with
150 mm

We set Dirichlet-type boundary conditions along the main rivers, which act as
major water sources or sinks in the area. The hydraulic head of the river
nodes was set to the respective elevation of the river nodes (compare also
the light blue dots in Fig.

Comparison of observed and simulated depth from surface to
groundwater level; observation data courtesy of TLUG, based on regionalized
observations of groundwater head measurements:

Absolute differences in depth from surface to groundwater level
between scenarios S1 and S2. Light gray lines show faults in the domain of
S1; dark gray lines show rivers. The dashed rectangle displays the frame for
Fig.

Figure

Detail of flow paths near faults and bottleneck structure; pathlines
colored by elevation; area shown in Fig.

The difference between the two scenarios S1 and S2, i.e., the influence of
a heterogeneous vs. a homogeneous parametrization, can be observed in
Fig.

The parametrization of the faults and its influence on the scenarios become
even more important when observing flow patterns:
Fig.

While the simulation generally yields reasonable results for the hydraulic
head and its distribution within the model domain, and the comparison to
observed values of the depth to the groundwater table shows generally similar
values, the model still offers large potential for improvement from
a hydrogeological point of view (compare also
Sect.

Dealing with complex subsurface structures, e.g., sedimentary basins such as
the Thuringian Syncline, is still a challenge for process modeling, e.g.,
fluid dynamics of subsurface systems. Even though sophisticated professional
tools for structural modeling exist (e.g., GOCAD, PETREL), the translation of
structural to accurate numerical models remains a demanding task due to
different (i) classic, and even (ii) new aspects. (i) The
complex-physics–simple-structure vs. complex-structure–simple-physics
paradigm is still valid

With GO2OGS, we provide a complete workflow to utilize complex
hydrogeological information of a geological structure model for the setup of
numerical subsurface simulation models
(Fig.

It has been shown that output of composite geological modeling tools (here GOCAD) may include elements of bad quality; these meshes cannot directly be used in FEM simulation tools.

Corresponding conversion algorithms have been developed that generate mesh domains of good quality and arbitrary spatial resolution. The conversion tools use GOCAD SGrid input ➀ and write VTK framework files (➁ to ➃).

The output of ➁ to ➃ is provided in an open data
format, making the presented workflow and tools versatile and potentially
useful for a multitude of applications, e.g., aiding in fundamental research

Converted results from ➁ could be directly used for
numerical modeling, here for simulation of flow in small-scale
sedimentary structures

As a proof of concept for ➃, we applied the developed workflow on a regional-scale catchment (Unstrut River till the Oldisleben gauge). Therefore, the VTU output of ➃ was used as an input for the OpenGeoSys numerical modeling toolbox. Employing two different scenarios, we could show that the parametrization of faults in the study area is important for the regional and local flow regime. In particular, hydraulic head and particle flow paths are altered by faults, which may affect, e.g., groundwater age or contaminant spreading. Further investigations are needed to reduce the uncertainty of the parametrization and the model.

With the presentation and implementation of a general workflow for the setup
of high-resolution numerical models based on structural information, this
technical paper constitutes a prerequisite for further investigations. With
the workflow, we are now able to tackle scientific questions related to
a better understanding of the dynamics of the Thuringian Syncline's
subsurface systems with accurate numerical subsurface models. Among others,
relevant scientific question are the following.

What is the influence of tilted fractures or the anisotropy of
fractures on the communication between layers? Can isotope
measurements

How will distributed heterogeneous recharge (e.g., from the SWAT
model) or distributed hydrological models

How important is the representation of rivers, i.e., third-type boundary conditions to represent colmation of river sediments, or alternatively, what is the impact from the usage of the river network model output?

In continuation of

by employing automatic calibration methods to obtain better parameter estimates
(hydraulic conductivity, porosity, leaching factor of rivers, recharge,
etc.), e.g., by utilizing PEST

in combination with the former, a parameter sensitivity study will help to answer the question of what the driving components of the groundwater system are, which in turn may yield information for future measurement campaigns (i.e., to identify regions with too few data);

approach numerical improvements, e.g., to include other mesh elements in the reconstruction, which offers the possibility of reducing the mesh size, i.e., the number of elements, while still incorporating relevant structural small-scale features through local grid refinement.

Besides the improvements of the numerical simulation, we see the possibility
of extending the workflow itself. These improvements may address the
following issues:

include more complex, non-vertical fault structures;

use meshes with different element types, e.g., prismatic or tetrahedral elements;

the latter will enable the option of adaptive meshing, especially in the vicinity of faults; and

introduce meshes with different element dimensions, e.g., using 2-D triangles as faults within 3-D cubes as the matrix structure.

Finally, as the underlying geological model is subject to an active and constant update process, the state-of-the-art information of a more detailed model may easily be incorporated into future numerical setups employing the GO2OGS workflow.

The source code of GO2OGS is available through the open-access
repository

For support, contact thomas.fischer@ufz.de, info@opengeosys.org

viaThomas Fischer and Dmitri Naumov developed the model code and performed the simulations together with Marc Walther. Sabine Sattler developed the geological GOCAD model. Thomas Fischer, Marc Walther, Sabine Sattler, and Olaf Kolditz prepared the manuscript.

This work was funded by the German Ministry of Education and Research (grant no. 03IS2091D). The authors appreciate the sincere cooperation with H. Huckriede, J. Geletneky, I. Zander (Thüringer Landesanstalt für Umwelt und Geologie), C. Kunkel (Friedrich Schiller University Jena), as well as C. Herold (Friedrich Schiller University Jena) and S. Attinger (Helmholtz Centre for Environmental Research, UFZ – Leipzig and Friedrich Schiller University Jena). The authors would like to thank the three anonymous reviewers for their comments that significantly increased the quality of the publication. The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association. Edited by: J. Kala