A technique to identify the predominant pore direction in a porous medium and application to reservoir rocks

The study of the pore space structure of a porous medium has been very much improved with the aid of microtomographic imaging and its analysis through image processing. In this paper, a technique to identify the predominant pore direction (PPD) in the pore space is introduced and according to that the pore space can be partition as vertical (V), horizontal (H) or diagonal (D). The PPD technique has been developed for 2D and 3D spaces based on microCT images of a porous medium and can be used to both, pore or grain spaces. An implementation of the PPD has been done in an in-house computer program using Phyton. A set of application tests for 2D and 3D PPD partitioning is given, being, respectively, i) a synthetic binary image and a binary image of a Berea sandstone rock sample and ii) a Berea sandstone and a Carbonate reservoir rock core samples, both provided by the MicroCT Images and Networks of Imperial College London database, and a 3D grain partitioning of a trabecular bone structure. Additionally, the V, H and D PPD effective pore networks of a pore space are determined; the porosities for the PPD subspaces and for their effective pore networks are computed and results are provided. Finally, a brief discussion about the implementation and computational cost for the 2D and 3D cases is provided. It is worth to mention that the findings of this study may help for better understanding of directional fluid flow and mechanical stress in porous medium.


CT
-Micro-computed tomography I n×m -2D binary image of size n × m I n×m×k -3D binary image of size n × m × k L( i ) -2D Euclidean length along the ray for the angle i , i is the ray number N -Total number of rays launched from the pore pixel (2D) or voxel (3D) L( i , j ) -3D Euclidean length along the ray for polar angle i and azimuthal angle j , i the polar ray number and j the azimuthal ray number N -Number of polar rays N -Number of azimuthal rays p 0 -Seed pore pixel

Introduction
A porous medium is essentially treated as a physical structure divided into a void space, referred as pores, and a solid matrix, named grains. In many applications of porous media, the main interest is sometimes related to the pores, particularly when the research focus is on fluid flow and mass transport, while in others situations the focus is on the grains, especially when the structural behavior under load is the relevant aspect of study. These occurrences are observed, for instance, in porous media like reservoir rocks, soil or trabecular bone. Currently, most of the research on porous media microscale is based on imaging, particularly microCT ( CT) imaging of samples as shown in the following papers Sidorenko et al (2021), Sakugawa et al (2020), Irayani et al (2019), Andrä et al (2013a), Andrä et al (2013b). In addition, for simplicity, the 2D or 3D grayscale image is converted to binary, with the black pixels or voxels representing the pores and the white ones the grains, such that the color attribute of the pores is 0 and that of the grains is 255 (or 1). Notice that sometimes it is convenient to work with the complement by inverting it, turning white (255) to represent the pores and black (0) the grains.
When studying the fluid flow through a porous medium, the effective pores define the route through where the fluid may flow. In other words, the connected pores that have an inlet at one face and an outlet at another face are the ones that allow the fluid to move through the region or volume of interest (ROI or VOI). There are a couple of properties associated with the pores that have to be considered, such as the pore thickness, the pore connectivity and their degree of tortuosity, for instance. Additionally, it is also important to know what is the preferential direction of the pores in the ROI/VOI, what percentage of the pores are aligned horizontally, vertically or diagonally. The pore space main directions influence how the fluid will naturally or forced flow under the action of pressure or gravity, for instance.
The trabecular bone is a spongy-like structure, with the trabeculae as the solid matrix (grains) and the marrow filling the pore space. The predominant direction of the trabeculae is an important characteristic of the structure as it acts as load-bearing paths to spread out applied stress Roque and Alberich-Bayarri (2015). In soil physics, permeability and consolidation are the most variable soil properties. These parameters may vary with depth even in case of homogeneous soil layers. However, due to the anisotropy, the coefficient of consolidation and the coefficient of permeability in the horizontal direction are typically different from their values in the vertical direction Laskar and Pal (2018). The vertical and horizontal pore directions play a relevant role for water flow and nutrient transport to different ground depths as well as for the water retention curve.
In the field of carbonate reservoir rocks, porosity and permeability are quite important rock properties; they are, respectively, essentially responsible for the storage capacity and flow of hydrocarbons in the reservoir. Additionally, the pore connectivity, tortuosity and pore predominant directions are determinant for the fluid and electrical current to flow through the rock Tiab and Donaldson (2016).
The fluid flow in a porous medium depends on the pore space. On the other hand, the pore network influences the directional permeability, and thus, the predominant direction of the pore network constrains the fluid flow. For instance, air permeability is an important property of fabrics as it can be a decisive parameter for the functionality of the fabric. It has been shown that the pore direction in woven fabrics affect the air permeability Havlov (2014).
The special core analysis laboratory (SCAL) methods used to estimate porous medium quantities are costly and time-consuming, but more recently digital rock analysis Kalam (2012) has significantly minimized these downside effects. Digital rock physics (DRP) Berg et al (2017), Al-Marzouqi (2018) is a breakthrough technology to digital oilfield that merges high-resolution imaging, image processing techniques and numerical simulations on the pore scale based on porous medium core samples.
In this paper, it will be presented and discussed a technique to identify the predominant pore direction in a porous medium and through it do a partition of the pore space according to the pore orientation. The layout of the paper is as follows: First, Sect. 2 presents the theoretical framework with the 2D (2.1) and 3D (2.2) general equations defining the predominant pore directions and providing the partition of the pore space into vertical, horizontal and diagonal pores; then, Sect. 3 presents the results of the application of the techniques to 2D and 3D different case studies. Section 4 provides a brief discussion on the necessary calculations and computational costs to compute the 2D and 3D pore predominant directions. Finally, Sect. 5 concludes the paper.

Theoretical framework
For the purpose of the techniques that will be presented here, it will be considered 2D and 3D binary images of the porous medium. Let us consider initially a 2D binary image I of a porous medium and assume that I is an image of size n × m pixels that a pixel p ∈ I is represented by a square and has coordinates p = (x, y) set at the center of the pixel, following the standard reference frame adopted in image processing, where the origin pixel (0, 0) is at the top left corner of the image. Here, it is adopted that in a binary image of a porous medium, the solid matrix (grains) is represented by black pixels with a 0 color attribute and the voids (pores) are white pixels with a 255 color attribute. Likewise, 3D binary images will have the same attributes as for the pixels, with the voxel coordinates set at the center of the voxel.

2D directional partition of pores
In this, it will be introduced the foundations for the identification of a 2D pore orientation in a porous medium. A pore will be classified according to their horizontal (H), vertical (V) or diagonal (D) predominant direction, and the pore (or grain) space will be partitioned into a subset of pores according to their predominant orientation. For that, given a single pore pixel, say p 0 = (x 0 , y 0 ) ∈ I , the problem consists in how to identify what is its predominant direction. Figure 1 gives an illustration of the three major pore predominant directions defined in this study.
To determine which is the predominant direction of a pore, according to the partitions of the plane given in Fig. 1, the star length distribution principle Smit et al (1997) will be applied from a pore pixel p 0 = (x 0 , y 0 ) with k = 16 rays (straight lines) being launched in every = ∕6 interval, performing a total of N = 192 rays. According to the star length distribution principle, to identify the pore pixel predominant direction we look for the first grain (black pixel) to be intercepted by the ray. If no pixel is found, we take the pore pixel that has reached the boundary of the image I n×m . Then, the Euclidean distance along the ray L( i ), i = 0, … , 191 is computed from the seed pixel p 0 to the first grain one or to the boundary, as Fig. 2 illustrates.
To obtain the 2D pore predominant direction, the following expressions are calculated for each partition: V e r t i c a l p a r t i t i o n V : ∕3 ≤ < 2 ∕3 a n d 4 ∕3 ≤ < 5 ∕3.
where i = 2 i N corresponds to the angle of each ray launched from the seed pore pixel p 0 . The predominant direction of p 0 is obtained from that means the direction with the largest value defines the predominant pore direction (PPD) of p 0 .
Remind that H , V and D are mathematical partitions such that they satisfy H ∩ V ∩ D = � and H ∪ V ∪ D = PS , where PS is the pore space of the sample. (1)

3D directional partition of the pores
The 3D pore direction partitions are a natural extension of the 2D case. Figure 3 illustrates the 3D directional partitions. Let us consider the usual image processing spherical coordinates and that from a pore voxel v 0 a set of rays are launched according to (r, , ) , where r are the rays departing from v 0 to the nearest grain or to the boundary of the VOI, is the polar angle, and is the azimuthal angle. The Euclidean distance along a ray will be given by To identify the pore voxel predominant direction, let us consider that 16 rays are launched from a pore seed voxel v 0 within each interval = ∕6 and = ∕6 , leading to N = 24 polar rays and N = 48 azimuthal rays, providing a total of N = 1152 radial rays from each pore voxel. The Euclidean distance along a ray will be given by L( i , j ) , 0 ≤ i ≤ , i = 1, … , 24 and 0 ≤ j < 2 , j = 1, … , 48 . Notice that now the total number of rays is 6 times greater than for the 2D case, which causes an increase in the computational cost associated with the size of the problem, in other words, to the sample volume, its porosity and mathematical operations.
The 3D directional partitions are given by: The predominant direction of a pore voxel v 0 is obtained from in other words, the direction with the largest value defines the predominant pore direction (PPD) of voxel v 0 .

PPD applications
In this section, application tests of the PPD identification and directional partitioning of the pores are presented. First, the PPDs for a 2D synthetic binary image are evaluated and then for a 2D binary image of a Berea sandstone rock sample. These are two simple cases just to illustrate the feasibility of the technique. For the 3D case, three application tests are considered. First, a Berea sandstone and a Carbonate reservoir rock core samples are investigated, both were obtained from the MicroCT Images and Networks of Imperial College London database and, finally, an application considering the grain partitioning of a distal-radius trabecular bone structure. These are well-known type of porous medium structures; in particular, the rock core samples are typical of a sandstone and a carbonate reservoir.
The 2D and 3D algorithms for directional partitioning of a pore space were implemented in an in-house software using Phyton running on a Windows Ⓡ OS. In addition, the computation of the number of pixels (2D) or voxels (3D) in the sample, the pore space porosity, the partition porosities, the effective pore networks and their porosities is also provided by the software.

2D synthetic and Berea image samples
As a simple and effective test for the PPD technique, a 2D synthetic binary image was created with 6 different geometrical pictures that include straight horizontal, vertical and diagonal lines, a curved line, two solid rectangles and two solid ellipses. Figure 4 shows, on the left, a 2D synthetic binary image and on the right the identification of its PPD partitions: horizontal in red, vertical in blue and diagonal in green. The straight lines show correctness of the 2D PPD. As the curved line is a digital curve, that is, a discrete line formed by pixels, the 2D PPD H, V and D partitions are present along the same line. The same occurs for the solid pictures. Remind that these are digital pictures and the 2D PPD is defined for an angular interval. Figure 5 shows the application of the 2D PPD equations given in 2.1 for a 2D binary image of a Berea sandstone sample seen on the left and its PPD partitions on the right.
The 2D PPD technique can be applied for thin sections Saxena and Mavko (2016); Saxena et al (2017), and a 3D reconstruction of a sample can be made; nevertheless, this approach is not accurate as it does not take into account the rays in the z-direction as indicated in Sect. 2.2.

3D PPD applications
Let us consider the partitioning based on the 3D PPD Eqs. 5, 6, 7 and 8 applied to the Berea sandstone and Carbonate C1 reservoir rock core samples to show their pore predominant direction subspaces, and to the distal-radius trabecular bone to show its grain predominant direction subspaces.

Berea sandstone
The first sample is a Berea sandstone that has volume size 400 × 400 × 400 voxels and 5.345 m resolution. Figure 6 shows the 3D full pore space in grayscale; red, blue and green are the horizontal, vertical and diagonal PPD subspaces, respectively, and finally the full pore space exhibiting the union of the three PPD subspace partitions. Just the horizontal, vertical and diagonal PPD partitions subspaces are displayed in Fig. 7. Figure 8 shows the effective pore networks of the full sample and the H', V' and D' PPD partitions obtained from the effective pore network.
Notice that the H', V' and D' PPDs obtained from the full effective pore network shown in Fig. 8 are not necessarily effective pores in the respective directions. Figure 9 shows the effective pore network from each one of these directional partitions. The vertical PPD partition shown in Fig. 8 does not have any effective pore network. Table 1 gives the Berea sandstone sample porosity , the PPD porosities {H,V,D} , the effective pore space porosity e , the PPD porosities {H � ,V � ,D � } obtained from

Carbonate C1
The second application is for the Carbonate C1 sample from the database. Its volume size is 400 × 400 × 400 voxels and 2.85 m resolution. Figure 10 shows the 3D full pore space, the horizontal, vertical and diagonal PPD pore subspaces and finally the full pore space exhibiting the union of the three pore subspace partitions. Figure 11 shows the effective pore network of the full sample and the H', V' and D' PPD obtained from the effective pore network. Figure 12 shows the effective pore networks from each one of these partitions. Table 2 gives the Carbonate C1 sample porosity , the PPD porosities {H,V,D} , the effective pore space porosity e , the PPD porosities {H � ,V � ,D � } obtained from the effective pore space and PPD porosities {eH,eV,eD} of the effective directional pore networks. Notice that 92% of the pore space are effective pores. From the H', V' and D' partition, only 89%, 34% and 32% of the pores form the fundamental directional effective networks.

Trabecular bone
Let us consider now a 3D distal-radius trabecular bone image obtained by CT from an ex vivo sample, which has a volume of 257 × 257 × 240 voxels and resolution of 34 m. This sample corresponds to an osteoporotic trabecular bone structure; additional details about this sample, named 269, can be found in Roque et al. (2013). Figure 13 shows the 3D full grain space of the distalradius trabecular bone, the horizontal, vertical and diagonal predominant grain direction (PGD) partitions and the full grain space exhibiting the union of the three partitions. Figure 14 shows the effective grain network of the sample and its H', V' and D' PGD partitions. Figure 15 shows the effective fundamental directional grain network from each one of these partitions. Table 3 gives the trabecular bone volume fraction ( = BV TV , BV is bone volume and TV is total volume of the sample), the PGD TBVF {H,V,D} , the effective TBVF e , the PGD TBFV {H � ,V � ,D � } obtained from the effective TBFV and PGD TBFV {eH,eV,eD} of the effective directional trabecular networks. Notice that 98% of the trabecular bone space are effective trabeculae. From the H', V' and D' partitions, only 10%, 59% and 2% of the trabeculae form the fundamental directional effective networks.
In this section, a set of two 2D and three 3D application tests of the PPD have shown how the technique can be used to identify the pore orientation on 2D or 3D microCT images, respectively, and partition the pore space of the   sample (reservoir rocks for pore and trabecular bone for grain) into horizontal, vertical and diagonal directional subspaces. In addition, the effective pore networks were identified and their PPDs partitions were obtained with their porosities calculated and the results are summarized in Tables 1, 2 and 3.

Discussion
The predominant direction of a pore in a 2D or 3D space is dependent on the reference frame orientation and on the length of the rays. The image processing frame follows the right-hand rule, as such a pore direction changes due to frame rotations, and a change in the pore structure itself may change the length of the rays, which provoke orientation changes as well. A pore orientation is not affected by a simple change in the sample porosity as long as the pore structure where it belongs to is not modified. The PPD partitions (see Figs. 1 and 3) can be easily adjustable, leading to narrower or wider pore orientations. The effective pore networks of a sample show the set of pores that are connected to each other from an inlet face to a distinct outlet face. They form the pore network that allows a fluid to flow throughout the sample. The H, V and D effective fundamental directional pore networks give an indication on how the pores get directionally connected to each other. This induces a natural direction choice for fluid flow under the action of a strength. From Table 1, it can be seen that the effective horizontal porosity of Berea sandstone is higher than the vertical and diagonal ones, which gives an indication that individually there will be more facility to horizontal flow than in the other two directions. Similarly, in the Carbonate C1 case ( Table 2) the H effective fundamental directional porosity is higher than the other two directions, which indicates that a fluid would flow more naturally in the horizontal direction instead.
In the last case, if we look at the trabecular structure, the vertical effective fundamental directional trabecular bone volume fraction is higher than in the other two directions. According to previous studies Roque et al. (2013), this result was actually expected as the trabeculae are primarily aligned along the direction where they suffer more frequent stress; in such sample, it corresponds to the vertical direction (z-direction, distal direction) to provide more strength to the bone structure.
In Figs. 9 and 12, it can be seen that the effective fundamental directional pore network and in 15 the effective fundamental directional trabecular bone network were obtained from the respective H', V ' and D' partitions The H',V' and D' partitions shown in Figs. 8,11 and 14 were obtained from the effective pore (or grain) network of the respective samples. However, in the process of identifying these partitions some of the effective pore (or grain) network loses connectivity between previous connected pores that had distinct predominant directions, turning to directional partitions where the pores (or grains) are no longer effective in the corresponding direction. That is why the effective fundamental directional pore (or grain) networks shown in Figs. 9, 11 and 15 have a much less porosity (or TBVF) as shown in Tables 1, 2 and 3.

Algorithm implementation and computational cost
The main difficulty to implement the algorithms to compute the 2D and 3D PPDs is to find out the first black (grain) pixel/voxel that is intercepted by the ray or, otherwise, the white (pore) pixel/voxel that are at the boundary of the region/volume of the sample. Suppose a 2D image of size n × m with P 2D pore pixels. From each pore pixel, a total of 192 rays is launched (see Sect. 2.1), and thus, there will be 192 2D n m rays launched for the sample, where 2D is the 2D sample porosity, and for each one of the pores, there is the need to find out the first black pixel (grain) that is intercepted by the ray or the pixel at the boundary along the ray direction. For the example of 2D Berea sandstone given in Fig. 5, whose size is 400 × 400 pixels and porosity 2D = 0.2012 , the cost is ≈ 6.18 × 10 6 searches.
For the 3D image of size n × m × k , with P 3D pore voxels, a total of 1152 rays are launched from each pore voxel (see Sect. 2.2); therefore, there will be 1152 3D n m k rays launched for the volume sample, where 3D is the 3D sample porosity, and for each one of the pores, there is a need to find out the first grain voxel that is intercepted by the ray or the voxel at the boundary along the ray direction.
If we just consider the total number of rays launched from each pore voxel as the size of the problem, the order of magnitude of the problem from 2D to 3D grows 6k 3D 2D times. The 2D search algorithm for the first black pixel intercepted by the ray is relatively simple; however, for the 3D case, it is more complex, which turns the computational cost much higher yet, demanding more computational resources for large 3D samples. The computational cost may impose a certain limitation for the PPD determination when the size of the samples and/or when the number of rays launched from the source pore pixel (2D) or pore voxel (3D) is increased.
The PPD and effective pore network identification algorithms for the 2D and 3D cases were implemented in an in-house software using Phyton. The applications presented in Sect. 3.2 were run in a PC with i7-8565U CPU, 1.80GHz, quad-core and 16Gb of RAM, running under Windows operating system. Table 4 gives the number of rays launched for each sample, the CPU time spent to run the 3D H, V, D PPD for the Berea sandstone and Carbonate C1 samples and the PGD for the distal-radius trabecular bone sample. It is also given the CPU time spent to find the effective pore (EF) or grain (EG) networks for these samples.

Summary and conclusion
In summary, in this paper we have presented and discussed: -A technique and its computational implementation to identify the predominant pore (or grain) direction (PPD) as horizontal, vertical or diagonal, in a porous medium. -The partitioning of a 2D or 3D pore (or grain) space of a porous medium according to the PPD. -The algorithms for 2D and 3D PPD partitioning were implemented in Phyton, and a set of application tests were done: First, it was considered a single 2D synthetic image which included vertical, horizontal and diagonal digital straight lines, a curved line and solid pictures, and then a Berea sandstone sample. These have shown the accuracy of the technique. Secondly, a 3D Berea sandstone and Carbonate reservoir rock core image samples, that are available at the MicroCT Images and Networks of Imperial College London database, were considered as they are very typical of reservoir rocks. In addition, a 3D application was done considering the grain partitioning of an osteoporotic distal-radius trabecular bone structure. -The effective pore network of the 3D pore space was obtained, and the PPD and partition of the effective pore network were computed for each sample. -The 2D and 3D pore space porosities effective porosities, H, V and D PPD subspace porosities and their respective PPD effective pore network porosities were also obtained for the reservoir rocks and, similarly, for the grain volume fraction of the distal-radius trabecular bone, whose results can be found in Tables 1, 2 and 3.
The technique presented in this paper is novel providing an additional tool to improve the digital analysis of porous medium properties with respect to their pore orientations, as such it needs to be further explored. Currently, work is in progress to investigate how the PPD influences the pore network connectivity, its tortuosity and the permeability of reservoir rock core samples. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creat iveco mmons. org/ licen ses/ by/4. 0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.