GO2OGS 1.0: a versatile workﬂow to integrate complex geological information with fault data into numerical simulation models

. We offer a versatile workﬂow 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 scientiﬁc questions or engineering tasks often involves multidisciplinary approaches. Conversion workﬂows 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


Problem description and motivation
To gain understanding of phenomena in our environment, it is common practice to develop models . Since real-world problems show a high complexity, it is necessary for various disciplines to cooperate in order to achieve a proper description of different aspects of a model (Wycisk et al., 2009;Refsgaard et al., 2012;Sharpe et al., 2002). However, the various scientific disciplines use varying types of models to investigate diverse scientific questions (Laniak et al., 2013). A multitude of techniques, workflows, and tools is applied, which often has the result that data conversion between separate tools or data sources is necessary (Sharpe et al., 2002;Gichamo et al., 2012;Walther et al., 2012a;Wojda and Brouyère, 2013). Understandably, it is desirable to provide conversion workflows that are effortlessly repeatable for an arbitrary number of applications, while the employed conversion algorithms are reproducible and understandable. Such workflows, i.e., to integrate various data sources from different sources, are prominent in many disciplines, and have recently been shown in publications like Wu et al. (2005), Hardebol and Bertotti (2013), Shen et al. (2014) and Ragettli et al. (2015).
In the following, we present a workflow to convert data between two types of models, i.e., from a geological structure model into a numerical simulation model. Specifically,  (Jarvis et al., 2008), Germany overview map from Wikimedia Commons (http://commons.wikimedia.org/). we propose a workflow to convert geological structure models, created by the GOCAD (Geological Object Computer Aided Design 1 ) software in the SGrid format, into simulation models usable by the OpenGeoSys finite element simulator in the VTK (Visualization Toolkit) framework format. Geological structure models are used, for instance, to gain better understanding of the development and formation of geological units. Based on the characteristics of the geological units, it is possible to create a hydrogeological simulation model to study the fluid flow in the subsurface. Thus, a geological structure model can aid in obtaining a valid simulation model, by providing the parametrization in a heterogeneous domain. A valid simulation model can help to understand various phenomena and makes predictions of the future possible. Even if a simulation model includes only a few parameters, it can be a challenging task to set up the model appropriately and solve it analytically. In most cases of simulation models, analytical solutions are not available at all. To solve such problems, numerical methods need to be 1 http://www.pdgm.com/products/gocad/ (last access: 16 December 2014) utilized; in this case, we use OpenGeoSys. Similarly to Park et al. (2014), who show an approach to integrating Petrel data into OpenGeoSys (PET2OGS), we developed the GO2OGS workflow for the transition from GOCAD to an open-source format and finally OpenGeoSys.

Modeling tools and site description
GOCAD is a sophisticated tool to model complex geological structures; see for example Zanchi et al. (2009). It can be used as a pre-processor to display hydrogeological simulation models to acquire and assess information on the heterogeneity of aquifer systems.
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. 1 and 2, "GOCAD Model"). The Thuringian Syncline is located in central Germany and has an extension of approximately 150 km from northwest to southeast and 80 km from southwest to northeast. It is surrounded by the Harz and the Kyffhäuser mountains in the north, the Thuringian Forest in   Table 1. the south and the Thuringian Highlands in the southeast. On the basis of the geological structure model, various scientific aspects should be investigated. We selected two different hydrogeological examples to study the applicability of our workflow. Setup A: Eolian, Fluvial, Lacustrin and Sandflat models. Among other research questions, scientists of the Influins project are interested in the course of flow paths in distinct sedimentary depositions. (Kunkel et al., 2013) used GOCAD to create a number of small-scale geological structure models. In the following, we used one of these setups as an example.
Setup B: 3-D geological overview model of the Thuringian Syncline. In contrast to Setup A, a second setup (B) features a more complex structure, with a heterogeneous domain of layered structures, as well as faults and outcropping layers. A complete description of the 3-D geological model of the Thuringian Syncline, which is the basis for the GOCAD model used in this study, and the geological background, will be publicly available in the report TLUG (2015). In the following, we will restrain relevant information for our workflow to acquire a finite element (FE) mesh and set up a simulation model.
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 (Kolditz et al., 2012a) is a well-established, open-source, platform-independent simulation package for coupled THM/C processes (thermal, hydraulic, mechanical and chemical processes). The numerical simulation model utilizes the finite element method (FEM) and has been applied in various instances in the environmental sciences (Sun et al., 2011;Walther et al., 2012bWalther et al., , 2013Nagel et al., 2013;Maxwell et al., 2013;Kolditz et al., 2015).
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.

Previous work and the literature
As a link between geological and visualization sciences, Par-aViewGeo 2 offers interfaces to many file formats used in the geological community. Amongst others, it is able to read GO-CAD ASCII files. However, we were unable to utilize this for our purposes, as not all features of the GOCAD SGrid format, e.g., the so-called split nodes, important for integrating the faults, are supported by the reading algorithms. We could also not continue the work of the ParaViewGeo GO-CAD SGrid reader, as ParaViewGeo, with the last contribution in January 2012, is based on a deprecated ParaView version (3.12.1), while the current ParaView version is 4.2 (last access: 17 November 2014). Maier et al. (2004) use various scripts to convert GOCAD input data to set up a MODFLOW model. Similarly, Luo et al. (2012) describe a way to set up a FEFLOW simulation model using a collection of external scripts and tools. In the study of Ni and Chen (2014), a converting algorithm from GOCAD SGrid models to FLAC3D simulation models is shown. Zehner et al. (2015) describe and compare three GOCAD internal workflows to generate a mesh suitable for numerical simulations of electromagnetic fields. The workflows use external open-source meshing tools (e.g., GMSH, Geuzaine andRemacle, 2009 3 , or TetGen, Si, 2015 4 ), or an interface to the Ansys 5 modeling environment. Recently, Qu et al. (2015) offered a method to generate volumetric gridded fault zones by adding elements in the respective areas to the existing grids. This method, although shown to be applicable on a multitude of setups, will require an existing base mesh of good quality, requires the use of additional software, and is not open-sourced. Unfortunately, none of these tools satisfies every criterium mentioned before, i.e., publicly available algorithms, well-documented workflows, platform independency, as well as the ability to include faults in an unstructured mesh, which renders them unusable for our purposes.
Numerical simulation studies on the hydrogeology in the Thuringian Syncline have been conducted by Rödiger (2005), who created a hydrogeological simulation model for the eastern part of the Thuringian Syncline. Rödiger shows awareness of the existence of faults, but states that the used numerical modeling software (VISUAL MODFLOW) does not support incorporation of the faults. Furthermore, Zech (2013) investigated transport processes in the subsurface within the Thuringian Syncline. Zech used a workflow provided by Zehner (2011) to obtain two vertical 2-D cross sections from GOCAD for numerical simulations with OpenGeoSys. In this case, the data preparation was solely done within GOCAD. Unfortunately, we were unable to base GO2OGS on this, as the reduced complexity of the 2-D setups did not incorporate the necessary features into a full 3-D domain.
As far as the authors are aware, similar groundwater flow simulation codes, e.g., GMS 6 , or FEFLOW, do not have a streamlined interface that satisfies the aforementioned requirements to make use of the GOCAD SGrid format.

Structure of the paper
The conversion workflow to transform a geological structure mesh (SGrid) into a usable simulation mesh (VTU) is sketched in Fig. 2, which simultaneously shows the structure of the paper. In the following, we use the notation , , and to refer to subfigures of Fig. 2 and steps of the workflow, respectively.
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 7 , and is available to the community via https://github.com/envinf/ GO2OGS. In the paper, we will use the style GO2OGS: algorithm to refer to a specific algorithm.
After 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 .

Background
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 ⊂ R 3 can 6 http://www.aquaveo.com/supported-file-types (last access: 16 December 2014) 7 For support, contact thomas.fischer@ufz.de, or info@opengeosys.org be written as where a ij are coefficients, u is a potential and f describes the sources and sinks of a quantity. Furthermore, conditions on the boundary = ∂ have to be specified: where Eq. (1b) is denoted first-type or Dirichlet boundary conditions and Eq. (1c) describes flow through the surface normal of the boundary, i.e., second-type or Neumann boundary conditions. Unfortunately, even for this simple linear problem, an analytical solution does not exist. To solve the equation, a numerical method can be used computing an approximation to the solution.
In the following, we want to apply the finite element method (FEM, Zienkiewicz et al., 2000), which has already been incorporated into the OpenGeoSys open-source simulation toolbox. In order to apply the FEM, we need a partition of the domain , where the problem Eq. (1) is defined; i.e., we need to generate a mesh where the equation can be parametrized.

Meshing and element quality
We used volume elements for the partition of a polyhedral approximation of the problem domain into simpler parts. This partitioning is referred to as meshing.
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 (see Johannsen et al., 2006;Johannsen, 2006;Graf and Degener, 2011).
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. Knupp (2007) gives an overview of the relevance of quality criteria, among others, for FE methods. The OpenGeoSys DataExplorer , a data exploration framework as part of the OpenGeoSys project, provides several quality measures. The OpenGeoSys DataExplorer can be used, for example, to calculate the edge aspect ratio r ea for an element e (face or cell) For reconstruction, one could employ TetGen (Si, 2015), a mesh generator that uses tetrahedra to construct a 3-D domain. In order to mesh the domain with TetGen, one needs to provide piece-wise linear complexes along all outer and internal boundaries. The acquisition of these complexes is especially challenging in areas close to faults, where split nodes constitute discontinuities. Additionally, TetGen will generate an extraordinary number of elements in converging layers (compare Fig. 6 middle and right) to satisfy element quality criteria. This will lead to too many elements to finish the FE simulations and acquire results within reasonable amounts of time. To encounter these challenges, we follow an alternative way for reconstructing the mesh domain by utilizing a hexahedra resampling method. TS2 For the hexahedra resampling based reconstruction, we use the implementation provided by Algorithm 2.

Reconstruction using hexahedra resampling
In the first step, employing GO2OGS: resample, a structured grid within the bounding box of the original unstructured grid y is generated. Depending on the needs of the simulation, it is possible to adjust the resolution in the x, y, or z dimensions. Properties like the affiliation to a geological unit are transmitted from unstructured cells to the cells of the structured grid. In order to determine the property value of a cell c s in the structured grid we look for the cell c u of the unstructured grid that contains the geometric central point of c s . The cell c s is assigned the same property value as is assigned to c u . Cells within the bounding box, but not located in the unstructured grid, are marked as invalid. They are preserved to maintain the structure of the original grid and to simplify the indexing.
In the second step, cells intersecting faults surfaces are marked using GO2OGS: IntegrateFaultFacesets. As a consequence, a fault which is represented as an areal element in the GOCAD SGrid model becomes a volume element in the structured grid.
In the last step, in GO2OGS: removeInvalidCells.py, cells marked as invalid in the first step are removed. The resulting unstructured grid is stored in the open data format VTU z. We called this grid "VTU + " in Figure 2 to signify the improvement of the element quality due to the recon-struction algorithm. Fig. 8 compares the GOCAD mesh y with the resampled mesh z; notice, how the GOCAD mesh not only shows degenerated elements along the fault, but also the zero-thickness elements at the top of the domain.
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. r ea 1), and most often achieved by increasing the number of elements for the representation of heterogeneous features. Yet, this is contrary to issue (2), which at first depends on the system resources of the computing machine (i.e. whether the simulation can be run at all), but eventually influences the runtime of the FE model. Furthermore, issue (3) refers to the fact that in order to omit loosing information of the topology of the geological units, an appropriate z resolution depending on the information density of the underlying geological model is required. Exemplary for the representation of the geological features with respect to the resolution, Fig. 9a shows the original GO-CAD mesh in comparison with two resolutions after reconstruction ( Fig. 9b and c). Finally, (4) is a necessary consideration that will become important after the FE simulation has finished, i.e. whether the investigated processes have been rendered correctly, or that numerical stability criteria were not violated (e.g. Peclet, Courant numbers).

Application example
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 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 project Influins and future work in the study area.
note the remarks at the end of the manuscript. following where i is the ith straight line segment and |e| is the number of segments of e. As a rule of thumb, elements with r ea 1 are unsuitable for the application of the FEM (Zienkiewicz et al., 2000).

Data structures and visualization
OpenGeoSys can read and write the VTU file format of the well-documented Visualization Toolkit (VTK, Schroeder et al., 2006). The documentation for this open data format, as well as the source code of the VTK project, is publicly available from the project website 8 . We use ParaView (Ayachit, 2015) to visualize and evaluate the particular converting steps. ParaView is a widely used open-source visualization toolbox, also in hydrogeological sciences Helbig et al., 2014;Walther et al., 2014). ParaView uses the VTK framework to apply various filters in a pipeline combining a large variety of visualization algorithms.

Converting GOCAD SGrid data to an open data format
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 .

Read original GOCAD data
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 8 www.vtk.org (last access: 16 December 2014) as well as the definition and implementation of data structures used within the algorithm can be found in GO2OGS: GocadSGridReader. Selected intermediate results of the algorithm are depicted in Fig. 4. Executing steps 1-3 of the algorithm results in a structured grid (see Fig. 4a and b). After integrating the split nodes (step 7), the mesh becomes an unstructured grid (Fig. 4d). How the split nodes are stored in the GOCAD SGrid format is documented in Bourke, P. (2008). For the integration of split nodes into the open data format mesh (step 7 of Algorithm 1), we generate new nodes and update the topology of the associated elements. Figure 5 exemplarily visualizes the integration of a split node, and why the GOCAD created meshes may show such bad element quality: notice the creation of a void within the mesh at node (1,1), or the degeneration of element (1,0), after integration of split nodes (1,1').
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 Kunkel et al. (2013); results not shown here). 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. 4d) led to corrupted mesh elements in the original GOCAD SGrid model, with, e.g., non-planar element faces, or near-zero volume elements. Figure 6 visualizes a discontinuous stratigraphic unit in three levels of detail; at the finest level (Fig. 6, right), one can recognize an example for a corrupted mesh element in the Rotliegendes (lower, brown layer).
To evaluate the overall quality of the mesh of Setup B, we calculated the edge ratio criterion (compare Eq. 2) for every cell read from the original GOCAD SGrid (Fig. 7, gray bars). About 60 % of the mesh elements show an aspect ratio r ea ≤ 10 −2 ; as mentioned in Sect. 2.1, elements with r ea 1 are not suitable for FE simulations. Therefore, we will have to employ further conversion steps on Setup B to correct and provide a mesh that can be used in a FE model.

Re-construction of a complex geological 3-D model
To increase the mesh quality of the original GOCAD SGrid mesh, we will use the method of reconstructing the mesh, converting to (see Fig. 2). The method of reconstruction offers the ability to specify element dimensions prior to meshing and, thus, positively influence mesh quality. The aim of the reconstruction is to achieve the closest possible approximation of the geological structures of the original GO-CAD SGrid model, preferably using a small number of highquality cells.
For reconstruction, one could employ TetGen (Si, 2015), a mesh generator that uses tetrahedra to construct a 3-D domain. In order to mesh the domain with TetGen, one needs to provide piece-wise linear complexes along all outer and in- Figure 6. Mesh elements at non-continuous geological units, vertical cross section A-B (see Fig. 1) through the GOCAD model at different magnification levels, vertical exaggeration 20×; legend given in Fig. 3. ternal boundaries. The acquisition of these complexes is especially challenging in areas close to faults, where split nodes constitute discontinuities. Additionally, TetGen will generate an extraordinary number of elements in converging layers (compare Fig. 6, middle and right) to satisfy element quality criteria. This will lead to too many elements to finish the FE simulations and acquire results within reasonable amounts of time. To encounter these challenges, we follow an alternative way for reconstructing the mesh domain by utilizing a hexahedra resampling method.

Reconstruction using hexahedra resampling
For the hexahedra resampling-based reconstruction, we use the implementation provided by Algorithm 2.
In the first step, employing GO2OGS: resample, a structured grid within the bounding box of the original unstructured grid is generated. Depending on the needs of the simulation, it is possible to adjust the resolution in the x, y, or z dimensions. Properties like the affiliation with a geological unit are transmitted from unstructured cells to the cells of the structured grid. In order to determine the property value of a cell c s in the structured grid, we look for the cell c u of the unstructured grid that contains the geometric central point of c s . The cell c s is assigned the same property value as is assigned to c u . Cells within the bounding box, but not located in the unstructured grid, are marked as invalid.
Geosci. Model Dev., 8, 3681-3694 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 gauge Oldisleben (compare Fig. 1), is part of the Thuringian Syncline. Following Eq. (1), and neglecting storage for a long-term mean, the groundwater flow in porous media can be modeled by the equation with the conditions on the boundary: Here, [κ(x)] (L 2 ) is the tensor of permeability of the porous medium at the point (L) denotes the hydraulic head, and [q] (T −1 ) describes sources or sinks of the fluid. For a solution of Eq. (3a) it is necessary to set the material parameters for each element of the mesh and to prescribe Dirichlet (see Eq. 3b) or Neumann (see Eq. 3c) conditions on the boundaries.
For the model at hand, i.e. setup B, we employed Algorithm 2 on y and were able to produce a reconstructed mesh z considering all issues raised in Sect. 2.3. The generated mesh consists of 45.7 mio. cells with a homogeneous cell size of 250 m × 250 m × 10.56 m (x × y × z). This yields an edge aspect ratio r ea ≈ 0.04, i.e. a much better element quality than the original GOCAD mesh (see green bar in Fig. 7). Figure 9c visualizes the level of detail of the heterogeneous geological information for the chosen resolution.

Delineation of a catchment domain
To define the border of the numerical model, we had to cut the geological model of the Thuringian Syncline z with the catchment boundary of the Unstrut river. Therefore, we established the following workflow: (a) the catchment boundary was derived from the SRTM 90 m digital elevation model (Shuttle Radar Topography Mission, Jarvis et al., 2008) employing ArcGIS watershed delineation algorithms, (b) the resulting shape file was used to flag the properties of the mesh elements outside of the catchment boundary (GO2OGS: and subsequently the respective elements were removed (GO2OGS: removeMeshElements). The result was the simulation model {.

Parametrization of the simulation model
With specified boundaries, we could parametrize Eq. (3), and apply boundary conditions to the domain. For the proof of concept, and out of simplicity, the numerical model was set up and parameterized based on plain, but reasonable assumptions. Doubtlessly, there is potential for improving the model setup and, through calibration techniques the results, to approach more realistic projections. Yet, this is beyond the scope of this work, but we will yield recommendations in an outlook for future modeling endeavours.

Fluid and medium properties
We assumed common constant values for dynamic viscosity µ = 1 × 10 −3 Pa · s, and density ρ = 1 × 10 3 kg · m −3 of water. As our model features a steady state simulation, the only additional parameter group that had to be defined was intrinsic permeability k F, L (with superscript F for faults and L for layers). Within the material groups, the parameters were set to a homogeneous value. We defined two distinct scenarios (S1, S2): S1 is the result { of the presented workflow with a heterogeneous domain, with commonly cited values for the inherent geological units (Seidel, 2003). In scenario S1, we treated the faults to be impermeable (i.e. low permeability κ F κ L ). In contrast to S1, scenario S2 is defined through a homogeneous value for all material groups providing a setup for assessing the importance of heterogeneities and faults within the chosen catchment on regional scale. Table 1 gives an overview for the used values, and scenarios of the parameter permeability.
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, TS3 ), we neglected this material group in the setup. The respective elements were removed following a similar approach as described in Sect. 3.1 (GO2OGS: removeMeshElements). The same is probably true for material groups 10-12 (i.e. Zechstein and Upper Rotliegendes), but removing the associated elements resulted in a non-continuous mesh, i.e. holes in the domain; therefore, we did not remove these elements.

Boundary conditions
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 Please note the remarks at the end of the manuscript. They are preserved to maintain the structure of the original grid and to simplify the indexing.
In the second step, cells intersecting fault surfaces are marked using GO2OGS: IntegrateFaultFacesets. As a consequence, a fault that is represented as an areal element in the GOCAD SGrid model becomes a volume element in the structured grid.
In the last step, in GO2OGS: removeInvalidCells.py, 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 + " in Fig. 2 to signify the improvement of the element quality due to the reconstruction algorithm. Fig. 8 compares the GOCAD mesh with the resampled mesh ; notice how the GOCAD mesh not only shows degenerated elements along the fault, but also the zero-thickness elements at the top of the domain.
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., r ea 1), and is most often achieved by increasing the number of elements for the representation of heterogeneous features. However, this is contrary to issue (2), which at first depends on the system resources of the computing machine (i.e., whether the simulation can be run at all), but eventually influences the runtime of the FE model. Furthermore, issue (3) refers to the fact that in order to to prevent loss of information on the topology of the geological units, an appropriate z resolu- tion depending on the information density of the underlying geological model is required. Exemplary for the representation of the geological features with respect to the resolution, Fig. 9a shows the original GOCAD mesh in comparison with two resolutions after reconstruction (Fig. 9b and c). Finally, (4) is a necessary consideration that will become important after the FE simulation has finished, i.e., whether the investigated processes have been rendered correctly, or that numerical stability criteria were not violated (e.g., Peclet numbers, Courant numbers).

Application example
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. 1), is part of the Thuringian Syncline. Following Eq. (1), and neglecting storage for a long-term mean, the groundwater flow in porous media can be modeled by the equation with the conditions on the boundary: Here, [κ(x)] (L 2 ) is the tensor of permeability of the porous medium at the point is dynamic viscosity, [h] (L) denotes the hydraulic head, and [q] (T −1 ) describes sources or sinks of the fluid. For a solution of Eq. (3a) it is necessary to set the material parameters for each element of the mesh and to prescribe Dirichlet (see Eq. 3b) or Neumann (see Eq. 3c) conditions on the boundaries.
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. 2.3. The generated mesh consists of 45.7 million cells with a homogeneous cell size of 250 m×250 m×10.56 m (x×y×z). This yields an edge aspect ratio r ea ≈ 0.04, i.e., a much better element quality than the original GOCAD mesh (see green bar in Fig. 7). Figure 9c visualizes the level of detail of the heterogeneous geological information for the chosen resolution.

Delineation of a catchment domain
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 m digital elevation model (Shuttle Radar Topography Mission, Jarvis et al., 2008) employing ArcGIS watershed delineation algorithms; (b) the resulting shape file was used to flag the properties of the mesh elements outside of the catchment boundary (GO2OGS: ResetPropertiesInPolygonalRegion), and subsequently the respective elements were removed (GO2OGS: removeMeshElements). The result was the simulation model .

Parametrization of the simulation model
With specified boundaries, we could parametrize Eq.
(3) and apply boundary conditions to the domain. For the proof of concept, and for simplicity, the numerical model was set up and parametrized based on plain but reasonable assumptions. Doubtless, there is potential for improving the model setup and, through calibration techniques, the results, to approach more realistic projections. However, this is beyond the scope of this work, but we will yield recommendations in an outlook for future modeling endeavors.

Fluid and medium properties
We assumed common constant values for dynamic viscosity µ = 1 × 10 −3 Pa · s and density ρ = 1 × 10 3 kg · m −3 of water. As our model features a steady-state simulation, the only additional parameter group that had to be defined was the intrinsic permeability k F, L (with superscripts F for faults and L for layers). Within the material groups, the parameters were set to a homogeneous value. We defined two distinct scenarios (S1, S2): S1 is the result of the presented workflow with a heterogeneous domain, with commonly cited values for the inherent geological units (Seidel, 2003). In scenario S1, we treated the faults as impermeable (i.e., low permeability κ F κ L ). In contrast to S1, scenario S2 is defined through a homogeneous value for all material groups providing a setup for assessing the importance of heterogeneities and faults within the chosen catchment on a regional scale. Table 1 gives an overview of the used values and scenarios of the parameter permeability. material group in the setup. The respective elements were removed following a similar approach as described in Sect. 3.1 (GO2OGS: removeMeshElements). The same is probably true for material groups 10-12 (i.e., Zechstein and Upper Rotliegendes), but removing the associated elements resulted in a non-continuous mesh, i.e., holes in the domain; therefore, we did not remove these elements.

Boundary conditions
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. (3c), i.e., ∂ n h(x) = 0, x ∈ , at the vertical borders and the bottom of the chosen catchment domain. A homogeneous groundwater recharge was assumed on the top of the model with 150 mm · a −1 , which lies within a common bandwidth for the area (Seidel, 2003;Rödiger, 2005). We introduced a top layer (see material group 13 in Table 1) with a high permeability to represent infiltrating capabilities of the soil, which are enforced by weathering processes even on initially low permeable units that strike the ground surface. The top layer comprises all cells that own nodes that are exposed to the upper boundary of the mesh, resembling less than 1 % of the total thickness of the modeling domain.
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. 1 for the positions of river nodes). Figure 10 shows the simulated hydraulic head in the model domain of the Unstrut catchment for scenario S1 (i.e., hetero-  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. 13. geneous setup). The model provides a reasonable distribution of the hydraulic head; larger values can be observed in the higher southwestern and northern mountainous regions, while smaller values are calculated for the lower northeastern area near the catchment outlet. The rivers, as boundary conditions, alternate between losing and gaining sections. A comparison of the observed and simulated depth from the surface to the groundwater table (Fig. 11) affirms that model S1 is generally capable of representing the regional groundwater level distribution in the domain; topographical effects are recognizable through larger distances in elevated areas and smaller distances in lower regions, whereas additionally, groundwater levels reach the ground surface in the vicinity of rivers.

Results
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. 12, which shows the absolute difference of hydraulic heads from S1 and S2 at the top of layer 2. Although, for the many areas of the domain, the differences are mostly small (< 2 m), simulation results vary in the northern part of the catchment, where the rivers as boundary conditions dominate the flow regime, and in the southern part, where faults significantly alter the hydraulic head.
The parametrization of the faults and its influence on the scenarios become even more important when observing flow patterns: Fig. 13 shows a particle tracking within an area where two faults reduce the available flow width (see Fig. 12 for the location of the detailed view). In S1, Fig. 13a, the low-permeable faults act as barriers, limiting the pathways of the particles and directing the flow through a "bottleneck" in lower layers (yellow, green, blue, and purple pathlines). As the Top layer (material group 13) was interpreted as an upper stratum being exposed to weathering processes with increased permeability, the faults do not hinder the flow paths and particles may pass (red pathlines). Contrarily, in the homogeneous case S2, Fig. 13b, flow paths are not limited by faults and show an exchange between the rivers, as the overall flow regime is primarily governed by the rivers' boundary conditions.

Discussion
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. 4.1). For example, groundwater levels in the upper elevations (northern areas) are currently governed by the first-type river boundary conditions; utilizing a third-type boundary condition might reflect a more realistic situation. Also, the varying influent and effluent conditions of the streams that are results of the simulations should be compared with observations, and, e.g., the near-river hydrogeology adapted. Although the simulated particle flow paths most likely do not represent the local flow regime correctly, more importantly, the comparison of the scenarios underlines the relevance for the usage of available distributed geological information, and the influence of a heterogeneous parametrization on groundwater flow simulations in the study area.

Conclusions
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 (Fletcher, 1991), whether this concerns, e.g., data acquisition and model parameterization through data scarcity (Ritzema et al., 2010;Gräbe et al., 2013) or model validity (Nguyen and de Kok, 2007). Also, (ii) modern geophysical exploration methods provide increasingly detailed information and high data accuracy, so that adequate numerical representation is required (Thorleifson et al., 2010;Schmelzbach et al., 2011;Van Dam, 2012;Lucia et al., 2015). Thanks to the increasing availability of supercomputing facilities, e.g., Peta-Flop architectures, solving numerical models with 10 9 degreesof-freedom became possible, offering a new dimension of hyper-resolution subsurface modeling (Wang et al., 2014). With respect to numerical modeling, we consider meshing as a non-trivial, complex task that needs special attention in order to provide a sound basis for FE simulations.
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. 2). Specifically, important achievements of this work are the following.
-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 (McKenna et al., 2003;Pozdniakov et al., 2012), offering easier and straightforward integration of subsurface heterogeneous information (Walther et al., 2012b), or increasing prediction quality (Sutanudjaja et al., 2011;de Hoyos et al., 2012). Respective readers and conversion methods can easily be written following the VTK framework documentation.
-Converted results from could be directly used for numerical modeling, here for simulation of flow in smallscale sedimentary structures (Kunkel et al., 2013).
-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.

Outlook
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 (Matter et al., 2005;Tian et al., 2015) or simulation of groundwater age (Goode, 1996) help to identify the flow paths in the field?
-How important is the representation of rivers, i.e., thirdtype 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 Zech (2013), can a regional-scale, heterogeneous, 3-D subsurface model help to understand thermohaline effects?
In order to improve the prediction capabilities and quality of the regional-scale setup, one may instrumentalize one of the following approaches: -by employing automatic calibration methods to obtain better parameter estimates (hydraulic conductivity, porosity, leaching factor of rivers, recharge, etc.), e.g., by utilizing PEST (Gallagher and Doherty, 2007) as in Li et al. (2009) and Sun et al. (2011); -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 T. Fischer et al.: GO2OGS: a versatile workflow to integrate complex geological information -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.

Code availability
The source code of GO2OGS is available through the openaccess repository 9 via https://github.com/envinf/GO2OGS.