Monte Carlo simulation for uncertainty estimation on structural data in implicit 3-D geological modeling , a guide for disturbance distribution selection and parameterization

Three-dimensional (3-D) geological structural modeling aims to determine geological information in a 3D space using structural data (foliations and interfaces) and topological rules as inputs. This is necessary in any project in which the properties of the subsurface matters; they express our understanding of geometries in depth. For that reason, 3-D geological models have a wide range of practical applications including but not restricted to civil engineering, the oil and gas industry, the mining industry, and water management. These models, however, are fraught with uncertainties originating from the inherent flaws of the modeling engines (working hypotheses, interpolator’s parameterization) and the inherent lack of knowledge in areas where there are no observations combined with input uncertainty (observational, conceptual and technical errors). Because 3-D geological models are often used for impactful decision-making it is critical that all 3-D geological models provide accurate estimates of uncertainty. This paper’s focus is set on the effect of structural input data measurement uncertainty propagation in implicit 3-D geological modeling. This aim is achieved using Monte Carlo simulation for uncertainty estimation (MCUE), a stochastic method which samples from predefined disturbance probability distributions that represent the uncertainty of the original input data set. MCUE is used to produce hundreds to thousands of altered unique data sets. The altered data sets are used as inputs to produce a range of plausible 3-D models. The plausible models are then combined into a single probabilistic model as a means to propagate uncertainty from the input data to the final model. In this paper, several improved methods for MCUE are proposed. The methods pertain to distribution selection for input uncertainty, sample analysis and statistical consistency of the sampled distribution. Pole vector sampling is proposed as a more rigorous alternative than dip vector sampling for planar features and the use of a Bayesian approach to disturbance distribution parameterization is suggested. The influence of incorrect disturbance distributions is discussed and propositions are made and evaluated on synthetic and realistic cases to address the sighted issues. The distribution of the errors of the observed data (i.e., scedasticity) is shown to affect the quality of prior distributions for MCUE. Results demonstrate that the proposed workflows improve the reliability of uncertainty estimation and diminish the occurrence of artifacts.


Introduction
Three-dimensional (3-D) geological models are important tools for decision-making in geoscience as they represent the current state of our knowledge regarding the architecture of the subsurface.As such they are used in various domains of application such as mining (Cammack, 2016;Dominy et al., 2002), oil and gas (Nordahl and Ringrose, 2008), infrastructure engineering (Aldiss et al., 2012), water supply management (Prada et al., 2016), geothermal power plants (Moeck, 2014), waste disposal (Ennis-King and Paterson, 2002), natural hazard management (Delgado Marchal et al., 2015), hydrogeology (Jairo, 2013) and archaeology (Vos et al., 2015).By definition, all models contain uncertainty, being simpli-fications of the natural world (Bardossy and Fodor, 2001) linked to errors about their inputs (data and working hypotheses), processing (model building) and output formatting (discretization, simplification).Reason dictates that these models should incorporate an estimate of their uncertainty as an aid to risk-aware decision-making.
Monte Carlo simulation for uncertainty propagation (MCUE) has been a widely used uncertainty propagation method in implicit 3-D geological modeling during the last decade (Wellmann and Regenauer-Lieb, 2012;Lindsay et al., 2012;Jessell et al., 2014a;de la Varga and Wellmann, 2016).A similar approach was introduced to geoscience with the generalized likelihood uncertainty estimation (GLUE; Beven and Binley, 1992), which is a non-predictive (Camacho et al., 2015) implementation of Bayesian Monte Carlo (BMC).MCUE (Fig. 1) simulates input data uncertainty propagation by producing many plausible models through perturbation of the initial input data; the output models are then merged and/or compared to estimate uncertainty.This can be achieved by replacing each original data input with a probability distribution function (PDF) thought to best represent its uncertainty called a disturbance distribution.Essentially, a disturbance distribution quantifies the degree of confidence that one has in the input data used for the modeling such as the location of a stratigraphic horizon or the dip of a fault.In the context of MCUE, uncertainty in the input data mainly arises from a number of sources of uncertainty, including but not restricted to device basic measurement error, operator error, local variability, simplification radius, miscalibration, rounding errors, (re)projection issues and external perturbations.In the case of a standard geological compass used to acquire a foliation on an outcrop, device basic measurement error refers to error in lab and under perfect conditions (this information is typically provided by the manufacturer).
operator error refers to human-related issues that affect the process of the measurement such as trembling or misinterpreting features (mistaking joints or crenulation for horizons for example).
local variability refers to the difficulty of picking up the trend of the stratigraphy appropriately because of significant variability at the scale of the outcrop (usually due to cleavage or crenulation).
simplification radius refers to the uncertainty that is introduced when several measurements made in the same area are combined into a single one.
external perturbations refer to artificial or natural phenomena that have a detrimental effect on precision and accuracy such as holding high-magnetic-mass items close to the compass (smartphone, car, metallic structures) when making a measurement or the magnetization of the outcrop itself.All these sources of uncertainty may be abstracted to individual random variables, which are all added to form a more general uncertainty variable that disturbance distributions are expected to represent.The disturbance distributions are then sampled to generate many plausible alternate models in a process called perturbation.Plausible models form a suite of 3-D geological models that are consistent with the original data set.That is, the degree of uncertainty associated with the original data set allows these models to be plausible.In layman's terms, the perturbation step is designed to simulate the effect of uncertainty by testing "what-if" scenarios.The variability in the plausible model suite is then used as a proxy for model uncertainty.Several metrics have been used to express the model uncertainty in MCUE, including information entropy (Shannon, 1948;Wellmann, 2013;Wellmann and Regenauer-Lieb, 2012), stratigraphic variability (Lindsay et al., 2012) and kriging error.The case for reliable uncertainty estimation in 3-D geological modeling has been made repeatedly and this paper aims to further improve several points of MCUE methods at the preprocessing steps (Fig. 1).More specifically, we aim to improve (i) the selection of the PDFs used to represent uncertainties related to the original data inputs and (ii) the parameterization of said PDFs.Section 2 reviews the fundamentals of MCUE methods while Sect. 3 addresses PDF selection and parameterization.Lastly, Sect. 4 expands further into the details of disturbance distribution sampling.

MCUE method
Recently developed MCUE-based techniques for uncertainty estimation in 3-D geological modeling require the user to define the disturbance distribution for each input datum, based on some form of prior knowledge.That is necessary because MCUE is a one-step analysis as opposed to a sequential one: all inputs are perturbed once and simultaneously to generate one of the possible models that will be merged or compared with the others.MCUE is vulnerable to erroneous assumptions about the disturbance distribution in terms of structure (what is the optimal type of disturbance distribution) and magnitude (the dispersion parameters) of the uncertainty of the input data.However, it is possible to post-process the results of an MCUE simulation to compare them to other forms of prior knowledge and update accordingly (Wellmann et al., 2014a).
The MCUE approach is usually applied to geometric modeling engines (Wellmann and Regenauer-Lieb, 2012;Lindsay et al., 2013;Jessell et al., 2010Jessell et al., , 2014a)), although it can be applied to dynamic or kinematic modeling engines (Wang et al., 2016;Wellmann et al., 2016).This choice is motivated by critical differences between the three approaches, at both the conceptual and practical level (Aug, 2004).More specifically, explicit geometric engines require full expert knowledge while implicit ones are based on observed field data, variographic analysis and topological constraints (Jessell et al., 2014a).Geometric modeling engines interpolate features from sparse structural data and topological assumptions (Aug et al., 2005;Jessell et al., 2014a); they require prior knowledge of topology and are computationally affordable (Lajaunie et al., 1997;Calcagno et al., 2008).Dynamic modeling engines require knowledge of initial geometry, physical properties and boundary conditions; the modeling process is computationally expensive.Kinematic modeling engines require knowledge of initial geometry and kinematic history (Jessell, 1981); the modeling process is computationally inexpensive.The implicit geometric approach is preferred for MCUE because knowledge of initial conditions is nearly im-possible to achieve, and perfect knowledge of current conditions defeats the purpose of estimating any uncertainty.
Implicit geometric modeling engines use mainly three types of inputs: interfaces (3-D points), foliations (3-D vectors) and topological relationships between geological units and faults (stratigraphic column and fault age relationships).Drill holes and other structural inputs such as fold axes and fold axial planes can also be used (Maxelon and Mancktelow, 2005).Each data input is assigned to a geological unit and the model is then built according to predefined topological rules.The implicit geometric 3-D modeling package GeoModeller distributed by Intrepid Geophysics was used as a test platform for this study.The use of this specific software is motivated by its open use of co-kriging (Appendix C), which is a robust (Matheron, 1970;Isaaks and Srivastava, 1989;Lajaunie, 1990) geostatistical interpolator to generate the models (Calcagno et al., 2008;FitzGerald et al., 2009).In addition, GeoModeller allows uncertainty to be safely propagated provided that the variogram is correct (Chilès et al., 2004;Aug, 2004) as the co-kriging interpolator then quantifies the its intrinsic uncertainty.Nevertheless, MCUE is not inherently limited by the choice of the interpolator and, therefore, may be used with any implicit modeling engine.In the next section, a series of improvements are proposed to address the disturbance distribution problem.

Distribution types and their parameters
Often, the disturbance distribution used to estimate input uncertainty is the same (same type and same parameterization) for all observations of the same nature (Wellmann et al., 2010;Wellmann and Regenauer-Lieb, 2012;Lindsay et al., 2012Lindsay et al., , 2013)).Disturbance distribution parameters are defined arbitrarily (Lindsay et al., 2012;Wellmann and Regenauer-Lieb, 2012) in most cases.Additionally, uniform distributions have been regularly used as disturbance distribution and expressed as a plus minus range over the location of interfaces (Wellmann et al., 2010;Wellmann, 2013) or the dip and dip direction (Lindsay et al., 2012(Lindsay et al., , 2013;;Jessell et al., 2014a).Here, propositions are made about the type of disturbance distributions that should be used for MCUE, how to parameterize them and associated possible pitfalls.

Standard distributions for MCUE
The structural data collected to build the model are impacted by many random sources of uncertainty (Fig. 1) such as measurement, sampling and observation errors (Bardossy and Fodor, 2001;Nearing et al., 2016).Additionally, the uncertainty tied to each measurement is considered to be independent of the others.However, that is not to say that there is no dependence over the measured values themselves.For example, dip measurements along a fault line are expected to be spatially correlated though each measurement is an indepen-dent trial in terms of its measurement error.Consequently, MCUE may sample from disturbance distributions independently from one another.Under these conditions, the central limit theorem (CLT) holds true for these data (Sivia and Skilling, 2006;Gnedenko and Kolmogorov, 1954) if the variance of each source of uncertainty is always defined.Uncertainty would then be better represented by disturbance distributions that are consistent with the CLT, namely the normal distribution for locations (Cartesian scalar data) and the von Mises-Fisher (vMF) distribution for orientations (spherical vector data; Davis, 2003).However, MCUE does not a priori forbid the use of any kind of distribution.The normal distribution is the canonical CLT distribution (i.e., the distribution towards which the sum of random variables tends) defined as where ε and σ are the arithmetic mean and standard deviation, respectively.Note that the normal distribution is conjugate to itself or to Student's t distribution depending on which parameters are known a priori.That is, a normal prior distribution gives a normal or Student posterior distribution in the Bayesian framework given that the likelihood function is normal itself.The vMF distribution (Fig. 2) is the CLT distribution for spherical data; it is the hyperspherical counterpart to the normal distribution (Fisher et al., 1987) and is used under the same general assumptions for unit vectors on the pdimensional unit hypersphere S (p−1)1 .The most important property of the vMF distribution is the axial symmetry of the data around the mean direction.The vMF distribution is also the maximum entropy distribution for spherical data and is conjugate to itself given that the likelihood function is vMF distributed (Mardia and El-Atoum, 1976).These properties make the vMF distribution suitable for uncertainty analysis of spherical data (Hornik and Grün, 2013).Sampling from the vMF distribution is described in Appendix A. The general probability density of the vMF distribution for S p−1 is expressed as follows where γ T is the transposed mean direction vector and κ is the concentration.|| || denotes the Euclidean norm.κ is analogous to the inverse of σ for the normal distribution.High κ values denote distributions with low variance (Fig. 2), ultimately leading to a p-dimensional hyperspherical Dirac distribution and κ = 0 means complete randomness (equivalent to a p-dimensional hyperspherical uniform distribution).
Bear in mind that κ impacts the shape of the vMF distribution exponentially (Fig. 2).Therefore, confidence intervals are not linearly correlated to κ.For example, the 95% half aperture confidence interval for κ = 1 is 150 where I v (κ) is the modified Bessel function of the first kind at order v, and p the dimensionality of S (p = 3 for S 2 ).

Disturbance distribution parameterization
Regardless of which type of disturbance distribution is chosen, it is inappropriate to use the same distribution with the exact same parameters for each measurement in many cases, including but not restricted to cases in which some data inputs are actually a statistic -such as the mean -that is derived from a sample instead of an actual individual occurrence (Moffat, 1988); cases in which inputs (at the same location) are samples themselves (Kolmogorov, 1950); and cases in which the magnitude of the uncertainty of measurements may be impacted by the value of the measurement itself (Moffat, 1982).Statistics derived from samples (e.g., mean, median) or the actual sample are expected to lead to less dispersed disturbance distributions compared to single observations (Patel and Read, 1996;Bewoor and Kulkarni, 2009;Bucher, 2012;Sivia and Skilling, 2006;Davis, 2003).Therefore, making multiple measurements at each location of interest is recommended and field operators should not group, dismiss or otherwise alter these kinds of data.For example, a 15 m long limestone outcrop is expected to yield numerous structural measurements that should be fed to MCUE so that disturbance distributions are parameterized more precisely.However, because structural data inputs are sparse and often scarce, a Bayesian approach to disturbance distribution parameterization is proposed.More specifically, a prior disturbance distribution is updated by measurements over a CLT compatible likelihood function to generate a predictive posterior disturbance distribution.The following demonstration applies to both the normal and the vMF distributions.
The uncertainty about an input structural datum (location or orientation) can be described by a distribution G where µ true and ϑ true are the true mean and dispersion of the population, respectively.Measured data at a single location are a n-sized sample X = {x 1 , . .., x n } of G.The disturbance distribution that should be used for MCUE must take into account prior knowledge about ϑ true and the observed data X.This is achieved through a simple application of Bayes' theorem: where µ and ϑ are the expression of prior knowledge about µ true and ϑ true , respectively.The dispersion ϑ true is expected to be a deterministic function estimated via rigorous metrological studies, the methodology of which is beyond the scope of this paper.Thus, Eq. ( 5) simplifies to The prior distribution function p (µ) expresses prior belief about µ.In this case, p (µ) is defined as Jeffreys improper prior (Sivia and Skilling, 2006) for locations to express a complete lack of knowledge about p (µ): and, for the same reason, as a uniform spherical distribution for orientations: The likelihood distribution p(X|µ) expresses the probability of observing X given µ and is obtained under the assumption of independence by computing the joint density function for X: The posterior predictive distribution p x|X expresses the theoretical distribution of a new observation given X regardless of µ; it is the target disturbance distribution to be sampled for MCUE and is given by where x is the element to be sampled.To illustrate the rationale for the usage of Eq. ( 10), one may consider the following example.
When a single measurement is made at a specific location, the sample size is 1 (X = {x 1 }) and the observed average is equivalent to the measurement itself (µ ≡ x 1 ).Assuming that the dispersion function is deterministic then ϑ true is known.Obviously, the posterior distribution will be p (µ = x 1 |X = x 1 ) = p (x|x 1 , ϑ true ).One might think that p (µ = x 1 |X = x 1 ) is the target disturbance distribution that should be used for the perturbation step.However, doing so would lead to systematic underestimation of the effect of ϑ true because p (µ|X) only quantifies the knowledge of µ in regard to x 1 .Indeed p (µ|X) tells about all the possible G values that x 1 might be sampled from; it is a distribution of the average and, ultimately, how far µ is from µ true is (and will remain) unknown (Fig. 3).That is, sampling directly from p (µ|X) ignores the fact that µ = (µ − µ true ) 2 is unknown.Such a procedure would introduce an undesired unknown bias to the perturbation step.To account for this, p (µ|X) is compounded to itself to obtain p (x|p (µ|X) , ϑ true ), which is equivalent to (10) and practically amounts to a double sampling of p (µ|X).Consequently, regardless of the quality of the prior knowledge about ϑ true , sampling from the posterior predictive distribution is better than sampling from the posterior distribution.
For a normal distribution, the posterior predictive distribution p N is where µ 0 and σ are the sample mean and prior standard deviation, respectively.Note that µ 0 is data contribution only while σ is prior reliable knowledge obtained via metrological analysis.For the vMF distribution the posterior predictive distribution is given by where µ 0 and κ and are the mean direction vector of the sample and prior concentration (Appendix B) of the observed sample.Note that Eqs.11 and 12 can be applied to data recorded as a mean value provided that the size of the sample is known.In Eq. ( 12) µ 0 is given by where where φ is the colatitude, θ is the longitude and R is the resultant length of the observed sample.In Eq. ( 12) R is given by From Eqs. ( 11) and ( 12) it appears that sampling from prior distributions directly will lead to systematic underestimation of dispersion because σ 2 ≥ σ 2 n and κ ≤ κR 1+R .In turn, this bias will narrow the range of models explored by MCUE and will make the final results look less uncertain than they should be.This is highly important because disturbance distribution sampling in MCUE is a one-step process and incorrect disturbance distributions are not to be updated or refined at any point.Therefore, accurate parameterization of a disturbance distribution at the beginning of the process is crucial to ensure accurate sampling.Bayesian schemes exist to validate models based on some external observations or assumptions (Fig. 1) that are used to build likelihood functions (de la Varga and Wellmann, 2016).However, these schemes are known to not yield good results when incorrect informative priors are used (Freni and Mannina, 2010;Morita et al., 2010).Incorrect informative priors have low dispersion (high precision, "self-confident") and high bias (low accuracy, "off target").This results in an inability of standard Bayesian schemes to update these priors regardless of the strength of the evidence.For example, a posterior distribution extracted from a single foliation measured on an outcrop cannot be used as a disturbance distribution.Indeed, it is too narrow and may be heavily biased (Fig. 3).To avoid this detrimental effect, one should instead sample from the posterior predictive distributions (Eqs. 11,12) for more accurate results about uncertainty.

Measurement scedasticity
Scedasticity is defined as the distribution of the error about measured or estimated elements of a random variable of interest (Levenbach, 1973).It expresses the relationship between the measured values and their uncertainty.In the case in which uncertainty is constant across the variable space the variable is homoscedastic (Fig. 4a); such behavior is commonly assumed in gravity surveys (Middlemiss et al., 2016).When uncertainty is not constant throughout the variable space, the variable is called heteroscedastic (Fig. 4b,  c).Note that heteroscedastic cases include both structured (Fig. 4b) and unstructured (Fig. 4c) relationships between the measured values and their respective errors.Structured heteroscedastic variables show a clear relationship (e.g., correlation, cyclicality) between the variable and its uncertainty while unstructured ones do not.Structured heteroscedastic behavior is observed electrical resistivity tomography (Perrone et al., 2014), magnetotellurics (Thiel et al., 2016;Rawat et al., 2014), airborne gravity and magnetics (Kamm et al., 2015), and controlled-source electromagnetic (Myer et al., 2011) surveys.It is usually possible to transform a structured heteroscedastic variable to a space where it becomes homoscedastic (commonly the log space), perform analysis and transform back to the original space.Unstructured heteroscedastic behavior is common in seismic surveys and impacts inversions (Kragh and Christie, 2002;Quirein et al., 2000;Eiken et al., 2005).The heteroscedastic case essentially allows for any level of correlation between the measured values and their uncertainty or error to be possible (Fig. 5).
The failure to account for scedasticity often implies the assumption of homoscedasticity as this assumption allows for a wider range of statistical methods to be applied.With heteroscedastic data, the results of methods that depend on the assumption of homoscedasticity, such as least-squares methods (Fig. 4), give results of much decreased quality (Eubank and Thomas, 1993) and this may lead to the validation of incorrect hypotheses.Scedasticity analysis from raw data without prior knowledge is challenging (Zheng et al., 2012) and this topic of research is still being investigated (Dosne et al., 2016).If there is no option for an appropriate transform, it is advisable to perform an empirical analysis of scedasticity beforehand.This is usually achieved through experimental assessment of uncertainty under various conditions (metrological study) of measurement and over the entire range of measured values (Allmendinger et al., 2017;Cawood et al., 2017;Novakova and Pavlis, 2017).The results of such analysis can then be used to define the prior dispersion (ϑ in 5) more accurately as a function of the measurement instead of a constant.
Each data input is expected to carry its own parameterization for disturbance distribution depending on the nature of the input (single measurement, sample, central statistic).Additionally, the parameters of the disturbance distributions are better defined when scedasticity is accounted for.It is worth mentioning that both the normal distribution and the von Mises-Fisher distribution have a complete range of analytical or approximated solutions for both posterior and posterior predictive distributions (Rodrigues et al., 2000;Bagchi and Figure 5. Distribution of errors for the cases described in Fig. 3. Homoscedastic case shows constant uncertainty and no relationship of uncertainty to the data.The structured heteroscedastic case has a linear relationship of uncertainty to the data.The unstructured heteroscedastic case demonstrates no obvious relationship of uncertainty to the data and is not constant.Guttman, 1988;Bagchi, 1987).In the next section, disturbance distribution sampling for spherical data (orientations) is discussed.

Sampling of orientation data for planar features
In the geoscience, the orientation of planar features such as faults and bedding is described by foliations.These foliations can be recorded in the form of dip vectors using the dip and dip-direction system.This system is equivalent to a reversed right-hand rule spherical coordinates system.The following covers sampling strategies for such spherical data and demonstrates their impact on MCUE results.

Artificial heteroscedasticity
Recent research using MCUE (Lindsay et al., 2012(Lindsay et al., , 2013;;Wellmann and Regenauer-Lieb, 2012; de la Varga and Wellmann, 2016) uses dip and dip-direction values independently (as two scalars) from one another.The dip and dip-direction system is a practical standard for field operators to record and make sense of orientation data.However, it is highly inappropriate for statistics.Geoscientists generally perform statistical analysis on stereographic projections of the dip vectors to the planes.Because stereographic projection involves the transform of dip vectors to pole vectors (normal vector to the plane), it gives a sound representation of the underlying prior uncertainty distribution.The pole transform step is essential to avoid variance distortion (Fisher et al., 1987) as shown in Fig. 5.The distortion will increase as the dip of the plane diverges from π 4 and is maximal for degenerate cases of the dip-direction system such as horizontal and vertical planes (Fig. 6).In the case of an uncertain horizontal plane, dip vectors distribute themselves directly below and about the equator of S 2 , following a girdle-like distribution (Fig. 7a, b).Consequently, the resultant length is null and the spherical variance S 2 s (Eq.16) equals unity as the barycenter of all dip vectors is located at the center of S 2 .
Naive interpretation of S 2 s may lead one to misinterpret uncertainty to be infinite S 2 s = 1 and the plane's orientation to be uniformly random where it might be very well constrained in reality.That is so because S 2 s is a scalar quantity used to represent dispersion for samples of spherical unit vectors.Therefore, it is expected that S 2 s is ambiguous in some cases.The opposite effect occurs for (sub)vertical planes where S 2 s will appear to be lower than expected.In Fig. 7, the effect of dip vector sampling and pole vector sampling is demonstrated for theoretical cases.Here, the blue clusters are the direct result of pole vector sampling and always describe the plane's behavior accurately in terms of pole vectors.They have constant point density and are isotropic; parameterization is easy and reliable for distributions such as von Mises-Fisher (Fig. 7b, d, f) or bounded uniform (Fig. 7a, c, e).Green clusters are the result of pole vector sampling (blue) converted back to dip vector and they describe the plane's behavior accurately in terms of dip vectors.These clusters have varying shapes and may not be modeled satisfactorily by any existing spherical distribution for all possible cases.Red clusters are the direct result of dip vector sampling and fail to describe the behavior of the plane accurately.Therefore, accurate sampling based on dip vectors (green) is nearly impossible to achieve without increasing the number of parameters of the distributions to take into account the aforementioned effects (i.e., adding a set of functions to compensate for scedasticity errors as well as boundary effects).For example, in a scenario in which dip vectors are used directly to estimate a sample's spherical variance or sample over a disturbance distribution, one may attempt to define separate values for dispersion of dip and dip direction (Lindsay et al., 2012) in order to compensate for scedastic incoherence.A horizontal plane's uncertainty is then obtained by setting circular variance as null over the dip direction and as any real positive value over the dip.In addition, some form of boundary control or polarity correction of the dip is necessary to remove incorrect occurrences.Conversely, poles to planes carry information about polarity implicitly (e.g., the Cartesian pole of a horizontal plane is [0, 0, 1] while its reversed counterpart is [0, 0, −1]).Note that this method still does not solve the scedasticity issue entirely, especially for high uncertainty values about the dip and vertical dips.Similarly, if dip vectors are used directly, near-vertical planes display uniform random behavior of dip direction (Fig. 7e, f) instead of the expected "bow tie" pattern.This pattern is impossible to model accurately using CLT spherical distributions as they are unimodal and symmetric.
The use of distributions in MCUE makes it very sensitive to scedasticity over inputs.The uncertainty of a dip vector which is quantified by any dispersion parameter similar to S 2 s will show non-systematic heteroscedasticity because of variance distortion.A plane dipping at any angle would show increased heteroscedasticity of its uncertainty as the dispersion parameter used to parameterize the underlying distribution increases.Note that uncertain planes show increased heteroscedasticity as their dip diverges from a 45 • dip (Fig. 7c, d).Boundary effects also play a role for horizontal and vertical limit cases (Figs. 6, 7) as standard dip angles are constrained to 0 − π 2 .That of course considerably lowers the quality of any subsequent procedure that relies on accurate propagation of uncertainty as these planes are expected to have the lowest uncertainty in terms of dip direction.These impracticalities make it generally better to work on pole vectors rather than dip vectors.Pole vectors mostly eliminate the need for variance correction and allow coherent sampling over a plane's orientation (Fisher et al., 1987).The pole vector transform is widely used in structural geology (Phillips, 1960;Wallace, 1951;Lisle and Leyshon, 2004) through stereographic projection.Therefore, data collected as dip vectors using the dip and dip-direction system (green clusters, Fig. 7) must be transformed to poles (blue clusters, Fig. 7) for accurate estimation of spherical variance.Disturbance distributions should then be defined and sampled based on the pole vectors (mean of blue clusters, Fig. 7), as described in Sect.3.2, instead of the mean dip vector (mean of green clusters, Fig. 7) to avoid distortion (red clusters, Fig. 7).The sample can then be converted back to dip vectors if required.

Impact of pole vector sampling versus dip vector sampling
The impact of pole versus dip vector sampling on the results of MCUE is evaluated on a simple synthetic model and on a realistic synthetic model.The simple model is a standard symmetric graben with four horizontal units; it has been chosen for its simplicity and is commonly used as a test case (Wellmann et al., 2014b;de la Varga and Wellmann, 2016;Chilès et al., 2004) in MCUE for proof of concepts.The realistic model is a modification of a real demonstration case that is part of the GeoModeller package based on a location near Mansfield, Victoria, Australia.It features a Carboniferous sedimentary basin oriented NW-SE that is in a faulted contact (Mansfield Fault) on its SW edge to a Silurian-Devonian set of older, folded basins.Outcropping units are almost all of the siliceous detritic type ranging from mildly deformed sandstones to siltstones and shales; the basement is made of Ordovician-Cambrian serpentinized sandstone.The original data for the Mansfield model were not altered in any way, instead data based on the Mansfield geological map (Cayley et al., 2006) geophysical map (Haydon et al., 2006) and airborne geophysical survey (Wynne and Bacchin, 2009;Richardson, 2003) were added to refine it.The graben model is built using orientations and interfaces only, with three interfaces and three foliations per unit and one interface and one foliation per fault (Figs. 8, 9c).The Mansfield model is built with 281 interface points and 176 foliations over six units and three faults (Figs. 10, 11c).For both models, perturbation is performed as described in Sect.3.For the graben model, units' interfaces are isotropically perturbed over a normal distribution with the mean centered on the original data point and a standard deviation of 25 m.The orientations of the faults are perturbed over a von Mises-Fisher distribution with the original data as the mean vector and concentration of 100 (p95 ∼ ±10 • ) following the recommended pole vector procedure described in Sect. 4.1 (Figs. 9a,11a) or the dip vector one (Figs. 9b,11b).For the Mansfield model, all interfaces and orientations (both for units and faults) are perturbed using the parameterization given for the graben model.The perturbation parameters for orientations were chosen to be compatible with metrological data.That is, values for the dispersion of the spherical disturbance distributions used for the foliations were estimated on the basis of the variability in plane measurements observed by other authors (Nelson et al., 1987;Stigsson, 2016;Allmendiger et al., 2017;Cawood et al., 2017;Novakova and Pavlis, 2017) in a variety of settings and for different types of devices.Perturbation parameters for interfaces were designed to meet observed GPS uncertainty (Jennings et al., 2010) and observed experimental interface variability in previous authors' works (Courrioux et al., 2015;Lark et al., 2014Lark et al., , 2013)).More specifically it was assumed that the observed end variability in the interfaces' locations in their models can be transposed to the presented cases.This is of course an approximation in the absence of specific metrological studies.
The influence of dip vector (Figs. 9b,11b) versus pole vector (Figs.9a, 11a) sampling of orientations is very noticeable over the output information entropy uncertainty models.Information entropy is a concept derived from Boltzmann equations (Shannon, 1948) that is used to measure chaos in categorical systems.Because of this, it is possible to use information entropy as an index of uncertainty in categorical systems.Dip vector sampling appears to add a layer of artificial "noise" on top of the uncertainty models.The noise prevents expected structures of the starting model (Figs.9c, 11c) to be easily distinguishable.In cases in which the orientation data are more vulnerable to improper sampling error (away from 45 • dips) important structures such as the nearvertical faults in the graben model (Fig. 9) or the circled areas in Fig. 11 may completely disappear.It also appears that areas where low uncertainty would be expected (orange unit in Fig. 11) are the loci of excess uncertainty.These observations support the assertion that pole vector sampling should be favored to improve uncertainty propagation in MCUE.

Discussion
Generally, CLT distributions are valid choices as prior uncertainty distributions (and disturbance distributions) because they describe the behavior of uncertainty well.However, there may be scenarios in which alternatives can offer a better solution.More specifically, the uniform or the Laplace distribution may better describe location uncertainty than the normal distribution.The uniform distribution indicates a lack of constraints as to the prior uncertainty distribution; it is a valid choice when there is little knowledge about data dispersion.The Laplace distribution is suitable if the measured data abide by the first law of errors instead of the second (Wilson, 1923).For example, to model the uncertainty on the thickness of a geological unit along a drill core, one might observe that the uncertainty of the location of the top and bottom interfaces of the unit is best represented by an expo-nential distribution.In this instance, the Laplace distribution would be a suitable option to model the thickness' uncertainty.Under similar circumstances, a spherical exponential distribution could be swapped with the vMF distribution.The Kent2 distribution is also a good candidate to describe orientation uncertainty when the pole vectors of measured orientations appear to be anisotropically distributed on S 2 (Kent and Hamelryck, 2005).
In this paper, it is explicitly assumed that the dispersion of prior uncertainty distributions is a deterministic function.Note that this does not necessarily make this function a constant and it might depend on the observed data.The dispersion function of field measurements (using a compass) of structural data would be expected to be nearly constant.Conversely, the dispersion function of interpreted measurements (using geophysics) would be expected to be dependent on the sensitivity of the intermediary method.Additionally, dispersion functions may be probabilistic as well as deterministic (Bucher, 2012).Determinism is a strong assumption when no metrological study was conducted beforehand to assess its plausibility.Such metrological studies involve experimental testing of devices and procedures in order to estimate precision, accuracy, bias, scedasticity or drift about measured data.These estimates can then be compiled into a dispersion function that can be used as an input parameter for other purposes, including prior uncertainty distributions for MCUE.Probabilistic dispersion functions apply non-negligible uncertainty to the dispersion function for prior uncertainty distributions.Uncertainty about dispersion makes the proposed workflow for disturbance distribution parameterization inadequate.Indeed, Eq. ( 5) may not be simplified into Eq.( 6) anymore and the following statements (Eqs.7 to 14) would then ignore the probabilistic nature of the dispersion function.Both the normal and the vMF distributions have analytical solutions or good approximations for such cases; the authors recommend the readers refer to relevant works (Gelman et al., 2014;Bagchi and Guttman, 1988) if required.Note that there is significant metrological work about borehole data (Nelson et al., 1987;Stigsson, 2016) as opposed to usual structural data such as foliations, fold planes, fold axes or interfaces.
Although the authors make the case for scedasticity analysis in MCUE, it is left open in this paper.Scedasticity is essentially an untouched subject in geological 3-D modeling and it was pointed out to make the geological 3-D modeling community aware of this fact and its potentially nefarious influence on MCUE outputs.However, standard metrological studies can determine scedasticity and include it in a dispersion function to be a parameter of the prior uncertainty distributions (Bewoor and Kulkarni, 2009;Bucher, 2012).The evidence brought at the theoretical and practical levels allows us to strongly advocate for the use of pole vectors over dip vectors.In fact, dip vector sampling shows poor performance away from 45 • dip planes, induces artificial heteroscedasticity and requires specific polarity indicators.This especially applies to MCUE methods in which Bayesian post-analysis is performed for the probabilistic model that results from basic propagation of uncertainty (de la Varga and Wellmann, 2016).In this respect, dip vector sampling leads to incorrect highly informative prior distributions, which is catastrophic for any Bayesian methods (Morita et al., 2010).Nonetheless, it is worth mentioning that the arguments in Sect.4.1 only apply to dip vectors of a plane and should not be extended to actual vector structural data such as fold axes or lineations.This is because these data represent linear features (lineations, fold axes, other planar features intersections) for which the concept of a pole does not apply.
Good prior knowledge about input uncertainty is critical to the propagation of uncertainty in general.This, in turn, makes metrological work mandatory to any form of modeling that relies on actual measured data.Note that it is acceptable to use preexisting metrological studies to define the priors (Allmendinger et al., 2017;Cawood et al., 2017;Novakova and Pavlis, 2017) provided that the measurement device and procedure used are similar to those of the studies.To gather multiple observations per site is strongly recommended as this practice sharply increases the quality of the disturbance distributions.From a practical point of view this would require field operators to perform several measurements on the same outcrop.If that is not possible one may group measurements of clustered outcrops together provided that the scale of the modeled area compared to that of the cluster allows it.The authors recommend not grouping clusters that are spread out more than 3 orders of magnitude below the model size (e.g., for a 10 km×10 km model, clusters of radius higher than 5 m shall not be grouped).Note that more refined structural data upscaling methods have been recently proposed to address this specific issue (Carmichael and Ailleres, 2016).However, there is another major source of uncertainty that stems from the necessarily imperfect modeling engine itself.Implicit geometric modeling engines (in this case, GeoModeller) use interpolation to draw the contact surfaces of geological units.Therefore, the parameterization of the interpolator may impact results.The co-kriging interpolator (Appendix C) used in this paper relies on (uncertain) variographic analysis (Appendix C, Eq.C2) and is natively able to express its own uncertainty (Appendix C, Eq. 29).Therefore, these sources of uncertainty are expected to be propagated along the input uncertainty as a hyperparameter in Eq. ( 5).

Conclusions
Propagation of uncertainty is the process through which different kinds and sources of uncertainties about the same phenomenon are combined into a single final estimate.MCUE methods seek to achieve propagation of uncertainty using Monte Carlo-based systems in which input uncertainty is simulated through the sampling of probability distributions called a disturbance distribution.Disturbance distributions are the distributions that normally best represent the uncertainty about the input data in the context of uncertainty propagation in geological 3-D modeling.
This paper discusses the importance of disturbance distribution selection and proposes a simple procedure for better disturbance distribution parameterization and a pole-vectorbased sampling routine for spherical data (orientations) used to represent the geometry of planar features.Pole-vectorbased sampling for spherical data and Bayesian disturbance distribution parameterization are proved -either through demonstration or through experiment -to be valid and practical choices for MCUE applied to implicit 3-D geological modeling.Namely, the normal and the vMF distributions are shown to be the best candidates for disturbance distributions for location and orientation, respectively.A Bayesian approach to disturbance distribution parameterization is shown to avoid underestimation of input data dispersion, which is important as such underestimation artificially decreases the output uncertainty of the 3-D geological models.Such underestimation may give a false sense of confidence and lead to poor decision-making.Pole vector sampling is evidenced to be the best alternative because it is guaranteed not to distort the disturbance distribution shape or generate artifacts in the output uncertainty models the way dip vector sampling does.
The proposed framework and methods are compatible with previous MCUE work on 3-D geological modeling and can be easily added to existing implementations to improve their accuracy.As MCUE is applicable to all fields in which 3-D geological models are needed, so is the proposed framework.The primary domains of application are the mining and oil and gas industries at the exploration, development and production steps.In addition, numerous secondary domains of potential application are available to this work, such as civil engineering and fundamental research.
Code and data availability.Both the Mansfield and graben Ge-oModeller models (including the perturbed data sets and series of plausible models) showcased in the present study are available online openly at https://doi.org/10.5281/zenodo.848225(Pakyuz-Charrier and Intrepid Geophysics, 2017) and https://doi.org/10.5281/zenodo.854730(Pakyuz-Charrier, 2017), respectively.Instructions on how to use the GeoModeller API can be found at http://www.intrepid-geophysics.com/ig/index.php?page=geomodeller-api (Intrepid Geophysics, 2018).where d and a are the lag distance and the range, respectively.The theoretical semi-variogram is fit to an empirical semi-variogram (Matheron, 1970).In practical cases, the empirical to theoretical semi-variogram fit is never perfect and is mostly parametric.The probability that the potential value estimated at a point x is comprised between t and t (Aug, 2004;Chilès et al., 2004) is given by where σ CK (x) is the co-kriging standard deviation and N is the normal cumulative distribution function.Equation (C3) can be used as an uncertainty estimator for the interpolator (Eq.C1) if and only if both t and t can be defined adequately as equivalent to the top and bottom of a formation.If it happens to be the case these probabilities can be combined to the final uncertainty model at the merging step (Fig. 1).However, such definition is not always possible.Note that

Figure 3 .
Figure 3.Effect of bias over predictive posterior and posterior distributions.

Figure 4 .
Figure 4. Synthetic examples of different levels of scedasticity of measurements of the same variable.(a) Homoscedastic case, (b) structured heteroscedastic case and (c) unstructured heteroscedastic case.Note how the least-square polynomial residual score (R 2 ) is heavily impacted by scedasticity.

Figure 6 .
Figure 6.Distortion of the maximum likelihood estimation (MLE) of concentration or spherical variance of 100 spherical unit vector samples of a size of 1000 individuals drawn from a von Mises-Fisher distribution with κ = 100.A pole-based estimate is always consistent with the data while dip-based ones either over-or underestimate it.

Figure 7 .
Figure 7. Effect of sampling over dip vectors or pole vectors on bounded uniform spherical distribution at a range = 10 • (a, c, e) and von Mises-Fisher distribution at κ = 100.0(b, d, f) for uncertain horizontal planes (a, b), 45 • dip planes (c, d) and vertical planes (e, f).Correct (pole perturbed) dip vectors are green, incorrect (dip perturbed) dip vectors are red and blue vectors are the poles.See Sect.4.1 for details.

Figure 8 .
Figure 8. Structural data for the graben model and modeled surfaces for units and faults.Spheres represent interfaces and cones represent pole vectors.

Figure 9 .
Figure 9.Effect of pole (a) versus dip (b) perturbation for a graben model (c); orientations are perturbed over a vMF distribution with κ = 100.

Figure 10 .
Figure 10.Structural data for the Mansfield model and modeled surfaces for units and faults.Spheres are interfaces and cones are orientations.

Figure 11 .
Figure 11.Effect of pole (a) versus dip (b) perturbation on a cross section of the Mansfield model (c); orientations are perturbed over a von Mises-Fisher distribution with κ = 100.