Florian Wellmann

Guillaume Caumon

3-D Structural Geological Models: Concepts, Methods, and Uncertainties

Keywords: Concepts, methods, and uncertainties 3-D geometric modeling, spatial interpolation, uncertainty quantification, geological knowledge, joint inversion, complexity

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction and preliminary considerations

"Everything should be as simple as possible, but not simpler" Albert Einstein (maybe...)

The geological space in the subsurface contributes in essential parts to the anthroposphere, as it has a direct influence on mankind and is equally affected by our technological society. Geological processes create significant heterogeneities at multiple scales below the surface of the Earth. These heterogeneities are essential sources of wealth for mankind (e.g., natural resource deposits) but they are also a main driver of natural hazards such as earthquakes and volcanic eruptions. On the other hand, we influence and alter the subsurface geological space, for example due to resource and water extraction, or engineering and geotechnical applications. Therefore, the evaluation of the spatial distribution of properties in the subsurface plays an essential role for a wide range of applications, from scientific investigations, to geotechnical aspects [START_REF] Culshaw | From concept towards reality: developing the attributed 3D geological model of the shallow subsurface[END_REF], hydrogeological investigations [START_REF] Anderson | Hydrogeologic facies models to delineate large-scale spatial trends in glacial and glaciofluvial sediments[END_REF], raw materials and hydrocarbons [START_REF] Ringrose | Reservoir Model Design[END_REF], and finally policy and public hazard information.

Since the first geological map by William Smith more than 200 years ago, geoscientists have put significant efforts in understanding and visualizing how rocks are organized below the Earth surface, both for practical reasons and out of scientific curiosity to understand the world below our feet. This knowledge is now commonly encapsulated in structural geological models, often simply referred to as geological models or geomodels in this context, as a representation of geometric elements such as rock unit boundaries, faults, horizons and intrusions on a scale of meters to kilometers. The goal of this paper is to present recent advances in the field of geological modeling, with interpolation concepts and algorithms that allow us to represent what we do know-and approaches to analyse and quantify uncertainties to clarify what we don't.

Geological models are intimately linked to geophysics, as they can be seen as a spatial representation of specific aspects of geological knowledge. In this context, they often provide the basis for a spatial parametrization of geophysical forward models, included in geophysical processing and inversion frameworks. At the same time, geophysical observations and interpretations are often an important input to geological models, and geophysical measurements are commonly used to validate geological models, in combination with petrophysical models.

In our experience, the geometric viewpoint of geological models can provide a link between geological and geophysical observations with additional geological and physical information and knowledge, as many relevant aspects share the concept of discrete interfaces in the subsurface as a unifying element, if only in different forms and appearances. To further explain this idea, we will first describe concepts of interface detection in geophysics and geology and review related work. Subsequently, we will describe the current state of geological modeling methods, specifically with the focus on the potential integration into joint geological-geophysical investigations. Even in the best possible combination of information, uncertainties remain-and a substantial part of recent research addresses the investigation and quantification of uncertainties in geological models. The most relevant aspects are described below, followed by a discussion of the concepts and the major challenges, and an outlook to interesting paths for future work.

On geological realism considered in structural geological models

The essence of many subsurface investigations is to obtain an estimate of the true physical property value or class at any location in a region of the subsurface. As we yet have no possibility to directly measure this value exhaustively in space, we need to combine all available measurements, observations, and knowledge to obtain an informed estimate for it.

Of course, the spatial distribution of properties in the subsurface is not random-it is rather the result of the highly complex and long geological history that finally led to the state of the system that we observe in its present form. Unfortunately, we cannot model the entire evolution of the Earth in such a detail that we would be able to obtain a completely realistic picture of the spatial distribution of all relevant properties. Instead, we aim to formulate models that capture essential aspects of this evolution, which we deem relevant for the purpose and scale of a specific investigation. This specific consideration of geological knowledge can be considered as geologic realism, and an essential question for each subsurface investigation is: which type of geological knowledge and how much geologic realism is required for a specific modeling aim and purpose?

In this manuscript, we focus on geologic realism with respect to dominant structural geometric features. The motivation for the consideration of these features is best described with an image of a geological structure, provided in the following section.

A natural rock structure as a motivating example

We take here an example of a folded carbonate layer in the Pyrenees (Aragon, Spain), presented in Fig. 1a. The (simplified) geological story on the area associates sedimentary and tectonic processes [START_REF] Alonso | Forelimb deformation in some natural examples of fault-propagation folds[END_REF]: during the Upper Cretaceous, sandstone (Marboré Formation) was deposited on top of the Pyrenean basement rocks, and was then overlaid during Paleocene by dolomite and massive limestone beds. Layered carbonate beds were deposited during the lower Eocene, then a sudden rise of the relative sea level led to the deposition of deep-marine turbidites of the Hecho Group. This sea-level rise is explained by a fast subsidence related to the Pyrenean compression which provoked a flexure of the lithosphere on either sides of the axial (collision) zone. The North/South shortening then affected these sediments, which accommodated deformation by folding and faulting.

This history is inferred from observations by applying basic stratigraphic principles established in the seventeenth century [START_REF] Steno | De Solido Intra Solidium Naturaliter Contento Dissertationis Prodromus[END_REF]: more recent strata are deposited roughly horizontally on top of older strata; variants stem from later geological events. Geological concepts have much evolved since then, but these essential principles are still applicable in general.

Clearly, we cannot capture the entire complex history quantitatively, including the formation of the basement, the details of erosion and sediment transport and marine life which lead to the formation of sediments, global and regional effects like sea-level rises and subsidence, and finally the complexities of the Pyrenean orogeny. However, what we choose to consider when we talk about structural geological models are surfaces in 3-D space, which capture the essence of geological events in this complicated and long history and which are often manifested at discrete points in time 1 . These surfaces are, for example, related to changes in the depositional system, leading to different rock types that we observe today, or related to deformation events like faulting or ductile deformation 2If we now consider the realistic case that we can not observe the geological domain in an outcrop in the Pyrenees, but that it may be several meters to kilometers below the surface, then the next question is how this geometric representation can be obtained from only a limited amount of information that we typically have (Fig. 1b): boreholes may provide direct observations, but only at very limited points in space, and geophysical measurements provide additional indirect information. How can we "join the dots", considering all of this partial information, and equally use the aspect of geologic realism that we described above?

The fundamental goal for the theory of geological modeling is to obtain methods that allow a reasonable description of geometric and topological features in full 3-D space (Fig. 1c). As all models, this description is bound to approximate the truth within some tolerance. For example, Figure 1: Example of complex geological structures as seen on outcrop (a): this fault-related fold involves juxtaposition of sharp boundaries between turbidites, carbonates, and sandstone. As typical observations are incomplete and sparse (b), a challenge is to decide on which boundaries to represent, and to determine their geometry consistently with conceptual knowledge and the associated interpretations. Location: Aragues del Puerto, Aragon, Spain.

details on the right of Fig. 1a showing some internal structures, joints and details about the fault core zone, may be consciously ignored and deemed impractical to represent explicitly at the scale of interest. The construction of a geological model is limited by the ability of the selected methods to represent the level of geological realism that is considered on a specific scale. The geological models presented here often act as the framework models for models of finer-scale variability, for example using geostatistical approaches (e.g. [START_REF] Pyrcz | Geostatistical Reservoir Engineering[END_REF][START_REF] Ringrose | Reservoir Model Design[END_REF]. Geological models are therefore often hierarchical models, combining multiple scales. Even though we limit our exposure here to structural models, we will mention this link at appropriate positions in the paper and highlight its relevance in the discussion. Also, as our information and knowledge is limited, the description may be locally inaccurate. This is illustrated by some differences between the reality on Fig. 1a and the model on Fig. 1c: Surfaces are locally smoother in the model, and a small fault visible in the outcrop is missing in the model. As a combined result of limited precision, conceptual and observation limitations, each geometric representation will be uncertain. An increasing amount of recent research is going into the evaluation of these related uncertainties.

A note on model complexity

An important aspect of geologic realism is the question of complexity of the investigated geological setting, and specifically, how much of this complexity should be represented in a subsequent geological model. We refer here with geological complexity purely to the complexity related to the geometric setting, in line with the purpose of this paper, (see also Wellmann et al., 2010a;[START_REF] Jessell | Next Generation Three-Dimensional Geologic Modeling and Inversion[END_REF][START_REF] Pellerin | Elements for measuring the complexity of 3D structural models: Connectivity and geometry[END_REF], for more details on this consideration). Typical levels of complexity that we will consider are sketched in Fig. 2, and concern the number of geological features of interest (e.g., number of geological interfaces or contacts between interfaces), and the possible relationships between these features. The simplest possible structural elements in geological models are continuous sub-horizontal interfaces (Fig. 2a): at each location (x, y) in the model domain, this interface has exactly one corresponding z value as a function of this location, z = f (x, y). We refer to these methods here as map-based 3 . This type of structure can be modeled with a wide range of existing interpolation techniques, including geostatistical approaches. The situation is slightly more complex for multiple conformal layers (Fig. 2a, right), as the relationship between layers has to be considered (e.g. layers should not cross, constraints on thickness variations), but techniques for these elements also have long been available (e.g. [START_REF] Abrahamsen | Bayesian Kriging for Seismic Depth Conversion of a Multi-Layer Reservoir[END_REF]. More details are provided in Sec. 3.2.1. The modeling situation is more complex when faults and fault networks are considered (Fig. 2b). First, faults add additional geometric elements that need to be modeled. Furthermore, geological continuity across faults often has to be taken into account, as faults act as discontinuities in the geometric sense and lead to a significant increase in topological complexity (e.g. [START_REF] Pellerin | Elements for measuring the complexity of 3D structural models: Connectivity and geometry[END_REF]Thiele et al., 2016a). And finally, reverse faults (Fig. 2b, right) also increase the complexity for a geological layer, as they lead to a layer doubling, resulting in the interface potentially existing twice at a single location (x, y) in space, making map-based interpolation methods unsuitable.

A high level of complexity in both geology and topology is also quickly reached when more complex poly-deformed terranes and the interaction between geological sequences is taken into account (Fig. 2c), for example due to the consideration of several cycles of sedimentary sequences and unconformities. Other examples are intrusions and dykes in a region. Furthermore, overturned folds lead to interface doubling, similar to reverse faults. The same is true for dome structures (Fig. 2c,right).

Overall, faults, fault networks and deformation or intrusion settings quickly add a level of complexity that many conventional interpolation algorithms can not consider and thus call for dedicated full 3-D geological modeling techniques. These methods will be further discussed in Sec. 3.1.3 for explicit approaches, and Sec. 3.2.3 for implicit approaches.

This consideration of the required level of geological complexity is very important, as it is strongly linked to the selection of a suitable modeling algorithm (discussed in Sec. [START_REF]-D Modelling With Geoscientific Information Systems: The Problem[END_REF], and the possibility to analyze model uncertainties (Sec. 4). It is also an aspect where most geophysical subsurface investigations diverge from more geologically-oriented approaches. This distinction is explained in more detail in the following section.

Two viewpoints on interface detection

The problem of interface detection has occupied geophysicists and geologists for a long time. Arguably, two different viewpoints on this problem have emerged in this context: an approach to detect interfaces on the basis of geophysical data and measurements, described as the data-driven approach below, and an approach on the basis of geological observations and considerations, referred to as a geomodel-driven approach in the following. This contrast in the view on the same problem potentially derives from the academic background and language, and both approaches take different assumptions, but also share many commonalities. The essential aspects of both viewpoints are presented in Fig. 3 and we will develop our interpretation on these viewpoints in the following and especially illustrate how they interact. For clarity, encircled numbers are added in the following text (e.g. 3 ) to refer to positions in Fig. 3. We aim to provide an overview of different aspects related to geological modeling and the link to geophysics, but also hope to provide a framework which can be helpful to place the great body of existing work, or new approaches in this field, into context, and to detect potentially interesting links to related work.

Both viewpoints initiate from the same object of interest: the physical reality 1 of the subsurface. As we cannot directly measure this reality in a comprehensive way, we combine models of measurements of different types with mathematical and numerical models to finally arrive at an approximation of reality, 2 . Both paths are linked in multiple ways and share several commonalities. They are both based on fundamental principles 3 , for example physical laws, and also on geological principles. In the following description, we aim to specifically examine the role of geological knowledge 4 and highlight the role of joint inversion methods 5 as a way to combine geological and geophysical information in a principled way.

The data-driven approach

We first describe the data-driven approach to obtain an approximation of reality on the basis of geophysical measurements 5 . The commonly expressed aim is then to obtain a subsurface model 2 of interfaces or directly property distributions from this data. A great body of work has been performed in this field, and multiple methods are now available and used in scientific investigations and practical applications with specific considerations for each geophysical measurement type. We will limit the following description to a general overview and provide some references to relevant literature with more details.

All geophysical data-driven approaches logically start with geophysical measurements 6 . These measurement may have been taken for a specific purpose (as often the case in targeted exploration, for example), or as a general regional measurement (as, for example, magnetic and gravity stateor country-wide surveys). The measurements themselves are often derived from basic physical principles, but depending on the origin and type of measurement, geological knowledge 4 may have been considered for the measurements themselves, for example in the design of a survey. A typical example is the decision about a suitable trace for a 2-D seismic or gravitational survey, perpendicular to the strike of the investigated geological structures (e.g. [START_REF] Yilmaz | Seismic data analysis: Processing, inversion, and interpretation of seismic data[END_REF][START_REF] Telford | Applied Geophysics[END_REF]. At the same time, geophysical measurements can be used to derive first insights into the subsurface. An experienced applied geophysicist may be able to directly interpret measured signals, for example raw gravity or magnetic measurements to obtain a first insight into subsurface structures. Geophysical measurements themselves can already contribute to a general geological knowledge of a specific region.

Still, the main contribution of geophysical measurements is usually obtained after dedicated data processing 7 . Examples are the various steps in seismic processing from raw data processing and stacking, over migration, to time-depth conversion to obtain a reflection image (e.g. [START_REF] Yilmaz | Seismic data analysis: Processing, inversion, and interpretation of seismic data[END_REF], or the various corrections applied to potential-field data (e.g. [START_REF] Telford | Applied Geophysics[END_REF][START_REF] Jacoby | Gravity Interpretation[END_REF].

Geological knowledge may also enter the processing workflow at this stage, for example in the form of bounds for velocity models for migration in seismic processing. An example for the use of geological reasoning in the analysis of potential-field data are also the various signal processing methods, which are often based on physical principles (for example considering the upward continuation of wavelengths of gravity anomalies), but may also directly consider geological interfaces as the main object of interest, for example when using wavelet-based methods (e.g. [START_REF] Hornby | Analysis of potential field data in the wavelet domain[END_REF].

In seismic processing and imaging [START_REF] Yilmaz | Seismic data analysis: Processing, inversion, and interpretation of seismic data[END_REF], the large redundancy present in reflection The circled numbers and letters refer to important aspects and links, explained in more detail in the text.

seismic data makes it possible to estimate a velocity model of the subsurface and to produce seismic images with limited geological information. In this case, the common approach is to use smoothness constraints in the regularization method; whereas it may be a suitable approach to consider smoothly varying velocities under the Born approximation, many aspects of geology can not be considered in this way, because major and significant changes between subsurface areas with distinctively different properties appear on many sharp boundaries at multiple scales, as obvious in outcrops and rock-faces in nature (Figure 1a). From processed data, two paths can be taken: the first is a direct interpretation of the processed data 8 , and the second a further combination with additional mathematical methods in an inversion framework 5 . In the first case, the aim is to produce a geophysical image, for example a profile of reflection seismic, or a spatially interpolated map of processed potential-field measurements at the surface, also then often called a geophysical map (e.g. [START_REF] Jacoby | Gravity Interpretation[END_REF][START_REF] Dentith | Geophysics for the mineral exploration geoscientist[END_REF]. In this step, the geophysical measurements are necessarily combined with geological knowledge in order to derive a meaningful interpretation (e.g. [START_REF] Saltus | Unique geologic insights from "non-unique" gravity and magnetic interpretation[END_REF]. In a minimal form, the concept of interfaces is used in the interpretation step, as obvious in the detection of interfaces on reflection seismic images (e.g. [START_REF] Yilmaz | Seismic data analysis: Processing, inversion, and interpretation of seismic data[END_REF], or lineaments on gravity maps (e.g. [START_REF] Dentith | Geophysics for the mineral exploration geoscientist[END_REF]. However, there is also a significant contribution of more expressive forms of geological knowledge in this step. For example, in the case of seismic interpretation, [START_REF] Bond | What do you think this is? "Conceptual uncertainty" in geoscience interpretation[END_REF] and [START_REF] Bond | What makes an expert effective at interpreting seismic images?[END_REF] showed the significant influence of prior geological knowledge about a region on the interpretation of structures from seismic images. Even in the case of relatively well-constrained geological knowledge, uncertainties may exist in the location of the number, the geometry and the connectivity of faults.

In addition to geological concepts, geological models often inform the interpretation step of geophysical data. The most obvious link is in the form of geophysical forward models, which are based on property distributions assigned to regions in a geological model, for example a model of gravity response at the surface from defined mass distributions in the subsurface (e.g. [START_REF] Jacoby | Gravity Interpretation[END_REF].

Finally, interpreted images or geophysical maps are already often an end-point in the data-driven workflow and taken as a suitable approximation of reality 2 for a specific purpose. At the same time, the additional understanding during the intellectual process of geophysical image and map interpretation also feeds back into the general geological knowledge about a specific region.

We now consider the second path from processed geophysical measurements 4 to final models 6 using inversion approaches 5 . We aim here to specifically outline the link of various inversion approaches to geological knowledge 3 , specifically with the link to the geometric mode of interface detection.

Typical inputs to inversion approaches are based on either directly using the processed geophysical measurements 7 , or already results from the imaging and interpretation step 8 . A direct combination exists in the inversion of interfaces from gravity and magnetic data (e.g. [START_REF] Silva | Potential-field inversion: Choosing the appropriate technique to solve a geologic problem[END_REF][START_REF] Jacoby | Gravity Interpretation[END_REF]. The non-uniqueness of the inverse problem to determine density distributions from these measurements has been the content of early research in the field of geophysics (e.g. [START_REF] Parker | Best Bounds on Density and Detph from Gravity Data[END_REF], and subsequent work on suitable constraints for inversion (e.g. [START_REF] Li | 3-D inversion of gravity data[END_REF]. Geological knowledge 4 , potentially also derived directly from geological measurements, has, for example, been considered in the form of dip measurements (e.g. [START_REF] Li | Incorporating geological dip information into geophysical inversions[END_REF][START_REF] Lelièvre | Constraining geophysical inversions with geologic information[END_REF]. Many approaches also exist to invert density or velocity distributions or interface positions from a geological starting model (e.g. [START_REF] Bosch | Lithologic tomography: From plural geophysical data to lithology estimation[END_REF][START_REF] Bosch | Lithologic tomography: an application to geophysical data from the Cadomian belt of northern Brittany, France[END_REF][START_REF] Fullagar | Constrained inversion of geologic surfacespushing the boundaries[END_REF][START_REF] Guillen | Geological modelling from field data and geological knowledge: Part II. Modelling validation using gravity and magnetic data inversion: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] Joly | Gravity inversion, AMS and geochronological investigations of syntectonic granitic plutons in the southern part of the Variscan French Massif Central[END_REF][START_REF] Schmidt | Hybrid modelling of gravity, gravity gradients and magnetic fields[END_REF][START_REF] Gradmann | Integrated geophysical modelling of a lateral transition zone in the lithospheric mantle under Norway and Sweden[END_REF][START_REF] Autin | Deep structure of the Argentine margin inferred from 3D gravity and temperature modelling, Colorado Basin[END_REF][START_REF] Haase | A 3D regional crustal model of the NE Atlantic based on seismic and gravity data[END_REF][START_REF] Zheglova | Multiple level-set joint inversion of traveltime and gravity data with application to ore delineation: A synthetic studyMultiple level-set joint inversion[END_REF]. The direct link to the geological model itself 9 is here very obvious.

Similar approaches exist for the inversion of electric and electromagnetic data, as well as seismic data (e.g. [START_REF] Guiziou | 3-D seismic reflection tomography on top of the GOCAD depth modeler[END_REF][START_REF] Clapp | Incorporating geologic information into reflection tomography[END_REF][START_REF] Hauser | A probabilistic seismic model for the European Arctic[END_REF][START_REF] Hauser | Probabilistic inversion of airborne electromagnetic data under spatial constraints[END_REF] In fact, the relevance of boundaries has also long been realised as a central aspect of geophysical imaging techniques and a lot of work has gone in this direction since the pioneering work of [START_REF] Gjøystdal | Computer Representation of Complex 3-D Geological Structures using a new "Sold Modeling[END_REF]. However, several studies explicitly suggested that geological insight can help seismic interpretations (e.g., [START_REF] Rankey | That's why it's called interpretation: Impact of horizon uncertainty on seismic attribute analysis[END_REF][START_REF] Osypov | Model-uncertainty quantification in seismic tomography: method and applications[END_REF] and that, even when using modern seismic acquisition and processing practices, there is still ambiguity in velocity models that can lead to significant uncertainties in the true position of events in images. Very interesting is the also increasing research on full waveform inversion of seismic data (e.g. [START_REF] Kamei | Full waveform inversion of repeating seismic events to estimate time-lapse velocity changes[END_REF][START_REF] Zhu | A Bayesian approach to estimate uncertainty for full-waveform inversion using a priori information from depth migration[END_REF]. These approaches can certainly be considered as dominantly data-driven, however non-uniquenesses still remain and the inversion results depends on an initial (or baseline) velocity model.

As a summary of the data-driven viewpoint based on geophysical measurements 6 to obtain an approximation of reality 2 , we see that a lot of work is aimed at a quantitative and, as much as possible, data-driven and mathematical analysis of the obtained measurement data. However, there is the fundamental problem that the measured geophysical signals provide only an indirect measure of the subsurface interfaces and properties, and an infinite-dimensional manifold of solutions exists [START_REF] Backus | Numerical Applications of a Formalism for Geophysical Inverse Problems[END_REF][START_REF] Parker | Best Bounds on Density and Detph from Gravity Data[END_REF]. From a mathematical point of view, this ambiguity leads to ill-posed problems, and these have to be addressed with suitable regularization methods (e.g. [START_REF] Tikhonov | On the solution of ill-posed problems and the method of regularization[END_REF][START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF]. In fact, the ill-posed nature of geophysical problems lead to many of the mathematical developments in this field (e.g. [START_REF] Kabanikhin | Definitions and examples of inverse and ill-posed problems[END_REF]. In order to obtain non-unique and meaningful solutions, many decisions have to be taken, beginning from suitable regularization methods, to processing steps and the consideration of additional information in the form of geological models or constraints. There exists, therefore, no completely objective result for the inversion of geophysical data.

Many (if not all) of the steps in the "data-driven" approach therefore contain at least some aspects of geological knowledge. Several authors explicitly state that geological considerations increase the effectiveness of geophysical inversion (e.g. [START_REF] Fullagar | Constrained inversion of geologic surfacespushing the boundaries[END_REF]. The fundamental question is therefore how this geological knowledge is obtained-and how relevant parts can be implemented at suitable steps in this "data-driven" viewpoint. These aspects are best explained in the framework of the "geomodel-driven" viewpoint, described below.

The geomodel-driven approach

A different point of view, anchored on centuries of work in the field of geological investigations, is to initiate the subsurface investigation from the viewpoint of main geological interfaces and the related primary information, as well as the associated geological concepts. This way of thinking is underlying the concept of geological maps and has subsequently been extended to full 3-D geological models. This viewpoint is at the core of understanding the relevant steps of geological model construction and the related uncertainties, and all essential aspects will be described in more detail in subsequent sections of this manuscript.

In this viewpoint, primary data are geological observations and measurements, for example in outcrops or wells 10 . Here, the link between observations and geological knowledge 11 becomes essential and deserves a specific note: geological observations can never be separated from geological knowledge. Even taking a measurement of a surface orientation in an outcrop requires the previous geological knowledge about the interpretation of relevance for a feature, or the detection of a fault plane. This realization is at the core of the hermeneutic aspect of geology (e.g. [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF] and described in more detail below.

The set of observations, combined with the purpose of the model and the geological knowledge, typically first enter into a conceptual model 12 . This conceptual model should consider all geological elements and available information, as well as aspects of geological knowledge with relevance for a specific study. This conceptual model is then the basis for the decision on an appropriate mathematical model for the geometrical interpolation in space 13 . Typical interpolation functions and the associated assumptions and limitations are described in detail in Sec. 3. We note here that this decision is also based on the question of geological realism to be considered and, most importantly, on the level of complexity.

Once the representation and associated interpolation methods are defined in the form of a mathematical model, we can assign geological observations as parameters to this model 14 . This step includes the definition of all relevant geological elements, as well as their interaction, for example the number of layer interfaces, and the structure of a fault network. This decision is clearly influenced by the conceptual model, geological measurements and the amount of available information, but, it is also often limited in practice by the specific mathematical model that is applied. The choice of a specific type of model representation and geological interpolation function determines which data can be used-and which can not.

The most commonly used geological observation type is interface points (e.g., [START_REF] Hardy | Multiquadric Equations of Topography and Other Irregular Surfaces[END_REF][START_REF] Dubrule | Comparing splines and kriging[END_REF][START_REF] Mallet | GOCAD: A Computer Aided Design Program for Geological Applications[END_REF], but the direct integration of orientation measurements is only possible in a few interpolation algorithms (e.g., [START_REF] Lajaunie | Foliation Fields and 3D Cartography in Geology: Principles of a Method Based on Potential Interpolation[END_REF][START_REF] Hillier | Three-Dimensional Modelling of Geological Surfaces Using Generalized Interpolation with Radial Basis Functions[END_REF][START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]. We specifically note that the set of geological input parameters can also contain drawings of interfaces in cross-sections, generated by a geological expert, or more subjective knowledge conveyed through interpretive interface points. Yet again, other interpolation methods can consider thickness information (for example kinematic algorithms, [START_REF] Jessell | Noddy" -An interactive Map creation Package[END_REF][START_REF] Wellmann | Pynoddy 1.0: An experimental platform for automated 3-D kinematic and potential field modelling[END_REF].

Next to direct geological observations and measurements, interpreted geophysical images and sections 8 are commonly used as an input to the parametrization of interpolation functions. Typical examples are reflection seismic data and geophysical maps of potential-field data. In many cases, this information enters directly into the interpolation step (e.g. interface picks in a seismic section), in other cases, it is included through the intermediate step of geological knowledge in a region, which is often informed by the consultation of all available geophysical images and interpretations and then considered in the definition of the conceptual model 12 and the choice of the mathematical model 13 , as well as the parametrization of the interpolation functions 14 .

In addition to the parameters based on geological and geophysical observations and measurements, each interpolation function requires specific parameters. These parameters are sometimes estimated on the basis of observations (e.g. parameters of covariance functions in geostatistical interpolations), but often based on a general expected behavior (e.g. smoothness of the interface). In addition, we often need to consider multiple interpolation functions to represent different geological sequences and structures, and then also have to define the interaction of these functions. More details on typically used mathematical interpolation functions and model parametrization are presented in Sec. 3.

Once interpolation function(s) and parameters are defined, we obtain a geometric representation of the defined geological features in space. The most common example today are full 3-D visualizations, either of surfaces as boundary representations, or volume models on discretized mesh structures (e.g. [START_REF] Caumon | Surfacebased 3D modeling of geological structures[END_REF], and these 3-D representations are often referred to as structural geological models or simply geological models4 . We note, however, that 2-D sections are also a common representation, also in the form of a geological map, as the intersection between a 3-D geometric model and a digital elevation model.

The geometric model 15 is then an important aspect of the required approximation of reality 2 . Other common aspects include property fields with additional rock property models, possibly also in combination with further small-scale models of spatial heterogeneity to represent variability on a scale below the considered structures (see Fig. 1).

In any case, given all aspects considered before, it is already obvious that geological models can contain several sources of uncertainty, starting from the initial geological observations, over potentially vague and subjective aspects of geological knowledge, to the choice of the mathematical interpolation approach, its geometric and topological structure, and any additional interpolation parameters. In addition, we also observe that the diverse set of input data and the required choices in the model construction step also lead to different types of uncertainty, for example uncertainties that are related to fixed choices made by the modeling expert in the definition of the conceptual model, different interpretations of features on a geophysical image, or measurement uncertainties in primary data. The treatment of these uncertainties is therefore closely linked to different steps in the geomodel-driven viewpoint and the respective choices in each step. The different sources and types of uncertainty are further discussed in Sec. 2.

It is clear from the description above that the geological model or map is, in essence, only a representation of one or several interpolation function(s) and the defined parameters, even though this aspect is often obscured by the fact that geological modeling methods are often hidden behind "black-box" implementations in commercial software packages. This understanding will, however, become essential for the treatment of diverse aspects of uncertainty, as many methods require the possibility for an automatic update of the geological model. From the previous description, it is evident that this level of automation is generally possible when the choices and settings in the mathematical model 13 and the parametrization 14 are made transparent.

This automation opens the way to an estimation of uncertainties in geological models. This is a key aspect of recent research in this field, and methods to estimate, quantify and visualize the effect of these uncertainties on the subsequent geological models are treated in Sec. 4. An important point to mention here is that this automation allows for a tight integration of the geological modeling step into inversion approaches. This feature is relevant as it allows, for example, the consideration of geological observations and general aspects of geological knowledge that can not directly be integrated into the geological (forward) modeling step. In addition, it provides the possibility for a tight integration with geophysical measurements, through the use of geological models in joint inverse approaches5 . The basis for this integration is typically a geophysical forward simulation on the basis of the spatial distribution of rock properties, obtained from the geological model.

Combining elements

The separation into two distinct approaches is clearly a simplification, but it captures the two main viewpoints on this problem and it is also, in our experience, the cause for a lot of misunderstanding between researchers trained in different fields. Both approaches have strengths and weaknesses. However, we specifically want to highlight here the links and the high connectivity between these different viewpoints, as described above and illustrated in Fig. 3.

Two important links exist between both viewpoints and can address the above concerns: (a) the central aspect of geological knowledge, and (b) integrated inverse approaches. Both aspects are directly related: the underlying geological reality is the combining element between all types of investigations, and our knowledge about it is therefore a central aspect, whereas integrated inverse approaches include the technical possibilities to obtain a model of the reality under the consideration of all available and relevant information.

A central problem with geological knowledge, though, it that it is very difficult to express in a quantitative and objective way. This aspect is fundamentally related to the fact that geological observations themselves can not be taken without a framework of knowledge, as already described above. We can instruct someone how to measure a planar feature with a geological compass, but if we place this person into an outcrop, he or she will be completely lost without the knowledge of how to identify a meaningful planar feature to measure 5 . This aspect is encapsulated in the characterization of geology as a hermeneutic science with a strong interpretative character [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF] and it leads to the fact that many geological observations have a subjective character. This is very different with many geophysical measurements, which can technically be taken (although certainly not interpreted) without any knowledge about the subsurface, purely on the basis of the measurement technique itself. For example, gravity measurements can be performed by experienced technicians. This is certainly a generalization, but overall a significant difference between geological and geophysical measurements and observations.

To re-state this aspect in our context of subsurface investigations: we cannot deny the ubiquity of interpretative and subjective aspects-but we can find ways to address this subjectivity. One approach addresses the elicitation of knowledge from experts as an attempt to reduce subjectivity in primary information (e.g. [START_REF] Curtis | The science of subjectivity[END_REF]. A complementary approach is to find ways to reduce the effect of this uncertainty through the use of all available information-and to achieve this aim, we need to combine both viewpoints in the best possible way.

In the end, we obtain then the approximation of reality 2 , informed by the geometric representation of the geological model and potentially optimized with additional geophysical information. As a result of the early integration of geological concepts and the tight link to geological knowledge, geomodels may include not only features directly visible on the geophysical data at hand, but also features that are expected on the basis of related prior geological information (e.g. subseismic fracture networks). This concept is also commonly known as "Shared Earth Model" in the petroleum industry [START_REF] Gawith | Seismic Validation of Reservoir Simulation Using a Shared Earth Model[END_REF], or "Common Earth Model" in the mining industry (McGaughey, 2006b). Concrete examples of this approach to real-world problems have also been proposed in other fields such as geodynamics (e.g., [START_REF] Ziesch | Geological Structure and Kinematics of Normal Faults in the Otway Basin, Australia, Based on Quantitative Analysis of 3-D Seismic Reflection Data[END_REF], seismology (e.g., [START_REF] Shaw | Unified Structural Representation of the southern California crust and upper mantle[END_REF] and geological engineering (e.g., [START_REF] Kaufmann | 3D Geological Modelling from Boreholes, Cross-Sections and Geological Maps, Application over Former Natural Gas Storages in Coal Mines[END_REF].

The previous exposure does not answer the question of the required level of complexity for a specific study. As all modeling assumptions, this aspect requires the consideration of the specific model purpose. The geomodel-driven approach, therefore, should certainly not be confused with simply adding all detail and complexity that is present in nature, to obtain a "realistic" picture. The aim of modeling is always an abstraction, and this consideration holds here in the same way as in other cases (see also discussions in [START_REF] Ringrose | Reservoir Model Design[END_REF]; [START_REF] Caumon | Geological Objects and Physical Parameter Fields in the Subsurface: A Review[END_REF] and throughout this paper).

Typical input data

Three-dimensional geological modeling generally starts with spatial points, lines or surfaces interpreted from available observations and measurements. The data used in geomodeling is very heterogeneous and includes direct observations of rocks and indirect geophysical measurements. The heterogeneity of data is in itself a challenge in geomodeling: gathering all relevant information can take significant time due to accessibility and numerical storage issues. Another challenge is to assess the reliability of this information, as several sources of uncertainty affect these geoscientific data: position of the measurement, accuracy of the measurement device, volume being investigated. The quantity of interest (e.g., magnetic susceptibility) may also not directly relate to the feature to be modeled. As direct accessibility is limited to the Earth surface and to boreholes and galleries, spatial representativeness of Earth data is also very important. These limitations imply varying degrees of processing and interpretation even before three-dimensional modeling starts. Fig. 4 show some typical types of earth data used in geomodeling, classified according to their volume of investigation and to the degree of interpretation they carry.

More precisely, geomodeling applications typically use several data types (e.g., Kelk, 1992;[START_REF] Kaufmann | 3D Geological Modelling from Boreholes, Cross-Sections and Geological Maps, Application over Former Natural Gas Storages in Coal Mines[END_REF][START_REF] Ringrose | Reservoir Model Design[END_REF]:

• Surface data include pictures and other remote sensing data (which may be combined into digital outcroop models), textual descriptions and drawings, georeferenced structural observations and measurements describing lithology, stratigraphy, unconformity, fault, stratigraphic orientation, fault striae, etc. These data are typically used directly as input for threedimensional modeling, but are often initially translated into maps and cross-sections.

Raw measurements

Interpreted • Geological maps and cross-sections result from two-dimensional interpretations filling the gaps between surface observations. The reliability of this information essentially depends on the distance to actual observations and on the interpreter's skills, which are seldom objectively documented. In practice, 3-D visualization and modeling can help improving these interpretations by allowing the interpreter to consistently map in space.

• Borehole data include both well cores and geophysical well logs. Cores have the same features as surface data, but may be more difficult to interpret, as no rocks are generally observed away from the borehole. Geophysical logs are typically images or result from a local geophysical measurement which reflect physical properties in the neighborhood of the well. Images logs may contain significant structural information to identify layer dips and fractures. Overall, borehole data tend to be sparse and may have representativeness problems, as their location is often geared towards natural resources. The borehole data are generally interpreted in terms of location or range of location for the various geological surfaces.

• Geophysical surveys may either be directly interpreted or processed to provide a three dimensional image. Reflection seismic data particular provide extremely useful information to identify reflectors and the associated discontinuities, which are then interpreted in terms of geological objects (e.g. [START_REF] Jacoby | Gravity Interpretation[END_REF]. Depending on the type of the geophysical method, on the quality and on the spatial coverage of the data, interpretations may bear a variable value in resolving geophysical ambiguities.

• Flow data include flow rates at wells, transient pressure/ height and fluid composition measurements. This type of data does not explicitly include geological information, but can indirectly tell a lot about how subsurface reservoir rocks are connected.

Finally, if we neglect interpretation uncertainties, most geomodeling input data end up as labeled points indicating the most likely location and/or orientation of a particular geological surface. Additionally, we often have access to interval data (indicating the presence or absence of a particular geological feature in some region of space, a common case in geological field campaigns, for example) and trend data (indicating the slope, the fold axis or the dip and dip direction of geological interfaces). Additional information may also be used as input data, for instance the fault displacement or the thickness of a particular stratigraphic layer [START_REF] Mallet | GOCAD: A Computer Aided Design Program for Geological Applications[END_REF], or may be considered as a way to (in)validate a generated model (e.g. [START_REF] De La Varga | Structural geologic modeling as an inference problem: A Bayesian perspective[END_REF].

After the general overview of the considered aspects of geological realism, the conventional workflow with links to geophysics, and the typically used input data and the link to geological knowledge, we now consider the very opposite: the lack of precise knowledge about the subsurface, and the sources and types of uncertainty in geological models.

Classification and origin of uncertainties in Geological Models

"[...] you cannot do inference without making assumptions"

David J. C. MacKay

Uncertainties have been an integral aspect of geological investigations since the early days of mining and exploration. This uncertainty is enshrined in common sayings between miners ("All is dark ahead of the pick"6 ), but it is just as well recognized for the oldest form of geological models, geological maps, where it has always been common to highlight uncertain areas and interface continuations, for example using dashed lines or different line thicknesses and shading (e.g. [START_REF] Greenly | Methods in geological surveying[END_REF].

A notable increase in research around uncertainties in geological maps and models occurred in relation to the development of digital mapping and modeling techniques. Early developments are represented in the work of [START_REF] Mann | Uncertainty in geology[END_REF], and the excellent review by [START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF] about uncertainties in geological mapping and modeling. In the latter paper, the authors describe the relevance of different types of geological knowledge (see also Sec. 1.2 and Fig. 3) and link these types of knowledge to different typical observations during geological mapping (see table 4 in [START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF]. These classification approaches were oriented along descriptive aspects, related to the process of data gathering and model building.

Similar considerations can be made for the typical workflow of geological model construction described above (Sec. 1.2.2). At each step in this modeling process, important decisions have to be taken, partly based on limited information and general aspects of geological knowledge, as well as measurements and observations, all subject to uncertainties. All of these decisions will manifest themselves in the final model (and, logically, also in the approximation of reality) in the form of model uncertainties.

In this section, we first describe sources of uncertainty related to each modeling step. We then provide an overview of uncertainties of typical input data. Finally, we describe approaches to classify uncertainties and propose a terminology based on uncertainty quantification concepts.

Uncertainty in the context of model construction

If we consider the workflow of model construction (Sec. 1.2.2, Fig. 3), then uncertainties enter at several levels: (a) at the definition of the conceptual model and the decision of the mathematical model, (b) at the definition of the structure of the mathematical model, (c) when selecting parameter values for the mathematical model, and finally (d) at the the exact value of the measurements. These steps and related uncertainties are presented in Fig. 5.

If we start from the model construction phase (Sec. 1.1), we first define a conceptual model, and here uncertainties are related to the lack of geological knowledge. This step is also heavily influenced by personal experience and belief (e.g. [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF][START_REF] Bardossy | Evaluation of Uncertainties and Risks in Geology: New Mathematical Approaches for their Handling[END_REF][START_REF] Bond | What do you think this is? "Conceptual uncertainty" in geoscience interpretation[END_REF]. On this basis, we decide on the mathematical model to perform the modeling step (described in more detail below, Sec. 3.2). We now need to decide on the structure of the model itself, for example how many layers to include and how many faults. This step is closely related to the conceptual model, but we emphasize that it is not necessarily the same: a typical example would be the case where the exact number of layers is uncertain, and we would like to evaluate it from data. Examples will be given below. At this stage, the model space contains all of the models that are possible on the basis of the mathematical model and the defined structure. Finally, we consider the data itself and we only keep the models which actually explain the observed data, given the measurement uncertainties. This step is commonly called conditioning to the data. However, even after this step, many models are possible that explain the data in the range of measurement uncertainties.

Several previous approaches attempted a classification of these sources of uncertainty into different types. In a general investigation of uncertainties in geology in the context of Nuclear Waste disposal, [START_REF] Mann | Uncertainty in geology[END_REF] promotes a classification of uncertainty types for geology following the previous (more general) work of [START_REF] Cox | Artifactual Uncertainty in Risk Analysis[END_REF]. His basic types are:

• Type 1: Uncertainty due to observation error, bias and imprecision in the measurement process; could be reduced with additional measurements or better instrumentation;

• Type 2: Inherent variability and stochasticity; could be quantified with repeated observations (e.g. this includes long-term stochastic processes like earthquakes);

• Type 3: Ignorance, lack of knowledge, use of imperfect models and the need for generalizations.

He then describes the general sources of uncertainty and assigns these to the different types (see Table 1 in [START_REF] Mann | Uncertainty in geology[END_REF]. He also makes the important point that a single source can be assigned to several types of uncertainty, which is already indicating that the separation of types is not completely clear.

In his understanding, Type 1 and Type 2 are quantitative and can be described with probability density functions (pdf's) in a probabilistic framework, whereas he considers Type 3 as qualitative uncertainty, which could potentially be described with Fuzzy set theory [START_REF] Zadeh | Fuzzy sets*[END_REF] and also mentions the possibility for model comparison for Type 3 uncertainties. A similar separation into three types of uncertainties has been proposed by [START_REF] Dubois | Fundamentals of Fuzzy Sets (The Handbooks of Fuzzy Sets)[END_REF], and [START_REF] Bardossy | Evaluation of Uncertainties and Risks in Geology: New Mathematical Approaches for their Handling[END_REF] provide an even more basic separation into uncertainties due to natural variability and uncertainties due to human imperfections and incompetency (to which they assign measurement errors, see below). The classifications above were defined for general geological investigations, and Wellmann et al. (2010b) proposed a more specific adaptation to the context of geological modeling. They assigned these uncertainty types in the logical structure of a model building process, where Type 1 uncertainties refer to input data (e.g. the position of a formation boundary), Type 2 to uncertainty in interpolation and extrapolation from known data points, and Type 3 to incomplete knowledge about the structural existence. [START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF] provide a classification oriented along "levels" of mapping campaigns, separated into Data acquisition, Primary interpretation, and Compound interpretation (see Table 4 in [START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF]. In this list, they do not separate quantitative from qualitative uncertainty estimates, but they discuss the fact that uncertainties can partly be quantified directly from measurement errors (e.g. positional error taken from GPS) or obtained from repeated measurements, but that also highlight that uncertainty related to interpretation is more difficult to quantify. Interestingly, they include in their level of "Compound interpretation" the complete essence of uncertainty quantification in the typical question: "How can I quantify the uncertainty associated with this sophisticated interpretive model that I have slowly built up through a long iterative process of data collection and individual primary interpretations?" [START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF]. A similar view to the one of [START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF] has been taken by [START_REF] Mccaffrey | Unlocking the spatial dimension: digital technologies and the future of geoscience fieldwork[END_REF] in a review of digital geological mapping methods. Both works express the hope that future work in the field of digital mapping would allow a better tracking of uncertainties through the entire modeling process.

The uncertainties at different levels of geological mapping in [START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF] are also reflected in the detailed description of uncertainties due to human shortcomings, incompetency or inadequate conditions in [START_REF] Bardossy | Evaluation of Uncertainties and Risks in Geology: New Mathematical Approaches for their Handling[END_REF], where they distinguish problems during sampling (representativeness, sampling patterns, sampling density), insufficient lab measurements, non-measurable properties and vague qualitative descriptions ("rare", "common", "frequent", etc.). The last sources of uncertainty in their list (6-10) refer are also especially relevant in our treatment of geological modeling, and this requires some more detailed thought. As discussed in Section1.1, geological models are simplifications of reality and these simplifications and generalizations are sources of uncertainty. This includes all aspects related to the general background knowledge (see Fig. 3). These uncertainties also fall into the group of "beliefs" by [START_REF] Zimmermann | An application-oriented view of modeling uncertainty[END_REF]. [START_REF] Bardossy | Evaluation of Uncertainties and Risks in Geology: New Mathematical Approaches for their Handling[END_REF] describe that checking these generic (conceptual) models is the most difficult task (in geological model construction). The next sources of uncertainty are related to the choice and the correct application of the mathematical model, and, finally, uncertainties related to the conclusions drawn from a specific geological modeling investigation.

In practical applications, it may furthermore well be the case that uncertainties are apparently increased when more data are considered. This is especially the case when some data or concepts are in conflict (e.g. [START_REF] Zimmermann | An application-oriented view of modeling uncertainty[END_REF]. However, this is also a strong indication that either the model itself is not suitable or that some information which is not relevant to the specific model is being used. This type of uncertainty is also referred to as disinformation in the hydrology literature (e.g. [START_REF] Beven | Facets of uncertainty: epistemic uncertainty, non-stationarity, likelihood, hypothesis testing, and communication[END_REF].

Uncertainty of data and measurements

As described in the previous section, uncertainties enter at all stages of model construction (Fig. 5).

Here, we provide a brief overview of typical sources of uncertainty related to input data to geological models as described in Sec. 1.3, suitable probability distributions, and methods of expert elicitation to pool judgment from multiple experts.

Uncertainties related to different data types

An immediate input to geological models are geological field observations and measurements, either from dedicated field campaigns, or derived from geological maps (Fig. 3). Sources of uncertainty in this context have partly been described in the literature (as described in the previous section, e.g. [START_REF] Mann | Uncertainty in geology[END_REF][START_REF] Bardossy | Evaluation of Uncertainties and Risks in Geology: New Mathematical Approaches for their Handling[END_REF]Bowden, 2004a;[START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF]Wellmann et al., 2010a).

Especially the estimation of uncertainties in geological mapping, evaluated in [START_REF] Jones | Digital field data acquisition: towards increased quantification of uncertainty during geological mapping[END_REF] and [START_REF] Mccaffrey | Unlocking the spatial dimension: digital technologies and the future of geoscience fieldwork[END_REF] has partly been described with uncertainties in subsequent 3-D geological models in mind.

LiDAR and UAV devices, as well as diverse types of remote sensing data, are now increasingly used to gather geological surface information provide an unprecedented amount of information that can be expected to constrain geological structural models [START_REF] Bilotti | Quantitative structural analysis with stereoscopic remote sensing imagery[END_REF][START_REF] Fernández | Three-dimensional reconstruction of geological surfaces: An example of growth strata and turbidite systems from the ainsa basin (pyrenees, spain)[END_REF]Caumon et al., 2013a;[START_REF] Bistacchi | Photogrammetric digital outcrop reconstruction, visualization with textured surfaces, and three-dimensional structural analysis and modeling: Innovative methodologies applied to fault-related dolomitization (Vajont Limestone, Southern Alps, Italy)[END_REF][START_REF] Vollgger | Mapping folds and fractures in basement and cover rocks using UAV photogrammetry, Cape Liptrap and Cape Paterson, Victoria, Australia[END_REF]. [START_REF] Cawood | LiDAR, UAV or compass-clinometer? Accuracy, coverage and the effects on structural models[END_REF] performed a detailed analysis comparing LiDAR, UAV and conventional field observations and measurements with respect to accuracy and the effect on subsequent geological models. Although limited to a specific case study, this is an important step forward and we can expect to see more work in this direction in the future. Another highly active field of research is the use of remote sensing data to decipher geological surface information. For an introduction to the topic, see [START_REF] Lary | Machine learning in geosciences and remote sensing[END_REF], the comprehensive work on Remote Predictive Mapping (RPM) by the Geological Survey of Canada (e.g. [START_REF] Harris | Remote Predictive Mapping (RPM): A not so New Paradigm for Mapping Canada's North[END_REF][START_REF] Schetselaar | Remote Predictive Mapping 1. Remote Predictive Mapping (RPM): A Strategy for Geological Mapping of Canada's North[END_REF][START_REF] Harris | Remote Predictive Mapping 3. Optical Remote Sensing -A Review for Remote Predictive Geological Mapping in Northern Canada[END_REF] and the interesting comparison by Cracknell andReading (2013, 2014) for the specific aspect of geological mapping on the basis of remote sensing data.

Two aspects are receiving repeated attention in the discussion of uncertainties in geological interface points for geological mapping and modeling: (a) data density, and (b) geological complexity. An example for the effect of data density of interface points on a geological model is provided in the work of [START_REF] Putz | Three-dimensional model and late stage warping of the Plattengneis Shear Zone in the Eastern Alps[END_REF] for a model of a shear zone in the Eastern Alps. In a more recent work, [START_REF] Carmichael | Method and analysis for the upscaling of structural data[END_REF] evaluate the effect of orientation data density on model construction and they provide a method for an automatic upscaling of orientation measurements for an optimized use in geological models. The second aspect of geological complexity is a lot more difficult to capture. Several studies have proposed measures of complexity (e.g. [START_REF] Ford | Evaluating geological complexity and complexity gradients as controls on copper mineralisation, Mt Isa Inlier[END_REF][START_REF] Zhao | Mapping complexity of spatial distribution of faults using fractal and multifractal models: vectoring towards exploration targets[END_REF], also partly with direct consideration of geological models (e.g. [START_REF] Lelliott | A structured approach to the measurement of uncertainty in 3D geological models[END_REF][START_REF] Pellerin | Elements for measuring the complexity of 3D structural models: Connectivity and geometry[END_REF]. Interestingly, [START_REF] Schweizer | Uncertainty assessment in 3-D geological models of increasing complexity[END_REF] analyzed the effect of increasing model complexity on the estimation of model uncertainty. Also important to note is that geological complexity is closely related to geological history. However, there is yet no unifying measure of geological complexity to be directly used in a subsequent uncertainty study. This is certainly an interesting aspect of future research (see Discussion, Sec. 5).

Other typical input data are borehole data, either in the form of well cores or geophysical well logs and the interpretation of geophysical wireline logs is a major topic in itself. At this point, we specifically want to mention quantitative approaches based on statistical rock physics [START_REF] Mavko | A rock physics strategy for quantifying uncertainty in common hydrocarbon indicators[END_REF]; [START_REF] Mukerji | Statistical rock physics: Combining rock physics, information theory, and geostatistics to reduce uncertainty in seismic reservoir characterization[END_REF], as they have been applied in automatic extraction of geological attributes from well logs (e.g. [START_REF] Eidsvik | Estimation of Geological Attributes from a Well Log: An Application of Hidden Markov Chains[END_REF], in recent studies also with an estimation of uncertainties (e.g. [START_REF] Grana | Quantitative log interpretation and uncertainty propagation of petrophysical properties and facies classification from rock-physics modeling and formation evaluation analysis[END_REF][START_REF] Grana | Joint facies and reservoir properties inversionJoiFRI[END_REF].

Probably the most dominantly used geophysical data in geological models are reflection seismic data, either classically in the form of 2-D seismic sections, but increasingly also in the form of 3-D seismic cubes [START_REF] Dorn | Modern 3-D seismic interpretation[END_REF][START_REF] Jackson | 3D Seismic Interpretation Techniques: Applications to Basin Analysis[END_REF] and it can be argued that reflection seismic data sets provide the closest form of a direct "imaging" of the subsurface structures [START_REF] Davies | 3D Seismic Technology: Are We Realising Its Full Potential?[END_REF]. Uncertainties in seismic data have been analyzed in numerous studies (e.g. [START_REF] Thore | Structural uncertainties: Determination, management, and applications[END_REF][START_REF] Glinsky | Integration of uncertain subsurface information into multiple reservoir simulation models[END_REF][START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF][START_REF] Weinzierl | Volumetric impact of fault perturbation in the first Fresnel zone[END_REF] and can mainly be assigned to processing steps (stacking, migration, time-depth-conversion) and the subsequent interpretation step.

A significant body of the early work on uncertainties in geological models has been linked to uncertainties in seismic processing, specifically time-depth-conversion (e.g. [START_REF] Abrahamsen | Stochastic Models for Seismic Depth Conversion of Geological Horizons[END_REF][START_REF] Abrahamsen | Bayesian Kriging for Seismic Depth Conversion of a Multi-Layer Reservoir[END_REF][START_REF] Samson | Quantifying the Impact of Structural Uncertainties on Gross-Rock Volume Estimates[END_REF]. The advantage of this source of uncertainty is that prior knowledge exists in the form of expected functional relationships between time and depth domain, as well as the expected seismic velocities.

The relevance of this type of uncertainty is also reflected in more recent work. [START_REF] Li | Assessing seismic uncertainty via geostatistical velocity-model perturbation and image registration: An application to subsalt imaging[END_REF] developed a workflow combining geostatistical modeling and seismic imaging to asses the effect of velocity model uncertainty on generated images. The authors generate multiple image realizations on the basis of randomized velocity models. Interesting is then the link to subsequent seismic interpretations. The authors describe the adaptation of an image registration approach (the Thirion demon algorithm [START_REF] Thirion | Image matching as a diffusion process: an analogy with Maxwell's demons[END_REF] to track a manual interpretation of one seismic image (performed by a seismic interpreter) to all other randomly generated image realizations. [START_REF] Osypov | Model-uncertainty quantification in seismic tomography: method and applications[END_REF] also highlight the impact of uncertainties in seismic velocity models by re-migrating horizon interpretations.

Recent work has also focused on the effect of interpretation uncertainties, mainly related to the interpretation of seismic measurements, but also for potential-field data. Uncertainties in seismic interpretation has received increasing attention in recent years (e.g. [START_REF] Bond | What do you think this is? "Conceptual uncertainty" in geoscience interpretation[END_REF][START_REF] Bond | What makes an expert effective at interpreting seismic images?[END_REF]Alcalde et al., 2017a,b) and it is well recognized that seismic interpretations are strongly linked to the experience and knowledge of the expert (an aspect that goes hand in hand with the hermeneutic interpretation of geology, see Frodeman, 1995, and discussion). The studies also found a significant influence of expertise in specific tectonic settings, as well as the specific techniques used in the process itself, on the subsequent seismic interpretation, partly leading to biased interpretation. Studies using borehole data revealed similar patterns [START_REF] Bond | Structural model creation: the impact of data type and creative space on geological reasoning and interpretation[END_REF].

As these interpretations are often a primary source of input data to geological models, their uncertainty can be considered as highly relevant to the subsequent geological model uncertainty. However, even for practices where reservoir properties are directly estimated from seismic data, studies show that these workflows are still affected by the uncertainty of the primary structural interpretation [START_REF] Rankey | That's why it's called interpretation: Impact of horizon uncertainty on seismic attribute analysis[END_REF]. Interpretation uncertainty is also not resolved yet by autopicking algorithms (e.g. [START_REF] Herron | Pitfalls in horizon autopicking[END_REF], although increasing development on the basis of machine learning algorithms will potentially lead to improved algorithms in the future. Nonetheless, to manage possibly different interpretations and get a sense about how they quantitatively honor the initial seismic data, several authors have proposed to consider alternative geological structural models and then either generate synthetic seismic images by convolution (e.g., [START_REF] Lallier | Relevance of the stochastic stratigraphic well correlation approach for the study of complex carbonate settings: Application to the Malampaya buildup (Offshore Palawan, Philippines)[END_REF][START_REF] Botter | Seismic characterisation of fault damage in 3D using mechanical and seismic modelling[END_REF], or forward model wave propagation to generate synthetic seismograms [START_REF] Irakarama | Appraising structural models using seismic data: Problem and challenges[END_REF]. Computing misfit functions remains, however a significant challenge as cycle skipping is likely to occur in the presence of large interpretation uncertainties owing for instance to poor illumination.

Gravity and magnetic measurements are commonly used geophysical measurements for geological models, both as input to geological models (derived from interpreted geophysical maps), as well as conditions for joint geological and geophysical inversion. A detailed treatment for the interpretation of gravity data and the related uncertainties is provided in [START_REF] Jacoby | Gravity Interpretation[END_REF]. Also here, recent methodological developments focus on diverse machine learning and optimization approaches with an additional assessment of uncertainties (e.g. [START_REF] Pallero | Gravity inversion and uncertainty assessment of basement relief via Particle Swarm Optimization[END_REF][START_REF] Pallero | 3D gravity inversion and uncertainty assessment of basement relief via Particle Swarm Optimization[END_REF]. Sources of uncertainty in specific case of aeromagnetic data are, for example, discussed in [START_REF] Holden | Identifying structural complexity in aeromagnetic data: An image analysis approach to greenfields gold exploration[END_REF] and [START_REF] Aitken | Semiautomated quantification of the influence of data richness on confidence in the geologic interpretation of aeromagnetic maps[END_REF]. These authors combine estimates of geological data quality, density and interpretability, as well as orientation uncertainty to derive a measure of confidence in aeromagnetic interpretations. They also express the possibility to use these maps to obtain weighting factors for subsequent geophysical inversions.

Probability distributions for various sources of uncertainties

The various sources of uncertainty related to different input data types need to be formalized in some form to be used in a subsequent modeling study. This aspect alone can be difficult, as geological information often contains subjective elements due to the interpretative character of geological observations (e.g. [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF][START_REF] Wood | Geological prior information and its applications to geoscientific problems[END_REF]. If uncertainties can be quantified, then the most common description is in the form of probability distributions, both when relating these uncertainties to outcomes in random experiments, as well as for a description of degrees of belief (e.g. [START_REF] Gelman | Bayesian data analysis[END_REF][START_REF] Tarantola | Popper, Bayes and the inverse problem[END_REF][START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF].

Probability distributions assigned to the input data and parameters have a direct effect on the subsequent step of uncertainty quantification. Ideally, these distributions should be derived from the primary input data (e.g., device-specific measurement uncertainties, uncertainties in time-depth conversion in seismic images or from the segmentation of wireline logs) or from repeated measurements (e.g., orientation measurements of planar geological features). In this case, a probability distribution can then directly be obtained on the basis of this information, possibly also from subsequent processing steps, such as probabilistic plane fits to structural lineaments (e.g., [START_REF] Seers | Probabilistic constraints on structural lineament best fit plane precision obtained through numerical analysis[END_REF].

For the cases where more subjective information about uncertainty is included, different types of distributions may be suitable. For a general overview of suitable distributions, see [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF], [START_REF] Gelman | Bayesian data analysis[END_REF] and [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]. An overview for the case of geological uncertainties is also provided in Wellmann et al. (2010a). Some sources of geological uncertainty can reasonably be described with a normal distribution, for example when a position of an interface can not exactly be determined due to a gradual transformation (e.g., in wireline logs). Another common case in geological mapping studies is that only outcrops are available where specific lithological units can be determined, but where the contact is not visible. In these cases, a uniform distribution could be applicable, when no other information about the contact is present. However, it should be noted that uniform distributions to restrict ranges or bounds should generally be avoided for high-dimensional cases, as uniform samples in a hypercube will be concentrated at the edges (e.g. [START_REF] Curtis | Prior information, sampling distributions, and the curse of dimensionality[END_REF]. When sufficient data are present and deemed representative, non-parametric distributions (i.e. histograms) may also be used, possibly after a distribution smoothing step.

For the specific case of orientation measurements, distributions should be assigned according to a suitable directional statistical distribution from the family of Fisher-Bingham distributions [START_REF] Fisher | Dispersion on a Sphere[END_REF][START_REF] Fisher | Statistical analysis of spherical data[END_REF][START_REF] Fisher | Approaches to Uncertainty in Spatial Data[END_REF]. The analogue to a normal distribution on a sphere is the von Mises-Fisher distribution, which is appropriate for orientations as spherical vector data [START_REF] Davis | Statistics and data analysis in geology, 3rd Edition[END_REF][START_REF] Pakyuz-Charrier | Monte Carlo simulation for uncertainty estimation on structural data in implicit 3-D geological modeling, a guide for disturbance distribution selection and parameterization[END_REF]. Uncertainties related to parameters of the applied mathematical interpolation model can potentially be derived when additional data is available (e.g. [START_REF] Aug | 3D geological modelling and uncertainty: The potential-field method[END_REF][START_REF] Gonçalves | A machine learning approach to the potential-field method for implicit modeling of geological structures[END_REF].

Expert elicitation

Geological prior information is often subjective (e.g., [START_REF] Wood | Geological prior information and its applications to geoscientific problems[END_REF]. This subjectivity should not be seen as an impediment to the use of quantitative methods, but rather as a motivation to develop suitable methods to reduce the effect of subjectivity, as appropriately described by [START_REF] Curtis | The science of subjectivity[END_REF]. One important aspect of subjectivity is the influence of prior experience and expectation, for example in the interpretation of seismic sections (e.g. [START_REF] Rankey | That's why it's called interpretation: Impact of horizon uncertainty on seismic attribute analysis[END_REF][START_REF] Bond | What makes an expert effective at interpreting seismic images?[END_REF]. Especially in the case of interpretations of uncertainty, cognitive biases such as over-confidence, anchoring, and availability bias can play an important role [START_REF] Tversky | Judgment under Uncertainty: Heuristics and Biases[END_REF].

Even considering the potential methods to consider uncertainties as input to geological models, it is important to stress that unfortunately not all knowledge can easily be described formally, let alone in a quantitative manner (See also discussion in Sect. 3.2.4. The difficulty to communicate knowledge is often expressed in the form of different stages of personal knowledge, following the pioneering work of [START_REF] Polanyi | Personal knowledge[END_REF]:

• Explicit Knowledge: knowledge, which is readily available and communicated, e.g. accessible in digital form in tables or data sets;

• Implicit Knowledge: in the mind of an expert (or, generally, any person) and knowledge that could be made explicit and communicated, given time and effort;

• Tacit knowledge: this is the difficult aspect of knowledge, which exists in the mind of a person, can not be easily communicated and made explicit; sometimes, this is what people would consider "gut feeling".

One method to make personal knowledge better available, and to address issues of bias and to reduce subjectivity are structured elicitation processes, going back to the work of [START_REF] Kidd | Knowledge Acquisition[END_REF] and [START_REF] Ford | Expert knowledge elicitation to improve formal and mental models[END_REF]. These methods have since been applied in geosciences (e.g. [START_REF] Curtis | Optimal elicitation of probabilistic information from experts[END_REF], also specifically with the aim to obtain probabilistic information from geological interpretations (e.g., [START_REF] Cooke | Experts in Uncertainty: Opinion and Subjective Probability in Science[END_REF][START_REF] Baddeley | An introduction to prior information derived from probabilistic judgements: elicitation of knowledge, cognitive bias and herding[END_REF][START_REF] Bowden | Optimal elicitation of probabilistic information from experts[END_REF][START_REF] Polson | Dynamics of uncertainty in geological interpretation[END_REF]. More recent interesting approaches also include applications of expert elicitation to estimate subjective probabilities of natural hazards [START_REF] Aspinall | A route to more tractable expert advice[END_REF][START_REF] Runge | An Interactive Tool for the Elicitation of Subjective Probabilities in Probabilistic Seismic-Hazard Analysis[END_REF], where also here the goal to include quantified uncertainties in the decision process is clearly stated.

Geological modeling methods

"It is not his possession of knowledge, of irrefutable truth, that makes the man of science, but his persistent and recklessly critical quest for truth"

Karl R. Popper

Three-dimensional geomodeling is about modeling the topology (connectivity), the geometry and the properties of a particular region of the Earth [START_REF] Mallet | Geomodeling[END_REF]. In general, this process does not require to reconstruct the history of the area through time, but rather to interpolate between sample data while following some essential rules and principles. This lack of consideration for geological processes can be considered as a weakness: what is the scientific value of interpolation as compared to process-based models that try to explicitly reproduce coupled processes at geological time scales?

A first element of answer lies in the extreme difficulty to find appropriate physical laws and boundary conditions in process-based models that run over geological time scales. Honoring spatial observations with these methods calls for solving an extremely ill-posed inverse problem. Interpolation appears, therefore, as a viable and computationally efficient alternative to obtain a consistent picture of the subsurface while honoring a large variety of observations (Sec. 1.3). More philosophically, geomodeling methods form a very good example of the "distinctive set of logical procedures" discussed by [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF] to describe geological reasoning. We fully agree with Frodeman that interpretation and historical considerations are essential to geological science in general and to geomodeling in particular.

In practice, geological history and interpretive concepts have been translated into geomodeling as: (1) ways to define and represent geological formations mathematically and numerically, and (2) workflows to build models honoring available information. In this section, we start by describing the typical input data used to create geomodels. We then we discuss possible methods to represent geological models, trying to establish the common mathematical approaches to the problem and how they relate to geological concepts. We then review the existing ways to interpolate structures from available data, before discussing validation strategies.

Geological models: essential characteristics and concepts

Geomodel representations

The goal of this section is to describe how complex geological domains can be described in mathematical and numerical terms. Geomodels in the literature are generally described in terms of processes and modeling workflows using various types of meshes, grids or meshless methods. To address subsurface uncertainty and discuss how geological concepts can be translated into mathematical terms, we think it useful to use a more abstract description, which we describe in the following. We will then review existing implementations in Sec. 3.2. In this section, we build on the same formalism as [START_REF] Sambridge | Transdimensional inference in the geosciences[END_REF] and further explain and elaborate the ideas expressed by [START_REF] Caumon | Geological Objects and Physical Parameter Fields in the Subsurface: A Review[END_REF].

Figure 6 shows various geomodels representing possible seismic wave velocity fields in the subsurface, and a corresponding set of one-dimensional profiles showing the true field (in red) and how the various representations approximate it (in black). These various options were all generated from sparse measurements and interpolated given the same structural geometry. They illustrate the basic principle used in numerical solid and physical modeling, which is to decompose the area of interest into a boundary representation or into a structured or unstructured mesh [START_REF] Mäntylä | An Introduction to Solid Modeling[END_REF].

In computational physics, all these geomodel representations can be described mathematically as

(x) = K k=1 m k ϕ k (x) (1)
to obtain a value of interest m at all possible locations x in a model domain (i.e. an approximation of reality, in the sense of Fig. 3).

In Figure 6, m is the velocity field, K is the number of nodes (or grid blocks) bearing the static scalar velocity values m k at the nodes (or grid blocks), and ϕ k denotes the interpolation functions on the grid. The functions ϕ k define the geometric structure of the parametrization. There is presently no general consensus about the type of functions ϕ k (x) to be used to best reflect what is known about the domain of interest. Most works choose a particular representation based on implementation comfort or on physical principles. The goal of this section and of Fig. 6 is to discuss about the adequacy of some choices to represent geological media.

For example, Figure 6(a) discretizes the space by regular orthogonal grid blocks which correspond to piecewise constant basis functions. It overlays a fine decomposition representing the reflection seismic amplitudes to a coarser mesh representing velocities. The structural geometries in all models were obtained by expert-based interpretations if the seismic image. This Cartesian grid representation is very simple and widely used to represent geomodels in may fields, for instance in mining applications (e.g. McGaughey, 2006a). Depending on the grid resolution, this representation can approximate layer boundaries and smear petrophysical contrasts. In Figure 6(b), the basis functions are stretched or squeezed vertically so that they align on layer boundaries. This corresponds to a classical modeling choice made in reservoir modeling, corresponding to corner-point stratigraphic grids [START_REF] Farmer | Geological Modelling and Reservoir Simulation[END_REF]. However, it is clearly visible that faults are approximated by stair-steps in this model representation.

Figure 6(c) shows another option where a linear tetrahedral mesh is used to better approximate fault geometry. The tetrahedra are conformable to the fault surfaces, but cross horizons. The velocity field is linearly interpolated between tetrahedral nodes.

Once the basis functions have been chosen, the unknown coefficients m k can be determined. Depending on the model purpose, m and m k may represent one or several petrophysical parameters such as electric conductivity, seismic velocity, Lamé parameters, rock facies, rock composition, etc. Geostatistics is an essential framework to model the values m k everywhere in space from available measurements. As very good reviews and textbooks on geostatistics already exist [START_REF] Hu | Multiple-point geostatistics for modeling subsurface heterogeneity: A comprehensive review: REVIEW OF MULTIPLE POI[END_REF][START_REF] Bosch | Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: A review[END_REF][START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty, 2nd Edition[END_REF][START_REF] Bourges | A practical review of geostatistical processing applied to geophysical data: Methods and applications: Geostatistical processing applied to geophysical data[END_REF]Pyrcz and Deutsch, 2014a), we only stress here the basic idea: all the unknown parameters m k are treated as correlated random variables. This general framework allows to describe spatial correlations for a given variable, and both co-located and spatial correlations between several petrophysical variables. Geostatistics also provides a framework to deal with the support (or volume) associated to each variable [START_REF] Journel | Mining Geostatistics[END_REF][START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty, 2nd Edition[END_REF]), but it is not widely used outside the field of mining geostatistics.

To describe complex domains, geomodeling generally follows the classical divide-and-conquer approach that is typically used in geological mapping. The domain is first subdivided into rock units that have some historical significance. Each unit may then be further subdivided into sub-regions that describe smaller features. Finally, the domain of interest is subdivided into a set of J regions {R 1 , . . . , R J }. The physical properties of interest can be represented within each unit (or sub-unit) R j by a set of basis functions ϕ j k (x), k = {1, . . . , K j } whose support covers and is completely included in R j . The descriptions of the values of interest everywhere are then given by: This type of representation is used for instance in the boundary element method, where an analytical petrophysical description is used within each subdomain. This type of approach has also been used a lot in geophysics for velocity modeling and ray tracing (e.g. [START_REF] Gjøystdal | Computer Representation of Complex 3-D Geological Structures using a new "Sold Modeling[END_REF][START_REF] Guiziou | 3-D seismic reflection tomography on top of the GOCAD depth modeler[END_REF] or for potential field modeling and inversion (e.g. [START_REF] Wang | Inversion for magnetic anomalies of arbitrary three-dimensional bodies[END_REF][START_REF] Jacoby | Gravity Interpretation[END_REF]. As borehole measurements generally highlight that heterogeneities within regions cannot be modeled analytically, a discrete representation of values is often used, as shown in the models of Fig. 6. In practice, the region boundary needs to be discretized on the model mesh (e.g. Fig. 6a), or the model mesh is designed to conform to the region boundary, so that the mesh is aligned to layers (Fig. 6b) or faults (Fig. 6c).

m(x) = J j=1 K k=1 m k ϕ j k (x). ( 2
) (a) Two horizon surfaces (b) Hard truncation (c) Soft truncation E E + E - H H + H - E H H + H - E + E H H + H - E +
As region boundaries often correspond to different rock types, it is common to extrapolate observed petrophysical values by treating region boundaries as discontinuities. This explains why sharp changes can be seen across horizons and faults in Figures 6(a),(b) and (d) and across faults in Figure 6(c).

It can also be useful to represent each region by a an indicator variable 1 j , equal to 1 inside R j and to 0 outside. Equation ( 1) is then modified into:

m(x) = J j=1 1 j (x) K k=1 m k ϕ k (x) . (3)
The indicator variable 1 j (x), is itself a solid model which can be represented by the boundary surfaces of the region j or some other solid modeling technique (e.g. [START_REF] Mäntylä | An Introduction to Solid Modeling[END_REF]. Fig. 7 illustrates this on the two-dimensional example of a horizon H eroded by an erosion curve E. Both curves separate the 2-D domain into four subdomains E + , E -, H + and H -described by four binary indicator functions 1 E + , 1 E -, 1 H + and 1 H -, respectively. These functions are equal to 1 in their respective domain and to zero elsewhere. The erosion of the horizon curve H corresponds to a Boolean operation stating that no point of H, H + and H -may exist in the region E + above E. Let us denote the non-eroded horizon H as a linear combination of I basis functions χ H i (Fig. 7a, right):

H = x ∈ R 3 | I i=1 c i χ H i (x). (4)
In this example, the basis functions χ H i may for instance correspond to cubic Hermite polynomials centered on the curve nodes.

A first option to implement erosion is to explicitly trim the horizon H by the erosion curve E. As in Eq. ( 2)), this amounts to reparametrizing H by I conformable basis functions χ H i (Fig. 7b):

H erod = x ∈ R 2 | I i=1 c i χ H i (x). ( 5
)
Alternatively, eroding the horizon H may be achieved by canceling the Horizon H outside of the half-space E -below the unconformity E. As in Eq. ( 3), this amounts to multiplying Eq. 4 by the indicator function 1 E -(Fig. 7c):

H erod = x ∈ R 2 | I i=1 1 E -c i χ H i (x). ( 6
)
The same rules can be used to account for all geological discontinuities, as will be further explained in Sec. 3.1.2 and Sec. 3.1.3. As in the above two-dimensional example, all the models shown in Figure 6 started from a boundary representation of geological interfaces. Each variant corresponds to different choices for approximating the geometry of these boundaries: faces separating the blocks of a Cartesian grid (a) or of a curvilinear grid (b), or faces of tetrahedra aligned on the faults (c and d). As compared to (2), this formulation does not necessarily require to use a mesh conformal to the region boundary. This makes it easier to manage geometric changes of the region boundaries. This is one of the principles used for instance in the Extended Finite Element Method (XFEM) to model the growth of a propagating fracture [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model: NON-PLANAR 3D CRACK GROWTH-PART I[END_REF].

Both Eqs. ( 2) and ( 3) correspond to a "cookie-cutter" strategy, widely used in subsurface modeling to model strata (and facies) in reservoir grids (e.g. [START_REF] Swanson | A New Geological Volume Computer Modeling System for Reservoir Description[END_REF][START_REF] Johnson | Putting Geology Into Reservoir Simulations: A Three-Dimensional Modeling Approach[END_REF][START_REF] Lemon | Building solid models from boreholes and user-defined crosssections[END_REF]Pyrcz and Deutsch, 2014a). It was also used early on for velocity modeling by structurally complex domains by [START_REF] Gjøystdal | Computer Representation of Complex 3-D Geological Structures using a new "Sold Modeling[END_REF]. In that case, the indicator functions for the J regions are determined by boundary representation and Boolean operations. In basin modeling, [START_REF] Mello | Techniques for including large deformations associated with salt and fault motion in basin modeling[END_REF] proposed to create a boundary representation of the regions R J by stitching geological interfaces. Overall, this nested modeling strategy is essential to reflect the hierarchy of heterogeneities encountered in the subsurface.

Another general concept is to use some intermediate variable(s) u(x) to represent some geometric aspects of subsurface models. In Eq. ( 1), this amounts to replacing the position variable x by one or several variables u(x), yielding

m(x) = m(u(x)) = K k=1 m k ψ k (u(x)), (7)
where u(x) can be defined by another linear combination of basis functions χ i :

u(x) = I i=1 u i χ i (x). ( 8
)
As in Russian dolls, these geometric basis functions χ i (x) may also be truncated or designed to be conformable to some region boundaries.

The velocity values in Figure 6(b) illustrate this approach, as the corner-point grid bears an indirect representation of the u(x) coordinates by the regular indexing of the grid blocks. This indexing was used to extrapolate the velocity values from sparse observations so that directions of velocity anisotropy followed stratigraphic orientations. Figure 6(d) describes a more general function u(x) = [u(x), v(x), t(x)] representing the stratigraphic geometry. t(x) is shown on the onedimensional profile. In this example, relatively coarse linear basis functions are used to represent the intermediate variable u(x) on a tetrahedral mesh. The finer regular piecewise constant basis functions represent properties m(u).

This type of approach is common in the modeling of stratigraphic domains, as they allow to align petrophysical features on the structural directions, see for example [START_REF] Swanson | A New Geological Volume Computer Modeling System for Reservoir Description[END_REF]; [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF]; [START_REF] Lomask | Flattening without picking[END_REF]; [START_REF] Wu | Building 3D subsurface models conforming to seismic structural and stratigraphic features[END_REF].

To summarize, the examples shown in Figure 6 illustrate the interest of combining different sets of basis functions to decouple the representation of model geometry and the representation of model properties. Geological concepts are introduced as specific basis functions, as truncation of some basis functions, or as composition of basis functions. The consequences of these model choices are important, because they impact the ability of the models to approximate the true subsurface features. As the number of parameters and the relationships between these parameters change from one parametrization to the next, these choices also have a strong impact on computational performance, and on the structure of model space in uncertainty quantification and subsurface inverse problems.

Stratigraphy and representations of geological time

So far, we have mainly described subsurface parametrization in a mathematical way, starting from the petrophysical parameter field of interest. We now proceed to discuss model choices by starting from essential geological concepts and how they translate into constraints or into the parametrizations presented in the previous section.

As discussed in the introduction (Sec. 1.1), geology puts a significant focus on chronological considerations to analyze and interpret observations. This temporal (or historical) reasoning can partly be integrated into geomodeling strategies.

As shown in Fig. 1, a part of the history of rock formation and subsequent deformation phases is classically captured by a stratigraphic column. A stratigraphic column (Fig. 8) describes the vertical stacking of units, the nature and the characteristics of the rocks in each unit, and the type of stratigraphic surfaces bounding each unit. Following the superposition principle [START_REF] Steno | De Solido Intra Solidium Naturaliter Contento Dissertationis Prodromus[END_REF], the vertical axis in a stratigraphic column has both a spatial and temporal meaning. This description, valid at a certain scale, has significantly been used in geological modeling, as it organizes geological knowledge in a systematic way (e.g. [START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] Zhu | Building 3D solid models of sedimentary stratigraphic systems from borehole data: An automatic method and case studies[END_REF][START_REF] Perrin | Shared Earth Modeling: Knowledge Driven Solutions for Building and Managing Subsurface 3D Geological Models[END_REF][START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF]. Some freely available kinematic models exist to explicitly model deposition and deformation phases, for instance Noddy [START_REF] Jessell | Structural geophysics: Integrated structural and geophysical modelling[END_REF], PyNoddy [START_REF] Wellmann | Pynoddy 1.0: An experimental platform for automated 3-D kinematic and potential field modelling[END_REF] and the online web app VisibleGeology. These kinematic models are simple and computationally efficient; they can generate a large diversity of three-dimensional models, but they cannot directly honor the diversity of observations described in Section 1.3.

In structural modeling, the geometry of strata is generally treated by spatial interpolation (Sec. 3.2), while the unconformities between strata are modeled by independent interpolations followed by truncation. For example, the basis functions for horizons or rock units are canceled above an erosion surface (Fig. 7). Conversely, the basis functions corresponding to horizons or rocks younger than a given baselap surface are canceled below that baselap surface. Geological time is, therefore, treated indirectly by selecting the type of truncation operation which reflects the stratigraphic history. Figure 8: The stratigraphic column: a summary of geological history recorded in the sediments. The conformity or truncation of beds within larger stratigraphic units can be related to tectonic and stratigraphic history. In this example, beds within the units correspond to chronostratigraphic surfaces.

More generally, cross-cutting relationships between rock units are essential to determine the relative ages of geological objects. For instance, a layer cut by a dyke or by an evaporitic intrusion is older than the intrusive unit. Again, this can be handled by independent interpolations of the boundary surfaces of layers and intrusions and truncation of older rocks by younger ones. This reasoning, corresponding to Eq. 3 and illustrated in Fig. 7c, and has been implemented on various numerical representations [START_REF] Gjøystdal | Computer Representation of Complex 3-D Geological Structures using a new "Sold Modeling[END_REF][START_REF] Johnson | Putting Geology Into Reservoir Simulations: A Three-Dimensional Modeling Approach[END_REF][START_REF] Swanson | A New Geological Volume Computer Modeling System for Reservoir Description[END_REF][START_REF] Zeng | Stratal slicing; Part 1, Realistic 3-D seismic model[END_REF][START_REF] Lemon | Building solid models from boreholes and user-defined crosssections[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] Zhu | Building 3D solid models of sedimentary stratigraphic systems from borehole data: An automatic method and case studies[END_REF][START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF]. Geological surfaces can also be created or edited so as to directly enforce the termination of surfaces on other surfaces following the prescribed relationships [START_REF] Mallet | GOCAD: A Computer Aided Design Program for Geological Applications[END_REF][START_REF] Mello | Techniques for including large deformations associated with salt and fault motion in basin modeling[END_REF][START_REF] Caumon | Constrained modifications of non-manifold b-reps[END_REF][START_REF] Brandel | Automatic Building of Structured Geological Models[END_REF][START_REF] Perrin | Knowledge-Driven Applications for Geological Modeling[END_REF], which corresponds to implementations of Eq. ( 2) illustrated in Fig. 7b.

Another way to treat time in structural modeling is by explicitly or implicitly representing the geological time t(x) as a function of space (Eq. 7). A classical way of doing so is to include layering information within stratigraphic columns. Layering essentially informs the interpolation between the bounding surfaces of each stratigraphic unit. Three main interpolation functions have been defined, depending on the type of horizons bounding the units of a stratigraphic column [START_REF] Swanson | A New Geological Volume Computer Modeling System for Reservoir Description[END_REF][START_REF] Johnson | Putting Geology Into Reservoir Simulations: A Three-Dimensional Modeling Approach[END_REF][START_REF] Zeng | Stratal slicing; Part 1, Realistic 3-D seismic model[END_REF][START_REF] Mallet | Geomodeling[END_REF][START_REF] Perrin | Shared Earth Modeling: Knowledge Driven Solutions for Building and Managing Subsurface 3D Geological Models[END_REF]: proportional between top and base (conformable: unit B in Fig. 8), eroded (conformable to base: unit D in Fig. 8), baselap (conformable to top, Fig. 8). A fourth style combining eroded and baselap also exists (unit D in Fig. 8), and can be modeled by splitting the considered unit by a conformable horizon. In simple cases such as described in Fig. 2a, geological time t(x) can be determined by simple linear operations advecting the key stratigraphic surfaces within the adjacent units (see further details and references in Sec. 3.2.1). These geometric styles are a simplification of more advanced genetic concepts (see for instance [START_REF] Catuneanu | Sequence Stratigraphy: Methodology and Nomenclature[END_REF]. A better integration of concepts is the topic of recent research. In particular, due to the difficulty of determining the exact age of sediments, uncertainty about stratigraphic history may be described using several possible stratigraphic columns built from the data at hand and from various depositional concepts [START_REF] Borgomano | Stratigraphic well correlations for 3-D static modeling of carbonate reservoirs[END_REF][START_REF] Lallier | Relevance of the stochastic stratigraphic well correlation approach for the study of complex carbonate settings: Application to the Malampaya buildup (Offshore Palawan, Philippines)[END_REF][START_REF] Edwards | Uncertainty management in stratigraphic well correlation and stratigraphic architectures: A training-based method[END_REF].

Tectonic and epigenic structures

The Earth crust is affected by deformations and circulation of fluids under high temperature and pressure. These processes generate tectonic structures (which correspond to various types of rock deformation) and epigenic processes (which concern the transformations of rocks due to fluid circulations).

Tectonic events can also be formalized in terms of intrusions, fold and fault characteristics. In simple terms, intrusions are due to rocks behaving like viscous fluids (magma and evaporitic rocks at geological time scales), folds accommodate elasto-visco-plastic rock deformation, whereas faults correspond to brittle deformation localized on a slip surface. The transition between these modes of deformation is not simple to establish, and may vary depending on the considered temporal or spatial scale of interest.

Intrusions can have very complex shapes and topologies. Indeed, fluid-like rocks may follow intricate paths in the subsurface, see for instance [START_REF] Hudec | Terra infirma: Understanding salt tectonics[END_REF]. As intrusions crosscuts surrounding rocks, they are considered younger, event though they may be coeval with the deposition of some surrounding rocks (e.g., [START_REF] Giles | Halokinetic sequence stratigraphy adjacent to the El Papalote diapir, northeastern Mexico[END_REF]. The geometry of intrusions is generally determined by geometric interpolation. Relations with surrounding structures are treated by surface stitching or truncation using the same operations as with baselap stratigraphic surfaces.

Folds are typically encountered in compressive settings to accommodate shortening. In this case, folds are characterized by a higher curvature in the shortening direction than in the orthogonal direction. Folds may also exist in extensive or strike-slip regimes, where they often appear ahead of faults. In general, folded horizons should be globally developable (i.e., may be flattened with minimal distortion). This type of principle can be used to check the likelihood of a particular model (e.g. [START_REF] Gratier | Compatibility constraints on folded and faulted strata and calculation of total displacement using computational restoration (UNFOLD program)[END_REF]. In the specific case where horizon data consist of contour lines, [START_REF] Thibert | A direct method for modeling and unfolding developable surfaces and its application to the Ventura Basin (California)[END_REF] also proposed a direct method that directly creates horizon contours lines.

A fault is generally considered as a zero-thickness slip surface, that juxtaposes rocks of different nature due to the accumulation of displacement. Faults often have complex growth histories. Stress can be relaxed by frictional slip on the fault surface, by propagation of the fault tip, and by coalescence with other faults [START_REF] Marchal | Space and time propagation processes of normal faults[END_REF]. The deciphering of tectonic structures associated to faults involve several kinematic models which describe how deformation localizes around a fault in various contexts (e.g. [START_REF] Cardozo | Optimized trishear inverse modeling[END_REF]. In our Pyrenean outcrop example, a fault-propagation mechanism is invoked by [START_REF] Alonso | Forelimb deformation in some natural examples of fault-propagation folds[END_REF] to explain present-day observations: the sub-horizontal thrust faults separating the thrust sheets root into a deep detachment level below the Marboré Formation. Each fault developed towards the south, and the corresponding shortening was accommodated by a fold ahead of the fault tip. The fault has then accumulated displacement and has grown to cut the folded structure. Although it is not directly visible, a thin level of evaporite probably lubricated the decollement and filled the gaps by viscous flow between more competent rocks.

Faults are characterized by their orientation, size, roughness and displacement. Some of these parameters can be observed locally, and are also often described in statistical terms (e.g. [START_REF] Kim | The relationship between displacement and length of faults: A review[END_REF]. Relationships between geometric parameters and causality with hydromechanical processes are complex because of mechanical and hydrodynamic heterogeneities. In particular, interactions between faults are complex, leading to complex topological and geometrical configurations. Indeed, a fault may splay at the tips [START_REF] Marchal | Space and time propagation processes of normal faults[END_REF], or may grow until it branches onto an older fault. If propagation continues, the older fault is displaced by the younger one. When two conjugate branching faults intersect, they can accumulate displacement alternatively, creating X-shaped branch lines [START_REF] Ferrill | Crossing conjugate normal faults[END_REF]. Fault may also be eroded by stratigraphic surfaces.

In structural modeling, faults are generally treated as three dimensional surfaces. Determining the geometry of fault surfaces calls for interpolating between observations. The amount of fault displacement is determined by interpolation from observation data (e.g. [START_REF] Mallet | GOCAD: A Computer Aided Design Program for Geological Applications[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF], by interactive techniques [START_REF] Caumon | Constrained modifications of non-manifold b-reps[END_REF] or using specific near-field fault models [START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF][START_REF] Laurent | A parametric method to model 3D displacements around faults with volumetric vector fields[END_REF][START_REF] Godefroy | A parametric fault displacement model to introduce kinematic control into modeling faults from sparse data[END_REF]. As compared to other geological surfaces which always terminate onto other surfaces, fault surfaces can have tip lines not associated to any other surface. This can be translated into fundamental topological rules for geomodeling systems [START_REF] Caumon | Constrained modifications of non-manifold b-reps[END_REF]. This also raises practical topological and geometrical challenges to appropriately represent discontinuous properties across these internal discontinuities, see for example Chapter 2 of [START_REF] Mallet | Geomodeling[END_REF] and [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model: NON-PLANAR 3D CRACK GROWTH-PART I[END_REF]; [START_REF] Cherpeau | Stochastic simulations of fault networks in 3D structural modeling[END_REF].

Rocks may be transformed after their formation, e.g. by dissolution, precipitation and other mineral transformations in the subsurface. The features in the resulting epigenic objects are controlled by thermo-hydro-chemical processes which may span very long periods of time. Spatial prediction of these features is generally very difficult. Therefore, epigenic objects are often modeled as an overprinting of the initial structural models, using truncation as in Eq. ( 3), often combined with a latent variable as in Eq. ( 7) as in pluri-Gaussian methods [START_REF] Galli | The pros and cons of the truncated Gaussian method[END_REF][START_REF] Armstrong | Plurigaussian Simulations in Geosciences[END_REF] and some orebody delineation methods (e.g. [START_REF] Martin | Iterative refinement of implicit boundary models for improved geological feature reproduction[END_REF]. Depending on the geological context, it can be interesting to identify the geological meaning of this latent variable, e.g., the distance to fractures or other geological interfaces [START_REF] Henrion | ODSIM: An Object-Distance Simulation method for conditioning complex natural structures[END_REF][START_REF] Rongier | Simulation of 3D karst conduits with an object-distance based method integrating geological knowledge[END_REF].

Modeling workflows

In Chapter 6 of [START_REF] Perrin | Shared Earth Modeling: Knowledge Driven Solutions for Building and Managing Subsurface 3D Geological Models[END_REF], it is recommended in general to create and assemble geological surfaces by increasing geological age. This is justified by the principle that "an older geological event cannot modify a younger one". This principle is sound, and used in kinematic modeling methods such as Noddy [START_REF] Jessell | Structural geophysics: Integrated structural and geophysical modelling[END_REF][START_REF] Wellmann | Pynoddy 1.0: An experimental platform for automated 3-D kinematic and potential field modelling[END_REF]. It also led to an interesting geomodeling approach which assembles geological boundaries by topological operations completely determined by chronological considerations [START_REF] Brandel | Automatic Building of Structured Geological Models[END_REF][START_REF] Perrin | Knowledge-Driven Applications for Geological Modeling[END_REF]. However, this type of approach raises challenges because relative chronologies between objects can be difficult to establish. Moreover, all available observation points characterize the present geometry and topology of the system and not the past state. this may explain why the most common approach in commercial geomodeling packages (e.g., SKUA-Gocad, Petrel and Geomodeller) is to start by first modeling geological discontinuities such as faults and intrusions, and then to focus on stratigraphic surfaces.

Other model building orders can also be deemed convenient depending on the data at hand. For instance, one may want to start with interfaces which are most sampled or most visible from available data. Starting by modeling objects which are deemed to have most impact on the modeling outcome may also be an option [START_REF] Ringrose | Reservoir Model Design[END_REF].

For better or for worse, it is therefore not always possible to relate the modeling order of objects to geological chronology. Consider for example a salt intrusion truncating the surrounding strata: one may interpolate both geometries independently, and then cut the sediments; alternatively, one may first model the geometry of the intrusion, then interpolate the geometry of the strata to end on the unconformity (see Fig. 7). These two strategies can probably reach the same final result, but they certainly have an impact in terms of how easy model updating can be achieved and on overall model building performance and memory cost.

In today's practice, this arbitrary modeling order often implies iterations, and model updating as construction proceeds (Caumon et al., 2013b). These iterations translate how uncertainty is iteratively reduced by letting experts incorporate their geological knowledge interactively to obtain the "best" possible model. Approaches which aim at sampling uncertainties (Section 2) rely on mathematical rules such as the truncations discussed above. In both cases, interpolation between observations points is an essential component to come up with realistic geometry of the considered geological surfaces.

Numerical representations and interpolation of geological structures

This section further delves into the possible ways to numerically describe and build build structural geometry from available data and described some possible choices and meanings for the basis functions χ i (x) representing the model geometry.

Determining the geometry of geological structures from sparse data has long been a concern for geologists. Numerical methods focus either on:

• Directly building geological interfaces (explicit structural modeling methods). In this case, each geological interface is defined by interpolating between the nodes of a two-dimensional orientable graph. As surfaces are embedded in three-dimensional space, surfaces do not necessarily define a closed volume, and modeling involves projections of the data points onto the surface.

• Building a scalar function whose equipotential surfaces are the geological boundaries of interest.

The scalar function may be seen for instance as as a signed distance to the interface or as a relative geological time function. This type of method is, in general better posed, as no data projections are involved.

Explicit methods: maps and extrusion approaches

The advent of computers in the second half of the twentieth century led to so-called "automatic contouring" techniques, where the goal was to identify the elevation (or depth) z(x, y) of a particular geological surface at the geographic location (x, y) from a set of I sample points (x i , y i , z i ) i={1,...,I} [START_REF] Walters | Contouring by Machine: A User's Guide[END_REF][START_REF] Hardy | Multiquadric Equations of Topography and Other Irregular Surfaces[END_REF][START_REF] Olea | Optimal Contour Mapping Using Universal Kriging[END_REF]. These methods can only address relatively simple geological configurations as described in Fig. 2a. However, they have been much used out of consideration for their simplicity and computational efficiency. This section summarizes the main variants of mapbased approaches and introduces some mathematical interpolation principles which are common to all geomodeling methods. Early on, the various formulations of the problem were either very practical (asking how to program a contour plotting device from the data) or more theoretical (finding a semi-analytical solution by combining a series of polynomial or Fourier basis functions). The latter case can be posed as a classical interpolation problem, see for example [START_REF] Hardy | Multiquadric Equations of Topography and Other Irregular Surfaces[END_REF], where the aim is to estimate the coefficients c i of

z(x, y) = I i=1 c i χ i (x, y), (9)
where χ k (x, y) is some basis function centered at location (x k , y k ), for instance

χ i (x, y) = [(x -x i ) 2 + (y -y i ) 2 ].
In the presence of I data points (x i , y i , z i ) i=1,...,I , identifying z(x i , y i ) with z i forms a linear system whose resolution provides the coefficients c i :

[χ j (x i , y i )] • [c j ] = [c i ]. ( 10
)
As discussed by [START_REF] Hardy | Multiquadric Equations of Topography and Other Irregular Surfaces[END_REF], the use of additional data such as slope information could be addressed by a least-squares system using the analytical derivatives of χ i (x, y) and/or by adding polynomial terms p l (x, y) to (9):

z(x, y) = I i=1 c i χ i (x, y) + L l=1 d l p l (x, y). (11)
Kriging is also significantly used in contouring [START_REF] Olea | Optimal Contour Mapping Using Universal Kriging[END_REF][START_REF] Dubrule | Comparing splines and kriging[END_REF]Lemon and Jones, 2003, e.g.,). A complete coverage of kriging is out of the scope of this paper, but we refer to [START_REF] Cressie | The origins of kriging[END_REF] for a historical perspective and to [START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty, 2nd Edition[END_REF] for a thorough description. Essentially, kriging is preceded by a spatial data analysis phase, which aims at statistically characterizing the variable of interest (here the depth z(x, y)), considered as a spatial random field. In this frame, each data point z i at location (x i , y i ) is seen as the local outcome of a spatial random process. In the universal kriging formalism, the depth at location x, y is given by a linear combination of the data values plus some polynomial trend function:

z(x, y) = I i=1 λ i (x, y)z i + L l=1 d l p l (x, y). ( 12
)
The weights of this combination are obtained by minimizing the estimation variance of z (i.e., the spread between the model z and the unknown true depth z), under unbiasedness condition (i.e., stating that the average error between the model and the truth should be centered on zero). To determine the kriging weights λ i , geostatistics uses a variogram model, which captures important features about spatial variability, assuming it only depends on distances. In the end, the weights are obtained by solving a linear system fully determined by the variogram model and by the geometry of the data. When all data points are used for the kriging estimation, some algebraic manipulations produce the dual kriging system of Eq. ( 12), which has exactly the same form as Eq. ( 11), χ i corresponding to the covariance model. One practical benefit of the geostatistical approach is in the process to determine the shape of the basis functions φ, which are inferred from the available data using variogram analysis [START_REF] Dubrule | Comparing splines and kriging[END_REF].

Other interpolation methods such as inverse distance interpolation and natural neighbor interpolation have also been used in geological modeling [START_REF] Lemon | Building solid models from boreholes and user-defined crosssections[END_REF][START_REF] Ming | GSIS A 3D geological multi-body modeling system from netty cross-sections with topology[END_REF]. Like kriging, both methods are linear estimators :

z(x, y) = I i=1 λ i (x, y)z i .
Inverse distance weight each data point as using its relative inverse distance to the location to be estimated. Natural neighbors [START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF] involve Voronoi diagrams compute the weights: λ i is given by the relative area of overlap the domain nearest to (x, y) than to all the data points and the domain nearest to data point i than to any other data point. This robust weighting scheme amounts to locally adapting the support of the underlying basis functions χ i depending on data configuration. This is interesting for geoscience applications [START_REF] Watson | The natural neighbor series manuals and source codes[END_REF], but needs specific care in extrapolation [START_REF] Bobach | Natural neighbor extrapolation using ghost points[END_REF]. More generally, both methods are easier to apply than kriging, as less parameters and no statistical inference is are involved. Conversely, the geostatistical approach provides a procedure to control the spatial structure of the solution and to account for other data types, such as slopes [START_REF] Hardy | Multiquadric Equations of Topography and Other Irregular Surfaces[END_REF].

Many important details matter in these map-based approaches, from the well-posedness of the system to the management of faults, which are generally treated as discontinuities [START_REF] Pouzet | Estimation of a surface with known discontinuities for automatic contouring purposes[END_REF][START_REF] Mallet | Automatic Contouring in Presence of Discontinuities[END_REF][START_REF] Marechal | Kriging Seismic Data in Presence of Faults[END_REF]. Also, once maps have been interpolated in three dimensions, information about unconformities can be used to truncate the various surfaces and obtain the final layer geometry. For all these operations, implementation vary depending on the chosen numerical representation: rectangular grids [START_REF] Pouzet | Estimation of a surface with known discontinuities for automatic contouring purposes[END_REF][START_REF] Mallet | Automatic Contouring in Presence of Discontinuities[END_REF][START_REF] Marechal | Kriging Seismic Data in Presence of Faults[END_REF][START_REF] Swanson | A New Geological Volume Computer Modeling System for Reservoir Description[END_REF][START_REF] Fremming | 3D Geological Model Construction Using a 3D Grid[END_REF], or triangulated surfaces [START_REF] Lemon | Building solid models from boreholes and user-defined crosssections[END_REF][START_REF] Kessler | The Capture and Dissemination of Integrated 3D Geospatial Knowledge at the British Geological Survey Using GSI3D Software and Methodology[END_REF][START_REF] Zhu | Building 3D solid models of sedimentary stratigraphic systems from borehole data: An automatic method and case studies[END_REF].

Another aspect relates to the inability of map-based methods to represent and visualize structures having several elevations for the same geographic coordinates, for instance ore bodies, overturned folds, or horizons affected by reverse faults.

When faults intersect in map view but not in vertical view, this limitation can be addressed by defining a direction field tangent to the fault surfaces, see [START_REF] Fremming | 3D Geological Model Construction Using a 3D Grid[END_REF] and Mallet (2002, Chap. 8). This approach, often called pillar model, is illustrated in Fig. 9b, where pillars are shown as red lines. The horizon surfaces are then modeled by interpolating their apparent depth (or layer thickness) along these pillars, by treating faults as discontinuities. All the map-based methods are then applicable by replacing the vertical direction by the pillar direction.

As visible when comparing Fig. 9a and Fig. 9b, this method is limited to representing fault that have approximately the same dips. An alternative (or complementary) approach consists in extrapolating faults so as to completely subdivide the domain of interest into closed fault blocks, then building building strata within a larger domain and trimming by the resulting fault blocks (Hoffman and Neave, 2007a, e.g.,).

As discussed by [START_REF] Fremming | 3D Geological Model Construction Using a 3D Grid[END_REF] and Mallet (2002, Chap. 8) the above pillar-based approaches (and the map-based methods, [START_REF] Lemon | Building solid models from boreholes and user-defined crosssections[END_REF]), are often used in a hierarchical way to adapt to hierarchy of stratigraphic series. A few key horizons are first created (e.g., from reflection seismic interpretations), then intermediate horizons are interpolated between these key horizons using available data (e.g., borehole data). In the context of stratigraphic uncertainty modeling, [START_REF] Abrahamsen | Bayesian Kriging for Seismic Depth Conversion of a Multi-Layer Reservoir[END_REF] also used this principle to propose a geostatistical framework for the coherent conversion of a series of stratigraphic surfaces from the reflection seismic time domain to the depth domain. These methods are often terms extrusion methods or 2.5D methods, as they use twodimensional algorithms to create three-dimensional surfaces. As horizons are all represented as depths values along each pillars, minimal and maximal layer thickness can be enforced during interpolation. Stratigraphic unconformities can be obtained by setting consecutive horizons to the same depth value.

Another benefit of these approaches is that they relatively easily allow to create volumetric grids supporting petrophysical and physical modeling [START_REF] Swanson | A New Geological Volume Computer Modeling System for Reservoir Description[END_REF][START_REF] Johnson | Putting Geology Into Reservoir Simulations: A Three-Dimensional Modeling Approach[END_REF][START_REF] Fremming | 3D Geological Model Construction Using a 3D Grid[END_REF]Hoffman and Neave, 2007a). This explains why map-based and pillar-based models, in spite of their representational limitations, are vastly used in reservoir and near surface modeling applications. Various implementation variants exist in software such as Petrel, Gocad, RMS or Jewel Suite.

Full 3-D modeling approaches: Explicit

Instead of describing the depth as a function of the geographic position, a three-dimensional surface mesh can be used to describe geological interfaces. Surface meshes consist mainly or triangulated surfaces (also called triangular irregular networks or TINs) [START_REF] Mallet | GOCAD: A Computer Aided Design Program for Geological Applications[END_REF][START_REF] Ming | GSIS A 3D geological multi-body modeling system from netty cross-sections with topology[END_REF] or parametric surfaces such as NURBS [START_REF] Auerbach | Computer-aided geometric design of geologic surfaces and bodies[END_REF][START_REF] Fisher | Three Dimensional Solid Modeling of Geo-Objects Using Non-Uniform Rational B-Splines (NURBS)[END_REF][START_REF] De Paor | Bézier Curves and Geological Design[END_REF][START_REF] De Kemp | Interpretive Tools for 3-D Structural Geological Modeling Part I: Bézier-Based Curves, Ribbons and Grip Frames[END_REF][START_REF] Sprague | Interpretive Tools for 3-D Structural Geological Modelling Part II: Surface Design from Sparse Spatial Data[END_REF]. As these approaches have already been reviewed by [START_REF] Sides | Geological modelling of mineral deposits for prediction in mining[END_REF]; [START_REF] Caumon | Surfacebased 3D modeling of geological structures[END_REF], we only give here a brief summary or the methods so as to discuss their strengths and limitations.

Triangulated surfaces consist of a set of three-dimensional nodes {x 1 , . . . , x K } connected by triangles. In the simplest case, each node bears a linear basis function χ k equal to 1 at the vertex k and to zero at all the other vertices. The continuum of points belonging to a triangle are, therefore, described by linear combination of the three basis functions centered on the triangle vertices. This type of mesh can approximate the shape of any free-form surface provided that the mesh density is sufficient.

A first approach to create triangulated surfaces consists in associating curves defined on parallel cross-sections [START_REF] Sides | Geological modelling of mineral deposits for prediction in mining[END_REF][START_REF] Mallet | Geomodeling[END_REF][START_REF] Caumon | Surfacebased 3D modeling of geological structures[END_REF]. These methods are applicable to cylindrical or conical structures, but need many cross-sections to address domains affected by nonparallel faults, as the topology of faulted surfaces change from one cross-section to the next. Surfaces can be also created between boundary lines and borehole data, for instance by Delaunay triangulation and mesh refinement, see [START_REF] Sides | Geological modelling of mineral deposits for prediction in mining[END_REF], Mallet (2002, Chap. 3).

As all of these methods connect points and lines by linear segments, they tend to create very angular geological surfaces. Smoothing methods allow to obtain more realistic shapes while honoring subsurface data. For example, Discrete Smooth Interpolation (DSI) method [START_REF] Mallet | Discrete smooth interpolation[END_REF][START_REF] Mallet | Discrete modeling for natural objects[END_REF][START_REF] Mallet | Geomodeling[END_REF] uses numerical optimization to minimize the surface roughness while honoring available data. The principle is to start from an initial surface, then to express all pieces of information as a linear relation applied on the coordinates of the surface nodes. Finally, each of these linear relations is solved in the least-squares sense. For example, the roughness criterion can be locally expressed by setting the coordinates of each node to the average of the coordinates of its neighbors (note that this requires a specific treatment of boundaries). Another criterion is defined to locally move the surface towards data points: A particular data point x D that the surface should honor must first be projected onto the nearest surface triangle (x A , x B , x C ). Denoting (a, b, c) the barycentric coordinates of the projection in the triangle, on can write that the triangle includes the data point

x D if x D = a • x A + b • x C + c • x C . ( 13
)
More simply, a particular surface node may be locked so that is position is not altered during interpolation, which make it possible to exactly honor some data points in they are included in the surface mesh. [START_REF] Mallet | Geomodeling[END_REF] presents several other linear criteria or inequality constraints that can be used to honor other types of information, for instance local surface orientation, contacts between surfaces or thicknesses. In the end, all these linear equations are weighted and assembled in a least squares system to find the surface shape7 . The roughness is essential to ensure the well-posedness of the DSI system, as there are in general fewer data points than surface nodes. This method, together with the aforementioned surface creation and cutting algorithms, forms the core of the Gocad geomodeling platform developed in the 1990's [START_REF] Mallet | GOCAD: A Computer Aided Design Program for Geological Applications[END_REF].

These methods to create and smooth triangulated surfaces are well-defined mathematically and can be applied in an automatic way to create the various surfaces of interest independently one from another. However, when the available data cannot be projected onto a plan without changing the neighborhood (e.g., in the case of an overturned fold as in Fig. 2b or of a complex intrusion as in Fig. 2c), the creation of the initial surfaces generally involves manual operations to reach the desired result. One then needs to decompose the data to create and then merge sub-planar patches, or to locally change of the directions of attraction between the surface and the points [START_REF] Caumon | Surfacebased 3D modeling of geological structures[END_REF].

Instead of linear basis functions, NURBS and Bézier parametric surfaces essentially use higher order basis functions interpolating a control network of surface nodes to describe the geometry of surfaces. This makes it possible to represent surface geometry with a smaller number of nodes than with piecewise linear surfaces. This also allows to design complex folds and overturned surfaces by interactively moving a few control nodes in space. de [START_REF] De Kemp | Interpretive Tools for 3-D Structural Geological Modeling Part I: Bézier-Based Curves, Ribbons and Grip Frames[END_REF]; [START_REF] Sprague | Interpretive Tools for 3-D Structural Geological Modelling Part II: Surface Design from Sparse Spatial Data[END_REF] use this feature for instance to interactively design complex fold geometries by combining orientation data and graphical interactions.In general, these parametric surfaces use a tensor product of parametric curves. As a result, surfaces essentially have a rectangular shape and adaptive refinement is more difficult than with triangulated surfaces. As show by [START_REF] Zhang | 3D Parametric Modeling of Complex Geological Structures for Geotechnical Engineering of Dam Foundation Based on T-Splines: 3D parametric geological modeling based on T-splines[END_REF], T-Splines open interesting perspectives to address these challenges and manage complex interactions (Fig. 2) in the design of parametric geological surfaces. Indeed, whether triangulated or parametric, the relations between surfaces created independently are generally inconsistent with the geological principles reviewed in Section 3.1. In particular, surfaces created independently tend to cross one another or to leave gaps within the domain, which violates essential geomodel validity conditions [START_REF] Caumon | Building and editing a sealed geological model[END_REF].

To avoid this, a first option is to carefully design the surface construction methods to avoid inconsistencies in the first place. Consider for instance surfaces created by associating parallel cross sections. In this case, specific sections need to added by the modeler to delimit fault tip lines and where ore bodies terminate [START_REF] Sides | Geological modelling of mineral deposits for prediction in mining[END_REF]. The branch lines between surfaces are then linearly approximated between parallel sections. Another approach is to compute the fault surfaces and the corresponding branch lines. Horizons can then be created by picking the horizon cutoff lines on the fault surfaces used as cross-sections. Each horizon can then be created automatically by triangulating the interior of the so-called fault polygons using the available data points.

As these approaches are tedious, another option is first to create explicit surfaces covering the whole domain from all the data at hand, then corrects the problems using automatic or semi-automatic geometrical and topological operations: cutting surfaces (e.g., horizons by faults), surface trimming (e.g., faults by fault tip line of faults by other faults), extrapolation (e.g., to enforce the branching of two neighboring faults), and constrained interpolation to update surface geometry while maintaining the topological relationships [START_REF] Caumon | Surfacebased 3D modeling of geological structures[END_REF], See review of). These methods borrow a lot to computer-aided design (CAD) and may have been applied early on by editing the mesh of triangulated surfaces [START_REF] Mallet | GOCAD: A Computer Aided Design Program for Geological Applications[END_REF] and on parametric surfaces by truncating the parametric surfaces [START_REF] Gjøystdal | Computer Representation of Complex 3-D Geological Structures using a new "Sold Modeling[END_REF] 8 . In the end, mending all surfaces along their common boundary generates a boundary representation (Fig. 9c).

All the above methods imply a significant level of user interaction and validation, which require practice and skills from the modeler, and involves subjective choices. Several authors have proposed avenues to automate (or at least facilitate) geomodel building and make it both faster and more reproducible. This includes for instance automatically removing horizon data in the neighborhood of faults to reduce the risk that a horizon point attracts the surface on the other side of the fault. These practically important heuristics are implemented in commercial software, but are seldom described in scientific papers.

Notably, another strategy to increase automation uses a completely different approach based on Voronoi diagrams. These methods are rooted in computational geometry approaches, which are often used for general-purpose surface reconstruction from point sets [START_REF] Amenta | Surface reconstruction by Voronoi filtering[END_REF]. In geoscience applications, the main idea is to reconstruct all the surfaces at once by aligning the faces of Voronoi polyhedra on geological interfaces. This is achieved , either by a careful processing of data points [START_REF] Courrioux | 3D volumetric modelling of Cadomian terranes (Northern Brittany, France): an automatic method using Vorono diagrams[END_REF], Delaunay triangulation of all data points followed by incremental corrections [START_REF] Beni | Developing an adaptive topological tessellation for 3D modeling in geosciences[END_REF] or by numerical optimization [START_REF] Hale | Seismic interpretation using global image segmentation[END_REF][START_REF] Merland | Voronoi Grids Conforming to 3D Structural Features[END_REF]. These methods automatically provide a reasonable approximation of geological interfaces, but are sensitive to data density and quality. They actually represent a first class of methods which consider volume data structures to create geological surfaces.

Full 3-D modeling approaches: Implicit

Implicit modeling techniques represent geological surfaces as iso-values (or level sets, or equipotentials) of a three-dimensional scalar field s(x). This idea is quite old in itself [START_REF] Mallet | Three-dimensional graphic display of disconnected bodies[END_REF][START_REF] Houlding | 3D Geoscience Modeling Computer Techniques for Geological Characterization[END_REF][START_REF] Lajaunie | Foliation Fields and 3D Cartography in Geology: Principles of a Method Based on Potential Interpolation[END_REF], but it only became practical during the 2000's as computer memory increased.

Considering a single geological surface (ore body, horizon, fault), the idea of implicit methods is to compute the scalar field s(x) so that all data points sampling the surface have the same scalar value. In this case, it is convenient to see this scalar field as the signed distance function to the surface, being negative when x lies below (or inside) the surface and positive when x lies above (or outside) the surface, and null on the surface.

The essential principles to interpolate the scalar field s(x) between data points are the same as with map-based interpolation (Sect 3.2.1). However, stating that s(x) should be equal to zero at all data points x i is not sufficient to form a well-posed system. Indeed, a trivial solution would be that s(x) = 0 everywhere. Additional conditions are therefore needed so that the gradient of s(x) is never null. Adding extra artificial data points with non-zero values, adding some orientation measurement or forcing the gradient norm to be unitary is therefore needed to come up with a viable solution. In the case of ore body or salt modeling, data consist of internal and external points, so these can be handled by inequality constraints during least-squares minimization (Frank et al., 2007, e.g.,) or choosing appropriate positive or negative values for the points away from the interface [START_REF] Martin | Iterative refinement of implicit boundary models for improved geological feature reproduction[END_REF]. In the case where all data points are on the surface being modeled, at least one orientation point is needed [START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF]Caumon et al., 2013a), or, equivalently, the points need to be duplicated and offset before interpolation given some initial normal estimation [START_REF] Carr | Reconstruction and representation of 3D objects with radial basis functions[END_REF].

Two main numerical methods have been proposed to create individual implicit surfaces from sparse data points:

• As in Eq. 11, meshless methods estimate the scalar field s(x) as a linear combination of I basis functions χ i and monoms p l (x):

s(x) = I i=1 c i χ i (x) + L l=1 d l p l (x) (14)
This general form includes radial basis function (RBF) interpolation [START_REF] Carr | Reconstruction and representation of 3D objects with radial basis functions[END_REF][START_REF] Cowan | Rapid Geological Modelling[END_REF][START_REF] Hillier | Three-Dimensional Modelling of Geological Surfaces Using Generalized Interpolation with Radial Basis Functions[END_REF] and dual kriging [START_REF] Lajaunie | Foliation Fields and 3D Cartography in Geology: Principles of a Method Based on Potential Interpolation[END_REF][START_REF] Chilès | Modelling the Geometry of Geological Units and its Uncertainty in 3D From Structural Data: The Potential-Field Method[END_REF][START_REF] Aug | 3D geological modelling and uncertainty: The potential-field method[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF], which have, as in the twodimensional case, a similar system structure. In both cases, terms can be added to Eq. ( 14) to also use orientation data, which correspond to gradient of s(x) [START_REF] Lajaunie | Foliation Fields and 3D Cartography in Geology: Principles of a Method Based on Potential Interpolation[END_REF][START_REF] Chilès | Modelling the Geometry of Geological Units and its Uncertainty in 3D From Structural Data: The Potential-Field Method[END_REF]. In general, the basis functions χ i (x) depend on the isotropic distance ||x -x i || between the position to evaluate x and the data point x i . As this may generate blobby geometries away from the data, [START_REF] Martin | Iterative refinement of implicit boundary models for improved geological feature reproduction[END_REF] propose to use locally anisotropic distances instead or isotropic distances. In their approach the local directions of anisotropy directions are computed iteratively by solving several interpolations.

These meshless methods are very convenient but they involve solving a large and dense linear system of equations of size (I +L) 2 . Optimization methods have been proposed [START_REF] Greengard | A fast algorithm for particle simulations[END_REF], e.g., the fast multipole approach originally introduced by). A recent alternative has also been proposed by [START_REF] Renaudeau | Implicit structural modeling with local meshless functions[END_REF] using an intermediate set of points which bear local moving least squares basis functions. The overall system becomes sparse, as the basis functions have a local support, but the continuity of the solution is maintained.

• Mesh-based methods have also been proposed to compute the scalar field s(x) [START_REF] Moyen | 3D-Parameterization of the 3D Geological Space -The GeoChron Model[END_REF][START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF][START_REF] Souche | Volume Based Modeling-Automated Construction of Complex Structural Models[END_REF]Caumon et al., 2013a;[START_REF] Laurent | Iterative Thickness Regularization of Stratigraphic Layers in Discrete Implicit Modeling[END_REF]. In this case, the domain is covered by a linear tetrahedral mesh conformable to discontinuities. The interpolation method uses numerical optimization combining linear data terms as Eq. ( 13) with a roughness term stating that the scalar field should vary smoothly. As in explicit methods, weights can be used to globally or locally change the relative contribution of the data and regularization terms. This can be used for instance to change from smooth fold hinges to kink-style hinges when axial surfaces are mapped in the domain (Caumon et al., 2013a).

A significant feature in implicit methods is that several conformable surfaces may be modeled by the same scalar field. This principle was first introduced for foliation fields [START_REF] Lajaunie | Foliation Fields and 3D Cartography in Geology: Principles of a Method Based on Potential Interpolation[END_REF] and then widely adopted for modeling of stratigraphic series, as they typically consist of subparallel surfaces. In this case, the scalar field can be seen as a relative geological time function (denoted as t(x) in Sect. 3.1). This approach is illustrated in Fig. 9d. Different ways to compute this function use the aforementioned numerical methods. As compared to the reconstruction of a single surface, a practical question is to decide about appropriate values for each horizons. Indeed, a wrong choice of values may generate local minima in the interpolated solution, which imply closed bedding surfaces that cannot be generated by stratigraphic processes. Heuristic approaches based on thickness considerations have been proposed to address this problem. (Caumon et al., 2013a;[START_REF] Collon-Drouaillet | 3D geomodelling combining implicit surfaces and Voronoi-based remeshing: A case study in the Lorraine Coal Basin (France)[END_REF][START_REF] Collon | 3D modeling from outcrop data in a salt tectonic context: Example from the Inceyol minibasin, Sivas Basin, Turkey[END_REF]. Alternatively, in the case where enough data are available for the problem to be well posed, the increment of the scalar field can be estimated as differences to some reference value [START_REF] Lajaunie | Foliation Fields and 3D Cartography in Geology: Principles of a Method Based on Potential Interpolation[END_REF][START_REF] Chilès | Modelling the Geometry of Geological Units and its Uncertainty in 3D From Structural Data: The Potential-Field Method[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF][START_REF] Renaudeau | Implicit structural modeling with local meshless functions[END_REF].

In the presence of unconformities, the domain is split into several conformable sequences, then Boolean operations are applied to honor the information present in the stratigraphic column [START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF]Caumon et al., 2013a;[START_REF] Souche | Volume Based Modeling-Automated Construction of Complex Structural Models[END_REF]. The treatment of faults varies between the meshless an mesh-based techniques. Both methods need to determine fault geometry explicitly or implicitly before computing the scalar field corresponding to faulted formations. Mesh-based methods may then generate a computational mesh conformal to the fault network, where faults are treated as boundaries. Both mesh-based methods and implicit methods may also create a scalar field for each fault block, followed by truncation by the fault surfaces. In the case of meshless methods, a notable option [START_REF] Chilès | Modelling the Geometry of Geological Units and its Uncertainty in 3D From Structural Data: The Potential-Field Method[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF] is to treat faults by adding two discontinuous drift terms to Eq. ( 14), each corresponding to the average offset induced by the fault. The basis function of each term is maximal at the fault center and decreases to zero at the fault tip and orthogonally away from the fault; it is activated only on one side of the fault set to zero on the other side.

Although we have focused here on the construction of implicit surfaces from sparse data points, we should also mention several approaches to create similar implicit volumes directly from reflection seismic images [START_REF] Lomask | Flattening without picking[END_REF][START_REF] Wu | Building 3D subsurface models conforming to seismic structural and stratigraphic features[END_REF][START_REF] Wu | Horizon volumes with interpreted constraints[END_REF].

In summary implicit modeling, has opened new ways to automatically construct geological models with a high degree of complexity (Fig. 2), as typical in multiply deformed terranes (e.g. [START_REF] Jessell | Next Generation Three-Dimensional Geologic Modeling and Inversion[END_REF] or around dome structures (e.g. Wellmann et al., 2010a). This automation has led to significant progresses, as will be further discussed in Section 4.

Validating and updating geomodels

Overall, model quality is often addressed by computing the mismatch between interpretation points and the model surfaces. A complementary way is to discard some data points and use them as crossvalidation data, but this approach is only applicable in the presence of dense data sets relatively to the geological complexity. Realism is somewhat more subjective to characterize. In addition to elementary coherency rules, it may be assessed visually by model scrutiny. Quantitative methods involve thickness map computations, and volumetric considerations. Restoration techniques and deformation analysis can also be used to check if the reconstructed geometry is compatible with a likely deformation history. Last but not least, computing model forecasts allows to compare them with observations (e.g., water heads or tracer data in an aquifer, arrival times or waveforms at seismometers, etc.). In the next section, we present the various existing approaches to capture and quantify uncertainty in subsurface models.

Methods for uncertainty analysis in geological models

We described above the various sources and types of uncertainty in the typical workflow of geological model construction (Sec. 2) and the range of mathematical interpolation methods to generate single representations (Sec. 3). In the following section, we examine methods to analyze uncertainties in these model representations. In essence, most approaches can be summarized as an extension of a single, best-fitting deterministic model (Fig. 10a) to multiple possible models (or realizations) (Fig. 10b), for example multiple realizations of interfaces and faults. This approach is equally well established in the field of geophysics (e.g. [START_REF] Tarantola | Popper, Bayes and the inverse problem[END_REF] and applied (geo-)statistics (e.g. [START_REF] Caers | Modeling Uncertainty in the Earth Sciences[END_REF][START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty, 2nd Edition[END_REF][START_REF] Pyrcz | Geostatistical Reservoir Engineering[END_REF] as a way to represent and analyze uncertainties in complex non-linear systems.

Generating multiple realizations involves to propagate geological uncertainties under a particular choice about what aspect will be treated in a probabilistic way. In Sec. 4.1), we will review the various choices which can be made and provide concrete examples to method which perturb the data and/ or the structural models away from the data. This generated set of models is in itself an important outcome of the analysis of uncertainties, but it raises challenges in terms of visualization and communication (Sec. 4.2). A second challenge is to reduce uncertainty (Sec. 4.3), which may be achieved by assessing the geological likelihood of the various models (Fig. 10d), or by incorporating indirect geophysical or flow measurements (Fig. 10e). Often, these two tasks involve an additional discretization and petrophysical modeling step (Fig. 10c). In this paper, we choose not to review the petrophysical modeling aspect and refer to the geostatistical resources [START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty, 2nd Edition[END_REF]Pyrcz and Deutsch, 2014a) and petrophysics modeling methods (e.g., [START_REF] Avseth | Quantitative seismic interpretation: Applying rock physics tools to reduce interpretation risk[END_REF]. We will, however, mention some aspects of petrophysical models when reviewing examples of geophysical inversions in Sec. 4.3.2.

In Sec. 3.1, we introduced the fundamental equations describing how geological models are represented in space by combining several basis functions in Eqs. (1)(2)[START_REF]-D Modelling With Geoscientific Information Systems: The Problem[END_REF](4)(5)(6)(7)(8). However, as discussed in the remainder of Sec. 3, geological knowledge involves many possible variants and compositions of these basic equations. In this section, we propose more abstract notations, which consider that a model contains both spatial and non-spatial parameters, in particular:

x : the spatial coordinates; Φ = {ϕ 1 , . . . , ϕ K } : a set of basis functions, and m = {m 1 , . . . , m K } : the associated coefficient vector representing the spatial discretization of the geological model;

α : additional parameters corresponding to the parameters of the interpolation function, such as the exponent of radial basis functions or the parameters of the chosen covariance function, the regularization weights, etc.;

κ : primary geological information corresponding to data point locations and values described in Sec. 1.3;

τ : topological description corresponding to relationships between geological surfaces, which represent conceptual information used to build or truncate the basis functions Φ, see Eqs. (2-3);

Ω : additional conceptual information or evidence, not necessarily formalized in mathematical terms or not directly incorporated in the geological modeling process, for example information about the style of deformation, the maximum burial depth of some sample points during the geological history, etc.

For generality, we will denote all the above parameters as θ in the following and clarify to which of them we refer to, if required: In many conventional approaches, the geological uncertainties which are considered are related to the model parameters m, while all the other aspects are kept constant (e.g., [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]. However, the notation for θ clearly suggests that any uncertainty in the parameters Φ, m, α, κ, β or Ω will contribute to uncertainties in the subsequent geological model. Under this framework, defining uncertainty amounts to choosing which aspects will be randomized for uncertainty quantification (UQ).

θ = f (x; Φ, m, α, κ, β, Ω) . ( 15

Uncertainty propagation

In practice, uncertainty quantification calls for sampling some model parameters from probability distributions. This approach, widely used in applied statistics, is commonly referred to as the Monte Carlo approach to forward uncertainty quantification or error propagation (e.g. [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF]. In general, Monte Carlo simulation performs independent sampling of model parameters. However, geological knowledge often suggests that parameters are related, so a concern in geomodeling applications is to reproduce the interactions between model parameters.

In the remainder of this section, we start by describing approaches which quantify uncertainties by focusing on uncertainty in input data κ, such as interface points (Sec. 4.1.1 and Fig. 11a). Then, we discuss methods which also address uncertainty about model parameters m and interpolation parameters α (Fig. 11b) for surface-based models (Sec. 4.1.2), explicit volume models (Sec. 4.1.3), and implicit models (Sec. 4.1.4). In Sec 4.1.5, we discuss stochastic models which sample the connectivity of geological structures between observations (Fig. 11c). Finally we present ways to address conceptual uncertainties addressing the uncertainties in fundamental geological modeling rules τ and conceptual geological knowledge Ω.

Sampling data uncertainty (κ)

As discussed above (Sec. 2.2.1), uncertainties may exist about the spatial data parameters κ used to build the geological model. Following the example in Fig. 5d, we consider here first uncertainties about the position of a particular surface at depth. For example, Fig. 11a (left) shows vertical boreholes where the transition between layers is not perfectly observed. A reason may be that wireline logs do not allow for an exact positioning of the horizon, because core is missing, or because the transition itself is in fact gradual. As in all measurement errors, we will assume that we can describe the uncertainty about the exact horizon depth with a probability distribution.

A common approach to evaluate how these uncertainties in input parameters manifest themselves in the physical space of the model is to draw a sample from each distribution to obtain a realization of the input data set κ. On the basis of this input data set, we can now generate the geological model m(x; κ) using one of the interpolation methods discussed in Sec. 3. As a single model realization does not provide any information about model uncertainties, a common approach is to sample n parameter sets κ(1) . . . κ(n) and to generate a geological model for each of these parameter sets (Fig. 11a,right).

This type of method has been applied both to data location and orientation measurements (Wellmann et al., 2010b;[START_REF] Wellmann | Uncertainties have a meaning: Information entropy as a quality measure for 3-D geological models[END_REF][START_REF] Lindsay | Locating and quantifying geological uncertainty in three-dimensional models: Analysis of the Gippsland Basin, southeastern Australia[END_REF]. In these cases, implicit interpolation approaches (Sec. 3.2.3) were used, as they allow to completely automate the model construction. In principle, all fully automatic interpolation functions can be used with this type of method. In these approaches, data values or locations are sampled independently from input probability distributions reflecting positioning or measurement errors. As discussed by [START_REF] Pakyuz-Charrier | Monte Carlo simulation for uncertainty estimation on structural data in implicit 3-D geological modeling, a guide for disturbance distribution selection and parameterization[END_REF], the sampling of orientation data should consider statistical orientation models such as the Fisher-Bingham distribution, instead of an independent sampling of strike and dip. Another question with independent data sampling is whether the resulting spatial continuity of the model is realistic. Indeed, not considering spatial correlation when sampling observations may yield large layer thickness variations or unrealistic surface undulations. In some cases, parameter dependence can be determined on the basis of the measurement itself and introduced in the form of hyper-parameters to the distributions. For example, uncertainties in time-depth-conversion will affect multiple seismic picks. The distributions for all of these points can then be described with respect to uncertainties in the time-depth conversion function. Similar examples are interface points determined from gravity inversion, or even multiple points on a single interface, included on the basis of expert judgment. In addition, Bayesian approaches have been developed to consider additional geological criteria (Fig. 10d) in order to obtain parameter correlations (see de la Varga and Wellmann , 2016, and Sec. 4.3).

Geometric model perturbation: surface-based methods (explicit surface parameters m and interpolation parameters α)

Instead of changing point positions and regenerating interpolations, geometric model perturbation essentially aims at adding geometric noise to an existing model (Fig. 11b). As geological structures are generally smooth, the noise is commonly correlated in space in the form of a covariance function (or variogram model). This type of method was originally introduced to assess petroleum accumu-lation uncertainties in oil and gas reservoirs [START_REF] Delfiner | Conditional Simulations : A New Monte Carlo Approach To Probabilistic Evaluation Of Hydrocarbon In Place[END_REF]. On a map-based gridded horizon representation (Sec. 3.2.1), the idea is simply to generate spatially correlated random fields for possible perturbations at each grid node. Efficient algorithms such as the turning bands or sequential Gaussian simulation may be used to generate these random fields while honoring discrete observations, see for example Chap. 8 of [START_REF] Goovaerts | Geostatistics for natural resources evaluation[END_REF] or Chap. 7 of [START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty, 2nd Edition[END_REF]. This simulation approach takes two sequential parameters as input: the unconditional probability distribution that the simulated random field should honor, and a model of spatial variability. Both parameters are generally assumed to be the same everywhere in space, and therefore a perturbation is generally simulated, instead of the horizon depth itself. Indeed, this separation allows the consideration of regional inclination and folding of the strata during interpolation, and then to choose appropriate statistical models to describe the residual between the interpolated model and the unknown truth. For example, choosing a Gaussian probability distribution of zero mean for the perturbation will lead to a set of model realizations which symmetrically deviate from the interpolated reference model. The spatial frequency of the perturbation is typically described by a variogram (or spatial covariance) model. Both the distribution and the variogram model have a strong impact on the produced realizations and should be calibrated from available data and prior geological knowledge. It is also possible to account for uncertainty in the variogram range or in the dispersion of the distribution. In all cases, the estimation variance produced by kriging provides local estimation errors, which depend on the chosen variogram model and on the data layout. This information may be used directly in its discretized form (Fig. 10c) for some applications, but most works prefer generating realizations which produce possible models.

An early but very complete example of these approaches is available in [START_REF] Abrahamsen | Stochastic Models for Seismic Depth Conversion of Geological Horizons[END_REF] and [START_REF] Abrahamsen | Bayesian Kriging for Seismic Depth Conversion of a Multi-Layer Reservoir[END_REF]: they combined uncertainties in seismic time-to-depth conversion by kriging multiple sub-parallel geological horizons in a "layer-cake" model using a Bayesian kriging approach [START_REF] Omre | Bayesian kriging-Merging observations and qualified guesses in kriging[END_REF]. The work of these authors denotes a step-change in the consideration of uncertainties, as they also explicitly consider the correlation of multiple layers through the seismic input data. They extended the work subsequently further to hydrocarbon-in-place estimates [START_REF] Abrahamsen | An Integrated Approach to Prediction of Hydrocarbon in Place and Recoverable Reserve With Uncertainty Measures[END_REF], gross volume estimates [START_REF] Abrahamsen | Estimation of Gross Rock Volume of Filled Geological Structures With Uncertainty Measures[END_REF] and the management of inequality constraints [START_REF] Abrahamsen | Surface prediction using rejection sampling to handle non-linear constraints[END_REF].

Another example of this approach is the probability-field method (e.g. [START_REF] Lecour | Modelling of stochastic faults and fault networks in a structural uncertainty study[END_REF][START_REF] Thore | Structural uncertainties: Determination, management, and applications[END_REF], which can be applied to geological interfaces and faults as surface elements, even in complex 3-D settings. A first idea in this case is to consider perturbation amplitude not along the vertical direction but orthogonally to each surface. The knowledge of the contacts between surfaces is used to define in which order the various structural surfaces are perturbed. For example, in the simple case of a horizon cut by a fault, the idea is to first perturb the fault geometry, then the horizon geometry to maintain the maintain the contact between both surfaces [START_REF] Lecour | Modelling of stochastic faults and fault networks in a structural uncertainty study[END_REF]. The second idea is to pre-compute the range of uncertainty around each surface by carefully propagating seismic imaging uncertainties and interpretation uncertainties. The uncertainty envelopes defined around each surface implicitly define a cumulative distribution function (cdf) along the surface normal (very similar to the representation in Fig. 11b. This cdf is then sampled by generating a spatially correlated random field between 0 and 1 using the probability field method [START_REF] Srivastava | Reservoir Characterization With Probability Field Simulation[END_REF]. As the envelope already reflects data uncertainties, this method provides an alternative to conventional conditional simulation procedures with the advantage that fast spectral methods can be implemented to simulate the probability fields (e.g. [START_REF] Pyrcz | Uncertainty in reservoir modeling[END_REF].

These methods have been widely applied (e.g. [START_REF] Samson | Quantifying the Impact of Structural Uncertainties on Gross-Rock Volume Estimates[END_REF][START_REF] Corre | Integrated Uncertainty Assessment For Project Evaluation and Risk Analysis[END_REF][START_REF] Charles | Experience with the Quantification of Subsurface Uncertainties[END_REF][START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF][START_REF] Irving | An uncertainty modelling workflow for structurally compartmentalized reservoirs[END_REF] to generate models sampling uncertainties both around horizons and faults. [START_REF] Lecour | Modelling of stochastic faults and fault networks in a structural uncertainty study[END_REF][START_REF] Thore | Structural uncertainties: Determination, management, and applications[END_REF][START_REF] Hollund | Havana-a fault modeling tool[END_REF][START_REF] Holden | Stochastic Structural Modeling[END_REF][START_REF] Rivenaes | A 3D stochastic model integrating depth, fault and property uncertainty for planning robust wells, Njord Field, offshore Norway[END_REF]. Some authors mention, thought, that the method can create geologically unrealistic surfaces (e.g. [START_REF] Srivastava | Probability Field Simulation: A Retrospective[END_REF][START_REF] Pyrcz | Geostatistical Reservoir Engineering[END_REF]). In the case of faults, recent improvements have been proposed by [START_REF] Røe | Flexible Simulation of Faults[END_REF][START_REF] Røe | An Uncertainty Model for Fault Shape and Location[END_REF], treating faults as surfaces in a rotated coordinate system and calculating horizon displacements using 3-D vector fields (see [START_REF] Hoffman | The fused fault block approach to fault network modelling[END_REF][START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF].

These methods are, in general, limited to relatively simple structural settings not affected by faults (Fig. 2a), or, in faulted domains (Fig. 2b), to cases where geometric uncertainties are small relative to fault spacing. Indeed, the sequential perturbation of many interfaces is very difficult to implement.

Geometric model perturbation: space warping methods

In geological contexts where many surfaces are present, the surface perturbation strategy relies on the careful consideration of the perturbation order to maintain model consistency. In the case of small perturbations, a volumetric approach can be applied instead. The principle is the same as for surface perturbation, but it uses a three-dimensional vector field instead of a scalar field applied in a prescribed direction. The main advantage is to consider all geomodel elements at once, treated as embedded in an elastic material. This is proposed for instance by [START_REF] Caumon | Elements for Stochastic Structural Perturbation of Stratigraphic Models[END_REF], who combine several two-dimensional perturbation fields generated independently on a set of faults to generate a global three-dimensional scalar field. As the input vector fields may cross, a divergencefree term is added to the interpolation to ensure that the perturbation is volumetrically consistent. This methodology was used for instance by [START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF] to create a population of structural models by direct deformation of a hexahedral reservoir grid.

In the context of local model updating, [START_REF] Tertois | Editing faults within tetrahedral volume models in real time[END_REF] propose a method to bound the variation of volume throughout deformation in the context of interactive model editing by experts. This approach was used by [START_REF] Caumon | Elements for Stochastic Structural Perturbation of Stratigraphic Models[END_REF] and [START_REF] Mallet | Solid Earth Modeling and Geometric Uncertainties[END_REF] to globally deform a model to reflect some small geometric fault uncertainties. More recently, [START_REF] Laurent | Interactive editing of 3D geological structures and tectonic history sketching via a rigid element method[END_REF] proposed to adapt an "as rigid as possible" method from computer graphics to also interpolate a 3-D displacement field between a set of fixed points. This type of method could certainly be extended to deform a reference model volumetrically. However, as large model distortions may occur in the case of large perturbations, these methods are only applicable for relatively mild uncertainties.

Geometric model perturbation: implicit methods (implicit m, κ and β)

In domains of high complexity (Fig. 2b andc), extending implicit interpolation methods (Sec. 3.2) has significantly helped accounting for uncertainties. This includes methods which perturb existing data or simulate synthetic data points (Wellmann et al., 2010a;[START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF][START_REF] Lindsay | Locating and quantifying geological uncertainty in three-dimensional models: Analysis of the Gippsland Basin, southeastern Australia[END_REF]Lindsay et al., , 2013b) ) and also methods which perturb the implicit scalar field away from the data [START_REF] Caumon | Elements for Stochastic Structural Perturbation of Stratigraphic Models[END_REF][START_REF] Caumon | Towards stochastic time-varying geological modeling[END_REF][START_REF] Mallet | Solid Earth Modeling and Geometric Uncertainties[END_REF][START_REF] Cherpeau | Stochastic structural modelling in sparse data situations[END_REF].

In the potential field method (see Sec. 3.2.3 and Lajaunie et al. (1997); [START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF]), the use of kriging algorithm provides, like in map-based approaches, a direct estimate of interpolation uncertainty through estimation variance. This approach has been implemented by [START_REF] Chilès | Modelling the Geometry of Geological Units and its Uncertainty in 3D From Structural Data: The Potential-Field Method[END_REF] and [START_REF] Aug | 3D geological modelling and uncertainty: The potential-field method[END_REF], but the calibration of the variogram model relies on the availability of many orientation measurements. However, the possibility to obtain a direct measure of uncertainty in the interpolation step is a promising aspect for future research.

The direct link to the initial geological input data also allows a direct consideration of effects of data density. This possibility has been used by [START_REF] Putz | Three-dimensional model and late stage warping of the Plattengneis Shear Zone in the Eastern Alps[END_REF] to evaluate the effect of reduced data density on a geological model, generated mostly from field observations in the Eastern Alps. The potential of this approach to consider uncertainties related to data positions in poly-deformed terranes has also been recognized by [START_REF] Maxelon | A workflow to facilitate three-dimensional geometrical modelling of complex poly-deformed geological units[END_REF], especially as many interface points have to be placed to constrain subsurface structures, where only very limited or no direct observation is available.

An important reason for the direct consideration of uncertainties in interface points is the possibility to treat high uncertainties related to the common requirement to include virtual data points to regulate model geometry, especially at depth when no abundant information is available. This is especially the case when creating models in regions where no high-quality seismic or borehole data is available, for example in regional-scale geological modeling [START_REF] Jessell | Next Generation Three-Dimensional Geologic Modeling and Inversion[END_REF]. These additional points are then commonly placed on the basis of expert knowledge, and therefore highly subjective and potentially also highly uncertain.

In mesh-based implicit methods, the kriging formalism is generally replaced by a discrete formulation of thin-plate spline energy or some other roughness criterion, see Sec. 3.2.3 and (Moyen et al., 2004;[START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF][START_REF] Souche | Volume Based Modeling-Automated Construction of Complex Structural Models[END_REF]Caumon et al., 2013a;[START_REF] Laurent | Iterative Thickness Regularization of Stratigraphic Layers in Discrete Implicit Modeling[END_REF]. Independently of how exactly a mesh-based implicit model is computed, a simple and direct way to perturb the model away from the observations is to add a spatially correlated random field (x) to the reference implicit scalar field s(x) [START_REF] Caumon | Elements for Stochastic Structural Perturbation of Stratigraphic Models[END_REF][START_REF] Caumon | Towards stochastic time-varying geological modeling[END_REF][START_REF] Mallet | Solid Earth Modeling and Geometric Uncertainties[END_REF][START_REF] Cherpeau | Stochastic simulations of fault networks in 3D structural modeling[END_REF][START_REF] Cherpeau | Stochastic structural modelling in sparse data situations[END_REF][START_REF] Aydin | Quantifying structural uncertainty on fault networks using a marked point process within a Bayesian framework[END_REF]. In this approach (Fig. 12), the random field values should be constrained to be zero at (certain) data locations. Away from the data, the random field must vary smoothly to avoid introducing kinks in the perturbed surface geometry. Scaling the random field (x) by 1/||∇s(x)|| approximates the horizon shift induced by the perturbation. In the presence of faults, [START_REF] Caumon | Elements for Stochastic Structural Perturbation of Stratigraphic Models[END_REF] suggest that computing the perturbation field in depositional space [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF] preserves the fault slip. [START_REF] Cherpeau | Stochastic structural modelling in sparse data situations[END_REF] used this idea to produce a relatively large continuous perturbation across faults (preserving the fault slip) with a relatively smaller discontinuous perturbation (changing the fault slip). Further approaches to geometrical perturbations of implicit stratigraphic models are also discussed by Mallet (2014, Chap. 9).

In all cases, implicit methods do not explicitly allow for controlling the connectivity of the level set surfaces. A consequence is that adding virtual data points or adding a random spatially correlated noise (x) may change the topology of the surfaces by creating closed iso-surfaces. This possibility may be a strength for some complex geological bodies such as ore bodies, as their connectivity can be uncertain; conversely, it may generate inconsistent models in the case of stratigraphic formations, as depositional processes can never generate closed surfaces.

Interesting extensions are also the recent developments of implicit methods for modeling of multistage fold geometries and overprinting deformation by [START_REF] Laurent | Implicit modeling of folds and overprinting deformation[END_REF], which allows an integration of these structural elements into uncertainty quantification and inverse frameworks [START_REF] Grose | Structural data constraints for implicit modeling of folds[END_REF][START_REF] Grose | Inversion of Structural Geology Data for Fold Geometry[END_REF]. In essence, this approach infers the parameters of Fourier series describing fold limb rotation and fold axis rotation angles from structural observations. The obtained posterior distribution are then used to estimate the uncertainty in poly-phased implicit fault models.

4.1.5. Topological uncertainties and data association (m, Φ, κ, β and τ )

The above model perturbation methods allow to sample uncertainty around a particular structural interpretation. Depending on the geological context, these perturbations may imply a few topological changes due to truncation by faults, unconformities or intrusions (see Sect. 3.1.1 and Fig 7). Indeed, the geometric perturbation followed by the truncation of a particular geological interface may change the number of connected components or the number of holes in this interface. In the example of Fig. 12, this could happen for instance when an implicit horizon goes above or below the branch line formed by X-shaped, λ-shaped or Y-shaped contacts between fault surfaces.

However, these changes remain relatively limited, in the sense that they do not significantly change the interpretation. More drastic topological changes imply considering variable numbers of geological interfaces, various ways to associate observations, and possibly also various relative ages (and truncation rules) for these interfaces.

In the context of fault uncertainty management in petroleum reservoirs, [START_REF] Munthe | Modelling Sub-seismic Fault Patterns using a Marked Point Process[END_REF]; [START_REF] Hollund | Havana-a fault modeling tool[END_REF]; [START_REF] Holden | Stochastic Structural Modeling[END_REF] pioneered by formalizing the main components of uncertain fault models. Their stochastic fault model is based on fault objects which are simulated using a marked point process. Prior information is translated using notions of fault families characterized by size, orientation and displacement distributions. Each fault is simulated by a stochastic point process which may include attraction or repulsion with previously simulated faults or fault families. Fault may also be truncated in the simulation process based on geometrical considerations. Last, each fault generates a discontinuous displacement field shifting the layers on either side of the fault. This approach has been mainly applied on extrusion-based 2.5D representations of subsurface models, see for example [START_REF] Rivenaes | A 3D stochastic model integrating depth, fault and property uncertainty for planning robust wells, Njord Field, offshore Norway[END_REF].

As the representativity of these models is limited, more recent work has focused on extending these methods to more general three-dimensional cases [START_REF] Maerten | Three-dimensional geomechanical modeling for constraint of subseismic fault simulation[END_REF][START_REF] Røe | Flexible Simulation of Faults[END_REF][START_REF] Cherpeau | Stochastic simulations of fault networks in 3D structural modeling[END_REF][START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF][START_REF] Cherpeau | Stochastic structural modelling in sparse data situations[END_REF][START_REF] Aydin | Quantifying structural uncertainty on fault networks using a marked point process within a Bayesian framework[END_REF]. Another recent area of research concerns the conditioning to observation data. Indeed, stochastic fault models may be seen as a way to associate spatial data (e.g., points where faults have been observed) consistently with prior knowledge (e.g., ideas about the fault family orientations, sizes). However, generating fault using a spatial point process until all observations are matched can be very inefficient. Therefore, [START_REF] Cherpeau | Stochastic simulations of fault networks in 3D structural modeling[END_REF]; [START_REF] Cherpeau | Stochastic structural modelling in sparse data situations[END_REF] have proposed a method where fault evidence are associated sequentially by sampling fault centers near the data and then trying to associate other fault data consistently with prior information. This problem has recently been formalized using a graph which can be sampled to generate stochastic fault networks [START_REF] Godefroy | Modlisation cinmatique et stochastique des failles partir de donnes parses pour lanalyse des incertitudes structurales[END_REF].

Similar methods have also been applied in stratigraphic settings to simulate the possible location of stratigraphic unconformities from sparse borehole data [START_REF] Lallier | Relevance of the stochastic stratigraphic well correlation approach for the study of complex carbonate settings: Application to the Malampaya buildup (Offshore Palawan, Philippines)[END_REF][START_REF] Lallier | Uncertainty assessment in the stratigraphic well correlation of a carbonate ramp: Method and application to the Beausset Basin, SE France[END_REF][START_REF] Wu | Simultaneous multiple well-seismic ties using flattened synthetic and real seismograms[END_REF][START_REF] Edwards | Uncertainty management in stratigraphic well correlation and stratigraphic architectures: A training-based method[END_REF]. As the location and number of gaps in the stratigraphic record vary, this amounts to changing the number of units and relationships between units in the stratigraphic column.

Alternative modeling approaches

In addition to the approaches mentioned above, several variants for stochastic model generation exist for alternative modeling approaches. Several authors used simple parametric models of lines and planes to evaluate uncertainties of faults, for example the propagation at depth related to uncertain observations and measurements at the surface using linear projection functions (e.g. [START_REF] Bistacchi | 3D fold and fault reconstruction with an uncertainty model: An example from an Alpine tunnel case study[END_REF], or to link fault and shear zone observations at the surface and in an underground lab [START_REF] Schneeberger | Methods and uncertainty estimations of 3-D structural modelling in crystalline rocks: a case study[END_REF]. In our notation, these approaches can be modeled by simple linear basis functions for linear regression, where position and gradient can be adjusted through the coefficient vector m.

Modeling approaches based on kinematic modeling concepts are difficult to describe with the notation in Eq.15, as they are based on modeling the interaction of multiple geological events with kinematic equations (e.g. [START_REF] Jessell | Noddy" -An interactive Map creation Package[END_REF][START_REF] Jessell | Next Generation Three-Dimensional Geologic Modeling and Inversion[END_REF]. However, this approach has also been used in a stochastic framework to estimate the effect of uncertain kinematic parameters (e.g. fault geometry, folding parameters) and timing of geological events (e.g. the order of faults) using the automatic modeling implementation pynoddy [START_REF] Jessell | Next Generation Three-Dimensional Geologic Modeling and Inversion[END_REF][START_REF] Wellmann | Pynoddy 1.0: An experimental platform for automated 3-D kinematic and potential field modelling[END_REF].

Dealing with multiple models: visualization and communication of uncertainties

There are generally two reasons for evaluating uncertainties in a study (e.g. [START_REF] Beven | Facets of uncertainty: epistemic uncertainty, non-stationarity, likelihood, hypothesis testing, and communication[END_REF]: (1) to answer a scientific question, improve understanding and evaluate different hypotheses, and (2) as a guidance for a subsequent decision making processes. In the second case, the communication of uncertainties becomes an elementary aspect. As these visualizations from an important basis for decisions (often by non-experts, and stakeholders with a varying background and expertise in the topic), the uncertainties associated with these models should be communicated. The analysis and representation of uncertainties in a spatial context is in its essence different to common approaches in uncertainty quantification approaches, where often the statistical distribution of one or several parameters is already sufficient as a measure of uncertainty (e.g. [START_REF] Gelman | Bayesian data analysis[END_REF][START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF]. In the spatial context, the main interest is often to determine the uncertainty of a model outcome (geological unit, porosity, etc.) at all locations in the investigated domain. The following description is therefore specifically relevant to these cases where the representation of the geological model is actually already a final aim, and it is not directly processed further to, for example, a process simulation, or an analysis which would lead to a dimension reduction (e.g. calculations of volume of a resource, geothermal heat in place, the outcome at a previous specified single location, etc.). This is to say, in all cases where the subsequent model use is not directly obvious or where the spatial structure itself will form part of a decision progress (e.g. in well planning, planning for next exploration steps or experimental design for additional measurements).

The visualization and communication of uncertainties is an active field of research (see, for example, the excellent overview by [START_REF] Spiegelhalter | Visualizing Uncertainty About the Future[END_REF] and, in the spatial context, it has long been driven by approaches to visualizations in a geographic context of maps [START_REF] Maceachren | Visualizing Uncertain Information[END_REF]Goodchild et al., 1994a,b;[START_REF] Leung | Visualization of fuzzy scenes and probability fields[END_REF][START_REF] Pang | Approaches to uncertainty visualization[END_REF][START_REF] Maceachren | Visualizing Geospatial Information Uncertainty: What We Know and What We Need to Know[END_REF]. Applications of these concepts have since been successfully applied in various geoscientific applications, notably in the analysis of geohazards (e.g. [START_REF] Pang | Visualizing Uncertainty in Natural Hazards[END_REF][START_REF] Kunz | Visualization of uncertainty in natural hazards assessments using an interactive cartographic information system[END_REF], with specific applications to flood prediction (e.g. [START_REF] Bruen | Visualizing flood forecasting uncertainty: some current European EPS platforms-COST731 working group 3[END_REF]Beven et al., 2014;[START_REF] Lim | Visualisation and evaluation of flood uncertainties based on ensemble modelling[END_REF][START_REF] Seipel | Color map design for visualization in flood risk assessment[END_REF], tsunami risk mapping (e.g. [START_REF] Goda | Uncertainty modeling and visualization for tsunami hazard and risk mapping: a case study for the 2011 Tohoku earthquake[END_REF], hurricane prediction (e.g. [START_REF] Cox | Visualizing uncertainty in predicted hurricane tracks[END_REF], seismic risks (e.g. [START_REF] Bostrom | Visualizing Seismic Risk and Uncertainty[END_REF], as well as in different analyzes of remote sensing scenes [START_REF] Van Der Wel | Visual exploration of uncertainty in remote-sensing classification[END_REF]; [START_REF] Comber | Spatial analysis of remote sensing image classification accuracy[END_REF] and in hydrogeological models (e.g. [START_REF] Benke | Visualisation of spatial uncertainty in hydrological modelling[END_REF]. But even though many approaches to visualize and communicate uncertainties have been developed in these related fields, most have been developed for the 2-D context of maps and sections and they are partly not directly applicable to the full 3-D setting of geological models considered here.

We can generally distinguish two approaches to visualize and communicate uncertainties from ensembles of generated models:

1. the visualization of multiple model realizations, and methods for morphing between different model realizations in a movie-form, and 2. the quantitative analysis and subsequent visualization of uncertainties using summary measures (e.g. variance, entropy, etc.).

These two aspects are described in the following.

Simultaneous representation of multiple models

An intuitive and direct visualization of spatial uncertainties is the representation of multiple possible model outcomes. This method has long been a standard method to represent uncertainties in geostatistical studies (e.g. [START_REF] Pyrcz | Uncertainty in reservoir modeling[END_REF][START_REF] Mariethoz | Multiple-point geostatistics : stochastic modeling with training images[END_REF] and to visualize multiple outcomes from geophysical inversions (e.g. [START_REF] Mosegaard | Monte Carlo sampling of solutions to inverse problems[END_REF][START_REF] Boschetti | Interactive inversion in geosciences[END_REF]. Also in the field of geological modeling, this type of representation is widely used (e.g. [START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF]Wellmann et al., 2010a;[START_REF] Lindsay | Locating and quantifying geological uncertainty in three-dimensional models: Analysis of the Gippsland Basin, southeastern Australia[END_REF]. A typical example is the presentation of multiple reservoir models, for example in [START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF], presented in Fig. 13a. These representations provide a good intuition for a possible variability in the generated models, but it quickly becomes difficult to interpret where exactly changes in the model occur and the representation is, for practical reasons, restricted to a limited number of figures.

A representation of multiple model outcomes in one figure provides a more direct representation of this variability. An example for multiple models in one section is given in Fig. 13b from [START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF], where the spatial variability for three different layers is clearly visible. This type of representation has also been attempted in 3-D, an example by [START_REF] Mallet | Solid Earth Modeling and Geometric Uncertainties[END_REF] is presented in Fig. 13c. However, it is obvious that this type of representation can quickly get confusing and difficult to interpret for increasingly complex 3-D structural settings.

A related possibility is the visualization of multiple models on a computer screen, either as a defined succession in the form of a movie, or in an interactive representation, where the user can iterate through a number of realizations (e.g. [START_REF] Srivastava | The Interactive Visualization of Spatial Uncertainty[END_REF][START_REF] Kunz | Visualization of uncertainty in natural hazards assessments using an interactive cartographic information system[END_REF]. This type of representation resolves some of the problems of a presentation in a fixed image, especially if it is possible to also rotate views in full 3-D, but it is more difficult to share or to integrate into a publication or report.

Approaches to quantify uncertainties based on a set of model realizations

Instead of the simultaneous representation, an additional approach is to use the set of generated models first to determine a spatial estimate of uncertainty, and then to determine a suitable method to visually present this uncertainty. Approaches to quantify uncertainty can be separated into methods directly using continuous measures (layer thickness at a location, depth to interface), and those using categorical values (geological unit, lithology class, etc.).

For continuous variables, a wide range of conventional statistical measures can directly be applied and visualized, for example showing mean values and standard deviations [START_REF] Potter | Visualizing Summary Statistics and Uncertainty[END_REF]). An 2010a) is presented in Fig. 14a, highlighting the uncertainty in layer depth related to the uncertain position of a fault. The requirement for realizations to follow a normal distribution for the analysis of values of mean and standard deviations can be circumvented through the use of percentiles and min/max surfaces (e.g. [START_REF] Potter | Visualization of Uncertainty without a Mean[END_REF][START_REF] Røe | An Uncertainty Model for Fault Shape and Location[END_REF]. A similar visualization can be attempted using volume rendering techniques, where uncertain areas are more transparent. An example is shown in Fig. 14b. Here, uncertainties of an interface position are shown and it is visible that these uncertainties increase when deviating from well locations. These types of analyzes are straight-forward to apply and implemented in many software packages, but their application is limited to 2.5-D structural settings of relatively low complexity (see Fig. 1), and can not easily be generalized to full 3-D settings (see Wellmann et al., 2010a, for an example).

An alternative to the consideration of continuous variables is the direct use of the categorical value of a geological unit or lithology in 3-D space. The following description of this type of spatial quantification of uncertainties follows Wellmann et al. (2010a), [START_REF] Wellmann | Information Theory for Correlation Analysis and Estimation of Uncertainty Reduction in Maps and Models[END_REF][START_REF] Wellmann | Uncertainties have a meaning: Information entropy as a quality measure for 3-D geological models[END_REF][START_REF] Wellmann | Information Theory for Correlation Analysis and Estimation of Uncertainty Reduction in Maps and Models[END_REF].

We consider here the subdivision of the model domain into defined subspaces, where m(x) represents a discrete geological class C (e.g., a lithology type or stratigraphic unit) 9 We can identify the class membership at each location x and, on this basis, define an indicator function for each model i in the ensemble of n generated models:

1 C (x, i) = 1 if x ∈ C in model i 0 if x / ∈ C in model i (16)
This indicator function is then used to estimate a probability P C (x) for each location in the model domain to belong to a class C on the basis of an ensemble of n generated models:

P C (x) = 1 n i 1 C (x, i) (17)
In this way, we obtain a multinominal probability field (see also Goodchild et al., 1994a), which can be used to visualize the probability for each separate possible class C (i.e. each geological unit) in 3-D representations (Fig. 14c). The variance σ C (x) of the categorical variable may also be used to map spatial uncertainty:

σ 2 C (x) = P C (x)(1 -P C (x)
). This analysis and representation is directly suitable when we are interested in a single possible outcome (e.g. a reservoir layer, or a mineable resource), but it is less useful as an overall measure of model uncertainty if more than two geological formations or lithology classes are considered, as we would then need multiple figures to present probabilities for each class.

For these cases, a suitable summary measure for a representation of uncertainties in the full 3-D space is required. Motivated by previous applications in on the definition of spatial entropy [START_REF] Batty | Spatial Entropy[END_REF] and the representation of uncertainties in maps [START_REF] Leung | Visualization of fuzzy scenes and probability fields[END_REF]Goodchild et al., 1994a), the use of cell entropy has been shown to be a useful method to quantify and visualize uncertainties about the full 3-D structural model [START_REF] Wellmann | Uncertainties have a meaning: Information entropy as a quality measure for 3-D geological models[END_REF]. It has similarly been applied to wireline log classifications [START_REF] Grana | Quantitative log interpretation and uncertainty propagation of petrophysical properties and facies classification from rock-physics modeling and formation evaluation analysis[END_REF] and cross-sections (e.g. [START_REF] Elfeki | Modelling Subsurface Heterogeneity by Coupled Markov Chains: Directional Dependency, Walther's Law and Entropy[END_REF]. Specifically, we use the Shannon entropy [START_REF] Shannon | A mathematical theory of communication[END_REF] to calculate an entropy value at each location x on the basis of the class probabilities P C :

H(x) = - C P C (x) log 2 P C (x) (18)
The advantage of this method is that we now obtain a scalar measure of uncertainty at each location x in space, and this aspect allows for a direct representation of model uncertainties in space.

Even though entropy values can be calculated in a continuous field, it is common practice to use a defined space-filling grid structure and to calculate entropy values for each cell for the purpose of visualization. A simple example for the use of cell entropy for 3-D geological models is presented in Fig. 14d (after [START_REF] Wellmann | Uncertainties have a meaning: Information entropy as a quality measure for 3-D geological models[END_REF]. The difference to the class probability representation (Fig. 14c) is here clearly visible: whereas in the latter case, uncertainties about a specific class outcome are shown, entropy provides a cumulative measure of uncertainty. Entropy as a measure of uncertainty has been applied successfully in geological modeling studies, both for the representation of model uncertainties (e.g. [START_REF] Bianchi | Integrating deterministic lithostratigraphic models in stochastic realizations of subsurface heterogeneity. Impact on predictions of lithology, hydraulic heads and groundwater fluxes[END_REF][START_REF] Schweizer | Uncertainty assessment in 3-D geological models of increasing complexity[END_REF], as well as the comparison of uncertainty reduction in the context of Bayesian inversion (e.g. [START_REF] Wang | A Segmentation Approach for Stochastic Geological Modeling Using Hidden Markov Random Fields[END_REF][START_REF] Wellmann | Uncertainty estimation for a geological model of the Sandstone greenstone belt, Western Australia -insights from integrated geological and geophysical inversion in a Bayesian inference framework[END_REF][START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF].

Apart from the use as a measure to visualize uncertainties, it should also be noted that the value itself has a clear quantitative meaning: the upper bound is defined by the number of classes k

H max = log 2 k (19)
and this value is obtained when all class outcomes are equally, corresponding to the case of maximum uncertainty. For more details on the underlying theoretical concepts, see [START_REF] Cover | Elements of Information Theory[END_REF].

Additional aspects about the quantitative spatial interpretation, the calculation of related fuzziness measures to quantify uncertainties of specific classes, as well as model average values are described in [START_REF] Wellmann | Uncertainties have a meaning: Information entropy as a quality measure for 3-D geological models[END_REF]. Furthermore, measures from information theory can be applied to determine the spatial correlation of uncertainties (see [START_REF] Wellmann | Information Theory for Correlation Analysis and Estimation of Uncertainty Reduction in Maps and Models[END_REF].

Another concept to quantify uncertainties in geological models on the basis of class membership is the measure of stratigraphic possibility, as the number of possible (lithology) classes L(x), and variability V (x) at each location, introduced by [START_REF] Lindsay | Locating and quantifying geological uncertainty in three-dimensional models: Analysis of the Gippsland Basin, southeastern Australia[END_REF]. In the notation used above, this measure can be calculated as:

V (x) = 1 -P Ĉ (x) ( 20
)
where Ĉ is the class with the highest estimated probability at position x. Both measures should be used for a representation of model uncertainty. For a more detailed description and examples, see [START_REF] Lindsay | Locating and quantifying geological uncertainty in three-dimensional models: Analysis of the Gippsland Basin, southeastern Australia[END_REF] and [START_REF] Lindsay | Inversion and Geodiversity: Searching Model Space for the Answers[END_REF].

Communication of uncertainties

With the methods described above, we obtain a scalar measure of uncertainty for the model domain, which can be visually represented in 2-D and 3-D views. However, the general communication of these uncertainties is, in itself, a challenging aspect. The question of uncertainty visualization is completely entwined with the question of how uncertainties are perceived by potential stakeholdersand therefore on psychological aspects (e.g. [START_REF] Spiegelhalter | Visualizing Uncertainty About the Future[END_REF]. Certainly, the choice for the best method to represent uncertainties depends on the specific type of data, the context, and the target audience (see [START_REF] Visschers | Probability Information in Risk Communication: A Review of the Research Literature[END_REF]. One possibility would be the adjacent representation of an uncertainty model next to a single model representation [START_REF] Viard | Adjacent versus coincident representations of geospatial uncertainty: Which promote better decisions?[END_REF]. But many studies showed that even the communication of probabilities to lay audiences can be difficult (e.g. [START_REF] Spiegelhalter | Visualizing Uncertainty About the Future[END_REF], and it will not be easier with the derived concepts of entropy or stratigraphic variability.

Significant research has gone into the coincident representation of uncertainties, i.e. the integrated representation of uncertainty together with primary data or a model, for a more direct communication of uncertainties. [START_REF] Maceachren | Visualizing Uncertain Information[END_REF] and [START_REF] Pang | Approaches to uncertainty visualization[END_REF] describe a variety of methods to represent uncertainties of (continuous and categorical) scalar, vector and tensor data. These approaches can be classified into intrinsic approaches (changing object appearance), and extrinsic methods (adding symbols or features to provide uncertainty information). Intrinsic methods include the use of hue and saturation of color, blurring, and focus to distinguish uncertain model areas (e.g. [START_REF] Maceachren | Visualizing Uncertain Information[END_REF][START_REF] Djurcilov | Volume Rendering Data with Uncertainty Information[END_REF][START_REF] Djurcilov | Visualizing scalar volumetric data with uncertainty[END_REF][START_REF] Maceachren | Visualizing Geospatial Information Uncertainty: What We Know and What We Need to Know[END_REF][START_REF] Seipel | Color map design for visualization in flood risk assessment[END_REF]. There have also been studies on the interactive visualization (e.g. [START_REF] Kunz | Visualization of uncertainty in natural hazards assessments using an interactive cartographic information system[END_REF] and animated visualizations (e.g. [START_REF] Bostrom | Visualizing Seismic Risk and Uncertainty[END_REF], even including vibrations in regions of uncertainty (e.g. [START_REF] Brown | Animated visual vibrations as an uncertainty visualisation technique[END_REF].

Finally, we would like to mention that there are also approaches to represent uncertainties in geological models, which are not based on the analysis of a model ensemble, as considered here. These methods include plots of spatial data density (e.g. maps of kriging variance or similar approaches representing the distance to data), possibly combined with data accuracy (e.g. [START_REF] Berg | Synopsis of Current Threedimensional Geological Mapping and Modeling in Geological Survey Organizations[END_REF].

Overall, it is important to keep in mind that a representation of uncertainty should be adjusted to the expected target audience. This includes the question if accuracy or uncertainty should be represented (e.g. [START_REF] Buttenfield | Representing Data Quality[END_REF], an intrinsic or extrinsic representation, and the specific type of uncertainty visualization. Overall, it can be said that, for a coincident representation of uncertainty, transparency has been found to be an intuitive measure, even though it is difficult to visualize [START_REF] Viard | Adjacent versus coincident representations of geospatial uncertainty: Which promote better decisions?[END_REF], but more research into the perception of uncertainties in geological models is still required.

Uncertainty reduction

We described above methods to generate randomized realizations of geological models. Depending on the suitability of the prior distribution and the model structure itself, a subset of these realizations could be violating basic geological principles or simply lead to a structure, which is not conforming to the expected geological setting. In this case, we aim to define geological criteria to evaluate the validity of a generated model realization (Fig. 10d). In a similar sense, geological models can be related to geophysical measurements through the additional use of rock physics models (Fig. 10e).

Examples are density values assigned to specific lithological classes, which allow a comparison to measured gravity anomalies, acoustic impedance for a comparison to seismic measurements, and permeability for a comparison to flow and transport measurements. In all of these cases, criteria can be defined to describe the validity of the generated model, and even to consider the geological model as part of a geophysical inverse approach.

A related task to the evaluation of model results with these additional geological and geophysical criteria is to sort the generated set of models into specific classes, or to pick out "extreme" cases, to evaluate their effect in subsequent process simulations or predictions. It is often straight-forward for a structural geologist to perform such tasks for a reasonable set of models. However, in order to scale this process up to a large set of realizations and to perform this process automatically, we require clearly defined methods for a quantitative comparison of models to other models, and of model results to additionally available geological and geophysical information, either in the form of measurements, observations, or even general aspects of geological knowledge.

We can here distinguish two different, but related, types of measures:

• Distance functions: measures of distance between model realizations;

• Mistfit functions: measures of distance from a model realization to an external additional observation or measurement.

Both aspects are related, as some measures could theoretically be applied in both contexts (e.g. a specific outcome at a given location could be used for model comparison, but also to define a misfit to an actual observation at this location). In the following, we first describe methods based on geologic criteria (Fig. 10d), before providing an overview of the link to geophysical criteria (Fig. 10e).

Model analysis and model comparison based on geologic criteria

Measures for model comparison can be directly based on the result of the interpolation function m(x). In addition, values at multiple locations are often combined in an additional function, for example to estimate summary measures or geometric aspects. [START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF] already provided a detailed list of misfit functions that can be applied to determine a distance metric between a model realization and diverse types of geological information, including spatial aspects of geological observations (point and volume information), as well as lineations and topological measures.

Especially the aspect of topology, as described by [START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF], is interesting as a measure of misfit, as topological relationships are often of practical interest [START_REF] Pouliot | Reasoning about geological space: Coupling 3D GeoModels and topological queries as an aid to spatial data selection[END_REF]. Further descriptions of topology in geology, as well as measures for the comparison of topological graphs are described in Thiele et al. (2016a,b). This research has also shown that topological graphs of geological structures quickly become very complex in full 3-D settings. But the abstraction to this level provides interesting possibilities to use this information as geological criteria to validate geological model realizations.

Additional geodiversity metrics are described by Lindsay et al. (2013a) with the aim of model comparison. The authors extend the analysis to formation depth and curvature, and additionally include measures that are calculated on discretized representations of a geological model, such as volume and neighborhood relationships, to evaluate the juxtaposition of stratigraphic units, as an element of topology. An interesting aspect is also that the authors combine all metrics and use a principal component analysis to determine outliers and barycenter examples in the set of generated models. [START_REF] Pellerin | Elements for measuring the complexity of 3D structural models: Connectivity and geometry[END_REF] describe measures of complexity for geomodel representations in boundary models, based on definitions from the field of geometric modeling and visualization [START_REF] Rossignac | Shape complexity[END_REF][START_REF] Sukumar | Towards understanding what makes 3D objects appear simple or complex[END_REF]. The authors distinguish structural model components, such as the number of surfaces and corners in a grid and local statistical neighborhood measures and evaluate the use of these measures to define a complexity space to classify models by their complexity. Several measures overlap with the definitions in [START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF] and Lindsay et al. (2013a), but are here focused on the use for models in triangulated boundary representations.

Model dissimilarity has also been investigated by [START_REF] Schweizer | Uncertainty assessment in 3-D geological models of increasing complexity[END_REF] in an assessment of model scenarios with increasing complexity. The authors also evaluated the use of the Jaccard similarity measure and the extension to consider (cell) probabilities with the normalized city-block distance.

In the application to a geological data set, the authors found that especially the city-block distance provides a suitable measure for model dissimilarity, especially in the combination with a spatial visualization of uncertainties using Shannon entropy, as described in Sec. 4.2.2. Related approaches have also been taken in the field of reservoir modeling and history matching (e.g. [START_REF] Scheidt | Representing Spatial Uncertainty Using Distances and Kernels[END_REF][START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF][START_REF] Park | History matching and uncertainty quantification of facies models with multiple geological interpretations[END_REF] with distance measures for production forecasts. In addition to these measures, an interesting path for future research is the investigation of more abstract methods to compress model space. For example, discrete cosine transforms (DCT) have successfully been used by [START_REF] Lochbühler | Summary statistics from training images as prior information in probabilistic inversion[END_REF] to compress information from multiple training images in sparse representations.

the fact that the effect of gravity can be formulated as a boundary value problem, and therefore all density distributions with the same boundary effect can not be distinguished (e.g. [START_REF] Jacoby | Gravity Interpretation[END_REF]. This leads to the problem that additional constraints and regularization parameters have to be defined in order to obtain a stable inversion result (e.g [START_REF] Telford | Applied Geophysics[END_REF][START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF][START_REF] Jacoby | Gravity Interpretation[END_REF]. Additional depth-weighting methods can be applied to avoid a concentration of mass near the surface (e.g. [START_REF] Li | 3-D inversion of gravity data[END_REF]. The consideration of structural geological information from geological measurements is also explored in the work of [START_REF] Lelièvre | A comprehensive study of including structural orientation information in geophysical inversions[END_REF]. An interesting recent discussion on the problem of non-uniqueness in the interpretation of potential-field data has been contributed by [START_REF] Saltus | Unique geologic insights from "non-unique" gravity and magnetic interpretation[END_REF]. The authors directly state the relevance of prior knowledge for a successful interpretation of potential-field data and highlight that typical potential-field studies incorporate additional constraints to obtain reasonable solutions out of the infinite universe of possibilities.

The most widely used method to consider geological information in the form of geological models is the use of these models to initialize a subsurface parameterization for a geophysical inversion or process simulation (see also [START_REF] Jessell | Next Generation Three-Dimensional Geologic Modeling and Inversion[END_REF]. All of these approaches are based on the concept that domains with similar properties can be defined in the subsurface, based on similar geological processes in the formation of a region. This concept is clearly related to the fundamental considerations for the use of structural models that we described at the beginning (Sec. 1.1), and directly integrated into the diverse interpolation formulations in Sec. 3.1 and the formulation of a physical value of interest m(x) (e.g. Eq. ( 1)).

In order to obtain a meaningful model parameterization for geophysical inversion, relevant rock properties have to be assigned according to the lithological classes given by the geological model. We do not aim to describe the underlying concepts of rock physics for this step, but refer to the relevant literature (e.g. [START_REF] Mukerji | Statistical rock physics: Combining rock physics, information theory, and geostatistics to reduce uncertainty in seismic reservoir characterization[END_REF][START_REF] Mavko | The Rock Physics Handbook, 2nd Edition. Tools for Seismic Analysis of Porous Media[END_REF].

We will first consider geophysical inverse approaches based on the adjustment of geometric bodies. These methods are also referred to as the geometric property mode (following [START_REF] Silva | Potential-field inversion: Choosing the appropriate technique to solve a geologic problem[END_REF]. Many examples for this approach exist again for gravity and magneitc inversion and they are often based on 2.5-D assumptions of polygonal objects in a section with an infinite lateral extent, going back to the developments of [START_REF] Talwani | Rapid gravity computations for two-dimensional bodies with application to the Mendocino submarine fracture zone[END_REF] and [START_REF] Talwani | Computation with the help of a digital computer of magnetic anomalies caused by bodies of arbitrary shape[END_REF]. This approach is, for example, implemented in the software GM-SYS [START_REF] Won | Computing the gravitational and magnetic anomalies due to a polygon: Algorithms and Fortran subroutines[END_REF] and IGMAS+ [START_REF] Götze | Application of three-dimensional interactive modeling in gravity and magnetics[END_REF] and widely used for interactive geophysical modeling and inversion. These approaches often implement map-and extrusion-based explicit modeling approaches (Sec. 3.2.1), which also makes a combination with probabilistic geophysical inversions directly possible (e.g. [START_REF] Hauser | A probabilistic seismic model for the European Arctic[END_REF][START_REF] Hauser | Probabilistic inversion of airborne electromagnetic data for basement conductors[END_REF]. However, as described before, map-and extrusion-based approaches are limited to geological settings with a relatively low complexity (Fig. 2a).

Mesh-based geophysical inverse approaches in the physical property mode [START_REF] Silva | Potential-field inversion: Choosing the appropriate technique to solve a geologic problem[END_REF]) are, instead, not limited to a specific complexity of the geological model, but require an initial discretization of the geological model into a grid structure (Fig. 10c). On this basis, geological models are often applied to obtain a subsurface parameterization for geophysical forward simulations and inversions. Still to date, many studies use geophysical forward simulations and then manually adjust the geological model based on a visual or quantitative comparison between forward model and geophysical measurements (e.g. [START_REF] Gradmann | Integrated geophysical modelling of a lateral transition zone in the lithospheric mantle under Norway and Sweden[END_REF][START_REF] Autin | Deep structure of the Argentine margin inferred from 3D gravity and temperature modelling, Colorado Basin[END_REF][START_REF] Haase | A 3D regional crustal model of the NE Atlantic based on seismic and gravity data[END_REF]. However, discretized geological models have also long been used as a basis for geophysical inversion routines (e.g. [START_REF] Bosch | Lithologic tomography: From plural geophysical data to lithology estimation[END_REF][START_REF] Bosch | Lithologic tomography: an application to geophysical data from the Cadomian belt of northern Brittany, France[END_REF][START_REF] Fullagar | Constrained inversion of geologic surfacespushing the boundaries[END_REF][START_REF] Guillen | Geological modelling from field data and geological knowledge: Part II. Modelling validation using gravity and magnetic data inversion: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] Moorkamp | A framework for 3-D joint inversion of MT, gravity and seismic refraction data[END_REF]. Main differences in the approaches are mostly in the type of geological criteria that can be defined, for example in the ability to fix specific known geometric objects, or to ensure neighborhood relationships between cells.

These concepts have been applied successfully in various gravity and magnetic studies [START_REF] Martelet | Integrated 3D geophysical and geological modelling of the Hercynian Suture Zone in the Champtoceaux area (south Brittany, France)[END_REF][START_REF] Joly | A multidisciplinary study of a syntectonic pluton close to a major lithospheric-scale fault -Relationships between the Montmarault granitic massif and the Sillon Houiller Fault in the Variscan French Massif Central: 1. Geochronology, mineral fabrics, and tectonic implications[END_REF][START_REF] Joly | A multidisciplinary study of a syntectonic pluton close to a major lithospheric-scale fault-Relationships between the Montmarault granitic massif and the Sillon Houiller Fault in the Variscan French Massif Central: 2. Gravity, aeromagnetic investigations, and 3-D geologic modeling[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF]Lindsay et al., 2013b;[START_REF] Wehr | A three-dimensional model of the Pyrenees and their foreland basins from geological and gravimetric data[END_REF], but it has to be noted that the geological model is here merely an initial step to parameterize the geophysical inversion. It is not only possible, but potentially expected, that the final property distribution in the inverted model does not conform to the geological constraints and input data that were considered in the initial geological model used to obtain the parameterization [START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF][START_REF] Lindsay | Inversion and Geodiversity: Searching Model Space for the Answers[END_REF][START_REF] Jessell | Next Generation Three-Dimensional Geologic Modeling and Inversion[END_REF].

The problem of this loss of the link to the primary information in the geological model has recently been addressed by [START_REF] Giraud | Uncertainty reduction through geologically conditioned petrophysical constraints in joint inver-sionConditioned petrophysical constraints[END_REF] combining uncertainty propagation with a probabilistic geological forward modeling method (Sec. 4.1) with joint geophysical inversions using a probability grid for each geological unit, instead of a fixed outcome in a cell (see also the calculation of probability grids in Sec. 4.2.2, Eqn. ( 17)). This probability value is then combined with petrophysical models to obtain starting models and constraints for joint geophysical inversion and results suggest that the inversion routine properly considers areas of higher geological certainty. Related to this approach is the work by [START_REF] Wellmann | Uncertainty estimation for a geological model of the Sandstone greenstone belt, Western Australia -insights from integrated geological and geophysical inversion in a Bayesian inference framework[END_REF], in which geological interface points of an implicit forward model are directly treated as parameters in a probabilistic inversion, and geological, as well as geophysical, likelihood functions are used to obtain posterior distributions of these geological input parameters. On the basis of the posterior distributions, multiple geometric model representations can then be obtained and further processed for uncertainty quantification and visualization.

As another recent example, Zheglova et al. (2018a) use implicit modeling of several geological interfaces coupled with petrophysical data for solving a joint gravity and seismic topography problem.

In the same spirit, [START_REF] Cherpeau | Method for Stochastic Inverse Modeling of Fault Geometry and Connectivity Using Flow Data[END_REF] propose to assess whether flow data can reduce uncertainties about fault locations and connectivity constrained by sparse geological cross-sections and prior distributions of fault parameters for fault families. They use a sequential stochastic fault model which associates fault observations based on prior information. Using the random numbers used in this algorithm as parameter vectors, they run a Monte Carlo Markov Chain (MCMC) algorithm to explore the parameter space. To include the possibility of having a different number of faults, they translate the prior probability of existence as modified distribution of fault size.

Other inverse models using a variable number of parameters are also proposed by [START_REF] Sambridge | Transdimensional inference in the geosciences[END_REF] and references therein. One of the key ideas is to use the number of model parameters as a variable in Bayesian inversion. This method has been used for instance for layer-cake geological models with a variable number of layers represented by a one-dimensional Voronoi diagram. An interesting aspect of these trans-dimensional inversion methods is that they tend to favor the models which can explain observations with the smallest number of parameters (parsimonious solution).

An additional interesting path combining geophysical inversion with geological knowledge is the concept of interactive inversion [START_REF] Boschetti | Interactive inversion in geosciences[END_REF][START_REF] Wijns | Inverse modelling in geology by interactive evolutionary computation[END_REF][START_REF] Wijns | Interactive geophysical inversion using qualitative geological constraints[END_REF], which aims at an explicit integration of expert knowledge in the inversion process. This is achieved through a system generating possible scenarios. These scenarios are visualized on a screen, and experts select and rank feasible solutions. In a next step, a genetic algorithm is used to propose new realizations from these selected models. In a sense, this approach is between classical deterministic inversion with "trial-and-error" adjustments of settings (and similar approaches to evaluate the model space) and more integrated approaches, where model differences are included in the form of likelihood functions (e.g. [START_REF] Wellmann | Uncertainty estimation for a geological model of the Sandstone greenstone belt, Western Australia -insights from integrated geological and geophysical inversion in a Bayesian inference framework[END_REF]. An advantage of this interactive approach is that realizations can be produced and evaluated which have not been expected before (in difference to the defined likelihood functions). A potential disadvantage is the limited explicit statement by the expert on the decision, and therefore a limited reproducibility of results. However, the combination of geophysical inversion with a supervised machine learning approach is an interesting path to integrate different forms of geological knowledge, including tacit expert knowledge, in geophysical inversions.

Discussion and conclusions

Geological modeling methods have been rapidly developing over the last years-from improved representations in both explicit and implicit approaches, to more complex structural interactions and fault and fold models, to a more comprehensive consideration of uncertainties. With this review, we aimed to provide an overview of these developments and recent applications in the diverse fields using 3-D geological modeling approaches. One recurring aspect is the link to geophysical data processing and inversion approaches, both because geophysical data and interpretations are an essential input to conventional geological modeling studies, but also because geological concepts and models are increasingly used in joint geophysical and geological inversions.

Even though we provided an extensive description of 3-D structural geological models, the underlying concepts, methods and uncertainties, we necessarily had to leave out several aspects in order to keep a focus on the most relevant topics, but also due to our own expertise. We have certainly omitted to describe some related work, as it is barely possible to cover the entire breadth of such an interdisciplinary topic. Also, several of our descriptions are certainly open to discussion, as no single best-fitting and all-purpose geological modeling and uncertainty analysis method exists. In the following, we will contribute to this discussion and, finally, conclude with some suggestions and potential paths for future research in this field that we envisage for the coming years.

Discussion

On the two viewpoints of model construction

We initiated the chapter with the description of the two common viewpoints on subsurface property and interface detection. The separation into two distinct approaches is clearly an over-simplification, but it captures the two main viewpoints on this problem and it is also, in our experience, the cause for a lot of misunderstanding between researchers trained in different fields.

Each approach in itself is open to criticism. The "data-driven approach", being parsimonious, may lead to surprises and errors. This occurs in particular when the geological features below the geophysical resolution exhibit a significant influence the model forecast (e.g. [START_REF] Julio | Impact of the en echelon fault connectivity on reservoir flow simulations[END_REF]. In addition, data-driven inversions can lead to subsurface realizations that are in clear conflict with geological concepts or observations (e.g. [START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF]. One can argue that this is not relevant if a reasonable prediction can still be obtained, but it would then mean that not all available information is considered, and also that some tacit assumptions hidden in an so-called objective model are inappropriate [START_REF] Journel | The abuse of principles in model building and the quest for objectivity[END_REF].

Conversely, a common criticism for the geomodel-based approach is that (a) it starts with subjective information which may induce bias in the model forecasts, (b) it contains interpretative elements which cannot be verified, and (c) it introduces a level of complexity that is not supported by data. Concerning subjectivity, we relate again to the various subjective choices that are also inherent in the "data-driven" approach and also state that subjectivity is not necessarily bad, as long as it is properly acknowledged and considered [START_REF] Curtis | The science of subjectivity[END_REF]. The use of all available a priori information is also strongly expressed by [START_REF] Tarantola | Popper, Bayes and the inverse problem[END_REF]. Interpretative elements are certainly part of any geological investigation, and we highlight here again the hermeneutic aspect of geology [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF]. This is an aspect which we simply can not avoid, and literally in the nature of complexity of geological investigations.

Another aspect that requires some consideration is the construction of geological models and the later usage of these models. It is a common statement that models are built for a purpose, and that this defined purposes not only is essential to define the scale of the model, but also the level of abstraction, etc. In the geoscientific literature, these approaches have recently been popularized under the general theme of "fit-for-purpose" models, or "purpose-driven" geological models [START_REF] Ringrose | Reservoir Model Design[END_REF].

In the context of geological modeling, these purpose-driven modeling approaches are typically applied in many engineering applications, for example in reservoir modeling, but also mineral exploration and geotechnical applications. In geological modeling, however, there are many other important usages, which can not so easily be linked to a defined purpose at the outset of the modeling study. This aspect can directly be seen when recalling that geological models are, at least partly, treated as the extension of a classical 2-D geological map into the third dimension, and it is evident in the expanding role of geological models is in the context of work of geological surveys (e.g. [START_REF] Kessler | The Capture and Dissemination of Integrated 3D Geospatial Knowledge at the British Geological Survey Using GSI3D Software and Methodology[END_REF][START_REF] Berg | Synopsis of Current Threedimensional Geological Mapping and Modeling in Geological Survey Organizations[END_REF]. In this sense, geological models can be understood as an effort to extend geological knowledge, commonly and previously captured in 2-D maps, now in 3-D geological models. In these cases, it is obvious that the models are a foundation for subsequent discussions and usages which can not be foreseen at the initial state of model construction.

Similar considerations hold in the field of hazard prediction, also a central aspect of geological surveys. Even though the model purpose is then better defined, it is not easily possible to provide a defined quantity of interest. For example, in the case of ground motion hazard predictions for Southern California [START_REF] Shaw | Unified Structural Representation of the southern California crust and upper mantle[END_REF], a model was created to test various earthquake scenarios given ruptures at specific locations to contribute to outreach and public policies. Also in this case, the potential interest in this model will be very different for people living at different locations.

In these cases, the geological model becomes relevant also as a central communication tool. This is important, as it does not only limit the possibility to evaluate a meaningful scale and level of detail at the beginning of the modeling campaign. Also, the presentation and visualization of the model now becomes a central aspect, combined with a representation of model quality and uncertainty.

On the choice of a suitable geomodeling method

The various methods described in Sec. 3 provide ways to produce 3D geometries from subsurface data and various degrees of conceptual knowledge.

The ability of mathematical interpolation (or approximation) to reach acceptable results has been much debated by geologists. Indeed, the ability of computers to automatically produce results acceptable by skilled structural geologists can be questionable. Early on, [START_REF] Walters | Contouring by Machine: A User's Guide[END_REF] argued that "to be effective, such a [contouring] program must allow the user to make judgments during the preparation or modification of the contoured map to enable the specialist to introduce his background". A reason for invoking expert input is that available observations never fully determine model features, which makes subsurface modeling an interpretation exercise. The interpretative aspect is, therefore, a feature to distinguish between different modeling methods. Several tools have been developed with the purpose of providing methods that allow to "draw" geological structures in 3-D, and aid in the use and combination of different data sets in this endeavor. These methods have been very successful, and widely employed, for example in the country-wide geological model of Britain [START_REF] Kessler | The Capture and Dissemination of Integrated 3D Geospatial Knowledge at the British Geological Survey Using GSI3D Software and Methodology[END_REF]. Alternatively, several methods aim at incorporating certain aspects of geology in the interpolation step itself, for example using locally variable anisotropy [START_REF] Martin | Iterative refinement of implicit boundary models for improved geological feature reproduction[END_REF] or prior knowledge of structural style (Caumon et al., 2013a). This also allows a significant interpretative contribution, but combines it directly with mathematical models that automatically consider certain (and perceived as important) aspects of geology. We posit that a lot of arguments for and against a certain geological modeling method are based on the perceived relevance of subjective (prior) knowledge and on how this knowledge is incorporated in the method (graphically by experts or automatically by algorithms and the associated parameters).

Overall, numerical structural modeling involves mainly two components: interpolation to fill the gaps between data points, and CAD operators to intersect, trim and weld the model elements so as to honor conceptual geological knowledge.

A difficulty at this point is to provide a comparison of the various numerical methods that does better than mentioning how "easy to use" or "flexible" they are. In practice, indeed, the level of one's familiarity with a particular software, some small but useful improvements to a particular im-plementation or simply different user interfaces for the same method can easily influence viewpoints. Therefore, before trying to summarize the strengths and the weaknesses of the various methods, we start by giving some evaluation criteria.

Expected level of geological complexity As already mentioned, the expressiveness of a method to simple of complex geological structures should be carefully considered complexity expected in the domain should be considered (Fig 2.

Computational and memory requirements The methods reviewed in this section do not have the same memory and computational requirements. As a rule of thumb the older methods are most efficient, whereas more recent methods are more demanding.

Ability to honor various data types Not all methods are equivalent in the type of data type they can honor.

Ability to incorporate geological knowledge

The various methods all use some sort of geological knowledge, if only by letting the modeler add interpretive data in the form of maps, cross sections or additional points. However, considering interactions between several objects (e.g., allowing for various folding styles) remains challenging with many three-dimensional methods.

Number of input parameters

An objective measure of a method's simplicity is the number of parameters it involves. Note that this criterion can go against the ability to incorporate geological knowledge. In this sense, a single mathematical criterion parametrized by three or four values is preferable to one thousand mouse picks if both approaches produce a comparable result.

Model updating A geological model is seldom engraved in stone. New data may come in, new conceptual models may arise, which both require model updates. Uncertainty quantification also needs to generate not one but several possible models. It may be worth asking how many operations and how much the computational efforts are needed to update a particular model.

Modeling purpose

As discussed in the Introduction, a geological model generally serves one or several purposes. The compatibility of a chosen modeling strategy to meet these purposes is important: the computation of a volume or the simulation of multiphase flow do not have the same implications.

With these criteria in mind, we now highlight some distinctions between explicit and implicit models, summarized in Table 1. We also refer to [START_REF] Jessell | Next Generation Three-Dimensional Geologic Modeling and Inversion[END_REF] and to [START_REF] Collon | 3D modeling from outcrop data in a salt tectonic context: Example from the Inceyol minibasin, Sivas Basin, Turkey[END_REF] for comparisons and discussions of the relative merits of the various numerical approaches.

Extrusion methods (depth maps and pillar approaches) are mature and have been driven by significant applications in reservoir modeling and near-surface applications (hydrogeology and geotechnical engineering). They are conceptually simple and the associated numerical representations are compact in terms of memory and computationally efficient. They are directly compatible with prismatic gridding methods which also use surface extrusion. In particular, they can create any stratigraphic surface in the volume of interest by simple vector operations. Their main limitation is the lack of flexibility to conform to multiple geological surfaces which may be oriented along arbitrary directions. This is particularly problematic in highly deformed structural settings. Additionally, there is no direct and robust way to account for off-surface structural orientation data.

Three-dimensional explicit methods have also reached a high level of maturity. They have virtually no representational limitation in approximating geological interfaces. As they are based on twodimensional surfaces, they also allow for a compact description of volumes. Interpolation acting on two-dimensional graphs is also relatively efficient. However, they generally involve the mastery of a relatively large set of geometrical and topological operations in order to be applied effectively. These operations (e.g., the intersection of triangulated or NURBS surfaces), although seemingly simple, are very difficult to implement in a robust way due to limited floating point accuracy. Model construction can be automated in relatively simple structural cases, but their application in complex domains generally involves significant manual editing. As data must be projected onto evolving surfaces, projections may cross discontinuities and generate local artifacts. Off-surface structural orientation data can be accounted for, assuming sedimentary beds are parallel along the projection line. Layer thickness may be locally honored using iterative projections [START_REF] Salles | Migration of a synclinal depocentre from turbidite growth strata: The Annot syncline, SE France[END_REF], but no approach so far has automatically used cleavage or axial directions data in explicit methods. Model updating is possible but delicate if one want to maintain a valid boundary representation [START_REF] Caumon | Constrained modifications of non-manifold b-reps[END_REF]. In terms of applications, explicit boundary representations may need to be simplified to create volumetric grids, involving stair-step approximations. Alternatively, they can create a frame for unstructured meshes, but the geostatistical modeling and process simulation are often more difficult on these meshes than on extruded regular grids.

Implicit methods are more recent, but they have raised significant enthusiasm, and have been implemented in several commercial software such as Leapfrog, Geomodeller, SKUA-GOCAD and Petrel (Volumetric-based modeling). Several features explain this growing interest. First, implicit methods consider volumes right away, which overall more robust and automatic than surface-based models. Overall, implicit methods leave a lot to interpolation; this means that complex shapes can be obtained automatically without need for expert input. As compared to surface-based methods, interpolation can also change the topology of the implicit surfaces. This may be a strength when the topology is unknown (e.g., for karst or ore body modeling), but may also be a problem when the topology should be controlled (for instance to prevent stratigraphic surfaces to form blobs), see [START_REF] Collon | 3D modeling from outcrop data in a salt tectonic context: Example from the Inceyol minibasin, Sivas Basin, Turkey[END_REF]. Another capability of implicit methods is to honor several types of structural data without any projection [START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF]Caumon et al., 2013a;[START_REF] Hillier | Three-Dimensional Modelling of Geological Surfaces Using Generalized Interpolation with Radial Basis Functions[END_REF]. A continuum of conformable stratigraphic surfaces can be generated, as in extrusion approaches, but without limitations in the type of fault network. Conversely, the management of thickness variations is complex to handle in implicit methods [START_REF] Laurent | Iterative Thickness Regularization of Stratigraphic Layers in Discrete Implicit Modeling[END_REF], as the various regularization terms used in interpolation tend to penalize simultaneously the variation of surface curvature and the variation of layer thickness. Even if recent work on anisotropy may provide interesting solutions for this [START_REF] Martin | Iterative refinement of implicit boundary models for improved geological feature reproduction[END_REF], current implicit methods can only represent a limited number of unconformities owing to computational and storage limitations. In terms of model updating, the automation of interpolation has significantly helped the development of new structural uncertainty modeling techniques (e.g., [START_REF] Cherpeau | Stochastic simulations of fault networks in 3D structural modeling[END_REF]Wellmann et al., 2010b;[START_REF] Lindsay | Locating and quantifying geological uncertainty in three-dimensional models: Analysis of the Gippsland Basin, southeastern Australia[END_REF][START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF].

In deterministic model updates, the lack of direct surface manipulation, can be a source of frustration. Indeed, the only option for local model updating is to add new data points and reinterpolate. In this case, methods that use basis functions with a local support are expected to perform more efficiently than implicit methods with a global support. As surface-based methods, implicit modeling techniques rely at some point on topological operations to finally create an image, a boundary representation or a mesh. These operations, acting on volumetric data structures, can be more difficult to implement than on surfaces.

In terms of application, a very interesting extension of implicit methods is the possibility to map every point of the subsurface into depositional space. this calls not only for the relative geological time but also two paleo-geographic coordinates [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF]; [START_REF] Moyen | 3D-Parameterization of the 3D Geological Space -The GeoChron Model[END_REF]; [START_REF] Mallet | Elements of Mathematical Sedimentary Geology: The GeoChron Model[END_REF]; [START_REF] Wu | Horizon volumes with interpreted constraints[END_REF]; [START_REF] Wu | Building 3D subsurface models conforming to seismic structural and stratigraphic features[END_REF]. This representation has also been used directly for geomechanical restoration of both preserved and eroded strata [START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF] and it is directly compatible with the extended finite element method [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model: NON-PLANAR 3D CRACK GROWTH-PART I[END_REF]. It has also been used to create reservoir simulation grids [START_REF] Gringarten | New Grids for Robust Reservoir Modeling[END_REF].

On the classification of geological uncertainties based on uncertainty quantification concepts

In Sec. 2.1, we described widely used classifications for uncertainties in geological models into different types [START_REF] Mann | Uncertainty in geology[END_REF][START_REF] Bardossy | Evaluation of Uncertainties and Risks in Geology: New Mathematical Approaches for their Handling[END_REF]. However, a difficulty in these classifications lies in the understanding of type 2 uncertainty ("inherent variability"). Indeed, it could be argued that there is no "inherent variability" in the natural world, as long as we stay above the scales and processes of quantum physics (e.g. leave out processes of radioactive decay, etc.)10 . [START_REF] Mann | Uncertainty in geology[END_REF] seems to expect a potential upscaling of Heisenberg's principle to larger scales and argues that some variables can not be measured without altering the rock. This is, however, not similar to the uncertainty in quantum physics, as new measurement types could be able to perform these accurate measurements, and it is thus only related to a lack of knowledge (see below). Consider the case of porosity distribution in a reservoir: we typically treat porosity as a stochastic variable and then use a random function concept to describe the spatial distribution (e.g. [START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty, 2nd Edition[END_REF][START_REF] Pyrcz | Geostatistical Reservoir Engineering[END_REF], and we could thus be inclined to classify porosity as a property underlying inherent variability or stochasticity, in the sense of Mann's "Type 2". However, using computer tomography (CT) methods, it is now possible to image porosity (almost) exactly (e.g. [START_REF] Wildenschild | Using X-ray computed tomography in hydrology: systems, resolutions, and limitations[END_REF]-it is thus evident that porosity is not underlying an inherent variability. Clearly, it is not feasible (at least with the currently available measurement techniques) to measure porosity exactly on the scale of an aquifer or a reservoir. We thus choose to treat porosity as a random variable, and then assign appropriate functions to model the distribution in space, for which we then can estimate parameters on the basis of observations11 . We consider this example of porosity as insightful, but similar considerations would hold for the random treatment of interface positions in geological models12 . Another practical avenue to classify uncertainties is to consider the viewpoint of uncertainty quantification. A common approach is then to consider two classes of uncertainties [START_REF] Kiureghian | Aleatory or epistemic? Does it matter?[END_REF]:

• Epistemic (systematic) uncertainties: correspond (etymologically) to uncertainties that can be described in the current state of knowledge. They covers things one could know in principle if enough information was gathered, but which cannot be observed, measured, or modeled accurately enough;

• Aleatory (statistical) uncertainties: relate the uncertainties that cannot be possibly reduced or explained other than with randomness.

As the line between aleatory and epistemic uncertainties is not very easy to draw, [START_REF] Kiureghian | Aleatory or epistemic? Does it matter?[END_REF] suggest to distinguish epistemic and aleatory uncertainties in the context of a specific model. Recognizing that, in principle, all uncertainties could be reduced if we would be able to measure exhaustively, or with a measurement device and principle which may be available in the future, or with more surveys that could in principle be made, we can consider all uncertainties as epistemic in nature. However, we separate those uncertainties that could be reduced in the near term from those that we consider as irreducible in the current state of a specific model. In a recent discussion on uncertainties in hydrological models, [START_REF] Beven | Facets of uncertainty: epistemic uncertainty, non-stationarity, likelihood, hypothesis testing, and communication[END_REF] classifies those uncertainties as aleatory, for which a stationary statistical variation can be defined, and for which, therefore, the full range of statistical methods is applicable. It can be argued that this separation also adds transparency to decision-making [START_REF] Kiureghian | Aleatory or epistemic? Does it matter?[END_REF], as we clearly describe our choices about which uncertainties to consider in which way. The aim is then to further move from epistemic to aleatory residual errors, a step that could be possible with further experimentation and investigation, at least for a part of the epistemic uncertainties.

In addition to the treatment of epistemic and aleatory uncertainties, there is always the risk of "unknown unknowns", an aspect that can never be avoided. It is therefore even more important to clearly communicate the "known unknowns" that have been omitted from the analysis, for example because of lack of knowledge, lack of computing power, or lack of time for a more detailed investigation (e.g. [START_REF] Beven | Facets of uncertainty: epistemic uncertainty, non-stationarity, likelihood, hypothesis testing, and communication[END_REF].

In the entire discussion so far, we have left out the uncertainty related to an applied belief system, an aspect that can be considered as an ontological uncertainty, for example the question whether probabilistic methods are suitable to represent beliefs about the nature of the residuals (e.g. [START_REF] Beven | Facets of uncertainty: epistemic uncertainty, non-stationarity, likelihood, hypothesis testing, and communication[END_REF]. This aspect is relevant, as it is related to the framework of knowledge, without which a definition of uncertainties themselves does not make sense (e.g. [START_REF] Nearing | Ensembles vs. information theory: supporting science under uncertainty[END_REF]), but we do not treat it in detail here and refer to the mentioned references for more information.

It remains the point that there is no single definition of types and sources of uncertainties in geological models. We are yet only at the beginning of a debate about useful definitions for subsurface models (e.g. [START_REF] Caers | Modeling Uncertainty in the Earth Sciences[END_REF], a process that has been ongoing in related fields for many decades (e.g. in hydrology [START_REF] Montanari | What do we mean by 'uncertainty'? The need for a consistent wording about uncertainty assessment in hydrology[END_REF][START_REF] Beven | Facets of uncertainty: epistemic uncertainty, non-stationarity, likelihood, hypothesis testing, and communication[END_REF][START_REF] Nearing | A philosophical basis for hydrological uncertainty[END_REF]. In our point of view, the consideration of aleatory and epistemic uncertainties forms an interesting framework, as this will provide a better link to concepts of uncertainty quantification in related fields. In addition, a reference to the modeling stage in which uncertainties are evaluated helps putting the work into reference (see Fig. 5). This combination will help to clarify uncertainties that have been considered in an investigation, and those that have been omitted-for reasons of deliberate choice, or because they are irreducible uncertainties in the current state [START_REF] Kiureghian | Aleatory or epistemic? Does it matter?[END_REF]. It should be made absolutely clear that the sentence "How do you know that you considered all uncertainties?" (as common criticism to the investigation of uncertainties) does not make any sense, but it is up to the modeler to be completely clear and transparent about the uncertainties that have been considered, and those that have not. Or, as [START_REF] Linde | On uncertainty quantification in hydrogeology and hydrogeophysics[END_REF] formulate very appropriately, "UQ in Earth Sciences can never be considered to be complete. Instead it should be viewed as a partial assessment that is valid for a given set of prior assumptions, hypotheses, and simplifications."

On the consideration of uncertainties in geological models

We presented in Sec. 4 several considerations and practical aspects for the quantification of uncertainties in 3-D geological models. The approaches mainly differ in the level of geological complexity that they allow to integrate (see Fig. 2), and the classes of uncertainties in the geological model construction step that they consider. However, in the description, we considered aspects of uncertainty that are related to modeling stages where, at least, the structure of the mathematical model has already been defined (see Fig. 5, c andd). In the description above, this refers to a choice of the type of basis functions Φ, and often also the topological description τ (i.e. the relationship between surfaces of geological interfaces and faults). But as we already pointed out in the section about geological uncertainties (Sec. 2), many authors stated the important relevance of uncertainties related to the initial conceptual model (Fig. 5a) and the structure of the mathematical model (Fig. 5b). These uncertainties are very difficult to capture in practice.

A first additional step is the consideration of uncertainties related to the structure of the mathematical model, for example: the number of geological interfaces and their stratigraphic and tectonic relationship, the number of faults and the topology of the fault network. The specific elements depend on the applied modeling method (see Sec. 3.2), but generally, these are all aspects that are related to a single mathematical interpolation function that has been selected for the modeling purpose. In the general model formulation of Eq. ( 15), this corresponds to a selection of the type of basis functions Φ, but possible differences in α, κ and τ . We will refer to differences on this level as different model scenarios in the following. Although several methods have proposed to address these uncertainties using sophisticated random sampling methods (e.g., [START_REF] Holden | Stochastic Structural Modeling[END_REF][START_REF] Cherpeau | Stochastic structural modelling in sparse data situations[END_REF][START_REF] Aydin | Quantifying structural uncertainty on fault networks using a marked point process within a Bayesian framework[END_REF][START_REF] Lallier | Uncertainty assessment in the stratigraphic well correlation of a carbonate ramp: Method and application to the Beausset Basin, SE France[END_REF], it is unclear at this point how they would compare to many scenarios generated by many experts.

It is evident that decisions at this stage can have an important impact on the subsequent modeling study, essentially also because these decisions are difficult to adjust in conventional geomodeling approaches. This aspect has brought several authors to argue that it is important to construct several "deterministic" models and that uncertainties in these models are more relevant than uncertainties in the stochastic model realizations (e.g. [START_REF] Bentley | Scenario-based reservoir modelling: the need for more determinism and less anchoring[END_REF][START_REF] Rowbotham | Bias in geophysical interpretation-the case for multiple deterministic scenarios[END_REF][START_REF] Ringrose | Reservoir Model Design[END_REF][START_REF] Bentley | Future directions in reservoir modelling: new tools and 'fit-forpurpose' workflows[END_REF].

The idea to use several model scenarios relates back to the visionary paper by Chamberlin (1897) and the concept of "multiple working hypotheses", in which he already foresees many aspects of human bias that have been formalized almost a century later by [START_REF] Tversky | Judgment under Uncertainty: Heuristics and Biases[END_REF]. Following this concept, one should always formulate multiple hypotheses, in order to avoid too much "affection" (which can be understood as a combination of several bias factors) to a specific hypothesis. In the context of geological modeling, we can interpret this as an imperative to attempt to construct multiple model scenarios, each valid as a representation of the geological concept.

When we start at the level of the structure of the mathematical model, then this approach is, in fact, feasible. For example, [START_REF] Corbel | Framework for multiple hypothesis testing improves the use of legacy data in structural geological modeling[END_REF] described a framework for the consideration of multiple geological hypotheses when only legacy data are available, with a related high uncertainty. [START_REF] Fichtner | The Collaborative Seismic Earth Model: Generation 1[END_REF] also recently proposed a Bayesian approach to reconcile several global or regional tomographic models created at different scales, under Gaussian assumptions. In fact, several authors also used different model scenarios as input for multiple stochastic model realizations on the basis of each of these scenarios. An interesting example is presented in the work of [START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF], where then model realizations from a large set of generated models are selected on the basis of distance measures (see also Sec. 4.3.1) for subsequent flow simulations. Related approaches have also been taken in the field of hydrogeology, although with simpler geological modeling methods, to evaluate the effect of uncertainties in the conceptual geological model on the predictive ability of a hydrogeological model (e.g. [START_REF] Refsgaard | A framework for dealing with uncertainty due to model structure error[END_REF][START_REF] Nilsson | Uncertainty in geological and hydrogeological data[END_REF][START_REF] Troldborg | The importance of alternative conceptual models for simulation of concentrations in a multi-aquifer system[END_REF][START_REF] Seifert | Assessment of hydrological model predictive ability given multiple conceptual geological models[END_REF].

One important question is how these different model scenarios are best combined with subsequent stochastic simulations in one uncertainty quantification framework. One possibility is the generation of a set of model realizations, based on different scenarios, and to apply suitable distance measures, as already mentioned above for the work of [START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF], but also used widely in the field of multiple point geostatistics (e.g. [START_REF] Demyanov | Geological feature selection in reservoir modelling and history matching with Multiple Kernel Learning[END_REF][START_REF] Demyanov | Uncertainty Quantification in Reservoir Prediction: Part 2-Handling Uncertainty in the Geological Scenario[END_REF][START_REF] Park | History matching and uncertainty quantification of facies models with multiple geological interpretations[END_REF]Jung et al., 2013, see also Sec. 4.3.1).

In the context of probabilistic inverse frameworks, transdimensional MC methods provide a potentially attractive path [START_REF] Sambridge | Transdimensional inference in the geosciences[END_REF], but they have so far only been applied to relatively simple geological model structures. An alternative would be the use of the possibility to perform structured model selection in a Bayesian framework for multiple model scenarios (see appendix, Sec. A.1). Similar approaches have been taken in hydrogeology (e.g. [START_REF] Schöniger | Model selection on solid ground: Rigorous comparison of nine ways to evaluate Bayesian model evidence[END_REF], and an application to geomodeling can be envisaged.

Several authors have also argued that probabilistic methods are not anymore suitable when considering epistemic uncertainties related to the conceptual model (e.g. [START_REF] Bardossy | Evaluation of Uncertainties and Risks in Geology: New Mathematical Approaches for their Handling[END_REF] and have promoted the use of extended probability methods, for example fuzzy logic [START_REF] Zadeh | Fuzzy sets*[END_REF][START_REF] Dubois | Fundamentals of Fuzzy Sets (The Handbooks of Fuzzy Sets)[END_REF]. These approaches are appealing, but they have so far not been used extensively in recent uncertainty studies of geological models.

Even though it is difficult to classify all approaches, we combined the discussed methods in an overview diagram in Fig. 15. One important divide exists between the consideration of uncertainties about the conceptual model and all subsequent steps. Once the modeling method has been decided, multiple approaches exist for several purposes and levels of geological complexity. Finally, we briefly described the link between geological models and geophysical inversions, based on several degrees of complexity in the structural models. Many developments in recent years work towards a tighter integration of information from both viewpoints, and we are confident that this will remain an active field of research, resulting in an even better combination of data and knowledge to investigate heterogeneities in the subsurface.

Modeling step and uncertainty class

Conclusions and outlook

After the above review of geological modeling methods, we now aim to outline several interesting perspectives for future research in the field, as well as some recommendations, which will hopefully lead to a better integration of geological and geophysical concepts. We identify these aspects by keywords:

Integration The consideration of all available information with geological knowledge remains a serious hurdle, not only because of methodological or fundamental limitations, but also because of the different culture and value frame between geologists, geophysicists, and engineers. Even though experts are needed in each sub-discipline, teaching the value of complementary perspectives on the subsurface of the Earth remains an educational challenge. Somehow, we would like to get away from strong feelings about different modeling approaches, and rather adopt a more fact-based consideration of the best possible method for a given task-which calls for explicitly defining what "best" means.

Communication Uncertainties make human beings uncomfortable; our brain (and our ego) is trained to develop convictions and to make choices to the point that choices are often implicit and not explicitly stated. Scientific papers are also more easily published when deemed "convincing". We strongly encourage geoscientists to more comprehensively consider uncertainties in all modeling endeavors, or at least to clearly acknowledge modeling assumptions.

Tools We encourage geoscientists to consider a suitable modeling algorithm and not just a modeling software. We are well aware that it is tempting to put a high level of trust in the result of a, potentially very expensive, modeling software. But the core decision should be based on the most suitable approach for a task at hand.

Visualization The visualization of spatial uncertainty in three dimensions remains a difficult and important challenge both to help making decisions and for communication purposes.

Geological realism Although geological modeling methods have made significant progresses incorporating geological knowledge, a significant number of geological concepts and observations

are not yet easily integrated in the available interpolation methods.

Machine learning

The advent of efficient machine learning methods opens very interesting perspectives in the field of joint geological modeling and geophysical inversion. Many interesting approaches have been around for a while and we expect interesting developments in the future, capitalizing on advances in formalizations of geological knowledge.

Efficiency In practice, the computational efficiency of modeling methods and inverse methods remains limited with regard to the ambition of creating reasonably accurate models explaining all observations. The development of efficient joint geological and geophysical inversions exploiting modern parallel computer architectures (and tomorrow, quantum computers) is an essential area for future advances.

Scale management All in all, it is not possible for a geological model to describe everything everywhere. The ability to navigate between scales depending on the considered physical processes is essential to represent geological details and uncertainties where they matter. The application of geological modeling methods in a local way and a sound homogeneization and de-homogeneization are essential components for the future of joint geophysical inversion.

We want to conclude this chapter with a remark by [START_REF] Boschetti | Interactive inversion in geosciences[END_REF] on the difficulty in communication between geologists and geophysicists: "Geologists may perceive geophysicists as being lost in an abstraction far removed from real geology; geophysicists often see geologists as hopelessly resistant to mathematical rigor". With the overview presented here, we hope to contribute to a better understanding of both sides in the aim to investigate the joint object of interest: the Earth below our feet.

In the following description and review of concepts, we will attempt to follow this logic as closely as possible, to delineate and explain the different effects of parameters, observations, and the hypothesis space in the context of different geological modeling and uncertainty quantification processes.

In the frame of this review, we do not provide an exhaustive overview of these distributions, but refer to relevant literature. Potential prior distributions for parameters θ as input to geological models have briefly been discussed in Sec. 2.2.2. Several aspects of geological likelihood functions are discussed in Sec. 4.3.1. An additional point that requires consideration is that we are using a model-based Bayesian approach. In all but the simplest geological (and, equally, geophysical) modeling and inversion scenarios, a non-trivial forward model is involved. In the case of geological modeling, it can, for example, be any of the model types described in section 3. When we consider the comparison to geophysical data, then geophysical forward simulations become an additional essential aspect. The requirements for combined geological and geophysical investigations using probabilistic methods are therefore somewhat different to many of the basic statistical models that are often used in introductory texts, as our generative models are generally much more complex. In addition, it is often not completely clear which role measurements and observations take (they are equally used as input parameters, but also as data in likelihood functions), and it is often impossible to even closely characterise the hypothesis space H, which would, for example, consist of the set of all possible model scenarios.

Following [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF], we therefore also consider the possibility to compare alternative assumptions or model scenarios H, also on the grounds of Bayes' theorem: P (H|D, I) = P (D|H, I) P (H|I) P (D|I)

where I denotes any additional unquestioned assumptions (e.g. continuity in space, fundamental geological and/or physical laws, etc.). In our context, different hypotheses can also be related to parameters of the mathematical model (also referred to as metaparameters or hyperparameters), see also the description in Sec.

For the common cases where decisions are the final aim a geological modeling investigation, it is directly obvious that the possibility to take uncertainties about assumptions and model states into account (Eq. 23), in combination with the adjusted parameters θ given observed data D (Eq. 21), can provide important insights. Constraint Defines a domain of non-admissible parameter values in optimization theory; more generally, describes a piece of geological information translated in numerical terms.

A.2. Glossary

Figure 2 :

 2 Figure 2: Complexity from the viewpoint of structural geological modeling: different geological settings and possible geomodeling approaches

Figure 3 :

 3 Figure 3: Two viewpoints on subsurface interface and property detection; blue boxes: measurements and observations, green boxes: processing and modeling steps, red boxes: intermediate representations, yellow boxes: fundamental principles and knowledge.The circled numbers and letters refer to important aspects and links, explained in more detail in the text.

Figure 4 :

 4 Figure 4: Typical Earth data used in geomodeling. All data relate to a range of spatial resolution (vertical axis), and imply varying degrees of processing and interpretation (horizontal axis). Only a few data types (in italic) provide exhaustive spatial coverage in three dimensions. Colors highlight the type of data for various levels of processing and interpretation.

  to represent geological setting, types of geological structures, complexity, scale, 2-D -2.5-D -3-D, etc. (b) Structure of mathematical model Structure of the chosen mathematical model: number and type of geological intefaces, fault network, etc. Parameters of the mathematical interpolation function: variogram model, RBF parameters, etc. (c) Parameters of mathematical model (d) Input points for interpolation Position of input points to interpolation function: interface position, orientation measurements, etc.

Figure 5 :

 5 Figure 5: Sources and types of uncertainty related to different modeling steps

  Basis functions aligned to boundaries (c) Linear tetrahedral mesh for better fault approximation (d) Tetrahedral mesh and stratigraphic coordinates

Figure 6 :

 6 Figure 6: Various representations of the p-wave velocity in a geological model. These models approximate the same structural geometry and show different seismic velocity models interpolated in space according to different representations. The right hand side column shows some one dimensional synthetic velocity profile (in red) and the model approximations (in black) corresponding to the left hand side models. See text for more details.

Figure 7 :

 7 Figure 7: Two possible ways to truncate a geological domain.

A

  Onlap unit (unconformable to base), turbidite lobes C Progradational unit (unconformable to base), deltaic deposits B Aggradational unit (conformable), carbonate buildup D Eroded unit (unconformable to top), alluvial deposits E Fractured basement, plutonic rocks S AB Baselap surface (e.g.

Figure 9 :

 9 Figure 9: Cross-sections of various numerical representations of geological structures. Red lines correspond to fault surfaces and black lines correspond to five stratigraphic surfaces (A, B, C, D and E). a: Reference continuous model. b: Pillar-based or extrusion methods (Sec. 3.2.1) use two-dimensional depth or thickness maps computed along a direction field tangent to the main faults. c: Boundary representations (Sect. 3.2.2) using piecewise linear surfaces with three dimensional coordinates to define sealed geological volumes. d: Implicit representations (Sec. 3.2.3) treat some surfaces -here, stratigraphic horizons-as iso-values of a three-dimensional scalar field (figured with the white dashed lines and the corresponding values).

Figure 10 :

 10 Figure10: The fundamental concept behind the majority of methods for uncertainty quantification in geological models is to go from one deterministic representation (a) to multiple realizations, for example multiple interfaces and faults (b). These realizations are eventually discretized (c). Uncertainty is then reduced by confronting these geological models to ancillary geological likelihood criteria (d) or geophysical observations (e).

  Figure 11: Typical elements of uncertainty in generation of model realizations

Figure 12 :

 12 Figure 12: Perturbing an implicit stratigraphic model: (a) Initial scalar field s(x) representing relative geological time; (b) Random field (x) used to perturb the stratigraphy; (c) Perturbed scalar field s (x) = s(x)+ (x); (d) View of two perturbed horizons; faults are displayed in semi-transparent red and the perturbed scalar field is displayed on the background planes.

Figure 13 :

 13 Figure 13: Visual representation of multiple model realizations: (a) model realizations in separate adjacent subfigures (Suzuki et al., 2008); (b) Simultaneous representation of multiple interfaces in a 2-D section (De la Varga et al., 2018), and (c) in a 3-D representation (Mallet and Tertois, 2010);

Figure 14 :

 14 Figure 14: Uncertainty visualization on the basis of a scalar quantity of uncertainty: (a) standard deviation of thickness in map view draped on mean surface (Wellmann et al., 2010a); (b) Ghost halo as a volumetric visualization of surface uncertainty (Viard, pers. comm.) (c) Calculated cell-based probability for a specific geological unit, and (d) cell entropy as combined measure of uncertainty (Wellmann and Regenauer-Lieb, 2012).

Figure 15 :

 15 Figure 15: Overview of methods to analyze uncertainties in structural models for different levels of geomodel complexity and different classes of uncertainties

  2.1 and Fig 5. The relevance of these different hypotheses has recently been highlighted for the case of seismic inversions[START_REF] Thore | Uncertainty in seismic inversion: What really matters?[END_REF], and we can equally consider them as highly relevant for geological modeling studies.Furthermore, we can use the estimated uncertainty about different models to then make predictions t, which are based on the observed Data D, the different model assumptions H and the background information I taking the probabilities of the assumptions or models into account (MacKay, 2003): P (t|D, I) = H P (t|D, H, I) P (H|D, I)

  2.5-D Representation of a three-dimensional model by extrusion (or sweeping) of a surface. This includes map-based representations of horizons and methods based on the extrapolation of a cross-section.Boundary representation Definition of a volume by one or several watertight surfaces.

Table 1 :

 1 Summary table comparing the various modeling approaches

	Method		Extrusion	Explicit 3D surface Implicit surface
	Representative-	Simple "2.5D" do-		High
	ness		mains
	Memory	re-	one 2D scalar field	three 2D scalar fields one 3D scalar field
	quirements
	Topology		Fixed and controlled	Fixed and controlled	Arbitrary	(emerges
			(open surface)		from interpolation)
	Point data		Yes, need projection	Yes, no projection in-
					volved
	Orientation		difficult (need projection, non linear	Easier, no projection
	data				involved
	Conformable	iterative, involves projections	direct
	surfaces
	Dealing with	Possible using Boolean operations or domain subdivision
	discontinuities

Discrete in the sense of very short, in a geological time frame.

Note that this explanation is mostly for sedimentary systems, but equal considerations and thoughts can lead to structural models of poly-deformed magmatic and metamorphic rocks (e.g.[START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF][START_REF] Maxelon | A workflow to facilitate three-dimensional geometrical modelling of complex poly-deformed geological units[END_REF], or even unconsolidated rocks and soil systems.

Similar in the level of complexity are extrusion (or "2.5-D") approaches which assume that geological volumes can be represented by extruding a surface along a direction

The term geological model is, of course, highly ambiguous and should only be used when the context and the relationship to geometric models is clear

To some extent, this is not true for all possible measurements. For example, it is possible to identify grain sizes on a core sample without knowing details about sedimentology. This is why geologists often have several terms to describe the same reality: some terms are mainly descriptive, while some others bear also an interpretation sense.

German original: "Vor der Hacke ist es duster"

As the form of the system depends on the initial surface, the least-squares DSI system is solved iteratively by updating the coefficients of Eq. 13 (and of the other data terms) between each iteration

Again, changing the mesh of a triangulated horizon surface to align on a fault surface amounts to make its basis functions conformal to the fault (Eq. 2), whereas truncating a parametric horizon surface depending on the side of the fault correspond to Eq. 3.

Note that this subdivision is not necessarily identical to the separation into regions R J as in Sec. 3.1.1, as multiple regions can belong to the same class C.

In fact, this statement has already been made by[START_REF] Kitts | Certainty and uncertainty in geology[END_REF], commenting on an early work of Mann[START_REF] Mann | Randomness in Nature[END_REF]

For a more detailed treatment of the aspect of estimating and choosing, we refer to the insightful work of[START_REF] Matheron | Estimating and Choosing: An Essay on Probability in Practice[END_REF].

Another interesting and insightful example about the "mind projection fallacy" of randomness is described in[START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF], with the common statistical example of drawing balls from an urn with replacement.

Acknowledgements

We are very thankful for the invitation by the editor Cédric Schmelzbach to present our viewpoint on the state of geological modeling and the related treatment of uncertainties in this issue. Comments and suggestions by Jürg Hauser and Eric de Kemp are greatly appreciated and helped significantly to clarify the content of this chapter. Furthermore, we owe great dept to the discussion about the presented topics with numerous colleagues and students over the last years and we are looking forward to many more discussions in the future! We would also like to thank the sponsors of the RING Consortium (http://www.ring.georessources.univ-lorraine.fr) for providing significant support to Guillaume Caumon and his team, and Emerson / Paradigm for providing the SKUA-GOCAD software used to prepare some figures.

Similar to [START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF], [START_REF] De La Varga | Structural geologic modeling as an inference problem: A Bayesian perspective[END_REF] evaluated the use of additional geological information to evaluate model realizations, but with a focus on the defintion of likelihood functions in a Bayesian inference framework (see appendix, Sec. A. 1), where this additional geological information is considered as additional data and encoded in suitable likelihood functions. This approach has subsequently been combined with an optimized implementation of the implicit cokriging interpolation [START_REF] Lajaunie | Foliation Fields and 3D Cartography in Geology: Principles of a Method Based on Potential Interpolation[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules: Recent Advances in Computational Geodynamics: Theory, Numerics and Applications[END_REF] in a probabilistic programming framework [START_REF] De La Varga | GemPy 1.0: open-source stochastic geological modeling and inversion[END_REF]. As this method provides the possibility to generate model realizations directly on basis of the geological input parameters (Fig. 11a), this combination enables a probabilistic inference of geological model parameters, in addition to the optimization of suitable model realizations. This approach also opens-up the path for a full probabilistic joint inversion of geological and geophysical data (e.g. [START_REF] Wellmann | Uncertainty estimation for a geological model of the Sandstone greenstone belt, Western Australia -insights from integrated geological and geophysical inversion in a Bayesian inference framework[END_REF].

A similar approach to combine different types of geological information was taken by [START_REF] Grose | Inversion of Structural Geology Data for Fold Geometry[END_REF]. The authors used a fold model to obtain posterior distributions for input parameters of an implicit fold modeling method [START_REF] Laurent | Implicit modeling of folds and overprinting deformation[END_REF]. In essence, the approach consists of a probabilistic optimization approach where parameters of a Fourier series describing fold limb rotation and fold axis rotation angles are inverted from structural observations. The obtained posterior distributions are then used to estimate the uncertainty in a subsequent modeling step, where spatial uncertainties are then visualised using Shannon entropy and an average angular distance to determine locations for uncertainty reduction. A similar approach, though not in a probabilistic framework, has been taken by [START_REF] Cardozo | Kinematic modeling of folding above listric propagating thrusts[END_REF] in a kinematic inverse modeling of folding above propagating listric thrusts using simulated annealing to determine a range of models that fit observed structures.

We also note that these misfit functions, described in a general term in many publications, also provide the possibility to be included in a Bayesian framework in the form of likelihood functions, which can either be formally defined on the basis of measurement uncertainties, or used in the form of summary measures using concepts of Approximate Bayesian Computation [START_REF] Beaumont | Approximate Bayesian Computation in Population Genetics[END_REF][START_REF] Marjoram | Markov chain Monte Carlo without likelihoods[END_REF] or generalized likelihood functions (e.g. [START_REF] Beven | GLUE: 20years on[END_REF]. A brief introduction into the probabilistic viewpoint on uncertainty quantification is provided in the appendix (Sec. A.1)

All these approaches show the relevance of determining suitable mesures to compare sets of model realizations and to find methods to integrate all available geological information in a modeling study. In this context, probabilistic approaches have recently shown promising results, and we expect further interesting work in this direction in the future.

Geological models and geophysical inversion

The combination of geophysical measurements and geological knowledge appears in many forms and contexts in geosciences. The fundamental problem of non-uniqueness in geophysical inverse problems is well known and has led to many developments to combine geophysical inversion with additional information since the early work of [START_REF] Backus | Numerical Applications of a Formalism for Geophysical Inverse Problems[END_REF], [START_REF] Parker | Best Bounds on Density and Detph from Gravity Data[END_REF] and [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF]. Many recent developments address this issue through the use of multiple geophysical data sets in joint geophysical inverse approaches, and these approaches are described in detail in the geophysical literature (e.g. Gallardo andMeju, 2003, 2011;[START_REF] Moorkamp | A framework for 3-D joint inversion of MT, gravity and seismic refraction data[END_REF]. We also mentioned important combinations of geological modeling and geophysical inversion in Sec. 4.1.2. Here, we consider more broadly approaches with important links between geological knowledge, encoded in 3-D geological models, and geophysical forward models and inversions.

Many authors have stated the requirement to combine geophysical and geological information. [START_REF] Fullagar | Constrained inversion of geologic surfacespushing the boundaries[END_REF] actually describe geophysical inversion as an extension of geological modeling, and [START_REF] Jessell | Towards an Integrated Inversion of Geoscientific data: what price of Geology?[END_REF] express the same view, from the viewpoint of geological modeling. Common approaches to include geological knowledge into geophysical inversions are the methods implemented in potential-field inversions. Non-uniqueness in the interpretation of gravity is already obvious from A. Appendix

A.1. The probabilistic viewpoint on uncertainty quantification

In the manuscript, we described approaches for uncertainty quantification without the direct link to a specific ontological system. However, we mentioned the use of suitable probability distributions in Sec. 2.2.2 and the application of these distributions in Monte Carlo simulations (Sec. 4.1) to estimate uncertainty propagation. Here, we briefly outline the Bayesian viewpoint for inverse approaches in the context of geological modeling.

In the Bayesian viewpoint, degrees of belief can be described using probabilities if they follow the Cox axioms [START_REF] Cox | Probability, Frequency and Reasonable Expectation[END_REF]. Probabilities can therefore be used to express assumptions, and equally to perform inferences given these assumptions. The important point is that these beliefs can then be communicated, as assumptions need then to be clearly expressed. Due to the rules of probability, we can then at least ensure that the same conclusions are drawn on the basis of the same assumptions (e.g. [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF].

This viewpoint of considering beliefs as probabilities is also known as the subjective interpretation of probability (e.g. [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF], but it has to be noted that this is neither an unwanted, nor a limiting aspect. Many authors even state that this way of thinking is essential in the context of geoscientific investigations (e.g. [START_REF] Curtis | The science of subjectivity[END_REF], and we share this viewpoint at least to the extent that we consider the aspect of geological knowlege as a central part in the context of geological modeling and geophysical inversions (see Fig. 3), and the relevance of subjectivity in the hermeneutic interpretation of geological reasoning (see [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF].

The concept of Bayesian statistics is described at length in many textbooks (e.g. [START_REF] Sivia | Data Analysis: A Bayesian Tutorial[END_REF][START_REF] Gelman | Bayesian data analysis[END_REF]. We use here the description in MacKay ( 2003) using Bayes' theorem to combine the unknown parameters θ with data D under the explicit consideration of the hypothesis space H:

or, written in a descriptive form:

For generality, we are here using P () equally to denote probability densities over continuous variables and probabilities of discrete variables.

Conceptual (geological) model Abstract representation of geological concepts that can be understood and discussed. Often conveys temporal and causal considerations.

Data-driven Approach that primarily trusts the data (and relies on a relatively simple conceptual model).

Explicit mode A geomodel where geological interfaces are represented as two-dimensional networks of nodes and basis functions interpolating between these nodes.

Geological interface a surface separating two geological units of different nature or age, for example: a fault, a horizon, an intrusive boundary.

Geomodel, Geological model Representation of geometric elements in the subsurface, also structural geological model.

Geologic realism Essential aspects of the complex evolution of a specific part of the Earth, considered for the purpose of a geoscientific study.

Geomodel-driven Approach that uses relatively complex conceptual models which may compensate for the lack of data.

Implicit model A geomodel where geological interfaces are represented as level sets (or iso-value surfaces) of one or several scalar fields.

Kriging A linear unbiased interpolation method which takes a model or spatial variability (the variogram) and a model of spatial average.

Mathematical (interpolation) model

A numerical framework to fill the gaps and estimate unsampled values from neighboring sample points. This involves a set of interpolation parameters α.

Sources of uncertainty

Relates to the origin of uncertainties, such as lack of knowledge, measurement error or ambiguity.

Stochastic modeling An approach sampling several random variables to produce a set of possible models (also termed realizations).

Structural geological model Representation of geometric elements in the subsurface, also geological model.

Topology In this paper, concerns the connectivity and adjacency between geological entities (rock units, interfaces and lines).

Types of uncertainty

Relates to what could reduce the uncertainties (better measures, better data, better concepts, irreducible).