A Geostatistical Approach for Grid-Independent Geomodeling in Complex Tectonic Environments

Abstract

It is of major importance in geological modeling to account for the geometry of the volume to be modeled. Several methods are available for introducing tectonic deformations or paleo-topographic surface shapes in the geological modeling process. This chapter proposes a synthetic overview of these methods based on geometric deformations or specific geostatistical models. The first approach consists of distorting the modeling grid with a more or less sophisticated unflattening algorithm. Other approaches consist of using a geostatistical algorithm that can take into account geometric deformations when populating the grid with properties. Two algorithms of this type are detailed: local geostatistics (LGS) and multiple-point statistics (MPS). With such algorithms, which can work with any type of grid, the flattening step can be skipped. Special attention is paid to the possibility of modeling the average value of a property, instead of a point value, with the three approaches. It is shown that, in such a case, it is better to perform the modeling in regular grids. The different methods are not exclusive to each other and can be combined, offering a wide range of modeling possibilities, assuming that the inference of technical parameters has been made properly.

1. Introduction

It is of major importance in geological modeling to account for the geometry of the volume to be modeled.

For example, in sedimentary environments, correlations are based on timelines identified in wells/drillholes from sequence stratigraphy analysis (maximum flooding surfaces, typically). To be realistic, numerical geological models must be built following such reference surfaces that correspond to the horizontal plane at the time of deposition. It is obvious that the initial shape of such reference surfaces can be altered by tectonic deformations, which, therefore, must be accounted for in the modeling process.

Another common example is the case of orebodies coming from lateritic alteration of soils (manganese or nickel orebodies, typically). The shape of the mineralization is strongly influenced by the shape of the paleo-topographic surface, which must be considered to properly model the spatial distribution of metal grade.

Several methods are available for introducing tectonic deformations or paleo-topographic surface shapes in the geological modeling process. This chapter proposes a synthetic overview of these methods based on geometric deformations or specific geostatistical models, with a focus on the pros and cons of each approach. Geometric approaches will be briefly presented, but special attention will be paid to the geostatistical models, which allow following and restoring deformations without distorting the modeling grid.

The conditions of application of each method will be studied, and the consequences of modeling grid characteristics and the change of support possibility will be detailed because of their important practical impact.

2. First approach: distorting the modeling grid

2.1 Basic principle

The most common approach for integrating deformations in a geological model consists of distorting the modeling grid.

In the simplest case, the shape of a user-defined reference surface is used as a guide for flattening the modeling grid, as illustrated in Figure 1.

Figure 1.

Geological modeling with flattening of a reference horizon.

The reference surface (the undulated surface in the upper left quarter of Figure 1) is flattened to become horizontal, and the original modeling grid is distorted accordingly. Then, the geological modeling is accomplished along directions parallel to the reference surface. The final unfolding operation sends the ensuing results back to the initial grid in which objects that were horizontal in the flattened grid are now folded and follow the reference surface shape.

It is a very efficient approach that is compatible with any type of original “structural” grid. It can be perfectly regular and cartesian, or cartesian and irregular, or it can be a Corner-Point Geometry (CPG) grid or even an unstructured one.

Several flattening options allow accounting for various geological correlation configurations corresponding to different depositional conditions, as illustrated in Figures 2 and 3 coming from [1].

Figure 2.

Flattening options parallel to a reference surface.

Figure 3.

Flattening options proportional to two reference surfaces.

It must be noted that CPG and unstructured grids made of cells with various shapes allow for more accurate modeling of complex deformations than regular rectangular grids. For that reason, they are often preferred by geologists. However, in many cases, the geological models are used as a basis for further technical operations, like mine planning or fluid flow modeling, which often require grids to be as regular as possible. This may create some difficulties and loss of accuracy when transferring the geological model from geologists to engineers.

2.2 Some advanced techniques

Flattening a grid with reference to one or two horizons is efficient, but standard algorithms are unable to handle complex geometric shapes that occur in compressive contexts. In such a case, the reference surfaces may be folded in a way that leads to defining several depths at a given location, as shown in Figure 4.

Figure 4.

Complex case for flattening.

In this kind of geometric configuration, which is common in compressive tectonic environments, standard flattening algorithms fail. Several sophisticated mathematical theories and algorithms have been developed to overcome these limitations. Many theoretical developments came from the domain of geosciences, like the Geochron Model [2] or automatic procedures detailed in [3]. Some others come from various technical domains and have been adapted to geosciences [4]. Most of these techniques are implemented in major integrated software platforms for the oil and gas industry. Many methods come from computer-aided design (CAD) and have been implemented in major integrated packages for the mining industry. Advanced algorithms recently implemented in commercial software can be found in [5, 6].

2.3 Change of support issues

The geometrical approach described above proves to be inappropriate when it is important to ensure a constant (or almost constant) cell volume in the whole grid. This constraint typically occurs when a change of support operation is needed, which is common in the mining industry. It is reminded that the “change of support” operation consists of deducing the histogram and variogram of the average value of a variable in a block from those of the variable at points. Usually, points correspond to measured data in composite samples, and blocks correspond to volumes of tens or hundreds of cubic meters.

Similar constraints occur when operations requiring limited and smooth spatial variations of cell volume and cell shape will follow geological modeling, such as flow simulations made with finite-element flow simulators. This is common in the oil and gas industry or in hydrogeology.

In such cases, sophisticated grids, highly distorted to perfectly follow the geometry of the geological layers in complex tectonic environments, may cause some trouble in further applications based on the geological model. In such grids, cells can be distorted. The high level of geological accuracy may be partially lost in subsequent steps, in particular in properties upscaling from high-resolution geological grids to reservoir grids used for flow simulations. Several methods have been developed to overcome this difficulty.

The purpose of this chapter is to present different geostatistical approaches that address this problem. They are based on advanced geostatistical algorithms used for filling grids with properties, which can account for geometric constraints, even if the modeling grid is perfectly regular and rectangular. These algorithms are compatible with any type of grid if they are used in point mode (modeling of the value of a variable at a given point), but focus will be put here on their relevance to modeling the average value of a variable on a volume.

3. Second approach: using local geostatistics

3.1 Accounting for local deformations with geostatistical algorithms

In classical geostatistical applications, one can perform estimations (kriging, for example) or simulations using a variogram model. It may have anisotropies, which are generally defined globally. Local geostatistics consists of considering spatial variations: in 3D, of the variogram parameters like the sill, the range, or the anisotropy ellipse orientation (see Figure 5) that can be of interest in some circumstances. This method is used for enhancing mapping in large-scale domains [7, 8].

Figure 5.

Local orientations of the variogram anisotropy ellipse.

Working with varying orientations of the anisotropy ellipse makes it possible to account for tectonic deformations in estimations and to reproduce deformations in simulation without having to distort the grid. An academic example is proposed in Figure 6 to illustrate this property. In this figure, two models are displayed. Both are simulations in the same regular rectangular grid. On the top panel of the figure, the first model is built using a strongly anisotropic variogram, with a global horizontal elongation oriented along the X-axis of the grid. On the bottom panel, the second model was built using local variations of the variogram anisotropy ellipse orientation, which were calculated from the undulated surface displayed below the grid. The same ratio between ranges as in the first model was used. It can be noticed that the simulation results are reproducing these undulations, even if the grid is purely cartesian with squared cells.

Figure 6.

Reproduction of deformations with LGS in a regular grid.

A zoom on another realization, shown in Figure 7, highlights the deformation in the regular grid.

Figure 7.

Reproduction of deformations with LGS in a regular grid (enlarged view).

3.2 Application with continuous and categorical variables

Local variogram parameters can be used with continuous or categorical variables, as shown in Figures 8 and 9.

Figure 8.

Reproduction of deformations with LGS in a regular grid (examples with facies proportion).

Figure 9.

Reproduction of deformations with LGS in a regular grid (example with facies).

In Figure 8, three examples are proposed, which all correspond to the expected local proportion of a given geological facies. The first image on the left is the reference, with no deformation. The other two images show the impact of the use of local orientation of the variogram anisotropy ellipse calculated from the surfaces below the displayed grids. One corresponds to undulations of small amplitude and high frequency, the other to undulations of high amplitude and low frequency.

In Figure 9, an example of an application to a geological model is displayed. Each color corresponds to geological facies, where the undulated surface used for calculating local variogram parameters has high amplitude and low frequency. The following comments can be made.

• First, since the method relies on local variogram parameters, it can be used only with variogram-based facies simulation methods. Typically, Sequential Indicator Simulation (SIS) and PluriGaussian Simulation (PGS), detailed in [9] and in [10], are compatible with this approach and available in commercial packages (if developers have included the use of local variogram parameters).

• For PGS, the simulation method used for generating the intermediate Gaussian random functions that will be a threshold must be compatible with local variogram parameters. Not all the simulation techniques for Gaussian random functions have this capacity, like the Turning Bands Simulation (TBS) method [10]. On the contrary, the simulation method based on the use of Stochastic Partial Differential Equations (SPDE) detailed in [11] does.

• It is critical to take care of the consistency between the proportions distribution and the facies distribution. Both SIS and PGS require a facies proportion model that can be global (stationary case) or local (non-stationary case). Local deformations of geological bodies in a stationary context (at the scale of the studied zone) are possible. In a highly non-stationary context, where local facies (or rock-types) proportions are used, the deformations should be accounted for in both the proportion model and the facies model to ensure consistency. Intermediate configurations can be encountered, which require user assumptions based on a geological analysis.

3.3 Inference of local variogram parameters

Local geostatistics can model tectonic deformations efficiently, but local variogram parameters must be well-calibrated to the geological context to ensure realistic results, usable in practice. For this purpose, local orientations of the anisotropy ellipse must be calculated using geology as a guide. The same two cases discussed for the grid flattening operation can be considered:

1. The anisotropies are parallel to a reference surface.

2. The anisotropies are guided by top and bottom surfaces. At the top of the domain, the anisotropy ellipse orientations are parallel to the top surface. At the bottom of the domain, they are parallel to the bottom surface.

From the gradients calculated at each point of a surface, it is possible to define a local rotation matrix that defines the local orientation of the anisotropy ellipse of the variogram. Results can be stored in a dedicated rotation object, which can be used later for kriging or simulation, as shown in Figure 10. The orientation displayed in the figure corresponds to a simple undulation along a single axis in the horizontal plane. The Z-axis component of the rotation vector is alternatively oriented upwards and downwards.

Figure 10.

Variations of rotation vector components.

When the local orientation of the variogram anisotropy ellipse is calculated from local gradients observed on geological horizons, this approach for accounting for deformations is automatically compatible with compressive tectonic contexts. Complex geometric configurations like the one shown in Figure 4 are manageable.

Special attention must be paid to the inference of the variogram that will be used in the geological model. In isotropic environments, there is no particular constraint, but some precautions must be taken in anisotropic environments. The variable orientation of the anisotropy ellipse is applied to a global variogram model obtained by fitting the experimental variogram to data that include tectonic deformations. In anisotropic contexts, this alters the experimental variogram shape and therefore the variogram fitting. It is recommended to calculate the experimental variogram on flattened data, even if further calculations are made directly in a non-flattened stratigraphic grid, to ensure consistency. Applying a flattening operation on a limited number of well data is fast, contrary to applying it on a high-resolution grid.

This approach is powerful, but it has some limitations. In particular, it is unable to account for discontinuities like faults. The method is well suited for continuous and smooth deformations, even when they are geometrically complex, but it cannot take into account the shift of a geobody along a discontinuity. There is no local parameter of variogram characterizing a fault throw, and therefore, the brutal displacement of geological layers or of mineralized bodies cannot be properly modeled with variogram-based methods.

Accounting for deformations using local geostatistics (LGS) is feasible only within folded domains, potentially confined by faults. However, it is not applicable in the presence of any faults characterized by significant displacement.

3.4 Modeling blocks

When considering points, as long as estimation or simulation proceeds on the basis of using uncertain variogram parameters, the modeling grid type (cartesian, CPG, etc..) does not really matter. When the goal is to determine point values at the centroid of cells, the only important parameters are the variogram, the centroid locations, and the locations of the surrounding data (if there are any). The shape and volume of cells around the centroids have no influence on the results. Therefore, this approach is naturally compatible with any type of grid, which can be of interest in the oil and gas industry, hydrogeology, geotechnics, and possibly other disciplines.

If a change of support must be considered, therefore, if the average value of a property in a block must be calculated from point data (e.g., the average grade in a block in the mining industry), there is a real issue. Use of geostatistics often involves employing several change-of-support models, like the discrete Gaussian model (DGM) detailed in [12], for characterizing the distribution of the average grade or accumulation in different sizes of blocks, knowing the characteristics of point data. The discrete Gaussian model can be seen basically as a block model, where the domain is partitioned into small blocks v. Then each sample point is considered as random within its block and conditional on its block value (here the multivariate value of the different elements). The point (multivariate) value does not depend on any other variable, whether they are values of other blocks or other points, even in the same block. These models are very slightly dependent on the orientation of the anisotropy ellipse, and it is possible to account for deformations with LGS, even for blocks, using a model like the DGM. Nevertheless, the following comments can be made.

• DGM is usually applied in the mining industry to determine the properties of additive continuous variables in panels and the selection of mining units (SMU), the size and shape of which are related to production constraints. Panels and SMU sizes are user-defined and constant in a production zone. The shape is usually regular and rectangular, as it is easier to produce regular blocks than irregular ones.

• By assuming that a block is discretized in a user-defined number of points, the DGM can efficiently calculate the grade or accumulation distribution for any block size. Therefore, from a theoretical point of view, it should be possible to rigorously calculate the average of any continuous additive variable in a grid with varying cell sizes, which could be of interest in disciplines other than mining. This would require applying the DGM on a cell-by-cell basis and to store in each cell the corresponding change of support coefficient along with the cell volume. In general, such options are not present in geostatistical and mining commercial software.

• As a consequence of the previous comments, it is recommended to consider a regular cartesian grid when using LGS for accounting for deformation in a block model.

4. Third approach: using MPS

4.1 Principle of MPS

The multiple-point statistics (MPS) technique is an efficient simulation method based on the analysis of an existing training image (TI) using a pattern of values (multiple points). It tends to reproduce the proportions (i.e., the relationships between the facies and the geobodies) from the TI. MPS can take into account conditioning data, local proportions, local homothety, local rotation, and trend variables.

This method does not use the variogram, which is a two-point statistic, but it does account for complex facies relationship characterized by the facies configuration observed on all the points of the pattern in the TI. It is potentially more powerful than variogram-based methods, but a representative training image is necessary. Such a TI is not always easy to find.

Many MPS implementations are dedicated to the simulation of categorical variables, like facies or rock types. In this chapter, all the illustrations are based on the DeeSse algorithm, detailed in [13], which is an advanced MPS implementation developed at the University of Neuchâtel (Switzerland) that is able to simulate both categorical and continuous variables and to manage multivariate scenarios.

The principle of the approach is to mimic a reference or training image (in 2D or in 3D). This image could be an analog if the geological environment is known, or an existing well-developed orebody or hydrocarbon field where many conditioning data are available. It can also be a simple theoretical image, which can be deformed through scaling and rotation operations included in the simulation process.

4.2 Option 1: working with a folded TI

Generating a simulation accounting for tectonic deformation is easy with MPS algorithms as long as the TI includes such deformations. It is an intrinsic property of this kind of algorithm that has been designed to reproduce complex relationships between facies. Tectonic deformations can be thought of as a particular geometrical relationship between geological facies.

It must be noted that the accuracy of the reproduction in simulations of complex relationships between facies or geobodies strongly depends on the complexity of the search pattern used for scanning the TI. The larger the search pattern, the more faithful the simulations are to the TI. The more complex the pattern, the more precise the characterization of the structures in the TI, and the better the result. Therefore, it is recommended to use complex patterns, with branches in many directions and points at different distances, to acquire a rich information set on the structures and to obtain accurate simulation results. The DeeSse algorithm is particularly interesting in this regard. Unlike standard MPS, DeeSse does not create a pattern database. Instead, each time it needs to estimate an unknown value, it directly scans the training image to find patterns that match the known values in the surrounding neighborhood. The matching patterns are not stored for future use; instead, they are identified anew each time an estimate is needed [13]. The local pattern for each simulated point is defined as the N closest informed points (data or already simulated cells), which, in general, allows sampling many directions and many distances. This number of N neighbors fully characterizes, in a very simple way, the accuracy of the simulation results. The larger it is, the better the result, but the longer the calculations.

Simulations with two different patterns are shown in Figure 11. It can be noted that the undulations are already correctly reproduced with a small pattern. It is obvious that if the deformations were more complicated than simple undulations, larger patterns would be required to get an acceptable result. The simulation with the larger pattern in the figure is reproducing the TI with a better accuracy than the first simulation: more details are visible.

Figure 11.

Reproducing folds in a TI for 2 patterns of varying accuracy.

4.3 Option 2: distortion from an unfolded TI

It is not always easy to get a folded TI representative of the geological environment for modeling. In such a case, one can use the image transformation possibilities offered by the MPS technique. Transformations based on rotation vectors, scaling factors, or trends can be included in the simulation process [13], resulting in a simulation that looks like a distorted TI.

In particular, the use of global or local rotation vectors allows for generating folded simulations, even if the TI is unfolded. A simple application of this method is shown in Figure 12. Starting from a TI with flat and horizontal geobodies, a result with tilted geobodies is obtained by applying a global rotation vector. In complex cases, local rotation vectors decomposed into three components, like for the LGS approach (Figure 5), should be considered. The way to get local rotation vectors is the same as in LGS, based on calculations from one or two reference surfaces.

Figure 12.

Using global rotation for generating distortions in a MPS simulation with an unfolded TI.

4.4 Grid types and change of support issue

The MPS technique does not manage the change-of-support and is not compatible with common change-of-support models like DGM. Therefore, the simulation results will correspond to the support of the TI. If the cells of the TI grid correspond to the average value of a given variable in the cell volume, then so will the simulation. Note that it makes sense only for continuous parameters and cells with the same volume; therefore, it is most applicable for regular grids.

Some research results have been obtained, permitting the management of change-of-support effects with this method [14], but they are rarely implemented in commercial software. At the present time, resolving this issue generally remains in the hands of the practitioner.

If just the point value of a variable at the cell center is considered for the TI and the simulation grid, then the algorithm will work with any kind of grid, and for either categorical or continuous variables, as in the case of the LGS approach.

Nevertheless, results must be handled with care. The TI grid and the simulation grid are often allowed to be different in the software implementations of the MPS method. It means that the grid resolution can be different for the TI and the simulation. Even for point simulations, working with grids of different resolutions will lead to distortions that could be unexpected, like in the heterogeneity level or the wavelength of geobody undulations.

4.5 How to find a representative TI

A representative and realistic TI is necessary for getting realistic results with the MPS technique. It must be noted that a folded TI is used for simultaneously characterizing the topological relationship between values of a given variable at different locations and the tectonic deformations. Therefore, the TI must be representative of two independent characteristics: the geological properties of the sedimentary deposit or of the mineralization, on the one hand, and the geometrical deformations on the other. Finding such a TI may be difficult.

In 2D, pictures taken from outcrops or satellite pictures can be very good TIs. In particular, they can provide good representations of deformations at different scales. A rigorous geological analysis of these images and of the available data on the domain to be modeled will govern how representative and realistic a given picture can be.

In 3D, it is more difficult to find appropriate TIs. A possible approach can be the use of an existing model built in a domain of the same tectonic trend. This is possible in the mining industry, where such reliable models built from considerable data are available from mines that have been produced for a long time. These models are good references for domains in the same neighborhood.

In the oil and gas industry, where the fields are generally at different scales than mineralizations in the mining industry, the same approach can be applied, assuming that the geologist is able to confirm that the reference and the field to be modeled are affected by similar deformations (two anticlines, for example).

Sometimes, it can be almost impossible to find an appropriate TI in 3D that is representative of all the characteristics. In such situations, the approach presented in Section 4.3 can be considered. It consists of splitting the complex problem into two simpler issues, as follows:

1. Finding a TI representative of the deposition conditions or of the mineralization process.

2. Characterizing the deformations.

Finding a “flattened” TI representative of a deposit or mineralization is simpler than finding a “full” TI that includes deformations. Existing models on mature fields or old mines are good examples. Conceptual models built from outcrops and/or seismic data can also be considered.

The second step consists of finding which deformations have to be applied in the simulation process. This is done by calculating rotation vectors from geological horizons mapped from well/drill hole data and seismic cubes or sections.

5. Conclusion

It has been shown in this chapter that there are many possible options to account for deformations in geological modeling. Each method is well suited to specific contexts, and all these approaches are more complementary than in competition.

Basically, two technical approaches can be considered:

1. Applying more or less sophisticated, unfolding algorithms to the modeling grid, and then filling the flattened grid and unfolding the grid filled with properties to restore the structural shape. This is a geometrical approach.

2. Modeling the deformations with some of the geostatistical algorithms used for filling the grid with properties directly in the structural grid. This corresponds to purely geostatistical approaches.

5.1 Combining approaches

Geometrical and geostatistical approaches can be combined. Instead of using the appropriate geostatistical algorithm in the structural grid, it can be used in the flattened grid. This allows for modeling multi-phased tectonic deformations. Different examples can be considered:

• Grid unfolding and use of LGS when populating the flattened grid.

• Grid unfolding and use of MPS when populating the flattened grid.

In both examples, a large-scale deformation can be managed by the unfolding operation, while secondary deformations at lower scales can be managed by the geostatistical algorithm used for populating the flattened grid. An example is shown in Figure 13.

Figure 13.

Combining grid flattening and LGS with local orientation of the anisotropy ellipse.

As the geostatistical algorithm is applied to the flattened grid, its parameters must be defined in this particular reference system, which has the following technical implications:

• If local rotation vectors are used in LGS or MPS, they must be calculated from geological horizons to which the same unfolding algorithm as for the grid has been applied.

• If the deformations are included in the TI for MPS, they must correspond to the low-scale deformations only. Here it is assumed that the large-scale deformations managed by the unfolding algorithm are not present in the TI. This can create some practical limitations for the use of this approach.

5.2 Discussion

The different approaches have pros and cons that are summarized in Table 1.

Method	Pros	Cons
Flattening	Flexible and general approach Can manage faults	Artifacts with complex geometry Geostat issues when cells shape varies Time and disk space consuming
LGS	Allows modeling deformations and heterogeneity variations Allows working directly in the structural grid Honors complex deformations in cartesian grids	Cannot manage faults Limited to variogram-based algorithms May require flattening of data for variogram modeling
MPS	No need for flattening Can manage facies and continuous variables Additional deformations can be generated during simulations	Requires an appropriate TI Issues with change of support Time consuming

Table 1.

Pros and cons of the different methods.

Flattening algorithms allow filtering of large-scale deformations and faults, which makes this kind of approach well suited to the modeling of large fields. On the other hand, this geometrical approach is not suited for managing change of support, and it is better to work at a point scale. Because of these two characteristics, the flattening approach is relevant for geological modeling in the oil and gas industry and can be used with modeling grids of any type.

In the mining industry, where the mineralizations are generally smaller than in oil or gas fields and where change of support is important and modeled, the use of regular grids is usually preferred. The purely geostatistical approaches are then well suited for modeling mining variables alone or in combination with a flattening operation.

It must be noted that the modeling of multi-phased deformations will require the sequential use of the two types of approaches. Any kind of grid type can be used if the modeling is done at point scale, whereas regular grids will be preferred if it is necessary to account for change-of-support effects.

5.3 Future directions

The geostatistical methods presented in this chapter can model spatial variations of more parameters than the orientation of the anisotropy ellipse or local rotation vectors. In particular, the degree of heterogeneity can also vary in space through local ranges (LGS approach or local scaling factors in MPS approach). Therefore, a wide range of geological configurations can be modeled with purely geostatistical methods, even in very complex cases. Complex deformations and spatial continuity variations can be considered simultaneously, which allows modeling complex environments, at different scales.

Nevertheless, it must be kept in mind that the spatial variations of the different variogram parameters must be calibrated to real data to get realistic models of the subsurface. This inference of the parameters is often difficult, due to lack of data, and such sophisticated models often represent the geologist’s ideas. Therefore, this kind of model will be mainly based on geological interpretations, and it is of primary importance to make an accurate and high-quality geological analysis before calculating the models. Using geological analogs (an old mine in a similar geological context, for example) can be very valuable if and only if the geological contexts of the zone to be modeled and the analog are really similar. Again, the quality of the geological analysis is critical. This issue must be kept in mind when using modeling results.

The use of secondary data, like geophysical data, can be considered as long as they are at the appropriate scale for helping inference of the modeling parameters. Some research work in this domain and appropriate software development should be useful to facilitate the use of these modeling approaches.

It is interesting to notice that these features are compatible with any type of grid when working at point scale and can easily account for change-of-support effects with usual support models like DGM when working with regular grids.

Acknowledgments

The author acknowledges Geovariances for having allowed and funded this publication, and for having provided access to the Isatis.neo software in which the presented methods are implemented, and which has been used for generating most of the Figures.

Special thanks are to Dr. Pedram Masoudi and François Geffroy, from Geovariances, and Dr. Emilie Chautru, from Mines Paris-PSL, for their useful comments and help.

Conflict of interest

The authors declare no conflict of interest.