Guillaume Rongier

Pauline Collon

Philippe Renard

Julien Straubhaar

Judith Sausse

Comparing connected structures in ensemble of random elds

Keywords:

Very different connectivity patterns may arise from using different simulation methods or sets of parameters, and therefore different flow properties. This paper proposes a systematic method to compare ensemble of categorical simulations from a static connectivity point of view. The differences of static connectivity cannot always be distinguished using two point statistics. In addition, multiplepoint histograms only provide a statistical comparison of patterns regardless of the connectivity. Thus, we propose to characterize the static connectivity from a set of 12 indicators based on the connected components of the realizations. Some indicators describe the spatial repartition of the connected components, others their global shape or their topology through the component skeletons. We also gather all the indicators into dissimilarity values to easily compare hundreds of realizations. Heat maps and multidimensional scaling then facilitate the dissimilarity analysis. The application to a synthetic case highlights the impact of the grid size on the connectivity and the indicators. Such impact disappears when comparing samples of the realizations with the same sizes. The method is then able to rank realizations from a referring model based on their static connectivity. This application also gives rise to more practical advices. The multidimensional scaling appears as a powerful visualization tool, but it also induces dissimilarity misrepresentations: it should always be interpreted cautiously with a look at the point position confidence. The heat map displays the real dissimilarities and is more appropriate for a detailed analysis. The comparison with a multiple-point histogram method shows the benefit of the connected components: the large-scale connectivity seems better characterized by our indicators, especially the skeleton indicators.

I

Connectivity is a key aspect of a geological study for its influence on fluid circulations. From a reservoir engineering perspective, it relates to geological structures with high and low permeabilities. But it also relates to the spatial distribution of these structures and the resulting inter-connections, which define the static connectivity. An incorrect connection can bias the results of the flow simulations [START_REF] Journel | Beyond covariance: the advent of multiple-point geostatistics[END_REF]Alabert, 1990, Gómez-Hernández and[START_REF] Gómez-Hernández | To be or not to be multi-Gaussian? A reflection on stochastic hydrogeology[END_REF]]. Reproducing the geological bodies together with their relations is so of prime importance [e.g., [START_REF] Deutsch | Fortran programs for calculating connectivity of threedimensional numerical models and for ranking multiple realizations[END_REF]Hewett, 1996, King and[START_REF] King | Flow Simulation of Geologic Models[END_REF].

Stochastic simulations aim at generating possible representations of the geological bodies with respect to the available data. Several methods exist, with a usual separation in two main categories:

• Pixel-based methods simulate one cell at a time, based on a prior model describing the structures of interest.

In sequential indicator simulation (SIS) [START_REF] Deutsch | GSLIB: Geostatistical Software Library and User's Guide[END_REF], the prior is a variogram built upon the two-point statistics of the data. Hard data conditioning with such method is easy. But the simulated structures do not look like geological bodies. This is especially true for bodies with curvilinear geometries such as channels, whose continuity is badly preserved. The plurigaussian simulation (PGS) [START_REF] Galli | The Pros and Cons of the Truncated Gaussian Method[END_REF]] limits this difficulty by accounting for the facies relationships. Multiple-point simulations (MPS) go a step further by borrowing multiple-point statistics not from the data but from an external representation of the expected geology, the training image (TI) [START_REF] Guardiano | Multivariate Geostatistics: Beyond Bivariate Moments[END_REF]].

• Object-based methods rely on the definition of geometric forms and their associated parameters. Each form represents a particular geological body [e.g., Viseur, 2001, Deutsch and[START_REF] Deutsch | FLUVSIM: a program for object-based stochastic modeling of fluvial depositional systems[END_REF]. The objects are then randomly placed in the domain of interest with parameters drawn in statistical laws. More recent approaches introduce some genetic aspects to improve the object organization [e.g., [START_REF] Lopez | Modélisation de réservoirs chenalisés méandriformes : une approche génétique et stochastique[END_REF][START_REF] Pyrcz | ALLUVSIM: A program for eventbased stochastic modeling of fluvial depositional systems[END_REF]. They provide more geologically consistent results. For instance channel continuity and relationships are better preserved than with pixel-based methods. But this is at the cost of the ease of parametrization. And object-based approaches have difficulty to condition the objects to data.

All these methods have advantages and drawbacks. This will influence the choice of a method and its parameter values when dealing with a case study. But few work aims at systematically analyzing the quality of a set of realizations regarding their static connectivity. The quality control often consists in comparing the histogram and variogram of several realizations with those of the data, or of the training image if any [e.g., [START_REF] Strebelle | Conditional Simulation of Complex Geological Structures Using Multiple-Point Statistics[END_REF][START_REF] Mariethoz | The Direct Sampling method to perform multiple-point geostatistical simulations[END_REF][START_REF] Tahmasebi | Multiple-point geostatistical modeling based on the cross-correlation functions[END_REF]. If more than the first twoorder statistics are necessary to simulate geological bodies [e.g., [START_REF] Guardiano | Multivariate Geostatistics: Beyond Bivariate Moments[END_REF]Srivastava, 1993, Journel, 2004], the same conclusion must apply when comparing realizations. Some authors propose to also use the higher-order statistics for quality analysis. [START_REF] Boisvert | Multiple Point Metrics to Assess Categorical Variable Models[END_REF] and [START_REF] Tan | Comparing Training-Image Based Algorithms Using an Analysis of Distance[END_REF] propose to analyze the multiple-point histogram. De Iaco and Maggio [2011] and De [START_REF] Iaco | On the use of different metrics for assessing complex pattern reproductions[END_REF] also explore the multiple-point statistics with high-order cumulants.

The purpose of most simulation methods is to reproduce statistics from a prior. Analyzing statistics highlights the method success in this reproduction, not if the realizations are geologically consistent. To do that, the statistical analysis is often completed by a visual evaluation of the global structures. The geological structures are compared to what is expected from the known geology, with a focus on the further use of the realizations. This use is often related to fluid circulations, and requires an assessment of the static connectivity, which is not directly imposed by the simulation methods contrary to the statistics. But a visual analysis remains subjective and limited to a few realizations, often in two-dimensions [e.g., [START_REF] Yin | A reservoir skeleton-based multiple point geostatistics method[END_REF][START_REF] Tahmasebi | Multiple-point geostatistical modeling based on the cross-correlation functions[END_REF].

Yet, some studies focus on analyzing the connectivity of the realization bodies. For instance [START_REF] Meerschman | A Practical Guide to Performing Multiple-Point Statistical Simulations with the Direct Sampling Algorithm[END_REF] use the connectivity function with the histogram and variogram to analyze the simulation parameter impact for the Direct Sampling MPS method [START_REF] Mariethoz | The Direct Sampling method to perform multiple-point geostatistical simulations[END_REF]. [START_REF] Deutsch | Fortran programs for calculating connectivity of threedimensional numerical models and for ranking multiple realizations[END_REF] uses directly the connected components determined from lithofacies, porosity and permeability models. He computes indicators such as the number of connected components or their sizes to rank the realizations. De [START_REF] Iaco | Validation Techniques for Geological Patterns Simulations Based on Variogram and Multiple-Point Statistics[END_REF] and De [START_REF] Iaco | On the use of different metrics for assessing complex pattern reproductions[END_REF] also use some measures related to the connected components, such as their number or their mean surface and volume. [START_REF] Comunian | 3D multiple-point statistics simulation using 2D training images[END_REF] rely on some of the previous indicators to analyze the quality of three-dimensional structures simulated from two-dimensional training-images. They also consider the equivalent hydraulic conductivity tensor as an indicator. However, this requires to have an idea of the hydraulic conductivities for the simulated facies.

Connected components enable to characterize the geometry and topology of the geological bodies, which is the purpose of the visual comparison of realizations. They also enable to study the static connectivity of the geological bodies, while being easy to compute. Contrary to a visual analysis of the realization, indicators from connected components are unbiased and can compare many realizations. Contrary to statistical or hydraulic property indicators, they focus on the sedimentary bodies by characterizing their connectivity and are more easy to apprehend. However, current methods based on the connected components are limited to few simple indicators, often analyzed independently.

This leads to the question of the result visualization to analyze more effectively the indicators. [START_REF] Scheidt | Representing Spatial Uncertainty Using Distances and Kernels[END_REF] and [START_REF] Tan | Comparing Training-Image Based Algorithms Using an Analysis of Distance[END_REF] both rely on the computation of dissimilarity values between the realizations. Those dissimilarities are computed based on the quality indicators measured on each realization. They are then visualized based on a Multi-Dimensional Scaling (MDS) [e.g., [START_REF] Torgerson | Multidimensional scaling: I. Theory and method[END_REF], Shepard, 1962a,b]. MDS represents the realizations as points, with the distance between the points as close as possible to the dissimilarities. The global analysis of the realization dissimilarities is so easier.

The present work aims at analyzing and discussing a set of indicators to quantify the quality of stochastic simulations from the viewpoint of static connectivity. This method performs on categorical three-dimensional images representing the facies constituting the geological bodies of interest. It can be applied on realizations from one or several stochastic simulation methods and/or parameter values. Conceptual images representing ideally the structures to simulate can also be considered. The chosen set of indicators relies on quantitative measurements on connected components and their skeletons (section 1). The indicators are used in dissimilarity computations to analyze the quality more directly (section 2). Several realizations obtained with different simulation methods (section 3.1) are then used to test the method and compare it to the multiple-point histograms (section 3), and discuss the results (section 4).

I -

The quality analysis fits in a stochastic process implying the simulation of many realizations in a grid. It further investigates the differences of static connectivity between these realizations.

About grids and grid cells

Many methods to simulate geological structures rely on a discretized representation of the domain of interest: a grid. The grid is a volumetric mesh composed of simple elements, hereinafter called cells. Many types of grid exist, with different cell types (e.g., tetrahedron or hexahedron). Most of the stochastic simulation methods rely on hexahedral grids, either regular or irregular. Irregular hexahedral grids help to be as conform as possible to the geological structures such as horizons and faults. The sedimentary bodies are then simulated within the parametric space of the grid [e.g., [START_REF] Shtuka | Petrophysical simulation within an object-based reservoir model[END_REF]. The parametric space mimics a deposition space to get rid of the deformation and faulting occurring after deposition and materialized in the grid geometry.

Consequently, the indicators are computed on hexahedral grids, both regular and irregular. Similarly to the simulation, the indicator computation is done in the parametric space of the grid. Thus, the indicators based on volumes or surfaces are rather computed using number of cells and number of faces. This avoids biases related to different grid geometries, which give different indicator values even if the objects are the same when transferred in the same grid.

Within a grid, the cells are connected one to another by their faces, their edges and/or their corners [START_REF] Renard | Connectivity metrics for subsurface flow and transport[END_REF]]. In the case of the hexahedral grids used for this work, one cell has three possible neighborhoods (figure 1):

1. One neighborhood composed of six face-connected cells.

2. One neighborhood composed of eighteen face-and edgeconnected cells.

3. One neighborhood composed of twenty-six face-, edgeand corner-connected cells.

This definition of the connectivity between a cell and its neighborhood can be extended to form connected components.

Central cell

Face-connected cell Edge-connected cell Corner-connected cell

Basic element: the connected component

The connected components result from the widening of the neighborhoods. They rely on the following definition of the connectivity between two cells: two cells belonging to the same facies are connected if a path of neighboring cells remaining within the same facies exists (figure 2). Applying this definition to all the cells of a facies gives the connected components of this facies.

This leads to a distinction between the geological objects, such as a channel or a crevasse splay, and the connected components. Indeed, the geological objects often tend to cross each others, giving one connected body where there is in fact several geological objects (figure 2). The range of possible shapes is larger for the connected components than for the individual object. This aspect complicates the comparison between images. But determining the connected components is far easier than trying to retrieve the geological objects. This is also close to the functioning of pixel-based methods, which do not try to reproduce geological objects but groups of cells, and therefore connected components. Here the nodes connected to only one segment -the nodes of degree 1 -are all along a grid border. Two nodes of degree three highlight the local disconnections between the channels at the bottom.

The connected component at the top has no node of degree higher than two, which shows the complete connectivity of all its cells, even locally.

Basic element: the skeleton

A curve-skeleton -simply called here skeleton -is a thin onedimensional representation of a three-dimensional shape. It is composed of nodes linked together by one or more segments (figure 3). The degree of a node is the number of segments connected to that node. Skeletons are often used to study some geometrical and topological features of a shape. Here the skeletons are those of the connected components. They enable to further characterize the global shape of the connected components, while giving more details about their topology than indicators directly computed on the components. Several methods exist to compute skeletons [e.g., [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Jain | Fundamentals of digital image processing[END_REF][START_REF] Brandt | Continuous skeleton computation by Voronoi diagram[END_REF]. The method considered for this work is based on slicing the grid along a given axis. The grid is subdivided into parallel slices of a given thickness. On each slice the connected components are computed and one node is assigned to each component. The nodes are then linked by computing the connected components over two adjacent slices. If two components from two slices form one connected component when the slices are combined, their nodes are linked. If they form several components, their nodes are not linked.

Indicators

The indicators studied in this paper focus on analyzing the connectivity of the geological bodies within a three-dimensional image. This static connectivity analysis is possible thanks to the connected components. All the indicators are quite simple and each one gives only partial information about the connectivity and its structure. But their combination provides a more detailed characterization.

Appendix A defines in detail all the indicators. Table 1 summarizes the indicator definition, by focusing on their relationship with the connected components. We distinguish three categories of indicators:

Global indicators: The global indicators characterize a facies

Table 1

The set of indicators with their definitions. The indicator definitions is given by a numerator and a denominator because the majority of the indicators comes from a ratio. Some ratios are computed on a single connected component, and their values are combined to obtain one indicator value per facies. We use the term "component" instead of "connected component" for the sake of simplicity. and not necessarily an individual connected component. Among them, the facies proportion is a classical indicator to compare realizations. Some others, such as the facies connection probability [START_REF] Renard | Connectivity metrics for subsurface flow and transport[END_REF], the connected component density or the traversing component proportion give an idea of the global connectivity.

Category Indicator

Symbol

Shape indicators:

Global measures such as facies proportions are not sufficient to characterize precisely the impact of the related facies on the flow [e.g., [START_REF] Western | Toward capturing hydrologically significant connectivity in spatial patterns[END_REF][START_REF] Mariethoz | Geological stochastic imaging for aquifer characterization[END_REF]. In particular, Oriani and Renard [2014] showed the influence of the connected component geometry -i.e., their shape -on the equivalent hydraulic conductivity, and therefore on the flow behavior. The shape indicators characterize the connected component shape through simple surface and volumetric measures. They all give one value per component.

The arithmetic mean of those values provides a value of the indicator for a given facies. This makes the indicator comparison easier.

Skeleton indicators:

The skeletons help to better character- ize the topology and global geometry of their connected components: their one-dimensional representation is easier to analyze. Here two indicators are introduced.

The inverse branch tortuosity characterizes the geometry of the skeleton. Its values for all the branches of all the skeletons related to a facies are averaged to obtain a single value for the facies. It completes the shape indicators in the characterization of the connected component shape. The node degree proportion depicts the topology of the skeletons. It helps to analyze the connectivity more precisely.

Q

The final purpose of this work is to easily and objectively compare several realizations. The indicators are thus computed on large sets of realizations, which may come from different methods and/or parameters. Then dissimilarity values based on the indicators help to compare the realizations.

In uence of different grid dimensions

Some cases imply to compare realizations on different grids, and the grids may have different dimensions. For instance in MPS, the training image is often larger than the simulation grid to maximize pattern repeatability. The grid dimensions influence the size of the traversing connected components, such as channels. This impacts in particular the connected component density and the number of component cells. When the grid size varies along the channel direction, the number of cells for the channels also varies. And even though the number of channels does not necessarily change, the grid volume does, impacting the density. These indicators highlight expected differences in such cases. Their direct use is then detrimental to the quality analysis.

We propose two workarounds to compensate for different grid sizes:

• Either sampling the images from the different grids so that all the samples have the same dimensions. The sample size are the largest dimensions common to all the grids. Each sample is randomly extracted and each image may be sampled several times to still catch the characteristics of the whole image.

• Or correcting the indicators of the difference between the grid dimensions. The smallest grid dimensions among all the grids form a hypothetical reference grid. The indicators are corrected to their expected value in such reference grid. Appendix B details this correction.

The sampling exempts from correcting the indicators, but it adds a step and requires the analysis of more images, which could slow down the process. If they are valid, the corrections should give similar results than the sampling in a more efficient process.

Indicator rescaling

The rescaling ensures that the differences between the ranges of indicator values will not affect the comparison. The histogram-based indicators -facies proportion, facies adjacency proportion and node-degree proportion -are not rescaled, to preserve their histogram behavior for the dissimilarity computation (section 2.3). Two methods can be used for rescaling: normalization and standardization.

The normalization method consists in rescaling linearly the indicators values between 0 and 1. The indicator I i is the i th indicator of the set previously defined. When computed for the facies f of the realization r, we will denote the computed indicator I r i f . The normalization is then obtained by rescaling it between its minimum and maximum values:

norm(I r i f ) = I r i f -m i f M i f -m i f (1)
with M i f the maximum value for the same indicator and facies among all the images and m i f the minimal value for the same indicator and facies among all the images. The standardization method consists in using reduced-centered indicator values. For an indicator i the standardized value for a facies f of a realization r is obtained using the following formula:

stand(I r i f ) = I r i f -µ i f σ i f (2)
with µ i f the mean for the same indicator and facies among all the images and σ i f the standard deviation for the same indicator and facies among all the images. Standardization is an interesting option to focus on the indicator variance. The normalization on the other hand decreases the influence of outliers and gives precise limits to the indicator values.

Dissimilarity calculation

The principle of comparing two images is to determine how dissimilar these images are. The indicators can be seen as coordinates of the compared images. These indicators are heterogeneous: they are either based on histograms or on continuous values. The computation of a dissimilarity value between two images requires a heterogeneous metric.

Following the example of [START_REF] Wilson | Improved heterogeneous distance functions[END_REF], two different metrics are combined into a heterogeneous Euclidean/ Jensen-Shannon metric. It uses the Jensen-Shannon distance, square root of the Jensen-Shannon divergence [START_REF] Rao | Differential metrics in probability spaces[END_REF][START_REF] Lin | Divergence measures based on the Shannon entropy[END_REF], for the histogram-based indicators -facies proportion, facies adjacency proportion and node-degree proportion -and the Euclidean distance for all the other indicators. The distance between two images r and s for a given indicator i of a given facies f is given by:

d(I r i f , I s i f ) = d JS (I r i f , I s i f ) if I r i
d JS (H r i , H s i ) = 1 2 n j=1     H r i j log    H r i j 1 2 (H r i j + H s i j )    + H s i j log    H s i j 1 2 (H r i j + H s i j )        (4)
with H r i and H s i the histograms of the indicator i for respectively the images r and s, n the number of classes for each histogram, H r i j and H s i j the proportions for the class j in the corresponding histograms. d E represents the Euclidean distance used with rescaled indicators:

d E (I r i f , I s i f ) = (resc(I r i f ) -resc(I s i f )) 2 (5)
with I r i f and I s i f the values of the indicator i for the facies f of respectively the images r and s and resc either norm (formula 1) or stand (formula 2). The final dissimilarity δ between two images r and s given their respective sets of indicators I r and I s is:

δ(I r , I s , ω, ν) = ω 1 d JS (I r 1 , I s 1 , ν) 2 + 12 i=2 n f =1 ω i ν f d(I r i f , I s i f ) 2
(6) with I r 1 and I s 1 the facies proportion histogram for the two images, I r i f and I s i f all the other indicator values depending on the indicator and the facies and n the number of facies. ω represents the set of weights ω i that control the impact of each indicator. ν represents the set of weights ν f that control the impact of each facies. Note that the facies proportion histograms are the only indicators with one result for all the facies. Thus the facies proportions are treated differently from all the other indicators. The Jensen-Shannon distance used in that case is slightly modified:

d JS (H r i , H s i , ν) = 1 2 n f =1 ν f     H r i f log    H r i f 1 2 (H r i f + H s i f )    + H s i f log    H s i f 1 2 (H r i f + H s i f )        (7)
The dissimilarity values computed by formula 6 between all the images constitute a non-negative symmetric matrix. This matrix has a zero diagonal corresponding to the dissimilarity between an image and itself. The dissimilarity matrix can be directly visualized with a heat map or treated by multidimensional scaling to get a more practical visualization.

Heat map

The heat map is a simple graphical representation of a matrix where the matrix values correspond to colors. In our case, the heat map is a two-dimensional representation. This colored representation highlights patterns in the dissimilarity matrix, either between realizations or between simulation methods. The main advantage of the heat map is to show the real dissimilarity values, contrary to the multidimensional scaling described in the next subsection.

The heat map also enables to classify the images and/or to apply clustering methods on it. A simple yet informative classification is the ranking according to the dissimilarities of the images toward one particular image. When using more advanced clustering methods, the matrix rows and columns are permuted to gather close values into the same cluster.

Multidimensional scaling

Multidimensional scaling (MDS) [e.g., [START_REF] Torgerson | Multidimensional scaling: I. Theory and method[END_REF][START_REF] Torgerson | Theory and methods of scaling[END_REF], see Cox and Cox, 1994 for a review] is a set of data visualization methods to explore dissimilarities between objects -represented by a dissimilarity matrix -through a dimensionality reduction: it aims at producing a configuration of the objects as optimal as possible in a lower dimensional representation.

Principle and method used

Finding the configuration of the images in a k dimensional representation consists in locating a set of points representing the objects in a k-dimensional Euclidean space -with k being at most equal to the number of images minus one. The point positioning is done so that the Euclidean distance d between two points matches as closely as possible the dissimilarities between the images:

d r,s = k i=1 (x r i -x si ) 2 (8)
with r and s two images, k the dimension number of the Euclidean space, x r i and x si the coordinates of respectively r and s in the i-th dimension. The number of dimension k for the MDS representation is an input parameter. When equal to the number of images minus one, the distances are normally equal to the dissimilarities. When k is lower, the MDS misrepresents more or less the dissimilarities. Several multidimensional scaling methods have been proposed [e.g., [START_REF] Cox | Multidimensional scaling[END_REF], depending on the type of dissimilarities and on the way to match the dissimilarities with the distances. The classical scaling [START_REF] Torgerson | Multidimensional scaling: I. Theory and method[END_REF][START_REF] Torgerson | Theory and methods of scaling[END_REF][START_REF] Gower | Some distance properties of latent root and vector methods used in multivariate analysis[END_REF] is the usual method for multidimensional scaling [START_REF] Scheidt | Representing Spatial Uncertainty Using Distances and Kernels[END_REF]Caers, 2009, Tan et al., 2014, e.g.,]. It assumes that the dissimilarities already are Euclidean distances. If this assumption can be relaxed to a metric assumption, i.e., the dissimilarities are distances, Euclidean or not, the classical scaling may further misrepresents dissimilarities based on a heterogeneous metric.

Here we use a different method: the Scaling by MAjorizing a COmplicated Function (SMACOF) [De [START_REF] Leeuw | Convergence of correction matrix algorithms for multidimensional scaling[END_REF][START_REF] Leeuw | Convergence of correction matrix algorithms for multidimensional scaling[END_REF][START_REF] Leeuw | Multidimensional scaling with restrictions on the configuration[END_REF]. Its goal is to get distances as close as possible from the dissimilarities using a majorization, i.e., the optimization of a given objective function called stress, through an iterative process. The stress derives from the squared difference between the dissimilarities and the distances. It is positively defined and equals to 0 only when the distances are equal to the dissimilarities. The optimization process corresponds to a minimization of the stress. The final stress value helps to assess the choice of the number of dimensions: the lower the stress is, the better is that choice.

Validation of the number of dimensions

Following the chosen number of dimensions for the representation, the point configuration matches more or less the dissimilarity values. Verifying that the dimension number is enough for a good match between the dissimilarities and the distances is so of prime importance. Two approaches allow testing the chosen dimension number:

The scree plot: It represents the stress of the SMACOF against the dimension number. The stress follows a globally convex decreasing function that tends toward 0 when the dimension number increases. A stress close or equal to zero means that the higher dimensions are unnecessary to represent the dissimilarities. The best number of dimensions is between the point with the highest flexion of the curve and the beginning of the sill at zero. The dimension value right after the point with the highest flexion is generally enough for a decent representation.

The Shepard diagram:

It represents the distances against the dissimilarities. The better the correlation, the better the choice of dimension number.

Two-dimensions are more practical for an analysis purpose. A three-dimensional representation remains a possibility if the improvement is significant enough from a two-dimensional representation.

Estimation of the point position con dence

The point position confidence is another way to assess the MDS ability to represent the dissimilarities. For each point r, an error e highlights the mismatch between the dissimilarities δ and the distances d with all the other points s:

e r = s |(aδ r,s + b) -d r,s | (9)
with a and b the linear regression coefficients found on the Shepard diagram. This measure gives a more local representation of the miss-representation than the scree plot or the Shepard diagram.

For visualization purpose, that error is then normalized, giving the confidence c for a given image r: c r = 1 -e r -e min e max -e min (10) with e max and e min respectively the greatest and the lowest error values amongst the errors of all the images. This confidence can then be attributed to its corresponding point in the MDS representation through the point transparency: the less transparent the point is, the best the dissimilarities related to this point with all the other points are represented.

E

The method, as described in the previous sections, consists in three steps:

1. Indicator computation.

2. Dissimilarity computation in a matrix.

Dataset

The dataset falls within the simulation of a channelized system. It contains several realizations representing the same sedimentary environment simulated with different methods. The analysis aims at highlighting the indicator ability to capture the differences of static connectivity between the realizations, and especially between the realizations from different methods.

As it concerns a sole case, it would be inappropriate to draw general conclusions on the simulation methods themselves.

The channelized system is composed of sandy channels with levees into a mudstone environment. A conceptual model, called the training image (TI) (figure 4, image at the top), provides an ideal representation of this system. The case study falls within the scope of a MPS study: several simulation methods are used to reproduce the sedimentary bodies observed in the training image. MPS performs better when the training image is larger than the realizations, to ensure enough pattern repeatability. It involves two grids: the first one for the training image (figure 4, image at the top) and the second one for the realizations (figure 4, images at the bottom).

The training grid contains two sets of images:

Analysis setting

The purpose here is to compare the realizations with the TI. It should lead to retain the method and associated parameters that reproduce at best the static connectivity of the TI for the studied case. The indicators used in this case study (table 2) rely on the face-connected components, because the faceconnectivity between cells is the most frequently used [START_REF] Renard | Connectivity metrics for subsurface flow and transport[END_REF]. All the indicators are equally considered (ω i = 1 for all i in formula 6). This avoids any subjective bias that could arise from favoring a given indicator. The mudstone environment is the resultant of the channels and levees placement. It has so no precise shape by itself and may blur the analysis. It gets a weight of 0 while channels and levees both get each a weight of 1 (ν mudst one = 0, ν channel = 1 and ν levee = 1). Channels and levees are considered equally important to reproduce, but this aspect is related to the case study and could be further discussed. The indicators are normalized to cancel the differences of different indicator ranges. Slices of 17 cells along the grid axis with the same orientation than the channels are used for the skeletonization. Several samples are also randomly extracted from the grids to evaluate the suitability of correcting the indicators when dealing with different grid sizes. The training grid having 500× 500 × 20 cells and the simulation grid 100 × 150 × 30 cells, the common largest dimensions for the samples are 100×150×20. The training grid is almost 10 times larger than the simulation grid. Therefore, 20 samples are extracted from the TI and each analog, whereas 2 samples are extracted from each DeeSse, IMPALA, OBS and SIS realization.

Visual inspection of the realizations

Looking at the connected components (figure 5) highlights some expectations for the dissimilarity analysis. Two aspects must be analyzed: the reproduction of the sedimentary body shapes and the reproduction of their connectivity, especially concerning the channels. In the studied case, the reproduction of the shape is pretty easy to analyze visually. The SIS realizations do not display any objects similar to channel/levee systems and are so far dissimilar from the TI. The OBS realizations look similar to the TI, which is what is expected considering that they come from the same method and parameters. DeeSse realizations have objects similar to channels, even if some continuity issues appear. They also seem to have an insufficient number of channels. IMPALA realizations have quite linear objects but which poorly reproduce channel and levee shapes.

Estimating the static connectivity in three-dimensional images is more difficult. The TI channels seem highly connected. The objects in the SIS realizations do not locally intersect like channels do and are far too connected. DeeSse realizations con- tain less objects and seem under-connected compared to the TI. The distinction between OBS and IMPALA realizations is difficult concerning the connectivity. Looking at the skeletons of the connected components (figure 6) corroborates those observations. DeeSse realizations are clearly under-connected compared with the other categories. SIS ones are over-connected. IMPALA realizations seem a bit more connected than OBS ones. The static connectivity within the training image is clearly heterogeneous.

Effect of different grid dimensions on the analysis

The TI, analogs and OBS realizations come from the same method with the same parameter values. The grid size is the only difference between all these images: the grid of the TI and analogs -the training grid -is larger than that of the OBS realizations -the simulation grid. This difference of grid dimensions directly impacts the connected component density and the number of connected component cells, which are corrected to take into account such difference (appendix B). But the realizations coming from the same method still differ when looking at the dissimilarities (figure 7, MDS representation from the original images). The OBS realizations within the simulation grid stand out from the TI and the analogs within the training grid. Such difference is absent from the samples, where all the images have the same size (figure 7, MDS representation from the image samples). The grid size seems to clearly impact the dissimilarity values.

However, both MDS representations (figure 7) have high stress values with two dimensions and can not be fully trusted. The heat maps (figure 8) clarify that situation. The heat map from the original images (figure 8, bottom left) appear nonhomogeneous. A red square symbolizes the significant dissimilarities between the TI and analogs on one side and the OBS realizations on the other side. The heat map from the samples if far more homogeneous, without red square. They confirm the impact of the grid size observable on the MDS representations.

Thus, correcting the connected component density and the number of connected component cells is not adequate, and other indicators are impacted by the grid dimensions. The TI, the analogs and the OBS realizations have all similar channel and levee proportions (figure 9). The channels and levees occupy the same volume inside the two grids. But the facies connection probabilities for both channels and levees differ between the realizations in the two grids (figure 9). The probability that two cells of the same facies belong to the same connected component is higher in the training grid than in the simulation grid. This is consistent with the difference of grid dimensions. When the grid dimension along the channel di- rection increases, the probability that two channels cross each other to form a single connected component increases too, especially here with sinuous channels. In such case, the grid size impacts the characteristics of the connected components and the associated indicator values.

Comparing samples appear to be essential with grids of different dimensions. And using samples reveals other aspects of the images. For instances, the different samples coming from the TI are highly dissimilar. This illustrates the non-stationarity of the TI concerning the connectivity: some areas contain only one connected component as the channels are all connected, whereas other areas contain more connected bodies.

Comparing the connectivity of the training image and of the realizations

The purpose is now to compare the training image to all the realizations. These realizations come from different methods, but all borrow their input from the training image and have to reproduce the sedimentary bodies of the training image. All the following analysis relies on the image samples and not on the original images to avoid any bias due to the difference of size between the training image and the realizations.

Analysis of the dissimilarities

The dissimilarities give a first insight on the relationships between the different realizations (figure 10). The training image samples fall within the OBS samples, highlighting the similarity of these images. The samples from the multiple-point methods, DeeSse and IMPALA, are close to the OBS samples, but they do not mix up much. All these images are so not completely similar. Furthermore, the DeeSse and IMPALA samples remain away from the TI samples. The SIS samples are clearly distinct from all the other samples, and are the most distant from the TI samples.

If the confidence of the two-dimensional MDS representation is not high, the heat map confirms those observations (figure 11). The first row shows the dissimilarities between the training image samples and the realization samples. The whitest samples, the OBS ones, are the closest to the TI. The reddest samples, the SIS ones, are the furthest from the TI samples. The DeeSse and IMPALA samples fall in between, and seem equally close to the TI samples. Globally, the differences between all the methods are significant, as highlighted on the MDS.

As observed in the previous section, the training image samples are dissimilar one from the other. It shows the heterogeneity of static connectivity within the training image. Concerning the realizations, the OBS realizations are also dissimilar one from the other, whereas the SIS realizations are all really close. Both DeeSse and IMPALA realizations are more spaced than the SIS ones, but not as much as the OBS realizations. All this tends toward a variable diversity concerning the static connectivity for the different methods in that case study. Going back to the indicators helps to further analyze such behavior.

Analysis of the indicators

The indicator values for the channels (figure 12 and14) and levees (figure 13 and 15) differ depending on the category. The differences are more or less clear depending on the indicator, whose behavior differs between the two sedimentary body types.

The OBS samples being similar to the training image samples appear also on the indicator values. These values are close -and for many indicators the closest -to the TI values for the channels. That trend is less obvious with the levees, with less close values. But the levee density is the only indicator to be really away from the TI values. All this confirms the close relationship between the training image and the OBS realizations concerning the static connectivity. It also confirms the visual observations. This is consistent with the use of the same method and parameters to simulate the training image and the OBS realizations.

Similarly, the significant dissimilarity between the SIS samples and the TI samples also appears on the indicator values. This is obvious on the traversing component proportion or the component density. The high component density means a higher number of connected components compared to the other samples. On the other side, the average number of component cells is quite low, meaning that most of these numerous components are small. The low traversing component proportion signifies that most of these components are not continuous enough to represent channels nor levees. Concerning channels, the significant difference between the SIS and TI samples for the shape indicators -number of component cells, box ratio, faces/cells ratio and sphericity -implies that the SIS components do not look like channels. This different shape also appears on the node degrees, with far higher node degree values than for the other categories, implying a less linear shape. Despite numerous and small components, the channel connection probability remains high. This means that these samples must contain one large component. This component must be traversing, as the traversing component proportion is not equal to zero. All those observations are consistent with the visual inspection of the realization, and confirm a significant difference of static connectivity between the training image and the SIS realizations. Many indicators also display a narrow range of values. This confirms the low variability between the SIS samples concerning the connectivity, as seen with the dissimilarities.

DeeSse and IMPALA samples have similarities with the SIS samples, especially more, and smaller, connected components than in the training image, as visible on the component density and the number of component cells. Similarly, the shape indicators show a significant difference between the TI samples and both the DeeSse and IMPALA samples. The higher sphericity implies in particular less linear shapes for the channels. Despite being equally dissimilar to the TI samples, the other indicators show significant differences between the DeeSse and IMPALA samples and the TI samples. The DeeSse samples have far lower channel and levee proportions. This impacts the facies adjacency, with channels and levees being more adjacent to the mudstone. But the most relevant difference between the DeeSse and IMPALA samples comes from the channel connection probability: the connection probability of the DeeSse samples is lower than that of the TI samples, whereas the connection probability of the IMPALA samples is higher than that of the TI samples. The IMPALA samples have a behavior sim-ilar to that of the SIS, with a few large component among smaller ones. But these component connectivity is not completely similar to that of the SIS. This is especially visible on the node degree proportions, with the IMPALA samples having an intermediary behavior between the TI and the SIS samples. On the other side, the DeeSse sample connectivity seems lower. The higher degree two proportion of the DeeSse samples implies few intersections between channels. The higher degree one proportion also implies more discontinuous components. Again, all of this is consistent with the visual observations: DeeSse channels are clearly identifiable but discontinuous, whereas IMPALA channels are less visible, with many intersections.

In this case, the indicators confirm what comes from the dissimilarities: the OBS realizations are the most similar to the training image from a static connectivity perspective. This is consistent with the visual observations, and with the use of the same method to simulate the TI and the OBS realizations. The next section endeavors to compare those results from what can be obtained with multiple-point histograms.

Comparison with multiple-point histograms

Multiple-point histograms or pattern histograms have made their way as indicators of a realization quality with MPS methods [e.g., [START_REF] Boisvert | Multiple Point Metrics to Assess Categorical Variable Models[END_REF][START_REF] Tan | Comparing Training-Image Based Algorithms Using an Analysis of Distance[END_REF]. We propose here to compare the results obtained with those histograms to the previous results. The histograms are based on a 3 × 3 × 3 pattern and are computed on three levels of multi-grids [START_REF] Tran | Improving variogram reproduction on dense simulation grids[END_REF], giving three histograms per image. The dissimilarity δ between two images r and s is adapted from the work of Tan et al. [2014]:

δ(H r , H s ) = 3 l=1 1 2 l D JS (H r l , H s l ) (11)
with H r and H s the sets of three histograms for each image, l the multi-grid level and D JS the Jensen-Shannon divergence, which is the squared Jensen-Shannon distance. A multi-grid level l of 1 corresponds to the finest level and here 3 is the coarsest level. The coarser levels characterize the large-scale behavior of the sedimentary bodies. But they induce a loss of information. This justifies the decreasing weights when the multi-grid level increases. Similarly to the work using multiplepoint histograms, the comparison is directly made on the original images, not on samples.

The observations about the category relationships made with the previous indicators (figure 10) remain valid on the MDS representation from the multiple-point histogram (figure 16). The training image falls within the OBS realizations. The DeeSse and IMPALA realizations are close from the OBS ones, but with a clear separation. They all remain separated from the TI. Again, the SIS realizations are far away from all the other images, including the TI. The main difference with the previous indicators comes from the variability within a category. This is especially noticeable with the SIS realizations, which seem to have a significant pattern variability.

The two-dimensional MDS representation is here again a poor representation of the dissimilarities, with a high stress. Only the dissimilarities with the training image are kept to directly study them and compare the ranking between different indicators (figure 17). Looking at all the connected component indicators -i.e., all the indicators described in table 1 -points out the conclusions coming from figure 10: the OBS realizations are the closest to the TI, the SIS ones the furthest, and the DeeSse and IMPALA realizations stand in between. Similar rankings come from the shape indicators -i.e., number of component cells, box ratio, faces/cells ratio and sphericityand the skeleton indicators -i.e., node degree proportions and inverse branch tortuosity.

The multiple-point histograms have also a similar ranking, with a clearer separation between the SIS realizations and the other realizations (figure 17, All the multi-grid levels). However, the dissimilarities between the training image and the DeeSse realizations vary significantly between the multi-grid levels. The largest multi-grid level even places the DeeSse realizations closer to the TI than the OBS realizations. This level characterizes the large-scale behavior of the sedimentary bodies. Such ranking is then particularly surprising due ot the presence of discontinuous bodies within the DeeSse realizations, but neither within the OBS ones nor within the TI. These continuity differences are confirmed by the skeletons, especially the higher proportion of node of degree one inside the grid for the DeeSse than for the OBS realizations.

D

The previous section highlights the ability of the method to distinguish realizations by focusing on the static connectivity through the connected components. This section discusses some aspects of the analysis process.

About the indicators

All the indicators proposed here rely more or less directly on the connected components. Some of them are classical, such as the facies proportion, but as highlighted on figure 9 the facies proportion is not enough to characterize the static connectivity. New indicators are introduced here compared to previous studies on connected components [Deutsch, 1998, De Iaco and[START_REF] Iaco | Validation Techniques for Geological Patterns Simulations Based on Variogram and Multiple-Point Statistics[END_REF]. Some indicators lead to better characterize the component organization, such as the traversing component proportion or the component density. Other indicators aim to better characterize the component shape, such as the sphericity. Using skeletons is also a new feature to compare realizations. The node degree proportion appears to give many details about the connectivity. The branch tortuosity has been less useful for the studied case, with a poorer discrimination of the realizations. This is due to the parameterization of the skeletonization, which favor the topology at the cost of the geometry of the skeletons. The use of multiple-point histograms as indicators in a method similar to [START_REF] Tan | Comparing Training-Image Based Algorithms Using an Analysis of Distance[END_REF] shows a ranking close to that with the connected component indicators. However, they do not characterize the realizations in the same way. The multiplepoint histograms of the finest multi-grid or multi-resolution level characterize in details the shape of the sedimentary bodies. The shape indicators are global measures over a whole connected component. As connected components can have variable shapes due to the sedimentary bodies intersections, being able to characterize more finely the component shape is an interesting asset. From this point of view, the multiple-point histograms could bring further information on the connected component shape.

However, the multiple-point histograms do not measure the static connectivity: they compare the patterns between the images, but not really the relationships between the patterns. The study of the coarsest multi-grid or multi-resolution levels attempts to look at the large scale behavior of the sedimentary bodies. But many details are lost in the process, what justifies the lower weights for these levels in the dissimilarity from multiple-point histograms [START_REF] Tan | Comparing Training-Image Based Algorithms Using an Analysis of Distance[END_REF]. And it still not characterizes the static connectivity. From this point of view, the skeletons describe more precisely the large-scale behavior of the components and their connectivity.

About indicator comparison

As stated in the previous section, a single indicator is not enough to fully characterize the static connectivity. Comparing several indicators lead to more relevant information about the realizations and how much they differ from the viewpoint of connectivity. Comparing realizations on grid of different dimensions leads to issues non-addressed by previous studies [Deutsch, 1998, De Iaco and[START_REF] Iaco | Validation Techniques for Geological Patterns Simulations Based on Variogram and Multiple-Point Statistics[END_REF]. A correction on the two most affected indicators is not sufficient to compensate for different grid dimensions. Sampling the images appears to be more efficient, and also helps to analyze the connectivity heterogeneity within the images. The question of the sampling representativeness remains to be explored.

Using a metric is very useful, because it gathers all the indicator values into one dissimilarity value and facilitates the comparison of the realizations and the analysis. [START_REF] Tan | Comparing Training-Image Based Algorithms Using an Analysis of Distance[END_REF] already used such process with multiple-point histograms. We have applied a similar principle to connected components, gathering many indicators into values easier to analyze. The introduction of a heterogeneous metric gives the opportunity to gather indicators of different types and further improves the method ability to characterize the realization static connectivity. At the end, the dissimilarities distinguish the realizations from different methods and parameter values, but also characterize the static connectivity variability between the realizations of a given method and parameter values.

Adding weights to the indicators in the metric computation means more flexibility for the user. Indeed, not all the indicators are significant to all the applications. For instance with a flow simulation purpose, the unit component proportion is not necessarily significant due to a fewer impact of the unit-volume component on the flow than channels. But such weights remain optional. In the case study, we did not discriminate the indicators with weights, because we wanted to study the information provided by all the indicators on the realizations. Studying the indicator values after the dissimilarities remains essential to better understand the static connectivity of the realizations.

About the skeletonization method

Skeletons enable to better characterize both the geometry and the topology of connected components. However, the skeletonization method influences both the geometry and the topol- ogy of the resulting skeletons. Among all the skeletonization methods, Cornea et al.

[2007] distinguish the thinning-based method as the method with the best control on the skeleton connectivity. This section aims at comparing the result of a thinning-based method with the method introduced in section 1.3 based on slicing the grid and computing the connected components, denoted as the slicing-based method. The slicingbased method used hereafter is the algorithm defined by Lee et al. [1994] and implemented in the geomodeling software Gocad by [START_REF] Barthélemy | Simulation of Anastomotic Karst Networks[END_REF].

The thinning-based method appears to perform better in two dimensions than the slicing-based method. But in three dimensions it tends to generate many small-scale loops (figure 18) which perturb both the topology and the geometry of the skeletons. The primary goal of the skeletons is to better characterize the large-scale topology -and possibly the geometry -of the connected components. The skeletons from the thinning-based method seem too perturbed to help in that characterization. The slicing-based method on the other side does not necessarily capture those small-scale elements due to the slice size. A large slice size may not capture the small components or all the component irregularities, but this is compensated in some way by the other indicators, in particular the shape indicators. Moreover, the thinning-based method tends to generate skeletons with many nodes, which are heavy to manipulate. The slicing-based method does not have the same issue when using quite high slice thicknesses. This aspect can be essential when dealing with several hundreds of images.

All this leads to favor the slicing-based method in this work. Some aspects still need to be explored, such as the impact of the slice size. But many more skeletonization methods exist, even if skeletonizing three-dimensional shapes is an open debate. Further work could be done to study other methods and the topology and geometry of resulting skeletons.

About MDS methods and accuracy

We rely on the Scaling by MAjorizing a COmplicated Function as multidimensional scaling method to represent the dissimilarities. The SMACOF significantly facilitates the dissimilarity analysis. However, the dimensionality reduction makes the MDS representations imprecise, and the distances between the points tend to differ from the dissimilarities.

Thus, the MDS is not a simple visualization tool and can impact the analysis. This can be illustrated by comparing the MDS representation from the classical scaling (figure 19) and that from the SMACOF (figure 20) to analyze the dissimilarities between the original realizations and the original training image (and not the samples). Normally, the TI should stand from the realizations (see figure 7). But the classical scaling puts the TI close from the OBS and IMPALA realizations. Only the point position confidence shows that the TI position is wrong on the representation. The SMACOF representation separates more clearly the TI from the other images. Moreover, if the global relationships between the realization categories are similar between the two representations, the relative position of the images can be significantly different. This is clear with the TI, but also with other images (table 3). This appears more largely on the Shepard diagram, with a better coefficient of determination r 2 for the SMACOF than for the classical scaling. The classical scaling tends here to decrease the dissimilarities. As a result, the realization ranking can differ between the dissimilarities and the MDS distance (table 3). Thus, analyzing the sole MDS representation can lead to erroneous interpretations.

The choice of the MDS method is significant, so as the choice of the number of dimensions. We have privileged two-dimensional MDS representations for the sake of visibility, but threedimensional representations would be worth testing. In any case, the MDS representation should always be cautiously studied and its misrepresentation of the dissimilarities should be kept in mind. From this point of view, the heat map facilitates the analysis of the real dissimilarity values. Analyzing a single row or column of the dissimilarity matrix -so comparing an image with all the others -is as easy to analyze as a MDS representation, but only on a subset of the images.

As the MDS facilitates the dissimilarity analysis, the dissimilarity simply makes the indicator analysis easier. After looking the MDS representation, it is essential to go back to the dissimilarity values to validate the observations. Similarly, studying the indicator values validates the observations and helps to further understand the difference of connectivity between the images.

Impact of the connected components on the ow

Facies heterogeneity shapes the fluid flow in subsurface. Thus, the petrophysical property simulation, directly correlated to the facies modeling, constitutes a preliminary step to further simulate flows. Being able to assess the static connectivity at the beginning of the workflow could constitute a real advantage in term of resources and time. It is also a way to ensure a better geological consistency, which in itself allows a better integration of field observations and measurements (seismic information, well data, etc.).

It would be interesting to apply the method on more detailed facies models than those of the case study. For instance, both channel and levee deposits are often heterogeneous, regardless of the sedimentary environment, with porous deposits more or less nested between flow barriers [e.g., [START_REF] Hubbard | Seismic geomorphology and sedimentology of a tidally influenced river deposit, Lower Cretaceous Athabasca oil sands, Alberta, Canada[END_REF][START_REF] Hansen | Genesis and character of thin-bedded turbidites associated with submarine channels[END_REF]. Such flow barriers take the shape of mudstone drapes along the channel margins, of margin failure deposits, of channel abandonment deposits, etc. They can have a significant impact of the fluid flow [e.g., [START_REF] Labourdette | Threedimensional modelling of stacked turbidite channels in West Africa: impact on dynamic reservoir simulations[END_REF][START_REF] Pranter | Analysis and modeling of intermediate-scale reservoir heterogeneity based on a fluvial pointbar outcrop analog, Williams Fork Formation, Piceance Basin, Colorado[END_REF][START_REF] Alpak | The impact of fine-scale turbidite channel architecture on deep-water reservoir performance[END_REF][START_REF] Issautier | Modelling of CO2 Injection in Fluvial Sedimentary Heterogeneous Reservoirs to Assess the Impact of Geological Heterogeneities on CO2 Storage Capacity and Performance[END_REF] and on the aquifer compartmentalization, with sometimes important consequences when they are ignored [see for instance [START_REF] Gainski | Turbidite reservoir compartmentalization and well targeting with 4D seismic and production data: Schiehallion Field[END_REF], in a oil exploitation context]. Our method could clearly help to distinguish between several images from their differences in static connectivity -such images including realizations from different methods and/or parameter values or referring models. From this perspective, the case study shows the method ability to identify the simulation methods that produce subsurface models consistent with the static connectivity of a referring model. Such approach is particularly adapted for fluvial and turbiditic channelized environments where channels tend to form high connectivity corridors, leading to channelized flow path. However, the static connectivity of a sedimentary body is not always representative of the flow behavior. For example, flow channeling can also emerge from non-channelized but highly heterogeneous bodies [e.g., Park et al., 2008, Fiori and[START_REF] Fiori | On Preferential Flow, Channeling and Connectivity in Heterogeneous Porous Formations[END_REF]. It highlights the dependence of the hydrodynamic connectivity on many parameters: the permeability contrasts between the different media, the internal heterogeneity of each media, etc. Flow simulations then require to assign the petrophysical properties to each facies -usually with geostatistical methods [e.g., [START_REF] Deutsch | GSLIB: Geostatistical Software Library and User's Guide[END_REF]. Our metrics obviously do not anticipate the results of such procedures and, thus, just measure the consistency of the facies simulations in term of static connectivity.

Depending on the studied environment, the reproduction of the static connectivity could be secondary and one could directly work on hydraulic connectivity through the correspond-ing properties. If reproducing the static connectivity does not guarantee to reproduce the exact hydraulic connectivity, it remains a step toward a better integration of geological information and knowledge in the physical description of the media. Our method provides a simple and objective basis for the comparison of large sets of realizations from this static connectivity point of view.

C

This work develops an analysis process to compare several categorical images in terms of static connectivity. The process relies on the computation of dissimilarities from a given set of indicators. The indicators rely on connected components. They cover various aspect of the realizations, from their facies proportions to the size of their components. Their analysis in the presented case highlights their ability to distinguish between realizations simulated with different methods.

The multidimensional scaling summarizes all this informa- tion through low-dimensional representations of the dissimilarities. It is a powerful and practical visualization tool to get a first idea of the relationships between the images. But it should be kept in mind that MDS representations provide only a partial view of the dissimilarities. Moreover, the representation and its quality are method-dependent. The MDS representation should so always be interpreted cautiously. Analyzing the dissimilarity matrix shall always be favored rather than analyzing the MDS representations. The heat map helps to have an easier check on the MDS interpretation. Then an analysis should be done on particular rows or column of the matrix depending on the studied cases. This gives access to simple but powerful visualization tools that work on the real dissimilarities.

A detailed analysis shows the consistency of the indicator behaviors. Indicators of the connected component shape are interesting for their simplicity. But their sensitivity to the grid size is a real issue. A correction is proposed here, but it does not remove completely the error, and analyzing image samples is a better choice. Multiple-point histograms also seem a valid option to analyze more locally the component shape.

Skeleton-based indicators appear to be promising. The node degree proportions provide a simple way to compare connectivity structures. Further work should be done on the skeletons to better characterize both the connected component topology and geometry.

Thus, the present work shows that analyzing the behavior of various indicators needs to be continued. More case studies should be considered. Channels are known to have a heterogeneous filling, with shale acting as a flow barrier. It could be interesting to model and work on the component formed by the reservoir facies rather than by the sedimentary objects they belong to. Starting with real data would also be beneficial.

If the proposed indicators focus on the static connectivity, an interesting perspective would be to explore their link with the dynamic connectivity.

A

This work was performed in the frame of the RING project at Université de Lorraine. We would like to thank the industrial and academic sponsors of the Gocad Research Consortium managed by ASGA for their support and Paradigm for providing the SKUA-GOCAD software and API. We also thank Sebastien Strebelle for its helpful remarks which influenced this research. We would like to thank the reviewers, including Aldo Fiori, for constructive comments which helped improve this paper.

A A D -

The following appendix describes more precisely all the indicators, starting with the indicators computed on the connected components and ending with those computed on the skeletons.

A.1 Global indicators

The global indicators do not characterize a particular connected component but an ensemble of connected components. They provide one value per facies.

A.1.1 Facies proportion

The facies proportion is of major importance in reservoir modeling considering its influence on porosity and permeability simulation. The proportion p of a facies f is defined as:

p f = n f n t (12)
with n f the number of cells of facies f and n t the total number of cells.

A.1.2 Facies adjacency proportion

The adjacency proportion between the facies quantifies the spatial relationships between those facies. The proportion p a of a facies f adjacent to a facies j is defined as:

p a f , j = n f , j n f ,t (13)
with n f , j the number of cells of facies f adjacent to the facies j and n f ,t the total number of cells of f adjacent to a facies different from f .

A.1.3 Facies connection probability

The connection probability Γ quantifies the global connectivity of a given facies f . It is calculated as the proportion of pairs of connected cells among all the pairs of cells of the considered facies [START_REF] Renard | Connectivity metrics for subsurface flow and transport[END_REF]:

Γ f = 1 n 2 f N f i=1 (n i f ) 2 (14)
with n f the total number of cells of facies f , N f the number of connected components of the facies f and n i f the number of cells of the connected component i associated to facies f .

A.1.4 Connected component density

The density ε of connected components for a facies f expresses the number of components per image cells:

ε f = N f n t ( 15
)
with N f the number of connected components of the facies f and n t the total number of cells.

A.1.5 Unit connected component proportion

The unit connected component proportion quantifies the number of connected components constituted by only one cell. The creation of such bodies may correspond to some erroneous and unwanted small scale noise. In that case, it is better to avoid them for the shape indicator computation (see A.2). The proportion p u of unit volume connected components of a facies f is:

p u f = N u f N f (16)
with N u f the number of unit volume connected components of the facies f and N f the number of connected components for the same facies.

A.1.6 Traversing connected component proportion

Connected components can have various behaviors against the grid: they can be in the middle of the grid without contact with a border (figure 2, connected component 3), along one border, along two adjacent borders (figure 2, connected component 2) or traversing the grid completely from one border to the opposite border (figure 2, connected component 1). These last components having a main impact on the flow, their reproduction is of prime importance. The proportion p c of traversing connected components of facies f is defined as the proportion of connected components that connect one border to the opposite one:

p c f = N c f N f -N u f ( 17
)
with N c f the number of traversing components for the facies f , N f the total number of components for the facies f and N u f the number of unit components for the facies f .

A.2 Shape indicators

Shape indicators all give one value per component. The arithmetic mean of those values provides a value of the indicator for a given facies. This makes the indicator comparison easier. The unit-volume connected components (see section A.1.5) give little information on the shape of the most important connected bodies considering flow circulation: the channels. But they can have a significant influence on the averaged shape indicators (section A.2). They are so ignored in the average computation.

A.2.1 Number of connected component cells

The number of cells n g f within a connected component g of a facies f represents the size of a connected component.

A.2.2 Box ratio

The box ratio is based on the axis-aligned bounding boxes of the components (figure A.1). It is related to the tortuosity of a component and to its orientation relative to the grid axis. The box ratio β of the connected component g of the facies f is expressed between 0 and 1 using the following formula:

β g f = n g f n b f (18)
with n g f the number of cells of the connected component g of the facies f and n b f the number of cells occupied by its axisaligned bounding box.

A.2.3 Faces/cells ratio

The faces/cells ratio is equivalent to the surface/volume ratio, which compares the surface area of a component with its volume. [START_REF] Deutsch | Fortran programs for calculating connectivity of threedimensional numerical models and for ranking multiple realizations[END_REF] uses it as a measure of the tortuosity of the components, but it is also affected by their roughness. The faces/cells ratio ζ of the connected component g of the facies f is expressed as:

ζ g f = m g f n g f (19)
with m g f the number of faces along the surface of the connected component g of the facies f and n g f its number of cells.

A.2.4 Sphericity

The sphericity φ compares the shape of a connected component g for a facies f with a sphere [START_REF] Wadell | Volume, shape, and roundness of quartz particles[END_REF]. It ranges between 0 and 1, where 1 corresponds to a sphere, and is expressed as:

φ g f = 36π (n g f ) 2 (m g f ) 3 (20)
with n g f the number of cells of the connected component g of the facies f and m g f its number of faces along its surface.

A.2.5 Skeleton indicators

Skeletons are one-dimensional representations of the connected components, on which measures of geometry and topology are done.

A.2.6 Node degree proportion

The node degree proportion is based on the number of segments connected to a given node (figure A.2). It is related to the topology of the skeleton and gives access to a deeper study of the connectivity of the structures of interest. The proportion p n of nodes of degree n over all the skeletons for a given facies f is expressed as:

p n f = o n f o f (21)
with o n f the number of nodes of degree n for all the skeletons of the connected components of the facies f and o f the total number of nodes of the skeletons of the connected components for the facies f . Here we separate into two different classes the nodes of degree one along a border of an image and those in the middle of the image. The first ones relate to component terminations due to the limited size of the image. The other ones may relate to a bad body reproduction, for instance with channels that should be continuous and without dead-ends. The proportions for the different degrees give a histogram of node degrees.

A.2.7 Inverse branch tortuosity

A branch is a part of a skeleton defined as an ensemble of segments that link nodes of degree 2 and delimited by two nodes of degree different from 2 (figure A.2). The inverse tortuosity t quantifies how twisted a branch b is, with values ranged between 0 and 1. This measure is related to the geometry of the skeleton. It is expressed for a branch b belonging to a skeleton s of a component of a given facies as:

A B I -

The following appendix describes the corrections applied to the connected component density and the number of component cells when the analysis implies different grid dimensions.

B.1 Correction for the connected component density

When some connected components of a facies are traversing (see A.1.6), reducing the grid dimension along the component If image 2 has no connection with a degree higher than 3, such as image 1, it displays some dead-ends highlighted by nodes with a degree one in the middle of the grid. This difference emphasize mis-reproductions of the channel connectivity. In image 1, l represents the curvilinear length of the branch and d the distance between its extremities. Those measures are used to compute the tortuosity, which can also differ between the images.

traversing extension can change the component size but not necessarily their number. For instance changing the grid size along the channel orientation does not change the number of channels. But it changes the grid volume and so does the density. Comparing two images of different sizes can affect the density whereas there is nothing wrong about the quality.

To overcome this aspect, the connected component density (formula 15) is replaced by:

ε f = N f ϑ t (23)
with N f the number of connected components of the facies f and ϑ t is a corrected number of grid cells corresponding to a mix between the dimensions of each compared images:

ϑ t = ϑ x × ϑ y × ϑ z (24)
with ϑ x , ϑ y and ϑ z respectively the corrected number of cells along the axis x, y or z. Those corrected numbers of cells are chosen depending on the smallest axis size among all the compared images and the size of the connected components.

For an ensemble S of images to be compared, we note m i the minimal size along the axis i over the images in S:

m i = min s∈S (n i s ) i ∈ x, y, z (25)
with n i s the number of cells along the axis i for the image s. Determining the corrected numbers of axis cells for a given image r of S depends on the size of its connected components:

ϑ i = m i if ∃ g ∈ r, n i g ≥ m i n i r otherwise i ∈ x, y, z (26)
with g a connected component of r, n i r the number of cells along the axis i for the image r and n i g the number of cells along the axis i for the component g.

For instance in (figure A.1), the component 1 of image 1 is larger than image 2 along the y axis. On the other hand, no component of that facies is larger than image 2 along the other axes. The corrected volumes for the image 1 ϑ t,1 and for the image 2 ϑ t,2 are so:

ϑ t,1 = n x 1 × n y 2 × n z 1 ϑ t,2 = n x 2 × n y 2 × n z 2 ( 27
)
with n x 1 , n y 1 and n z 1 the number of cells along the axis x, y and z for image 1 and n x 2 , n y 2 and n z 2 the number of cells along the axis x, y and z for image 2. Image 1 is so considered as truncated to limit the effect of the traversing bodies on the density.

B.2 Correction of the number of connected component cells

The objects are expected to have roughly the same range of sizes in all the images. But this expectation is no longer relevant in the case of traversing structures (see A.1.6), whose dimensions must vary following the size of the supporting grid.

As for the connected geobody density, the possible difference of image sizes imposes to rescale the number of cells of a connected component.

The rescaling is based on the estimation of the component bounding box size if this component was transfered within the grid of the smallest size. That bounding box size can be defined following a number of cells ϑ b :

ϑ b = ϑ x × ϑ y × ϑ z (28)
with ϑ x , ϑ y and ϑ z respectively the corrected numbers of cells of the bounding box along the axis x, y or z. Those corrected numbers of cells are following the same principles as to find the corrected number of cells for the density. Determining the corrected numbers of axis cells requires the same parameter m i (formula 25):

ϑ i = m i if n i b ≥ m i n i b otherwise i ∈ x, y, z (29)
with n i b the number of cells along the axis i for the bounding box of the component. Finally the corrected number of cells ϑ g f for the component g of the facies f is proportional to the ratio between the number of cells of its corrected bounding box ϑ b and the number of cells of its actual bounding box n b :

ϑ g f = n g f × ϑ b n b (30)
For instance the component 1 on image 1 of figure A.1 has a corrected number of cells ϑ 1 of:

ϑ 1 = n 1 × n x b 1 × n y 2 × n z b 1 n x b 1 × n y b 1 × n z b 1 = n 1 × n y 2 n y b 1 (31)

Figure 1 Figure 2

 12 Figure 1Possible neighborhoods for a given central cell in a regular grid (modified from[START_REF] Deutsch | Fortran programs for calculating connectivity of threedimensional numerical models and for ranking multiple realizations[END_REF]).

3 Figure 3

 33 Figure 3Example of skeletons for the connected components of the figure2. Here the nodes connected to only one segment -the nodes of degree 1 -are all along a grid border. Two nodes of degree three highlight the local disconnections between the channels at the bottom. The connected component at the top has no node of degree higher than two, which shows the complete connectivity of all its cells, even locally.

p

  Number of component cells Number of cells of the grid Sum of all the component values Facies adjacency proportions p a Number of component cells adjacent to a cell of a given other facies Number of cells of the facies adjacent to a cell of any other facies Sum of all the component values Facies connection probability Γ Squared number of component cells Squared number of cells of the facies Sum of all the component values non-unit-component values Box ratio β Number of cells of a component Number of cells of the axis-aligned bounding box of the component Average of the non-unit-component values Faces/cells ratio ζ Number of faces composing a component surface Number of cells of a component Average of the non-unit-component values Sphericity φ Surface area of a sphere Surface area of the connected component Average of the non-unit-component values Skeleton indicators Inverse branch tortuosity t Distance between the extremities of a branch Branch curvilinear length Average of all the branch values Node degree proportions p n Number of node connected to n segments Total number of node for all the skeletons Sum of all the skeleton values

3 .

 3 Dissimilarity visualization and analysis, especially with multidimensional scaling. The first two steps were implemented in a C++ plugin for the SKUA-GOCAD geomodeling software [Paradigm, 2015]. The last step was realized using the software environment for statistical computing R [R Core Team, 2012] with the addition of the R packages SMACOF [De Leeuw and Mair, 2009] and ggplot2 [Wickham, 2009].

Figure 4

 4 Figure 4 Training image and examples of realizations for each category.

Figure 5

 5 Figure 5 View of all the channel connected components within the TI and examples of realizations for each categories. The number in parentheses are the number of connected components of each image.

Figure 6

 6 Figure 6View of all the skeletons of the connected components for the TI and for a realization of each category.

Figure 7 Figure 8

 78 Figure 7 MDS representations of the dissimilarity matrices for the original images (with corrections of the indicators to cancel the effect of the grid dimensions) and samples (of same size). The scree plot for the original images only displays the stress values up to 10 dimensions on 200 possible. The scree plot for the samples only displays the stress values up to 10 dimensions on 2220 possible.

Figure 9

 9 Figure 9Box-plots comparing the facies proportions and facies connection probability for the TI, some TI analogs and the OBS realizations.

Figure 10

 10 Figure10MDS representation of the dissimilarities between the samples of the case study generated using SMACOF and validation graphs. The scree plot only displays the stress values up to 10 dimensions on 820 possible.

Figure 11

 11 Figure 11 Heat map representation of the dissimilarity matrix computed based on the samples of the case study.

Figure 12 Figure 13 Figure 14

 121314 Figure 12Box-plots comparing the range of indicators computed on the channels for the different categories, except the node degree proportions.

Figure 15

 15 Figure 15Mean node degree proportions of the levee skeletons for each category. The error bars display the minimum and maximum proportions. The first node degree 1 corresponds to the nodes of degree one along a grid border. The second node degree 1 corresponds to the nodes of degree one inside the grid.

Figure 16

 16 Figure 16 MDS representation of the dissimilarities between the images of the case study generated using SMACOF and validation graphs. The dissimilarities are based on the multiple-point (MP) histograms of the images. The scree plot only displays the stress values up to 10 dimensions on 400 possible.

Figure 17

 17 Figure 17 Box-plots comparing the realizations of each method with the TI. The dissimilarities depend on different indicators in each box-plot. Only the multiple-point histograms compare the original image directly, with one dissimilarity value with the TI per image. The other categories are based on image samples, and have several dissimilarity values with the TI per image. These values are averaged to obtain a single dissimilarity value with the TI.

Figure 18

 18 Figure 18Skeletons from the channels of the realization DeeSse 30 generated with two methods: the slicing method used in the case study and the thinning method.

Figure 19 Figure 20

 1920 Figure19MDS representation of the dissimilarities between the images of the case study generated using classical scaling and validation graphs. The scree plot only displays the eigenvalues up to 10 dimensions on 400 possible.

1 Figure A. 1

 11 Figure A.1 components of figure 2 represented with their bounding boxes (image 1). Image 2 is another image in a different grid.The two grids have different dimensions along the x axis, which can influence both the number of connected components and the number of grid cells. But they also have different dimensions along the y axis, along which the channels are oriented. Changes along the y axis influence the number of grid cells, but not the number of channels. In that case, the comparison of the density of the two images is biased. The same principle stands for the number of cells of each components. Differences between the number of cells for the component 1 of each image mainly come from the difference in grid dimension along y, not from real differences of connectivity. n x 1 and n y 1 represent the number of cells along the axis x and y for image 1, n x 2 and n y 2 the number of cells along the axis x and y for image 2 and n x b 1 and n y b 1 the number of cells along the axis x and y for the axis-aligned bounding box b of the connected component 1. m x and m y are the lower sizes between the two grids for the axis x and y (see appendix B for more details).

  s the curvilinear length of the branch b of the skeleton s and d b s the Euclidean distance between the two extremity nodes of the same branch. The values for all the branches of all the skeletons related to an image are averaged to obtain one indicator value for that image.

1 Node of degree 2 Node of degree 3 Figure A. 2

 1232 Figure A.2 Example of skeletons for the connected components of the figure A.1.If image 2 has no connection with a degree higher than 3, such as image 1, it displays some dead-ends highlighted by nodes with a degree one in the middle of the grid. This difference emphasize mis-reproductions of the channel connectivity. In image 1, l represents the curvilinear length of the branch and d the distance between its extremities. Those measures are used to compute the tortuosity, which can also differ between the images.

  The resampling step selects the first pattern with a distance lower than a given threshold. The training image is the TI and the set of parameters is given in table C.1) in the appendix. IMPALA: 100 MPS realizations simulated with the method IM- PALA [Straubhaar et al., 2011, 2013]. Contrary to the DeeSse, IMPALA still computes the conditional probabilities during the simulation. To improve the efficiency of this computation, the method stores the training image patterns in a list. The training image is scanned once at the beginning and the list is used instead during the simulation. The training image is the TI and the set of parameters is given in table C.2) in the appendix.

	OBS: 100 object-based realizations simulated with the same
	method and parameters used to simulate the TI (ap-
	pendix C, table C.3).
	SIS: 100 sequential indicator simulation realizations simulated
	using variograms based on the facies in the TI (appendix
	C, table C.4).
	Analog: 100 object-based realizations simulated with the same
	method and parameters used to simulate the TI (ap-
	pendix C, table C.3).
	The simulation grid contains four sets of images:
	DeeSse: 100 MPS realizations simulated with the DeeSse im-
	plementation [Straubhaar, 2011] of the direct sampling

TI: One object-based realization simulated using the objectbased method of the software Petrel [Schlumberger, 2015] (see appendix C, table C.3, for the simulation parameters). method [Mariethoz et al., 2010]. Contrary to more traditional MPS methods, the direct sampling bypasses the conditional probability computation and resamples randomly the training image. It relies on the compatibility -measured with a distance -between the conditioning data and the patterns scanned in the training image.

Table 2

 2 Set of indicators used for the case study. The indicator definitions are summarized in table 1 and more detailed descriptions are in appendix A.

	Category Indicator	Symbol Weight
		Facies proportion	p	1
		Facies adjacency proportion	p a	1
	Global	Facies connection probability	Γ	1
	indicators	Connected component density	ε	1
		Unit connected component proportion	p u	1
		Traversing connected component proportion p c	1
		Number of connected component cells	n	1
	Shape	Box ratio	β	1
	indicators	Faces/cells ratio	ζ	1
		Sphericity	φ	1
	Skeleton	Node degree proportions	p n	1
	indicators	Inverse branch tortuosity	t	1

Table 3

 3 Comparison of dissimilarities and distances between the TI and DeeSse realization 12 and 76 for to MDS methods.

	Compared images	Dissimilarity Classical scaling distance SMACOF distance
	TI -DeeSse 12	1.579	0.596	1.414
	TI -DeeSse 76	1.358	0.854	1.265
	DeeSse 12 -DeeSse 76	0.905	0.259	0.203

the numbers of cells of its axis-aligned bounding box along the axis x, y or z and n x 2 , n y 2 and n z 2 the numbers of cells of image 2 along the axis x, y or z. This reflects the impossibility to make that component fit into image 2. Its size needs so to be rescaled to be compared with objects of image 2.

A C S