loopUI-0.1: uncertainty indicators to support needs and practices in 3D geological modelling uncertainty quantiﬁcation

. To support the needs of practitioners regarding 3D geological modelling and uncertainty quantiﬁcation in the ﬁeld, in particular from the mining industry, we propose a Python package called loopUI-0.1 that provides a set of local and global indicators to measure uncertainty and features dissimilarities among an ensemble of voxet models. Results are presented of a survey launched among practitioners in the mineral industry, enquiring about their modelling and uncertainty quantiﬁcation practice and needs. It reveals that practitioners acknowledge the importance of uncertainty quantiﬁcation even if they do not 5 perform it. Four main factors preventing practitioners to perform uncertainty quantiﬁcation were identiﬁed: lack of data uncertainty quantiﬁcation, (computing) time requirement to generate one model, poor tracking of assumptions and interpretations, relative complexity of uncertainty quantiﬁcation. The paper reviews and proposes solutions to alleviate these issues. Elements of an answer to these problems are already provided in the special issue hosting this paper and more are expected to come.


Introduction
One objective of researchers who develop open-source 3D geological modelling algorithms (Loop, 2019;de la Varga et al., 2019) is to make them Findable, Accessible, Interoperable and Reusable (FAIR) for practitioners.As for any software, these algorithms should satisfy the needs and expectations of users (Franke and Von Hippel, 2003;Kujala, 2008).Thus, new developments should rely on a good understanding of modelling purposes, processes and limitations, following a philosophy of continuous improvement.In a general context, this becomes even more important given the increasing number of open-source algorithms in the fields of earth and planetary sciences (see Figure 1).However, to the best of our knowledge, the needs and uses of 3D geological modelling practitioners with respect to uncertainty quantification are only partially described in the literature, as it constitutes an emerging field that only recently gained traction in both academia and industry.
An essential purpose of modelling is to support decision makers by offering a simplified representation of nature that also provides a corresponding quantitative assessment of uncertainty, communicating what we know, what remains unknown and what is ambiguous (Ferré, 2017).Uncertainty quantification is essential, because it allows us to mitigate predictive uncertainty (Jessell et al., 2018) by expanding our knowledge, rejecting hypotheses (Wilcox, 2011) or falsifying scenarios (Raftery, 1993).
Questions related to 3D geological modelling and uncertainty quantification are not just limited to the minerals industry, but also concern the fields of CO2 sequestration (Mo et al., 2019), petroleum (Scheidt and Caers, 2009) and geothermal (Witter et al., 2019) energy resources as well as hydrogeology (Linde et al., 2017) or civil engineering in urban environments (Osenbrück et al., 2007;Tubau et al., 2017).Here, we are interested in the uses and practices of the minerals industry, that is dealing with both sedimentary basin and hard-rock and/or cratonic settings across regional to mine scales.
The three main pillars of uncertainty quantification are the characterization ::::::::::::: characterisation : of uncertainty sources, their propagation and mitigation throughout the modelling workflow (see Figure 2).The different sources of uncertainty, often overlooked, are related to measurement errors, interpretations, assumptions, modelling approximations and limited knowledge (sample size or unknown process).Measurement or data errors can be estimated by repetitive :::::::::: independent sampling or from instrument characteristics; they can be propagated through the modelling workflow by Monte Carlo data perturbation (Wellmann and Regenauer-Lieb, 2012;Lindsay et al., 2012;Pakyuz-Charrier et al., 2018).Combined with expert knowledge, initial dataset can be used to shape some assumptions and define plausible conceptual models; but despite the importance of conceptual uncertainty on predictions (Pirot et al., 2015), it is too often limited to the definition of a unique scenario (Ferré, 2017).From a perspective on algorithms, some assumptions such as how to set parameter ranges are needed and this can greatly impact the definition of geological parameters (Lajaunie et al., 1997) prior to running predictive numerical simulations.
While the main objective of uncertainty quantification and data integration might be to improve the confidence level of predictions for decision making, it usually involves the generation of model ensembles via Monte Carlo algorithm and it is rarely a straightforward step.Indeed, because of the high dimensionality and non-linearity of earth processes or the lack of data, history matching might prove difficult to achieve and predicted outcomes might present multiple modes (e.g.Suzuki and Caers, 2008;Sambridge, 2014;Pirot et al., 2017).In such cases, geological uncertainty analysis allows to improve our understanding of geological model (dis)similarities and how specific or shared features can be related to upstream parameters and downstream predictions.as voxel-wise entropy or cardinality (Lindsay et al., 2012) :::::::::::::::::::::::::::::::::::: (Lindsay et al., 2012;Schweizer et al., 2017), computed over an ensemble of geological voxets will inform about property field variability at specific locations (voxels) of the model mesh.Global indicators or summary metrics might be useful to identify how the statistics of specific patterns (e.g.fault or fracture network density, anisotropy, connectivity, etc.) evolve between different models and might also be a way to perform model or scenario selection (e.g.Pirot et al., 2019) or to reduce the dimensionality of the sampling space, from a high dimensional geological space to a low dimensional latent space (Lochbühler et al., 2013).Though some indicators have been used or developed for specific studies or softwares (Li et al., 2014), to the best of our knowledge, no independent uncertainty analysis tool applicable to both discrete and continuous property fields, combining local and global indicators is available to practitioners.
To investigate the uses and practices of the minerals industry regarding modelling and uncertainty quantification, we recently conducted a survey among numerical modelling practitioners from industry, government and academia in the sector of exploration and production of economic minerals.In this paper, first, we present the main results and interpretations from this survey, which questions are listed in Appendix A. Second, to answer :::: some :: of : these needs, we propose a set of indicators to quantify geological uncertainty over an ensemble of geological models, :::::::::: represented :: as :::::::::: regular-grid ::::::: voxets, ::: and : characterized by lithological units :: in ::: the :::::: discrete ::::::: domain and their underlying scalar-field derived from implicit modelling :: in ::: the ::::::::: continuous :::::: domain.
The various indicators are illustrated with a synthetic case derived from a simplified Precambrian dataset of the Hamersley, Western Australia.The Python code called loopUI-0.1 (Pirot, 2021) and the notebooks used to compute and illustrate these indicators are available at https://doi.org/10.5281/zenodo.5656151.5656151 2 Survey

Material and method
The survey was designed to be concise to encourage participation but with several open-ended questions to maximize :::::::: maximise our chance to learn about different uses and practices, as well as to minimize ::::::: minimise : induced bias whenever possible.

Results
This section summarises the answers provided by the survey respondents.Table 1 provides the general statistics of answers to the survey.The high answer rate (mostly over 80%) indicates that questions are meaningful for the respondents and indicates that the reader can have confidence in the presented results.The most answered question is Q1, about the modelling scale.
The least answered question is Q9, about data integration and upscaling .The second least answered question is Q8, about the modelling workflow used.All other questions have an answer rate above 80% on a global average.Note that global number of survey answers is smaller than the sum of each survey answers grouped by scale, as some survey answers were returned only once for multiple scales.
Due to existing overlap of collected answers on different questions, the results presented hereafter summarize and group the collected answers by theme (Output or objectives, Input data, Current modelling workflow and Limitations), as outlined in the survey (see Appendix).Each theme is treated by modelling scale (see Figure 3) when it involves different answers.
Q3-4.The main modelling objectives depend on the scale of investigation.At the largest scales, investigation and modelling are particularly useful to assist exploration and prospectivity mapping.At the Greenfields or Regional scale (>10km dimension, ˜1km resolution) , the main objective is to obtain a contextual and conceptual understanding of the regional geology, in particular regarding stratigraphy, topology and geochronology.At the Brownfields scale (1-10km dimension, ˜100m resolution), modelling aims to estimate deep structure beyond the depths reached by drilling, such as depth of a particular interface or the delineation of potential mineral system components (e.g.fluid-alteration pathways and ore deposition environments).
At the Mine scale (<1km dimension, ˜10m resolution), the objectives are to assist with near mine exploration, resource estimation, ore body localization, drill targeting, operation scheduling and efficient mining, by characterizing local structures and geometries of stratigraphy and mineralization.In addition to fulfilling these various objectives, 3D models are useful to improve hydrogeological characterization ::::::::::::: characterisation, to identify critical zones where more knowledge has to be gained , to allow for comparison between datasets and test internal consistency, scenarios and last but not least, for visualization and communication.
Q5. Data types used as input for geoscientific modelling are similar and complementary across scales.For instance, geological mapping, geophysical data (gravity, magnetics, electro-magnetics , seismic), and geochemistry are useful at all scales even though the related measurements might inform about a regional trend.Drillhole-derived data is usually much more abundant at the Mine scale and can be useful to understand the regional context if analysed with this alternative use in mind.For example, drill logs that only record the presence or absence of mineralisation are not useful for the regional context, however those that record rock properties, lithology and structure can be extrapolated to larger scales and integrated in regional models (e.g.Lindsay et al., 2020).
Q6. Regarding current modelling practices, although the survey respondents acknowledge the importance of uncertainty quantification and propagation into modelling to estimate the confidence around predictions, only about 10% perform quantitative uncertainty characterization :::::::::::: characterisation : and propagation from input data; about 20% do some qualitative characterization ::::::::::::: characterisation and a vast majority of 70% recognize that it is ignored (see Figure 4).Q7-9.Modelling steps, assumptions or geological interpretations are not recorded in one third of the cases.In most cases, assumptions and elements of the modelling process are described in metadata or in separate reports.In very few cases, specific procedures and tools are available to keep track of those.The respondents use a variety of software, platforms or programming tools to produce 3D geological models and further data-integration.The lack of a coherent workflow introduces the potential for uncertainty, error and loss of precision when forced to translate data formats across multiple software platforms, in addition to the time lost that could otherwise be dedicated to solving the problem under question or exploring alternative hypotheses.

Uncertainty indicators for categorical or continuous property fields
The estimation of prediction confidence in geosciences relies heavily on numerical simulations and requires generation of an ensemble of models.As indicated by the survey results, an important aspect of uncertainty quantification is its representation.
It helps identifying specific features or zones of interest.Moreover, respondents estimate that such visualization tools are missing.The uncertainty indicators presented hereafter provide a way to identify zones of greater or smaller uncertainty as well as (dis-)similarities between geo-model realizations.
Geo-models can be used to convey very different discrete or continuous properties.Discrete or categorical properties such as lithological formations or classification codes should be compared carefully.Indeed, while it is straightforward to state if two values are identical or different, additional information is needed to rank dissimilarities.Continuous properties such as potential fields or physical properties (e.g.porosity, conductivity, etc.) do not present this ambiguity to compare range of values.One can note, that depending on how physical properties are assigned during modelling, their value spectrum might be discrete.Global measures rely on the computation of summary statistics or on the identification of feature characteristics, independently from their locations.The dissimilarities between summary statistics or characteristics can be estimated via appropriate metrics such as the Wasserstein distance (Vallender, 1974) or the Jensen-Shannon divergence (Dagan et al., 1997) for instance.
The resulting global measures allow comparing models with voxets or meshes of different dimensions.However, the computation of summary statistics might be more time consuming than local measures.Their main advantage is that it allows to focus on pattern similarity between models which might be particularly useful to explore the selection of alternative scenarios analysis.
To illustrate the different indicators, inspired by a Precambrian geological setting, that is a simplified dataset from the Hamersley region, Western Australia, we generated synthetic ensembles of 10 model sets for three different scenarios.Each model set is composed of a lithocode voxet describing the lithological units (categorical variable), and of its underlying scalarfield voxet (continuous variable).The underlying scalar-field is obtained by composition of the different scalar-fields for each unconformable stratigraphic group.Scenario 1 considers all synthetic input data, while scenario 2 keeps only 50% of the data within a North-South limited band and scenario 3 retains input data with a 50% probability (see Figure 5).For each scenario, the positions and orientations of the input data are perturbed to provide 10 stochastic realizations.Positions are perturbed with a Gaussian error of zero mean and 3m standard deviation.Orientations are perturbed with a von Mises Fisher error of κ = 150 (corresponding to about ±5 deg of error).Each 3D model was generated using LoopStructural (Grose et al., 2021).

Cardinality
The Cardinality of a set, in Mathematics, is the number of elements of the set.Here, we define the Cardinality for a given voxel as the number of unique elements over the ensemble of models for the corresponding voxel (Lindsay et al., 2012).Computed

Entropy
Shannon's entropy (see Eq. 1) is a specific case of the Rényi entropy when its parameter α converges to 1 (Rényi, 1961).It has been applied to geo-models since a few decades already (Journel and Deutsch, 1993;Wellmann and Regenauer-Lieb, 2012) and is a finer way than cardinality to describe uncertainty over an ensemble of models, as it takes into account the histogram 190 proportions of the unique values encountered.Let us consider the categorical or discrete variable X that represent a voxel property and assume that it can take n distinct values among an ensemble of voxets.By denoting the probability of observing the i th possible value as p i , the entropy H of X is computed as follows.

Histogram dissimilarity
Given a pair of voxets V P and V Q , we measure the dissimilarities between their histograms by computing a symmetrized and smooth version of the Kullback-Leibler divergence (Kullback and Leibler, 1951) known as the Jensen-Shannon divergence or total divergence to the average (Dagan et al., 1997).Given two random variables P and Q, the Jensen-Shannon divergence is conmputed ::::::: computed : as JSD(P ||Q) = 1 2 KLD(P ||M ) + 1 2 KLD(Q||M ), where M = P +Q 2 and KLD is the Kullback-Leibler divergence.It requires P and Q to share the same support X and can be computed for continuous or discrete variables.Here, we assume that for our pair of voxets V P and V Q , the support of our random variables P and Q respectively, is discrete and of size n (possibly n bins for discretized continuous variables).
Denoting the support of P and Q by (x i ) i=1...n , p i = Prob(P = x i ), q i = Prob(Q = x i ) and m i = pi+qi 2 , the Jensen-Shannon divergence is computed as in Eq. 2.
Given two empirical semi-variograms γ1 and γ2 (see e.g. Figure 8), we propose to use a weighted l p norm as defined in Eq. 3.

Connectivity dissimilarity
The existence of preferential flow-paths or barriers in the subsurface often has a strong impact in many geo-applications.
Their characterization :::::::::::: characterisation : can improve the management of groundwater quality, the extraction of geothermal energy, and help mitigate the environmental impact related to either the production of non-and renewable resources from the subsurface or the sequestration of carbon dioxide and waste (e.g nuclear waste).Renard and Allard ( 2013 (Meerschman et al., 2013) and in hydrogeophysical application for model selection (Pirot et al., 2019).
Let us consider a binary spatial variable X ∈ 0, 1, and a distance h.Then, the lag-distance connectivity function τ (h) is defined as the probability that two h-distant points s and s + h whose value of X = 1 are connected.For a binary voxet, two voxels are connected if a path through the face of successive neighbour voxels with the same property exists.The lag-distance connectivity function (see Figure 9) can be written as: Now let us assume that the percolation threshold p produces a binary voxet characterized by the binary spatial variable X ∈ 0, 1.The global percolation metric Γ(p) (see Figure 10) is the proportion of the pairs of voxels that are connected amongst all the pairs of voxels for which X = 1: where N (X p ) is the number of distinct connected components formed by voxels of value X = 1 and p i = n i /n p is the proportion of voxels forming the i th distinct connected component, n i being the size in voxels of the i th connected component and 250 n p being the total number of voxels of value X = 1 in the voxet.Conversely, for the complementary set of voxels for which X = 0, we can compute Γ(p) c (see Figure 10).One can note that the two connectivity metrics are related as and Allard, 2013).
We propose two measures of connectivity dissimilarity between voxets, based either on τ (h) or Γ(p) and Γ(p) c .Let us denote by N c the number of considered class :::::: classes : of values.Note that for a categorical variable voxet, the classes are defined by the category values, while for continuous variable voxet, they can be obtained by thresholding with N c percolation thresholds p.Let us consider N lag the number of lags (defined similarly as for the experimental semi-variogram in the previous subsection), l p = 2 the distance norm and 1 τ = 1 − 1 Γ the indicator allowing to choose between τ or Γ connectivity.The we can compute the connectivity dissimilarity between two voxets as follows in Eq. 4:

Multiple-point histogram dissimilarity
Multiple-point histograms (MPH Boisvert et al., 2010) are based on pattern recognition and have been primarily used in the field of geostatistics (Guardiano and Srivastava, 1993) to quantify the quality of multiple-point statistics simulations.Patterns are delimited by a search window whose dimensionality matches the one of the dataset.One can count unique patterns, however the number of unique patterns might be relatively important, in particular for continuous property fields.In that case, it might require to restrain the analysis to the most frequent patterns (Meerschman et al., 2013).An alternative is to base the analysis on pattern cluster representatives (see Figure 11).Here, using an l 2 norm distance between patterns and k-means clustering (Pedregosa et al., 2011, scikit-learn implementation), we classify all patterns into N c = 10 clusters.Each cluster centroid or barycentre defines its representative.
In addition, voxets can be easily upscaled which allows MPH analysis of potentially large scale features with a small search window at high level of upscaling.Note that a given level l of upscaling, the size of the dataset is divided by 2 l along each dimension.Here, to avoid property values smoothing, we perform a stochastic upscaling, i.e. in a 2D case, the upscaled value of a 2×2 subset of pixels is achieved by a uniform random draw among the values of the 4 pixels.Cluster pattern identification is performed at the initial resolution level (l = 0) and at all possible upscaled levels.
For a given upscaling level, once k-means clustering of patterns has been performed on two voxets or datasets, distances between cluster representatives of two images can be computed: , where C i 1 is the i th cluster representative for voxet 1, C j 2 is the j th cluster representative for voxet 2 and w denotes the index of the windowsearch elements.
The cluster representatives between two datasets are paired by similarity (smallest distance), and re-orderd such that ∀i, To account for cluster size differences, the distance between paired cluster representative are weighted by proportion dissimilarities.It results in an MPH cluster based distance between voxets/datasets 1 & 2 defined as in Eq.5.
where p i 1 and p i 2 are the proportions of the paired clusters C i 1 and C i 2 with respects to voxets 1 & 2 respectively.One advantage of selecting cluster representative independently between two voxets is to lower computing requirements over large ensemble of voxets, performing the analysis for N v voxets instead of Nv(Nv−1) 2 pairs.However, performing k-means pattern clustering on two datasets might provide a more accurate and precise way to compute a distance between histograms with the same support of cluster representatives, allowing thus the use of Jensen-Shannon divergence for instance.One can note that we accounted for the size of the clusters, but we could also consider the density spread or concentration around cluster representatives.290

Wavelet decomposition coefficient dissimilarity
Wavelet decomposition is way to compress images.Each level of decomposition produces a series of coefficients.If computed for images to be compared, the dissimilarity of histogram of coefficients can be computed with the Jensen-Shannon divergence (Eq.2).Here, wavelet decomposition (Figure 12) is performed with the PyWavelets Python package (Lee et al., 2019) at all possible levels of decomposition and using the 'Haar' wavelet.Other wavelet could be used, however, tests have shown that 295 such dissimilarity measures are note very sensitive to the choice of the wavelet (Pirot et al., 2019).Then a wavelet-based dissimilarity measure between two voxet can be computed as in Eq. 6 by summing the Jensen-Shannon divergences computed for all pairs of approximation and decomposition coefficients at all possible levels: where C i,j 1 and C i,j 2 denote the distributions of the i th coefficients at upscaling level j for Voxet 1 and Voxet 2 respectively.

Topological dissimilarity
Thiele et al. (2016) give an overview of possible representations of the topology in the context of 3D geological modelling.
Different levels of complexity (e.g.1st or 2nd orders ...) can be used.Nonetheless, any topological indicator is a graph, that can take the form of an adjacency matrix.Thus to compute a topological distance between two 3D geological models (for instance as in Giraud et al., 2019), it seems natural to look at distances defined between graphs.Donnat and Holmes (2018) provide a comprehensive review of graph distances, used in the study of graph dynamics or temporal evolution.Though, here, in a geological context, we aim at comparing the topological diversity of an ensemble of geological models, we can use similar distances.Donnat and Holmes (2018) classify graph distances into three main categories as summarised below.
Low-scale distances capture local changes in the graph structure.The Hamming (structural) distance is the sum of absolute value of differences between two adjacency matrices, requires the same number of vertices (nodes) between the graphs -note that it is a specific case of the more general Graph Edit Distance.The Jaccard distance is defined as the difference between the union and intersection of two graphs.The Graph Edit Distance belongs to the NP-complete class of problems, and is not computed here.More info available in (Gao et al., 2010), at , or in Part IV Chapter 15 of the Encyplopediae of Distances (Deza and Deza, 2009, p. 301).Note that some packages and implementations exist to compute the graph-edit distance, but have not been tested here ( GMatch4py, graphkit-learn, other proposed heuristic ).
High level / spectral distances are global measures and reflect the smoothness of the overall graph structure changes by measuring dissimilarities in the eigenvalues of the graph Laplacian or its adjacency matrix.Some examples are the IM distance (Ipsen and Mikhailov, 2003), l p distances on eigenvalues, or the Kruglov distance on eigenvector coordinates (Shimada et al., 2016).
Meso-scale distances are a compromise or combination of low-scale and spectral distances: Hamming-IM (HIM) combination, Heat-Wave distance, polynomial distance, neighbourhood level distances, connectivity-based distances.
Here, we propose to build first order adjacency matrices (see Figure 13) from 2D or 3D voxet models.For continuous property fields, the voxet is discretize in N = 10 classes of values defined by N equi-percentile thresholds over the distribution of the combined voxets.We compute two topological distances: the structural Hamming distance and the Laplacian spectral distance (Shimada et al., 2016).Note that graphs characterizing geological model topology could be defined as attributed graph, to contain more information (edges properties such as age constraint, type of contact; vertices properties such as formation type, geophysical properties).
Thus more specific measures could be developed to take into account such characteristics.However, it would rely on the ability of geomodelling engines to provide these topology graphs with each model, and there is no guarantee that it would be meaningful for the inference of geochronology from geophysics.

Results: indicator comparison
Local measures of uncertainty (see Figure 14 & 15) and global indicators (see Figure 17 & 18) have been computed for 2D, 3D, categorical and continuous variable voxets (lithocode, scalar-field) for an ensemble of 30 model-sets.We also provide a comparison of the required computing time for the different indicators and highlight the contributing factors/parameters (see Table 2).
One can see from Figure 14   Figure 16 shows that histogram dissimilarities computed from lithocode or scalar-field voxets have a good correlation.Multidimensional scaling (MDS) plots reveal that dissimilarities are smallest within scenario 1 and then within scenario 3 while they are greatest within scenario 2. MDS plots show that scenarios 1 and 3 model-sets overlaps :::::: overlap : while scenario 2 model-sets are characterized by less similar histograms and thus scenario 2 sample cloud of points form a distinct cluster.
Figure 17 shows some correlation between histogram, semi-variogram, wavelet and structural Hamming based measures from the lithocode voxets.Figure 18 shows a good correlation, between histogram dissimilarity, wavelet based dissimilarity and structural Hamming distance from the scalar-field voxets, but not as strong as when computed from the lithocode voxets.

Discussion
Several factors might explain why a majority of practitioners do not consider input data uncertainty, but all are related to the limited resources available to practitioners (knowledge, algorithms, computing time, project timeframe or funding).One of  them is that data uncertainty is not quantified at the time of data acquisition or not available for some measurements, which is 350 the case when only one measurement or observation of geological data is made at a given location (e.g. for azimuth, dip, and lithology).Reasonable metadata standards may help to enforce error quantification, or at the very least provide some information about the nature of data collection as to infer where and what magnitude of error may be present, but even these are poorly or not recorded.Although repeated independent measurements would provide uncertainty estimates, procedures and limited time or budget resources are a hindrance.Sometimes, knowing the survey setup and the instrumentation characteristics, such as their precision and accuracy might avoid repeating field measures and allow for an estimation of measurement uncertainty.
Inversion algorithms used for geophysical data integration can also provide estimates of geophysical data errors.
A second reason is related to the required time and associated costs of modelling.Indeed, the process of data integration only uses a very limited amount of automation, thus the generation of a single model consumes already most of the practitioners' resources.In addition, the complexity of real world data often leads to a substantial number of parameters.Thus for high-360 dimensional problems, uncertainty propagation requires sufficiently large model ensembles to be representative, which might not be compatible with the limited resources available to the practitioners.
A third reason is due to the fact that assumptions, such as choices in geological interpretations, made during the modelling process are not always tracked.And when they are, they are often 'forgotten' at the next stage of the modelling workflow (Jessell et al., 2018).When these assumptions or justifications are recorded, they are described in metadata or in distinct reports.Consequently, conceptual uncertainty, which describes alternative yet plausible stratigraphy, tectonic and geodynamic settings is also ignored.
Another possible reason is that uncertainty is ignored out of convenience (Ferré, 2017) or by ignorance, lack of knowledge or education about the importance of uncertainty quantification.However, the formulation of the collected answers suggest that it is not the case for the surveyed practitioners, who rather acknowledge the importance and need for tools or algorithm to integrate uncertainty quantification in their modelling workflow.
While about 11% of the respondents indicate that they perform a quantitative uncertainty quantification, it is limited to aleatoric uncertainty, i.e. data measurement errors.However, the lack of spatial data samples contribute to epistemic uncertainty and our limited contextual knowledge adds up to conceptual uncertainty.In addition, it is generally expected that these sources of uncertainty have a bigger impact on predictive uncertainty (Pirot et al., 2015).Thus, in addition to develop tools to facilitate aleatoric uncertainty quantification for practitioners, accessible tools integrating epistemic and conceptual uncertainty quantification need to be developed and promoted in the minerals industry.
Restricting the survey to a few open questions encouraged participants to take the survey and to express their uses, needs and opinion with limited perception bias, however it did not allow to quantify the gathered answers, and often requires some interpretation.Nevertheless, it is a first step in acknowledging the practices and needs of 3D geological modellers in the minerals industry.Another limitation of the survey is that it does not look at the practice and needs of other fields like the petroleum industry (Scheidt et al., 2018), geothermal industry (Chen et al., 2015) or hydrogeology (Pirot et al., 2019), though they share similar scientific problems and also propose interesting solutions to deal with uncertainty quantification.techniques.:::: One ::::: could :::: also :::::::: consider summary metrics of lower dimensional model representations, that could be obtained from discrete-cosine transform (e.g.Ahmed et al., 1974), (kernel-) principal component analysis (e.g.Schölkopf et al., 1997) or (kernel-) auto-encoding (e.g.Laforgue et al., 2019) for instance.This could be particularly appealing as it could reduce indicator computing costs drastically, but might not allow to identify the specific features of interest in a representation space of lower dimensions.
Last but not least, one can note that we compared ensemble of models, whose underlying characteristics such as various lithological units derived from a shared stratigraphy and scalar-fields generated under similar assumptions are consistent.
However, we must warn that the various indicators are compatible with differences in property ranges or meaning (for categorical variables), and thus it is the responsibility of the user to ensure the coherence of the model ensemble used as input for the uncertainty computation.

Conclusions
The survey clearly shows that practitioners acknowledge the importance of uncertainty quantification; a majority recognize that they do not perform uncertainty quantification at all and all would like to do better.From this survey, we have identified geological model, efforts should be made to explore conceptual uncertainty (Laurent and Grose, 2020), as well as towards the implementation of systematic geological data uncertainty quantification, and the exploration of parametric and epistemic uncertainty (Pirot et al., 2020).It should be performed appropriately at all scales, across all geoscientific methods, such as the extraction of additional lithological data from drillhole databases (Joshi et al., 2021).To encourage uncertainty propagation among practitioners, accessible and compatible algorithms should be offered 1) to extract automatically geological data from open-databases (Jessell et al., 2021), 2) to quickly generate plausible geological models from a given dataset (Grose et al., 2021) in interaction with geophysical data integration (Giraud et al., (this issue) and at an appropriate resolution (Scalzo et al., 2021).This special issue (Ailleres, 2020) already provides elements of an answer to these problems and is expected to host future advances on these topics.Author contributions.The survey was initiated by Guillaume Pirot, developed jointly with Ranee Joshi and with the support of all other co-authors.All co-authors distributed the survey and collected answers.Answers were processed and analysed by Guillaume Pirot.The manuscript was mainly written by Guillaume Pirot, with contributions from all co-authors.
Competing interests.No competing interests are present.

Figure 1 .
Figure 1.Evolution of :::: yearly : open-source algorithm publications between 2000 and 2020; data from Web Of Knowledge.

Figure 3 .
Figure 3. Modelling scales and resolutions in the mineral industry with examples of scale-specific scientific enquiry; model illustrations of the Vihanti-Pyhäsalmi Area, Finland, adapted from Laine et al. (2015).

Figure 4 .
Figure 4. Consideration of input data uncertainty in modelling.

Figure 5 .
Figure 5. Example model set of lithocode and scalar fields for each of the three scenarios; the left column illustrates input data and model set for scenario 1 when all data is considered; the middle column illustrates input data and model set for scenario 2 when 50% of the data within a band is considered; the right column illustrates input data and model set for scenario 3 where each input data is decimated with a probability of 50%.

Figure 6 .
Figure 6.Horizontal sections of cardinality voxets for a categorical property field (first row) or similar indicators for a continuous property field (second to fourth row); the first row shows the cardinality over the ensemble of lithocode voxets, the second row displays the (max-min range) over the ensemble of scalar-field voxets, the third row presents the standard deviation over the ensemble of scalar-field voxets, the fourth row displays the averaged normalized range and standard deviation over the ensemble of scalar-field voxets; the left column refers to scenario 1, the middle column to scenario 2 and the right column to scenario 3.
fields, one need to discretize the continuous domain and integrate along the width of the bins(Marsh,   195    2013).Here, for the considered example, we choose 50 regular bins.Figure7displays Shannon's entropy for the lithocode voxets and the continuous entropy for the other ensemble of property fields: magnetic field, gravity field, density and magnetic susceptibility.

Figure 7 .
Figure 7. Horizontal sections entropy voxets; the first row shows Shannon's entropy over the ensemble of lithocode voxets, the second row displays the continuous entropy over the ensemble of scalar-field voxets; the left column refers to scenario 1, the middle column to scenario 2 and the right column to scenario 3.

Figure 8 .
Figure 8. Example of experimental semi-variogram for two lithocode voxets, computed for the lithocode value of 2; the two first columns show an horizontal section for each voxet; the last column shows the two corresponding experimental semi-variogram.
) have shown that connectivity cannot be captured by topological indicators such as the Euler characteristic, nor by one-point or two-point statistics (e.g. by histogram or semi-variogram analysis respectively).However, they have shown how a global percolation metric Γ(p) and a lag-distance connectivity function τ (h) are useful to characterize the connectivity of binary, categorical or continuous property fields.Connectivity indicators have also been used in multiple-point statistics applications to characterize the quality of stochastic simulations with respect to a training image

Figure 9 .
Figure 9. Illustration of the τ (h) connectivity function (right panel) computed on a 2D horizontal section from a binary 3D voxet (left panel).

Figure 10 .
Figure 10.Illustration of a scalar-field voxet horizontal section (top left) and its global percolation metrics Γ(p) and Γ(p) c :::: Γ c (p) : (top right) computed over the 2D sections; the bottom row displays horizontal sections of binary fields obtained by applying different percolation threshold p[%] to the scalar-field 3D voxet; blue areas, above the percolation threshold, contribute to Γ(p); yellow areas, below the percolation threshold, are used to compute Γ(p) c :::: Γ c (p).

Figure 11 .
Figure 11.Illustration of multiple-point histogram cluster representatives and sizes for two scalar-field horizontal sections, at the 3 rd upscaling level; the left column shows the two 2D voxets; for the other columns the first and second rows, respectively third and fourth rows, display the 10 cluster representatives and their size (number of counted patterns attached to the cluster representative) for the first voxet, respectively for the second voxet; their order reflect the best similarity between the cluster representatives for both voxets.

Figure 12 .
Figure 12.Illustration of a first level of haar-wavelet decomposition and the resulting coefficients histograms for two lithocode voxet horizontal sections.

Figure 13 .
Figure 13.Illustration of topological distances and adjacency matrices in the right column for two categorised voxets in the middle column, derived from two scalar-field voxet horizontal sections in the left column.
that Shannon's entropy and cardinality computed from lithocode voxets can have a good correlation.However, equivalent indicators continuous entropy and averaged normalized range and standard deviation computed from scalar-field voxet (Figure15) have very little in common.Indeed, the standard deviation or the range of values are sensitive to extremely different values, while the continuous entropy is sensitive to the proportion of categories of values.

Figure 14 .
Figure 14.Comparison of local uncertainty measures for an ensemble of 10 lithocode 3D voxets for scenario 1; 3D visualisation looking from the NW of the voxet, the top surface of the voxet an EW section at the northern face of the model looking from the south, a NS section on the western face of the voxet looking from the east.

Figure 15 .
Figure 15.Comparison of local uncertainty measures for an ensemble of 10 scalar-field 3D voxets for scenario 1; 3D visualisation looking from the NW of the voxet, the top surface of the voxet an EW section at the northern face of the model looking from the south, a NS section on the western face of the voxet looking from the east.

Figure 16 .
Figure 16.Comparison of model set histogram dissimilarities across the three scenarios; the left column displays a 2D MDS representation of histogram dissimilarities computed from lithocode voxets (top row) and from scalar-field voxets (bottom row); samples 0 to 9 belong to scenario 1, samples 10 to 19 belong to scenario 2 and samples 20 to 29 belong to scenario 3; the middle and right column subplots show histogram and density, joint density and cross-plot between histogram dissimilarities computed from lithocode voxets or from scalar-field voxets.

Table 1 .
Global number and rate of answers per question and detailed by modelling scale.
Marsh, 2013;Pirot et al., 2017)12;Pakyuz-Charrier et al., 2018)ty value spectrum depending on how physical properties are assigned or what input constraints are enforced during model construction.Local measures of uncertainty provide indicator voxets of the same dimensions than the voxets of a model ensemble.For a given voxel, uncertainty indicators are computed from the distributions of values taken by a given property at the corresponding voxel (same location) across the ensemble of model realizations.Such local indicators are very convenient for visualization:by sharing the same voxet as the model realizations, it is relatively easy to spot zones of low or high uncertainty.However, to be useful, it requires advanced modelling that integrates some spatial constraints, and possibly computationally expensive, in particular if they are produced by inversion algorithms.Here, we propose to compute Cardinality and Shannon's Entropy for discrete properties (e.g.Wellmann and Regenauer-Lieb, 2012;Pakyuz-Charrier et al., 2018), and similarly, range, standard deviation and continuous Entropy for continuous properties (e.g.Marsh, 2013;Pirot et al., 2017).

Table 2 .
Complexity and computing time for local and global measures of uncertainty using a single Intel(R) Core(TM) i7-8550 1.80GHz, based on an ensemble size of N = 10 geological models.