A fast FFT method for 3D pore-scale rock-typing of heterogeneous rock samples via Minkowski functionals and hydraulic attributes

The integration of numerical simulation and physical measurements, e.g. digital and conventional core analysis, requires the consideration of significant sample sizes when heterogeneous core samples are considered. In such case a hierarchical upscaling of properties may be achieved through a workflow of partitioning the sample into homogeneous regions followed by characterization of these homogeneous regions and upscaling of properties. Examples of such heterogeneities are e.g. fine laminations in core samples or different micro-porosity types as consequence of source rock components and diagenesis. In this work we utilize regional measures based on the Minkowski functionals as well as local saturation information derived through a morphological capillary drainage transform as a basis for such a classification/partitioning. An important consideration is the size of the measurement elements utilized, which could be considerable in the case of larger heterogeneities; in such case the calculation of the regional measures can be computationally very expensive. Here we introduce an FFT approach to calculate these measures locally, utilizing their additivity. The algorithms are compared against direct summation techniques and shift-overlap approaches for a selection of different averaging supports to illustrate their speed and practical applicability. We consider a range of artificial Boolean models to illustrate the effect of including hydraulic information on the resulting classifications scheme. This allows the determination of bias, since for these model systems local classes are known ab-initio. The classification framework is tested by comparing to the known initial micro-structure distribution and relative bias quantified in terms of choice of averaging elements (size and shape). Importantly, depending on the actual morphological transition between micro-type partitions, partitions including hydraulic attributes differ from pure morphological partitions with applications to electrofacies and hydraulic unit definitions.


Introduction
Rock type classification is an important topic in reservoir characterization and spans a large range of scales from seismic interpretation (meter scale [1]) to well logging (down to centimeter scale [2]) and further down in scale to core analysis. It forms part of an upscaling strategy involving the correlation of classes of different type, purpose, or origin; e.g. lithofacies vs. electrofacies (distinguishable classes based on well-logs). However, reservoir heterogeneity extends to even smaller scales.
The rapid development of (micro-) X-ray CT techniques exposes rock heterogeneity down to the order of microns, at the pore-scale, in 3D. While the calculation of petrophysical properties and correlations has been demonstrated for a range of properties including e.g. electrical conductivity, elastic moduli, mercury intrusion capillary pressure, and permeability [3]- [7], this is for homogeneous sandstones of limited clay content. For more complex carbonates cross-variograms between porosity and permeability and multiple imaging length scales [8] may be utilized. Often one may however want to utilize a hierarchical upscaling approach based on corescale rock-types to integrate several length scales and multiple physical properties, including for thinly laminated sandstones. While petrophysical rock-types may include such laminations and thus anisotropy caused, describing fluid flow leads to a homogenization problem which can be approached at the Darcy scale [9], or numerically through partial homogenization [10] by solving the Brinkmann equation. Either way, a classification technique for scale separation is required.
Classification techniques targeting petrophysical interpretation are not new. Archie developed a method of skeleton classification with a few petrophysical parameters in [11]. Skalinski and his colleagues used pore throat distributions as the main parameter for pore type classification with mercury injection pressure measurement data of many rock samples [12]. Characteristics of fractures and vugs of carbonate rocks were extracted from micro-CT images in [13]. Dernaika et al. [14] calculated porosity and mineralogy of core samples from advanced dual energy XCT imaging and linked to plug data after upscaling.
Schmitt and colleagues considered pore shapes (length, width, thickness, orientation) in reservoir rock from 3-D X-ray micro-CT images [15], limiting themselves to large pores due to the higher resolution requirements for deriving these measures. Alternately, one may characterize microstructure utilizing the additive morphological measures volume, surface area, mean and total curvature, namely the Minkowski functionals [16], which can be evaluated easily on a segmented tomogram [17], [18] and have been used in the past as sensitive descriptors of morphology [17]- [22]. The strong correlation of physical properties of Boolean models to morphological properties, the Minkowski functionals [23], [24], motivates the use in rock-typing approaches. This was done in [25] for a thinly laminated sandstone using a 1D sliding measuring window approach and in 3D using non-overlapping regions of the image. A highresolution field of morphological regional measure could not be computed due to computational cost. In [26] a support shift technique was utilized for the computations, which however still does not scale to large supports, as we require in this work. We define as support of a measurement a small (regional) volume within the sample, and over which an average is taken. A larger support is desirable to reduce the standard deviation of such a measure.
Another relationship often used in rock-typing approaches is given by mercury intrusion capillary pressure (MICP) curves [27]. In terms of micro-CT approaches this type of information can be gained either by direct imaging approaches [28], or by modeling the MICP process numerically [29], [30]. Capillary pressure heterogeneity was also inferred from regional saturation distributions in [31], while [32] considered non-wetting phase cluster distributions.
The resulting question we target in this work is whether the inclusion of hydraulic information -taken to be local fluid saturation -changes rock-types defined on the basis of the Minkowski functionals alone as basis for petrophysical rock-types due to their strong relationship to physical properties. We utilize the decomposition of the Minkowski functionals into histograms over vertex configurations to implement (fast) calculations of the Minkowski measures over large support via Fast Fourier Transforms (FFTs). The technique is compared against more direct implementation schemes on a set of artificial microstructures and applied to the classification of laminated Boolean composites. Classifications based on the Minkowski measures alone are compared with ones including a saturation attribute and resulting bias is discussed.

Methodology
In this section we review the basic techniques utilized in this work. We introduce the Minkowski functionals and measures as a full set of additive measures of 3D space, regional measures based on these definitions over a finite support, as well as resulting curvature density fields and algorithms to calculate them. Finally, for demonstration of the algorithms we introduce a set of Boolean composites featuring forming layered media.

Minkowski functionals
In this work we follow the notation of [26]. Consider a body Y in ℝ 3 with sufficiently smooth surface ∂ . Then the additive geometrical characteristics are volume V(Y), surface area S(Y), integral of mean curvature M(Y) and integral of total curvature K(Y) with surface integrals Here r1(s) and r2(s) are the maximum and minimum curvature radii in s. The above characteristics are related to the intrinsic volumes via On the segmented tomogram the intrinsic volumes can be determined by evaluating the histogram of the individual vertex configurations. Excellent illustrations of this process have been given in [17], [33]. Here we consider the grain phase as the phase of interest; all other phases are merged for this calculation and the characteristic function of the structure takes the value of 1 for the granular phase and 0 otherwise. A vertex made up of eight voxels can then be in 2 8 = 256 different configurations, where each configuration has a certain contribution to volume, surface area, mean, and total curvature. Lookup tables for these contributions has been given in [17], [33]. We use the 8-neighbourhood for 2D connectivity and 26-neighbourhood for 3D connectivity as given in [17]. For a recent review of the applications of integral geometry in porous media see [19].

Regional curvature measures
To define a measure over a finite volume we consider the average over regions of spherical and ellipsoidal shape. For a given position of this support in the actual tomographic images we then evaluate all vertices of the discrete tomogram, where the centers of the vertices are within the measurement window.
Consequently, the regional measures can then be derived from the local configuration histogram over the given spatial support. Evaluating the regional measures for a regular grid of support vectors will generate a set of morphological fields V(r), S(r), M(r), K(r). We report all measures as curvature densities, e.g. normalized by the averaging volume.

Curvature density fields by direct summation
The most direct approach for calculating the curvature density fields simply places the spatial support at a position-of-interest, evaluates the configuration histogram, and stores the result at the respective position, looping over the curvature measure grid. For dense strides of the support this leads to significant overlap of support volumes for adjacent evaluations; thus, large numbers of vertices are evaluated multiple times for large measurement support, e.g. large radii of a spherical support. This naïve approach therefore does not scale well with increasing support size, as it scales with support volume.

Curvature density fields by support-shift algorithms
A somewhat more efficient way recognizes this partial support overlap and only evaluates vertices anew which are added to the support volume after a shift on the lattice, e.g. stride along a lattice direction with a fixed step size; the same number of vertices will be removed and added from the histogram for each support shift when shifting in the same direction. This approach scales significantly better, and essentially scales with the surface "area" of the support for small strides; it takes direct advantage of the additivity of the Minkowski measures.

Curvature density fields by FFT methods
The additivity of the curvature measures allows a much more efficient approach to deriving regional measures. This approach utilizes a convolution of the support with the vertex configuration field, thus deriving the curvature density fields by a series of Fast Fourier Transforms (FFTs) and adding up the curvature contributions to each regional measure according to the lookup table of [17]. In this case the size of the support has essentially no influence on the speed of the algorithm. Very large support volumes and arbitrary support shapes can be utilized this way. [a] [b] Fig. 2. 2D slice through a layered Boolean model with two generating sphere processes, one for each layer.
[a] generated microstructure, [b] microstructure partially saturated by a nonwetting phase with invasion radius ri = 6.6. Sphere radii are r1 = 13 and r2 = 26. The microstructure has dimensions of 1200 x 1000 x 800 voxels after removal of boundary regions in arbitrary units. Blue denotes the invading non-wetting phase.
A comparison of the relative performance of the algorithms is given in Fig. 1 for a sub-volume of a segmented tomogram of Bentheimer sandstone. Clearly, using large averaging volumes becomes very costly with the non-FFT based algorithms applied here, while the FFT approach scales well, e.g. is independent of averaging volume. A significant speedup of the FFT method can be achieved by utilizing the rotational symmetry of the possible vertex configurations [17]. For the FFT method we consider both spherical and ellipsoidal measure supports. In practice the size and shape of the measurement support may be chosen through structural analysis, utilizing variograms [26], and considering that the resulting rock-types should be compact, e.g. do not have too many holes. Due to additivity and the availability of parallel FFT libraries for various architectures the proposed technique easily scales to largest datasets.

Boolean model
We generate layered synthetic micro-structures by controlling the generating grains of a Boolean (or Poisson particle) process according to an initial (known) layered distribution of two rock-types: after choosing a Table. 1 Structural details of dry sample in Figure 2. The generating shape for both layers is spheres. random location to place a particle the particle shape and size is decided by the corresponding pre-defined rocktype: Similar microstructure modeling approaches were presented in [10], [25], [30] for a more general 3D case, which was not required here for the demonstration of the technique. Placing a spherical particle of radius r at a random location results in a decrease of porosity. Particles can overlap into adjoining rock-types and the particle placement process is stopped when both rock-types achieve their target porosity. A 2D slice through a 3D realization of a Boolean layered system with two different sphere sizes of r1 = 13 voxels and r2 = 26 voxels is given in Fig. 2a. The structure parameters are summarized in Table 1.

Simulated non-wetting fluid invasion
We model the invasion of the synthetic rock with a nonwetting fluid by utilizing the capillary drainage transform [7], [24] with a fixed saturation chosen to provide additional contrast between the different rock-types. This is achieved by selecting an invasion radius below the critical radius of percolation for the larger-scale porosity, but above the critical radius of percolation for the smallerscale porosity. A 2D slice through the saturated medium at Sw = 0.38, corresponding to an invasion radius of ri = 6.6 is given in Fig. 2b. The critical radii for the two sphere composites making up the layers are rc1 = 4 and rc2 = 9. Note that the original samples were 1600 x 1200 x 1600 voxels. We selected the inner 1200 x 1000 x 800 region for further analysis in order to avoid boundary effects resulting from simulated capillary pressure fluid distributions.

Rock-type classification
We carry out the classification into different rock-types on the regional curvature measures and saturation field. This is achieved by defining the centroids of a Gaussian Mixture Model (GMM) in a supervised way by selecting the central zones of the layers as can be seen visually. The GMM centroids have cross-variance matrices of size 4 x 4 each (curvature measures only) or 5 x 5 (including saturation), which together with the respective means are derived upfront from statistics of unambiguous regions.

Results
The application of an ellipsoidal measurement support with half-axes of a=120, b=30, c=120 and resulting curvature and saturation fields is given in Fig. 3. The averaging support size and shape was chosen to avoid generating holes in the resulting classification and respecting the evident layering. The corresponding histograms of the regional measures are given in Fig. 4.
Recall that the porosity of the two rock-types is essentially the same. Consequently, there is no real contrast in the solid volume measure (Fig. 4a). The different color scales (Fig. 3b) are a consequence of different grain sizes, the larger structure showing the more extreme values. We observe a clear separation of rock-types in the surface area density as one would expect (Fig. 3c, Fig. 4b), with small spheres corresponding to the larger surface area contributions. The 2D slice illustrates the strong layering detected by surface area, while the histogram illustrates a good separation with a moderate area one could interpret as overlap region. An equally sharp contrast in layers is visible in Fig. 3d, corresponding to mean curvature. The smaller values (larger in magnitude) correspond to the small spheres. The corresponding histogram (Fig. 4c) shows an excellent separation with only a small area between peaks. Consider now the total curvature measure (Fig. 3e, Fig. 4d). The smaller sphere regions show a stronger connected behavior expressed by the more negative total curvature measure. Compared to this, the larger sphere regions exhibit a negative peak suggesting connectivity of both phases, but with some values becoming positive. Even larger measurement supports may lead to an even better separation of layers for this measure. Finally, the saturation field offers an approach to generate contrast between different rock-types by capillary pressure, here in a static sense; we consider this a hydraulic attribute for rock-typing. The measure is potentially attractive for low resolution image classification as such contrast can be generated either by simulation or by actual differential tomographic imaging (dry & wet imaging). From the 2D cross-section of the saturation map we observe that for the given invasion radius the rock-type with smaller features remains largely fully saturated with wetting fluid, with some exceptions close to the boundary, making the interface between rocktypes less consistent (Fig. 3f). The corresponding histogram shows a huge spike for Sw = 0 (truncated), and a somewhat distributed peak with a significant transition zone otherwise. The distribution of the curvature and saturation measures is further summarized in Fig. 5 through a 1D profile. We emphasize that the procedures introduced here are general for 3D; however, given the special nature of the laminated samples it is logical to also provide a corresponding 1D log.
Consider now the resulting rock-type classifications of the laminated microstructure for the case of the ellipsoidal measurement support of a=c=120 voxel and b=30 voxel via the GMM procedure. We illustrate the sharpness of the boundary both by taking 2D slices perpendicular to the bedding direction (Fig. 6) and by illustrating the resulting classification as 1D log (Fig. 7). For the case of using the Minkowski measures only as basis for the classification the recovered interface between the rock-types is sharp and coincides well with the original interface. Adding water saturation to the classification basis results in a somewhat smoother [a] [b] [c] [d] [e]  transition of the rock-types; furthermore, an apparent bias in the position of the boundary is introduced. We attribute this to an internal invasion boundary condition, e.g. a pore-scale saturation transition zone. This interpretation is supported by Fig. 6b and Fig. 6d, showing local misclassifications likely due to local heterogeneities in structure, allowing invasion of the boundary at lower capillary pressure. This local variation in saturation is also evident in Fig. 7. [a] [b] Fig. 7. 2D slices through the recovered 3D layer type distribution along Y dimension, overlapping original sample. The curve in (a) is recovered from the four Minkowski descriptors while (b) is recovered from Minkowski descriptors and fluid saturation. The blue color depicts the non-wetting phase saturation calculated by a MICP simulation.
We further quantify the relative boundary bias by scanning through each x-z column along the y-axis, calculating both the mean and standard deviation of the interface location for the case two different averaging element sizes and the Minkowski measures only versus including water saturation as additional hydraulic discriminator. The respective measures are reported in Table 2. The curvature measures alone recover the thickness of all layers almost exactly for both averaging ellipsoidal supports (b = 30 and b = 40, a = c = 4b). There is little variation in the location of the interface as well, e.g. no significant undulations are introduced, which is visible in the small standard deviation. Adding water saturation to the classification basis increases standard deviation significantly and introduces considerable bias towards the larger-scale rock-type.

Conclusions
In this work we introduced a fast FFT based algorithm for the calculation of regional curvature fields in 3D for the case of large measurement supports. The algorithm was demonstrated on a layered Boolean microstructure mimicking a laminated sandstone. It provides linear scaling with sample size and does not depend on the spatial support size. This allows the use of large spatial supports for rock-typing purposes utilizing the Minkowski functionals. We further illustrated the classification bias introduced by including hydraulic information into the classification basis. For the given model rock, the strongest discriminator in terms of interface sharpness was regional mean curvature, followed by regional surface area. Including hydraulic information into the classification basis resulted in a smoother boundary transition which overestimated the relative fraction of the larger scale porosity.
CHA thanks the Australian Research Council (ARC) for funding through grant DP160104995 and through an ARC Future Fellowship (FT120100216). HJ and CHA thank the Chinese Scholarship Council (CSC) for their support. This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. We thank reviewers Ryan Armstrong and Steffen Berg for a careful review, which further improved the quality of this manuscript.