Classification and quantification of pore shapes in sandstone reservoir rocks with 3-D X-ray micro-computed tomography

Introduction Conclusions References Tables Figures


Introduction
Natural and artificial materials are often characterized by their pore/grain shape or size distributions as determined by distinct analytical instruments.Several investigations have been conducted to classify pore types/shapes, and most of them are associated with 2-D quantification.However, a 3-D pore with irregular shapes cannot be appropriately characterized from 2-D image sections in pore-typing procedures, nor can the number of pores (Buller et al., 1990).The mass transport through porous media strongly depends on the structure of the pore network (Wiedenmann et al., 2013).Knowledge of pore morphologies is essential and they are one of the main parameters that control fluid flow underground.Petrophysical properties such as permeability, electrical conductivity and drainage capillary pressure are strongly influenced by throat sizes, which are constrictions of minimal cross-sectional area between pores (Buller et al., 1990).Many properties are controlled by pore/grain textures.Soete et al. (2015) demonstrated how seismic reflections in travertine systems are related to geobody boundaries, in which the seismic expression is a function of porosity, pore types, and shapes.The size and shape characterization of irregular particles is a key issue in many fields of science, which is often associated with large uncertainties (Bagheri et al., 2015).Moreover, to the best of our knowledge, the systematic 3-D pore shape quantification of sedimentary rocks based on sample size, pore network detachment (i.e., separation into individual pores), and distinct geometrical descriptor measure-Published by Copernicus Publications on behalf of the European Geosciences Union.
Mayka Schmitt et al.: Classification and quantification of pore shapes in sandstone reservoir rocks ments has not been comprehensively reported in the literature.
X-ray micro-computed tomography (x-ray µ-CT) is a robust technique that allows the three-dimensional (3-D) investigation of many materials and enables the quantification of relevant microstructural features.Although the application of µ-CT in the geosciences was previously introduced by Ketcham and Carlson (2001) and Ketcham (2005), many investigations that analyze irregular particle shapes are still performed with 2-D approaches (Buller et al., 1990;Zhang et al., 2016;Petrak et al., 2015).So far, µ-CT studies of 3-D features have mostly been related to the characterization of volcanic rocks (Eiríksson et al., 1994;Riley et al., 2003;Shea et al., 2010;Vonlanthen et al., 2015) or spheroidal objects (Robin and Charles, 2015), and these features are usually described by means of equivalent size or shape parameters, such as roundness or aspect ratio (Little et al., 2015;Saraji and Piri, 2015).Nevertheless, these intuitive descriptors provide more qualitative denotation and might obscure information regarding the original particle shape.Recently, McGrath et al. (2015) reported the difficulty in characterizing the 2-D and 3-D shapes of "free gold" particles in a flash flotation by means of circularity and sphericity.These authors concluded that these measurements cannot adequately describe the diverse shapes of irregularly shaped particles.Available high-resolution techniques might provide qualitative and quantitative pore structure information.However, many materials, including sedimentary rocks, lack information regarding the deep comprehension of irregular particle (pores or grains) shape analysis and a systematic way of classifying it, e.g., describing it to be similar to artificial object forms.A detailed pore shape characterization is necessary among many applications and can be used to infer the dominant mechanisms that act in a heterogeneous rock in response to its macro properties.Additionally, numerical methods and pore shape description are essential to optimize and to study parameter variations that influence laboratory measurements systematically.
In this study, the pore shapes of three sandstone rocks from distinct fields were analyzed and classified based on the xray µ-CT approach.This method involves counting and measuring detached pore objects directly from three-dimensional images.The 3-D images of the analyzed rocks were acquired at similar resolutions, and two main pore constituents were identified and defined: (i) main connected pore networks and (ii) disconnected pore ganglia.A watershed algorithm was applied to separate the main pore networks into individual pores to preserve the essential pore morphology, and the most suitable marker extent parameter (as it will be discussed in Sect.3.2.) was defined.Three subsamples with volumes of 1000 3 , 500 3 , and 250 3 voxels were extracted from each rock to investigate the effect of sample size on the pore shape classification.The advantage of using the proposed approach is that shape parameters can be calculated directly from the preserved pore textures with no need for equivalent volume (sphere or resistor) conversion, which creates inaccurate results for the shape classification.Additionally, the image visualization and analysis are completely automated and performed in the Avizo Fire software for an unlimited number of particles, which are geometrically described to facilitate the quantification analysis.

Geometrical aspects of 3-D particles
To describe and quantify a particle form (in this case, pores) in three dimensions, morphological parameters such as length, width, and thickness are required.These 3-D parameters must be perpendicular to each other but do not need to intersect at a common point (Blott and Pye, 2008).To perform the pore shape classification approach that is proposed and described in this work, we conventionally assigned the following practice for the geometrical descriptor of individual particles: (L) is the longest pore dimension, (l) is the longest pore dimension that is perpendicular to L, and (S) is the smaller pore dimension and perpendicular to both L and l.In practice, two methods were applied to measure L, l, and S from the 3-D irregular shaped objects, such as the pores that were found in the analyzed sedimentary rocks: (i) the bounding-box (BB) and (ii) the Feret caliper (FC) geometries.
In bounding-box analysis, the object coordination system follows the same as what is acquired for the 3-D sample volume, namely, L, l, and S; the large, medium, and small sizes from the boxing-based volume have axes of X, Y , and Z.In Feret caliper analysis, a maximum length (L) and a minimum width (S) caliper diameter that belong to the same plane are initially determined from all the object orientations, while l is the maximal diameter 90 • from the width.Hence, the Feret dimension of a 3-D particle is defined as the normal distance between two parallel tangent planes that touch the particle's surface.Because this value depends on the particle's orientation, a measurement for one single particle has little significance; therefore, measurements in 31 directions were taken from the image analysis software to determine the particle length (maximum length of Feret distribution) and width (minimum width of Feret distribution). Figure 1 shows the differences in the geometrical parameters L, l and M for bounding-box (a) and Feret caliper (b) methods, which are drawn in a three-dimensional pore particle.The shown particle was detached from the main pore network of rock S 1 (250 3 voxels) after applying the "Bin3" marker extent in the watershed algorithm and choosing the third highest 3-D length from the 25 highest 3-D volume particles.Table 1 depicts the values (in micrometers) of the grid-cell axes and 3-D geometrical parameters of the analyzed pore particle.The Feret caliper led to a value of L that is approximately 10 % higher than what is measured from the bounding-box calculation.
Another common dimensional feature that is used to characterize the thickness of 3-D particles is shown in Table 1.This feature is the equivalent diameter (EqD), which gives the analyzed object a corresponding spherical diameter size with equal voxel volume.Although this parameter is frequently used in the literature for 3-D pore structural description (Cnudde et al., 2011;Van Dalen and Koster, 2012), the error in the particles' thickness determination can be considerably high depending on the object's shape, especially for very irregular ones.In this study, the pore shapes were classified by considering the equivalent diameter as the thickness (l) descriptor in the Feret caliper analysis; the mismatches that were found in the results are discussed.

Theoretical background for pore shape classification and quantification
The relationship between the thickness and the length of a particle (l/L) has long been used to indicate its elongation (Lüttig, 1956), while the true flatness is best described by its width divided by thickness (S/ l).Regardless, the degree of sphericity or roundness, "equancy", has been used to describe equidimensional particles; Krumbein (1941) described equant forms as being spherical, with L = l = S.According to Blott and Pye (2008), a particle form can be qualitatively described in terms of its deviation from equancy; specifically, the equancy degree is defined by the combination of flatness and elongation.Figure 2a depicts five classes of equancy that are defined at 0.2 intervals for a qualitative classification of particle equancy plotted on a Zingg diagram (Blott and Pye, 2008).The methodology that is used in the present work to quantify sandstones' pore shapes is based on the espoused pore shape classification in Fig. 2b.This classification follows the original approach of Zingg's (1935) diagram, which was deeply discussed by Blott and Pye (2008) and applied recently by Soete et al. (2015) while studying travertines.As one can see in Fig. 2b, cubes are particles with comparable L, l, and S; rods are objects that are characterized by a much larger L and by the smallest S, which is similar to the thickness (l); and plates have equivalent L and S and, conversely, much smaller l.Blott and Pye (2008) discussed a classification by Sneed and Folk (1958), but the former authors consider that the original terminology that was proposed by Zingg (1935) is more appropriate for use in most natural sedimentary environments, providing a more even distribution of the form continuum.The so-called Zingg diagram was applied in the field of mineralogy to classify rocks by their shape.In mineralogy, two measurements differ if the ratio of the smaller one over the bigger one is larger than 0.6.However, the diagram thresholds may vary based on the application of distinct research fields because the technical cleanliness typically considers particles fibrous if their aspect ratio is 0.1 (Vecchio et al., 2011).Similar to Soete et al. (2015), who analyzed travertine rocks, we used the diagram in Fig. 2b to classify and validate pore shapes in sandstone samples by considering artificial objects as having shapes.Thus, when applying the proposed pore shape classification, individual pores are denominated according to the respective class numbers: class 1 (rod-like), class 2 (bladelike), class 3 (cuboid-like), class 4 (plate-like), and class 5 (cube-like).By using an identical 3-D subsample of volumes and images that were acquired roughly in the same resolution, the geometrical pore descriptors and shapes of three sandstone fields were systematically studied and quantified as follows.

Methodology
Pore shapes in sandstone reservoirs that were acquired from 3-D x-ray µ-CT images were classified and quantified based on artificial object similarity, as shown in Fig. 2b.After image filtering and pore-phase segmentation, the 3-D µ-CT data sets were analyzed by means of sample size, main pore network detachment, and particle descriptor measurements.A brief geological description of the analyzed sandstones and  Blott and Pye, 2008;Soete et al., 2015).the x-ray µ-CT parameters that were used in the analyses is provided as follows.
The Bentheimer sandstone is part of the lower Valanginium sequence, which is associated with the Lower Cretaceous (approximately 140-134 million years old) and it is comparable to the well-known Berea and Fontainebleau sandstone.This sandstone was segmented into three lithostratigraphical units according to Kemper (1968): lower Bentheimer sandstone, intermediate Romberger, and upper Bentheimer sandstone (also known as flasered sandstone).This rock has been petrographically classified as quartz sandstone with low clay and silt content (< 8 vol%).The grain size varies between very fine and coarse, although fine grains are predominant.In general, the components are well rounded and well sorted, excluding partial content of stable and unstable heavy minerals and feldspars.Stadtler (1998) provided a comprehensive and more detailed description of the S 1 -type sandstone.
The Obernkirchen sandstone is part of the upper Berriasium (also known as the Wealden sequence, approximately 145-150 million years old), which is associated with the Lower Cretaceous and forms the lower transition to the Jurassic.In principle, this rock is comparable to the S 1 rock type, despite containing better sorted grains and a higher degree of clay cementation.According to Hafner (1987), this sandstone is characterized by a homogeneous and fine-grained matrix structure, which is sometimes interrupted by thin clayey layers.Additionally, feldspar minerals indicate strong mechanical influences, which results in a higher degree of degradation than for the S 1 sandstone.The pore space is filled with different types of cementation, mostly carbonate and clay minerals (Börner et al., 2011).
The Flechtingen sandstone is part of the so-called Rotliegend sequence (approximately 302-257 million years old), which is associated with the Lower Permian.This rock greatly differs from S 1 and S 2 in terms of both petrophysical properties and the complexity of the pore network structures because of microfractures.According to Bärle (1996), this sandstone is generally characterized by a distinct heterogeneity throughout the entire sequence, which was caused by different depository environments.A fine layer of small-and medium-sized grains on a small (millimeter) scale is predominant.The red color of the rock is a result of fine deposited haematite, which was caused by the arid and desert-like environment during the Permian.Cementation minerals are often of calcitic and dolomitic nature, but baryte and clayey cementations are also very common.Clay minerals are often represented by chlorite and illite but are mostly represented by kaolinite.The resulting pore structures strongly vary on a small scale because of complex diagenetic phenomena, the mineralogy, and grading effects (Bärle, 1996;Paech et al., 2006).In summary, these different sandstones feature a broad range of pore structures because of their distinct diagenesis and tectonic history.Hence, we should be able to detect differences for the systematic and individual analysis of pore shapes in these three reservoir rocks by using the proposed methodological/statistical approach.

3-D x-ray µ-CT data acquisition and processing
The x-ray µ-CT analyses were conducted at the Leibniz Institute for Applied Geophysics (Germany) by using a nanotom 180 S instrument (tube characteristics: 180 kV, 500 µA), which was manufactured by GE Sensing & Inspection Technologies (product line of Phoenix X-ray, Wunstorf, Germany).The system is computer-controlled, and cylindrical samples with a diameter of 2-4 mm were fixed in the 360-degree rotator holder.The parameters that were used for the µ-CT data acquisition and image post-processing are shown in Table 2.The samples have roughly the same resolution, which was intentionally imposed to allow a comparative study among the analyzed sandstones.Image reconstruction was performed on a GPU(graphics processing unit)-based cluster (fourfold NVidia Tesla GPU), and no additional filters were applied during this stage.
The 3-D reconstructed volume was post-filtered, followed by a pore segmentation and separation process that was performed with the Avizo Fire 8.1.0software (Avizo, 2014), which allowed for data visualization and quantification.As shown in Table 2, specific filters were evaluated for each analyzed rock to improve pore noise reduction or the gain of the pore matrix/clay border contrast.The non-local means filter (Buades et al., 2005) showed the best results, although the additional median 3-D filter (Ohser and Schladitz, 2009) was applied to samples S 2 and S 3 to enhance the pore segmentation removing image noise without blurring pore edges.During the pore segmentation stage, each voxel was assigned to the object of interest (pore) based on a mixture of 2-D scanning electron microscopy (SEM) image registration and an automatized threshold shift of grayscale intensities by using the Otsu algorithm (Otsu, 1979).For each of the three rocks, 3-D volumes with dimensions of 1000 3 , 500 3 , and 250 3 voxels were cropped from the original 3-D data set.Segmented pore networks were extracted from these subsamples and prepared for further quantification.During this process, distinctively small and very small disconnected pore networks exist even before any pore detachment workflow is applied, in addition to the main pore networks.Hence, we define these residual pore networks as pore ganglia.

Identifying pore networks and residual pore ganglia from the 3-D pore structures
Limited voxel resolution in x-ray µ-CT analyses might cause imprecision when determining small-scale heterogeneities, such as pore throat sizes, to characterize rock samples (Beckingham et al., 2013).Therefore, the number of pore ganglia in a 3-D pore structure obeys the relationship between the detected pore volume (V p ) and its voxel resolution (R v ): the smaller this value, the higher the undefined gray level value in the image and more difficult it is to identify a voxel cluster as a pore with clear morphology.Consequently, a higher number of pore ganglia will be detected.To identify and characterize pore ganglia from a 3-D segmented pore volume, finding a lower cutoff limit for the ratio V p /(R v ) 3 , i.e., a pore volume size in voxels (N v ), which allows us to distinguish pore ganglia to be suppressed from the real pore system, is necessary.This examination is usually based on the distribution of pore sizes as obtained from topological analyses of the pore space (Dong and Blunt, 2009).According to Andrew et al. (2013) applied to determine pore ganglia: where V p and R v are in µm 3 and µm.This choice is based on the thickness of the individual pores from the 3-D data volumes before any pore detachment process is applied.As shown in Fig. 6a, the distinct descriptors (bounding-box, Feret caliper, and equivalent diameter) that were used to calculate the pore thickness of the analyzed rocks indicated that most of the smaller pores had sizes below 30 µm.Additionally, if one considers a defined pore cluster of voxels that was replaced by its digital analogues as being a sphere (N Veq ), then its equivalent size (D eq ) is given as So, if N Veq equals 2000 voxels, then D eq is approximately 16.This value, when multiplied by voxel resolutions of 1.75 and 2.25 µm, is equal to 27.3 µm (for samples S 1 and S 2 ) and 35.2 µm (for S 3 ), respectively, which matches perfectly within the observed pore ganglia as shown in Fig. 6a. Figure 3 depicts the 1000 3 voxel rendered volumes of the analyzed rocks, including the grayscale image (a), five largest (by volume) segmented pores (b), and detected pore ganglia (c), before any pore detachment process had been applied.In Fig. 3, the number of pore ganglia differed between the analyzed rocks, equalling 0.13 % for the Bentheimer sandstone (S 1 ), 0.31 % for the Obernkirchen sandstone (S 2 ), and 0.34 % for the Flechtingen sandstone (S 3 ).The relatively higher percentages for the S 2 and S 3 samples are explained by the rocks' porous nature.The two samples consist of both a defined connected pore network and smaller pores, some of which are not adequately distinguished between the gray levels of the analyzed resolution.For S 3 , this effect is notably more pronounced; therefore, instead of one main pore network (blue object for S 1 and S 2 ), S 3 has several networks that are located in the preferential direction perpendicular to the Z axis.Nevertheless, all the pore ganglia were suppressed from the 3-D renderings of data for an appropriate pore shape analysis and characterization before detaching the main pore networks to obtain individual pores.

Detachment of main pore networks by preserving pore morphology
The pore shape classification involves counting and measuring the 3-D geometrical descriptors (L, l, and S) of individual pores directly from the original segmented images.When quantifying a 3-D data set, authors frequently associate the segmented pore throats/bodies with an equivalent volume, usually cylinders/spheres, to perform a quantitative characterization of the porous media.However, the original pore morphology can be achieved by preserving the segmented pore phase and applying, e.g., a watershed algorithm (Ohser and Schladitz, 2009).The principle of the watershed algorithm is to compute watershed splitting lines on a segmented 3-D image that will detect surfaces and separate agglomerated particles, which are subtracted from the initial image.
When running a watershed algorithm, the most meaningful parameter is the "depth of valley" or "marker extent" as named in the Avizo Fire software, called "Bin" herein.In short, the smaller the chosen Bin is (lower levels of flooding), the more separated the main pore networks will be, and smaller individual pores will be generated and identified.Thus, without the Bin command, the 3-D pore space rendering comprises one or several connected main pore networks plus a high number of pore ganglia.Figure 4 shows the 3-D renderings and 2-D projections for 250 3 voxel data, with individual pores represented by distinct colors.To in-vestigate the algorithm effect on the residual connected pore networks, the smallest pore ganglia have been intentionally removed; also, the weight graduation 10, 3, and 1 in the Bin command were arbitrarily chosen.
The changes within the marker extent parameters in the watershed algorithm results are visualized in Fig. 4a for the 250 3 voxel volume of the Bentheimer rock (S 1 ).The evaluation of the Euler number (Vogel et al., 2010;Renard and Allard, 2012), which is used as a measure of the degree of fragmentation of pore networks, and the number of individual pores that is generated after each detachment step is shown in Fig. 5.In Fig. 4a, the main pore network (blue structure on the left side) is broken into a certain number of smaller pore networks that decrease in size with the Bin parameter number.Thus, the best marker extent (Bin) must be determined.Because the watershed algorithm depends on pore textures, an accurate investigation of distinct Bin parameters and 2-D/3-D visualizations was performed for each rock individually to reliably determine the detachment of complex pore structures, such as those found in the sedimentary rocks.Otherwise, the pore structures might be underseparated (Bin10) or over-separated (Bin1).Fortunately, the same behavior within the Bin parameter was observed for the Obernkirchen (S 2 ) and the Flechtingen (FL) sandstones, with Bin3 being the most faithful marker extent.Figure 4b  discussed based on distinct Bins.For the analyzed rocks, the "NoBin" 3-D renderings comprise few main pore networks and mostly pore ganglia.As soon as the detachment process is performed (Bin 10, 3, and 1), the pore ganglia that are suppressed from the 3-D volumes allow the emphasis of the main pore system morphology in the shape analysis.

Evaluating 3-D geometrical descriptors based on the pore thickness distributions
One of the most important steps to be established while performing pore shape analysis is the measurement of the 3-D geometrical descriptors L, l, and S from the segmented 3-D particle.In this procedure, a truthful quantification of individual pores in the original 3-D pore system morphologies should be provided.As mentioned previously, the boundingbox, Feret caliper, and equivalent diameter methods were used in this work.To calculate the pore ganglia and main pore network descriptors, this procedure was performed before applying the watershed algorithm; moreover, the descriptors were measured after each Bin 10, 3, and 1 process to quantify the detached pore networks.Figure 6 shows the thickness frequency number distribution in percentages for (a) the pore ganglia (NoBin) and (b) the detached pores (Bin3) for the 1000 3 voxel volumes of the analyzed sandstones.A volumeweighted size distribution is preferred for many applications (Van Dalen and Koster, 2012).However, only the numberweighted thickness distribution is provided because a direct comparison between the BB, FC, and EqD methods was intended.Table 3 shows the statistical results after a normal Gaussian fit to the rocks' data distributions, demonstrating the differences between the samples and the methods that were used to measure the pore thicknesses.
As shown in Table 3, the pore ganglia thicknesses are mostly below 20 µm, with an average of approximately 10 µm between the samples and the analyzed BB, FC, and EqD methods.Specifically, a high number of pore ganglia in S 3 , followed by S 2 and S 1 , can be attributed to the nature of much smaller pore structures that are found in such samples, as shown in Fig. 3a.The detached pore networks' thicknesses, on the other hand, showed considerable distinctions among the samples and methods: the sizes were larger for the FC method compared to the BB method, which indicates that the pore particles of the three analyzed rocks have distinct orientations from the acquired sample coordination system.The shift to smaller sizes, which is much more pronounced by the equivalent diameter method, indicates that the pores might be asymmetric and/or elongated; once calcu- lated, the particle size is reported as an equivalent spherical diameter, which is determined by the voxel number counting.Furthermore, the small differences that are observed among the BB, FC, and EqD analysis methods are discussed in the pore shape classification approach as follows.

Classification and quantification of the pore shapes in the analyzed sandstones
To validate the pore shape classification approach, distinct 3-D geometrical descriptors were investigated according to the diagrams in Fig. 2a, which allows qualitative particle description, and Fig. 2b, which allows pore shape classification    and subsequent quantification based on artificial object similarity.

Approach validation based on the 3-D visualization and geometrical descriptors
Figures 7, 8, and 9 depict 250 3 voxel volumes of rocks S 1 , S 2 , and S 3 , which show the results for five representative pores.These 3-D rendered pores were the five longest (by L), labeled from the 25 largest (by volume), and acquired after pore network detachment with Bin3 for each sample.Pore no. 3 (in red) in Fig. 7 is the same as in Fig. 1.The systematic choice of five equally categorized pores for each rock sample is appropriate for the discussion and validation of the results based on the 3-D visualizations.In Figs. 7, 8, and 9, the same colors that are attributed to the pore particles are depicted in the modified Zingg diagram (middle column) and equancy graphs (right column).Because of the unrealistically high S/ l ratios from the equivalent diameter descriptor, plotting this result within the 0.2 intervals of the equancy graph was impossible; therefore, the EqD descriptor is only discussed for the shape class diagram.Nevertheless, the equancy for the BB and FC descriptors was used to support the results of the shape class approach.By comparing the 3-D renderings in Figs. 7, 8, and 9, sample S 1 shows the largest pore volumes, followed by the S 2 and S 3 samples.Additionally, none of the five pores in S 1 are located in the "extremely non-equant" interval in the equancy graph for both the BB and FC descriptors, while at least three pores (no. 1, 4, and 5) in S 3 fell within this interval, indicating higher particle asymmetry for this rock.Although only qualitative evaluation can be drawn from the equancy diagrams (some particles fall within distinct intervals), the pores do deviate from being equant (regular shape) descending from S 1 to S 2 to S 3 .These results support good agreement with the proposed shape class approach for the three rocks when both the BB and FC descriptors were used.
However, when the bounding-box, Feret caliper, and equivalent diameter descriptors were investigated in the Zingg diagram for the pore shape classification and quantification, the results for the three rocks seemed to be more appropriate for the FC descriptor (compare e.g., pore 3 in Fig. 7).For this pore, the bounding-box descriptor indicates a cube-like shape; however, as shown in Fig. 1a, this cube shape is notably erroneous.Our assumption for the bounding-box mismatch is related to the main orientation of pore 3, which is neither perpendicular nor parallel to the sample z axis but is approximately within 45 • .In this case, the L, l, and S descriptors were inaccurately measured by the BB method, which will tend to overestimate the shapes as cubelike.The same behavior is observed for particles 4 and 5 in Fig. 8 and 2 and 3 in Fig. 9. Thus, the Feret caliper descriptor will be more adequate for materials that contain a reasonable number of pores with preferential directions, which match neither the normal nor parallel sample Z axis direction, providing enhanced pore morphological information over the bounding-box method.On the other hand, the given S/ l ratio is unrealistically high when evaluating pore shapes with the equivalent diameter descriptor (l), which proves that this descriptor must be handled with caution when characterizing very asymmetrical pore structures, such as those found in natural materials.Nevertheless, this parameter might be used to indicate the degree of particle asymmetry, as discussed in the following section.

Quantifying shapes of pore ganglia and detached pore networks
Selecting a region of interest (ROI) from the sample is required prior to any pore shape quantification or morphological analysis from a reconstructed 3-D image.To correlate many properties (e.g., porosity), the chosen ROI should be large enough to represent a sample's complexity and heterogeneity but small enough not to overwhelm the avail-able computing resources (Baker et al., 2012).This study's ROI volumes of 1000 3 , 500 3 , and 250 3 voxels are shown in Fig. 10, which notes the changes in the rod-(1), blade-(2), cuboid-(3), plate-(4), and cube-shaped (5) classes of the pore ganglia (NoBin) and detached pore networks after Bin3.In this figure, one can also observe the differences in the shape classes between the BB (blue), FC (red), and equivalent diameter (yellow) descriptors.No significant variations were observed from the 1000 3 to 250 3 volumes for each of the descriptor methods.For the pore ganglia's shape classes (dark colors), the changes were even less pronounced because of their much smaller size compared to the smallest ROI that was analyzed.For the three ROIs, most of the pore ganglia were shown to be cube-like, followed by plate-like (BB and FC methods) or cuboid-like forms (EqD method).Nevertheless, a slightly variation in the ROI size was observed for the detached pore networks' shape classes (bright colors) in S 1 and S 2 when the Feret caliper descriptor (bright red) was used.For these two rocks, the detached pores were predominantly plate-like in the 250 3 ROI, approximately equally plate-and cube-like in the 500 3 ROI, and mainly cube-like in the 1000 3 ROI.This change, which is notable for rocks with bigger pore structures, indicates that the pore shape quantification results might be snagged by excessively small sample sizes.Nevertheless, the systematic comparison of the analyzed ROIs, which is shown in Fig. 10, indicates reasonable reproducibility from the 1000 3 to 250 3 voxel volumes when performing pore shape analysis.Only the 1000 3 ROI is shown (Fig. 11) for the quantification and discussion of pore ganglia and detached pore network shapes.
Figure 11 shows the results in percentages and average values with standard deviations for each of the pore shapes (rods, blades, cuboids, plates, and cubes) in the analyzed rocks.From Fig. 11, one can see systematic changes between the pore shape classes (see values in percentage) within the evaluated Bin parameter and the BB, FC, and EqD descriptors.The NoBin results represent the pore ganglia shape's contribution (red color), whereas the detached pore network results after Bin10, 3, and 1 are shown in blue, green, and yellow, respectively.In the graphics, the pores tend to shift from plates to cubes from Bin10 to Bin1 for the detached pore networks in the BB and FC methods for the three reservoir rocks.The shift to cube shapes is expected because of the pore separation process, which obviously diminishes the particles' length to a value that is comparable to the width and thickness, shifting the S/ l and l/L ratios to the unit value.The results in Fig. 11 indicate the importance of choosing an appropriate Bin parameter, specific to each pore system, Bin3 being the best choice for the analyzed sandstones.A similar behavior was observed for the equivalent diameter descriptor results from Bin10 to Bin1, but with pore shapes that shifted from rod-like to cuboid-and cube-like.As already mentioned, the use of equivalent diameters to determine the thickness of the asymmetric pore structures in these sedimentary rocks creates unrealistic and mismatched results for the pore shape analysis.The errors increase with the pore morphology's asymmetry and divergence from equancy.Still, an interesting indication can be drawn from the EqD descriptor results for detached pore networks after Bin3: the rod-like shapes that were identified for the three sandstones, in addition to pore lengths (L) that were at least 3 times larger than the measured equivalent thickness (l), had l values that were around the same or 2 times smaller than the width (S).These extremely asymmetrical pore shape configurations can be visualized in the 3-D renderings in Figs. 7, 8, and 9, specifically, particles 1 (S 1 ), 3 (S 2 ), and 3 (S 3 ), which are the highest S/ l values.Thus, one could use the pore equivalent diameter in shape analysis to indicate the degree of pore heterogeneity/asymmetry based on rod-like shapes, which is equal to 6.70, 18.58, and 25.44 % for samples S 1 , S 2 , and S 3 .
As discussed in the literature (Anovitz and Cole, 2015;Chang et al., 2006;Okazaki et al., 2014), a number of pore form factors have been used to quantify pore shapes that might be correlated to the macroscopic properties of sedimentary rocks.Among these factors, flatness, elongation, and roundness are commonly used.Permeability, from a geological point of view, is still one of the most challenging properties to determine but is nonetheless a very important hydraulic property in solving accumulation and exploitation problems in the oil and gas industry.The permeability, k, depends on the parameters of the pore structure, including the porosity, Euler number (connectivity), and ge-ometrical factors, e.g., shapes.Pore shapes and quantities change during a reservoir's deposition (compare the distinct geological ages of the sandstones in Sect.2.1) because of several geological processes, such as burial diagenesis and compaction, which control the pressure solution, cementation, or grain rearrangement (Okazaki et al., 2014).Depending on the assumption or model of the pore microstructure's arrangement, many geometrical parameters might be correlated to predict permeability, such as for the Kozeny-Carman equation (Scheidegger, 1974;Tiab and Donaldson, 2004;Walsh and Brace, 1984), which is widely used by the petroleum industry.The simulation, estimation, or measurement of petrophysical properties from the investigated sandstones are beyond the scope of this paper, nonetheless, permeability, porosity, and specific surface area values previously measured in the laboratory (Halisch, 2013) are useful to observe reasonable reproducibility from the determination and quantification approach of rock pore shapes as presented in this work.For samples S 1 , S 2 , and S 3 , Klinkenberg permeability was equal to 555, 11, and 0.11 mD; specific surface area was equal to 0.347, 0.559, and 1.560 m 2 g −1 ; and mercury intrusion capillary pressure porosimetry was equal to 21.5, 16.1, 8.4 %.Additionally, Archimedes porosity ranged from 20 to 22, 16 to 17, and 9 to 12 %, while the bigger porosity in the low-field nuclear magnetic resonance measurements showed T 2 (in ms) around 100-1000, 10-100, and 0.1-100, for S 1 , S 2 , and S 3 .Thus, correlating the quantification results of pore shapes (observe mainly Bin3 and Feret caliper results in Fig. 11) of the analyzed sandstones, and considering that a cubelike shape imposes less fluid flux resistance compared to, e.g., a plate-like shape, one could assume that S 1 , which has the highest number of cube-like pores (58.80 %), is more permeable than S 2 (55.62 %) and S 3 (45.18%).Additionally, plate-like shapes are more common in S 3 (50.94%) compared to S 2 (42.01 %) and S 1 (39.49%).These findings match very well with the rocks' geological descriptions, ages, and Klinkenberg permeability measured in the laboratory; in particular, S 3 , comprised mainly of plate-like pore shapes, showed a permeability value 2 orders lower than S 1 and S 2 , and it is much older as well, thus more subject to grain compaction/deformation.The laboratorial results can confirm that the pore shape quantification approach obtained from the Feret caliper and Bin3 parameters faithfully represent the pore structure morphology of the analyzed sandstones.The organization of the pore system in these rocks is directly correlated to fundamental geological processes and evolution, which control underground hydraulic properties.

Conclusions and outlook
Three distinct reservoir rocks were systematically analyzed to characterize 3-D pore shapes based on x-ray µ-CT images.A modified Zingg diagram that associates particle shapes with similar artificial objects was created to validate distinct bounding-box, Feret caliper, and equivalent diameter descriptors.A detailed procedure was developed and described in this paper to quantify and draw faithful pore shape information concerning the real pore system.An important finding from the pore shape analysis was that the 3-D rock samples comprise main pore networks and several disconnected pore ganglia; the latter were suppressed from the pore volume samples to analyze the shape sensitivity of detached pore networks.The Euler number of these pore networks can be used to obtain a good quantification of the changes during the performed binning operations.
In the analyzed rocks, the most accurate particle descriptor was given by the Feret caliper because of the considerably high number of pores that were located neither perpendicular nor parallel to the Z axis of the acquired 3-D sample.Additionally, the most suitable marker extent parameter for performing the watershed algorithm in the Avizo Fire software was Bin3.Thus, most of the pore shapes were identified as being plate-and cube-like, equal to 39.49 and 58.80 % for S 1 , 42.01 and 55.62 % for S 2 , and 50.94 and 45.18 % for S 3 , respectively.Plate-like pores were evidently higher for the Flechtingen sandstone (S 3 ), which could be explained by the geological processes and history of this field compared to those of the S 1 and S 2 samples.
The pore thickness that was calculated using the equivalent diameter led to erroneous results in the Zingg shape classification for the three analyzed sandstones, which contain very asymmetric pore morphologies.However, the EqD descriptor might be used to indicate the degree of pore heterogeneity based on rod-like shapes, which were found to be equal to 6.70, 18.58, and 25.44 % for samples S 1 , S 2 , and S 3 , respectively.Additionally, the equivalent diameter was a faithful parameter for many investigations, e.g., obtaining pore size distributions from distinct techniques such as mercury intrusion capillary pressure when assuming capillary phenomena in parallel cylinders.Because of the growing worldwide interest in flow simulation in very complex pore structures, such as unconventional reservoirs, many techniques have been improved by allowing models to enter a pore shape parameter, e.g., in the gas adsorption method.In this context, the 3-D pore shape quantification that was described in this paper is essential for several studies from distinct research areas.This study is part of ongoing research to investigate the influence and impact of different pore shapes on the petrophysical properties of reservoir rocks, such as permeability, surface area, fractal dimension of the surface area, and complex electrical properties.Moreover, this research project will be extended to other types of reservoir rocks, e.g., carbonates and unconventional reservoir rocks from Brazilian fields, in the near future.

Fig
Figure 1.3-D geometrical descriptors for the bounding-box (a) and Feret caliper (b) methods with L, l, and M corresponding to D y , D z , and D x (a) and D max , D med , and D min (b).

Figure 2 .
Figure 2. Pore shape classification for equancy based on the S/L ratio (a); classes of pore-shape-like artificial objects (b) plotted on a Zingg diagram (modified fromBlott and Pye, 2008;Soete et al., 2015).

Figure 3 .
Figure 3. 1000 3 voxel renderings of the analyzed sandstones: grayscale image (a), five biggest (by volume) main pore networks (b), and residual pore ganglia that were excluded from the shape analyses (c).In (b) and (c), the labeled images show distinct colors for dissimilar pores.

Figure 4 .
Figure 4. 3-D and 2-D visualizations after the watershed algorithm is applied to the 250 3 voxel volume of S 1 (a), S 2 (b), and S 3 (c), with the labeled images showing distinct colors for dissimilar pores.

Figure 5 .
Figure 5. Quantification of the pore network fragmentation by using the Euler number (a) and number of segmented pores (b) during the binning operations for the three rocks and three different domain sizes, respectively.

Figure 6 .
Figure 6.Pore thickness (l) distribution of pore ganglia (a) and detached pore networks after Bin3 (b) as measured by the BB, FC, and EqD methods in the 1000 3 voxel volume of S 1 (yellow), S 2 (red), and S 3 (blue).

Figure 7 .
Figure 7. Pore shape results (S 1 ) for equancy and shape classes based on distinct descriptors for the five categorized pores in a 250 3 voxel volume.The lines shown in the graphics delimitate distinct particles shapes, as described in Fig. 2.

Figure 8 .
Figure 8. Pore shape results (S 2 ) of equancy and shape classes based on distinct descriptors for the five categorized pores in a 250 3 voxel volume.The lines shown in the graphics delimitate distinct particles shapes, as described in Fig. 2.

Figure 9 .
Figure 9. Pore shape results (S 3 ) of equancy and shape classes based on distinct descriptors for the five categorized pores in a 250 3 voxel volume.The lines shown in the graphics delimitate distinct particles shapes, as described in Fig. 2.

Figure 11 .
Figure 11.Pore shape quantification of samples S 1 , S 2 , and S 3 (1000 3 voxels) for the BB, FC, and EqD descriptors based on the proposed Zingg classification approach.

Table 1 .
3-D geometrical parameters of the detached particle from the main pore network of rock S 1 in the 250 3 voxel volume, as shown in

Table 2 .
3-D x-ray µ-CT data acquisition and processing parameters for the analyzed rocks.
a NLM: non-local means, b LN: local neighborhood.

Table 3 .
Thickness statistics for the pore ganglia and detached pore networks after a normal Gaussian fit to the distribution data, shown in Fig.6.SD denotes standard deviation.