Reliability of Algorithms Interpreting Topological and Geometric Properties of Porous Media for Pore Network Modelling

Pore network models (PNMs) offer a computationally efficient way to analyse transport in porous media. Their effectiveness depends on how well they represent the topology and geometry of real pore systems, for example as imaged by X-ray CT. The performance of two popular algorithms, maximum ball and watershed, is evaluated for three porous systems: an idealised medium with known pore throat properties and two rocks with different morphogenesis—carbonate and sandstone. It is demonstrated that while the extracted PNM simulates simple flow (permeability) with acceptable accuracy, their topological and geometric properties are significantly different. This suggests that such PNM may not serve more complex studies, such as reactive/convective transport of contaminants or bacteria, and further research is necessary to improve the interpretation of real pore spaces with networks. Linear topology–geometry relations are derived and presented to stimulate development of more realistic PNM.

1 3 1 Introduction 1 3 (Varloteaux et al. 2013a, b;de Jong 2006;Valvatne et al. 2004;Blunt et al. 2002). The use of the latter allows for the explicit modelling of wetting layers, which is difficult to achieve with the cylindrical geometry. Two general types of PNM exist: irregular pore network model (IPNM) and regular pore network model (RPNM). The most notable difference between them is that IPNMs do not have any pattern in the location of the pore bodies in the network, whereas RPNMs usually tessellate the volume of interest using repeating unit cells of a fixed geometrical shape such as a cube Matthews and Spearing 1992;Reeves and Celia 1996), a rhombic dodecahedron (Vogel and Roth 2001) or a truncated octahedron (Jivkov et al. 2013). The pores are then usually located at the centres of each unit cell and the interconnecting throats pass through the faces of the cell, although this is not always the case (Raoof and Hassanizadeh 2010;Xiong and Jivkov 2015;Xiong et al. 2016). The differences in the location of the pores predetermine the capabilities and limitations of each model type. IPNMs are usually constructed by mapping the "real" interior structure of the porous samples into an equivalent pore network.
The "real" microstructure is usually obtained through 3D imaging techniques such as micro-CT (Jivkov et al. 2013) and serial section focused ion beam/scanning electron microscopy. Such models are sample specific and may not necessarily be statistically representative of the medium as a whole (Ebrahimi et al. 2013;Dong et al. 2009;Rabbani et al. 2014;Van Marcke et al. 2010;Lindquist et al. 1996). Even though the volumes obtained by XCT are generally small enough to be computationally feasible for modelling transport process using any of the direct methods, IPNMs can provide valuable statistical information regarding the effect of geometry and topology of the porosity on transport properties for a given medium. This information could then be implemented into a RPNM, which is not bound by the volume restrictions of the 3D CT image. A RPNM constructed within a larger volume allows for capturing statistics from a number of imaged samples and for calculating transport over dimensions closer to the engineering length scale (Xiong et al. 2016).
The overall aim of this work is to provide the diverse microstructural information required for the construction of physically realistic RPNM for longer-scale simulations of transport. Critical to such constructions is the quality of the pore structure statistics extracted from XCT images. Previous comparisons between different pore network extraction algorithms suggest a number of discrepancies in between algorithm outputs and pose questions about the reliability of the data (Dong et al. 2008). Maximum ball (MB) (Dong et al. 2009) and the watershed algorithm (WA) (Rabbani et al. 2014) are widely used in current pore network studies (Chen et al. 2017;Li et al. 2017;Ju et al. 2017;Rabbani et al. 2017;Kelly et al. 2016;Li et al. 2018). However, to the best of our knowledge, WA and MB methods have not been compared comprehensively in the previous literature. Therefore, our main objective is to evaluate the performance of the two pore network extraction algorithms in terms of derivable statistical data by applying them to micro-CT images of different resolutions and to a synthetic pore space with known characteristics. This is intended to facilitate future selection of appropriate network extraction algorithm.
In the first instance, the two algorithms are tested on an idealised pore structure having known pore throat geometry statistics. Subsequently, the algorithms are applied to two benchmark materials-carbonate rock and sandstone. These rocks are selected due to the differences in their distinctive pore spaces resulting from the processes that formed them. Sandstones are sedimentary rocks composed of lithic fragments and detrital mineral grains from 63 to 2 mm in size sourced from other igneous, metamorphic and sedimentary rocks under erosion and deposited under a range of marine, fluvial, 1 3 lacustrine and volcano-sedimentary environments. On the other hand, carbonate rocks are the result of biological processes or chemical precipitation, which also determine the grain sizes although reworking processes may reduce clast/grain size. Carbonates also have a higher potential of intra-particle porosity due to the presence of voids within structures of shelly organisms that may have formed them (Nichols 2009).
In order to assess the effect of CT resolution on the extracted network statistics, the two rocks are imaged at different resolutions. For all studied cases (idealised and natural pore systems), single-phase permeability values are calculated by the constructed IPNM and by direct fluid mechanics simulations on the XCT images to demonstrate IPNM performance for simple flow simulations.

Idealised Porous Medium
A sample having pre-defined pore space statistics was first considered for evaluation of the two network extraction algorithms. It was created by starting with random distribution of 800 points within a predefined volume of 400 × 400 × 400 voxels. A Voronoi tessellation of the volume around the points was performed next using Voro++ (http:// math.lbl.gov/voro++/). Spherical pores were then assigned to all points with volumes proportional to the corresponding Voronoi cell volumes; the assigned pore radii were varied between 1 and 31 voxels. Further, cylindrical throats were introduced between pairs of pores in neighbouring Voronoi cells (cells with common faces). Throat crosssectional areas were made proportional to the areas of the cell faces they cross so that any intersections of throats in the physical space were avoided. Once the pore network was generated, additional checks were carried out in order to reduce any ambiguity in the interpretation of the features. In particular, the algorithm did not allow for throats with larger radii than the pores they connect. The last step of the generation of the idealised porous media was to convert the volume containing the artificial pore network into a voxelised representation. This was accomplished by checking whether the coordinates of every voxel from the volume lay within the volumes of any of the pores or throats from the originally generated network. If so, the voxel was marked as empty, or else it was marked as solid. The idealised network is shown in Fig. 1, and more information on voxelised representation and interpretation of the pore space could be found in Dong et al. (2009).

Benchmark 3D Image Data Sets
Two different types of sedimentary rocks were imaged at two different resolutions: a carbonate rock and Hollington sandstone. The CT scans and the segmentation of the carbonate rock samples have been reported previously by Dong et al. (2007). The Hollington sandstone samples were obtained from a quarry near Stoke-on-Trent, England, and imaged as part of this work. Examples of the two materials with scanning resolutions are shown in Fig. 2

Carbonate Data Sets
As described by Dong et al. (2007), the scans of the carbonate rock were conducted using Phoenix|x v|tome|ex micro-CT scanner. A total of 720 projections were collected as the sample was rotated 360°. Both scans were acquired at a voltage of 90 kV and a current of 100 μA. The voxel size of the fine-resolution scan was 2.85 microns (Carbonate H), and that of the coarse-resolution scan was 5.35 microns (Carbonate L). These voxel sizes translate into spatial resolutions (recognition of features) of approximately 8.6 and 16 microns, respectively. After image acquisition, data sets were reconstructed using the Sixtos software and then binarised.

Coarse-Resolution Scan
A sample of approximately 1 cm 3 (Sandstone L) was scanned using a Nikon Custom Bay CT system. The voxel size was 8.35 microns giving a spatial resolution around three times the voxel size, i.e. 25 microns. The specimen to source distance was 58.5 mm, and the specimen to detector distance was 1400 mm. During the XCT analysis, the specimen was rotated over a 360 o rotation range collecting 3142 projections (1 s/projection) per scan using an accelerating voltage of 60 kV and a beam current of 200 μA. After image acquisition, data sets were reconstructed using the Nikon Metris CT-Pro reconstruction software (version 2.2.4693.17506, 10 November 2012) and a voxel size of 8.35 microns was obtained.

Fine-Resolution Scan
The sandstone was crushed and pieces around 1 mm × 2 mm were put in a Kapton tube to ensure minimum X-ray absorption. A specimen, denoted hereon by Sandstone H, was placed on the Zeiss Versa XRM-520 CT system rotation stage with a specimen to source distance of 11.0 mm and a specimen to detector distance of 15.1 mm. During the CT scan, the specimen was rotated over a 360 o rotation collecting 1601 projections (1 s/projection) using an accelerating voltage of 80 kV and beam current of 87 μΑ, voxel size of 1.42 microns and a spatial resolution around 4.3 microns. After image acquisition, data sets were reconstructed using the Zeiss Reconstructor Scout-and-Scan reconstruction software (version 11.0.4779) and a voxel size of 1.42 microns was obtained.
After reconstruction, the data were loaded into Avizo standard 7.0 (Visualization Sciences Group, Bordeaux (VSG), France) for examination of the virtual slices and 3D volume renderings. For Sandstone H, a sub-volume of size 350 × 350 × 350 voxels was extracted from the interior of the sample to minimise any effects from the sample preparation. For Sandstone L, the extracted cubical sub-volume was of size 400 × 400 × 400 voxels. The samples underwent non-local means filtering, available in the Avizo software, prior to segmentation. The same software was used to segment the pore space from the solid phase. This was done manually by inspecting greyscale values along a longitudinal and transverse slice of the specimens to ensure accurate results. The resulting segmented volumes were exported for further processing as binary files.

Maximum Ball (MB) Pore Network Extraction Algorithm
The maximal ball algorithm, initially introduced by Silin et al. (2004) and Silin and Patzek (2006) to calculate dimensionless capillary pressures and then further developed by Dong and Blunt (2009), is used to identify the pore locations and generate equivalent pore networks. The throat identification algorithm, however, is based on a watershed segmentation of a distance map. The distance map is defined as the distance of each void voxel to the nearest solid wall. It is approximated as the distance of the void voxel centre to the centre of the nearest solid voxel minus half the voxel length. After computing the distance map, the next step is to find the location of pore centres as the local maximum of the distance map, which is obtained by generating a hierarchy of maximal balls inscribed in the pore space.
The maximal-ball hierarchy is generated by first assigning the balls to pore voxels uniformly throughout the pore space, one every two voxels in each direction. In this work, we did not allow balls less than three voxels in diameter, to filter out the high frequency roughness on the solid walls. Then, the balls that are fully overlapped are eliminated and the rest, called maximal balls, are used to generate the maximal-ball hierarchy. The pairs of maximal balls that partially overlap are identified, and the smaller ball is considered as a child of the bigger ball. The local maximum of the distance map is obtained as the collection of the maximal balls that are not a child of any other maximal ball, and a pore index is assigned to each of the local maxima.
To identify and characterise the throats, we use a watershed segmentation algorithm to obtain the exact boundaries of the individual pore bodies. To achieve this, the index of the maximal balls at the centres of the pores from the previous step is painted on the original image and then a region growing algorithm is used to grow these indices towards voxels with smaller radii until filling all the pore space. A throat is assigned wherever voxels from two pore bodies touch each other.
The last stage in the network extraction is computation of the throat and pore properties. The properties of the pores and throats computed in this work are compatible with the description used by network flow simulator codes (Valvatne and Blunt 2004;Oren and Bakke 2003). In summary, each pore and throat is assigned a volume, inscribed radius, a length and a shape factor which are discussed in detail by Dong and Blunt (2009).

Watershed Algorithm (WA) for Pore Network Extraction
The watershed method was initially pioneered by Thompson et al. (2006) and Sheppard et al. (2006) for analysing X-ray tomography images of rocks. Here, we provide a workflow for segmenting the porous space using the watershed algorithm (Rabbani et al. 2014).

Pre-processing
Here, the binary 3-D image is prepared for segmentation (majority transform, noise filtering and removing small objects). If the input image is noisy, the resulted image will 1 3 be over-segmented (Rabbani and Ayatollahi 2015). Majority and any other noise-filtering methods help to avoid small unconnected pores, which threaten the quality of results.

Processing
A city block distance transform is applied to the pore space of the rocks followed by median filtering to avoid detection of the prolate shapes (Wildenschild and Sheppard 2013). Now, the watershed algorithm, illustrated in Fig. 3, is applied to detect the ridge surfaces (Rabbani et al. 2016).
In three dimensions, the distance transform creates a cloud-like object within the porous space, in Fig. 3a. The denser parts of the cloud are located at the pixels farthest from the pore boundaries. The watershed algorithm initially places virtual fluid at the thicker parts of the clouds and gradually enlarges the fluid kernels to reach the outline, in Fig. 3b, c. Pixels, in which the filling fluids from two different pores touch, form the ridge surface, in Fig. 3d. Ridge surfaces are detected and their surface area measured. They represent the throats of the network. If these surfaces were removed, the remaining parts would be isolated pores. The pores are then labelled and their volumes evaluated (Rabbani et al. 2014).

Post-processing
In this step, we gather information to represent the extracted network. First, we find which pores are connected and the size of the connecting throat. The approach used by Rabbani et al. (Rabbani et al. 2016) is computationally cost-effective where they have dilated the binary image of the throats by a 1 × 1×1 or 2 × 2×2 cubic constructing element. Then, this dilated binary image is multiplied by the labelled image of throats and labelled image of pores. It Fig. 3 Processing steps of pore network extraction using watershed algorithm: a the realistic geometry of two connected pores; b the schematic distance map of the pore space; c growing fluid kernels within each pore; and d first contact between the fluid of two adjacent pores that has formed some pixels of ridge surface (orange line) 1 3 should be noted that this dilated image of throats has one for the pixels within and around the throats (ridge surfaces) and has zero for all other pixels. The result is an integer triplet: two different labels for adjacent pores and one label for the connecting throat-see Fig. 4c where a throat with label 1 couples two pores with labels 5 and 8. By analysing all triplets, the connections are delineated to construct a connectivity/incidence matrix describing which pores are connected via which throat.
In the processing step, we calculate the pore volumes and throat cross-sectional areas. The idealised pore is a sphere with the same volume as the real one, and the idealised throat is a cylinder with the same cross section area as the real throat. The idealised pores are centred at the centroids of the real pores. In addition, the length of a throat is obtained by subtracting the radius of adjacent pores from the distance between their centroids.
Finally, this network is tuned to match the porosity of the underlying image. For this, we can multiply the radius of all pores and throats by a correcting coefficient to obtain the pores and throats as shown in Fig. 4d.
The application of this algorithm to an isolated pore space in the Hollington sample is illustrated in Fig. 5.

Qualitative and Quantitative Assessment of Pore Space Interpretations
The images of the idealised and real materials, described in Sects. 2.1 and 2.2, respectively, are used to evaluate the performance of the two pore network extraction algorithms, described in Sect. 3, by the quality of the pore network statistics they provide. Figure 6 gives a qualitative insight into the way the two algorithms interpret the pore space of the idealised pore structure and the sandstone at the lower resolution.
It can be observed that MB tends to produce smaller pores in general, and particularly for sandstone, and in many cases interprets a pore space as a throat where WA interprets the same pore space as a pore. Further MB appears to miss a substantial fraction of the pore volume, while WA appears to shift the positions of the interpreted pores, possibly because of taking surrounding throats as part of the pore when calculating its volume and position. In general, both algorithms tend to overestimate the number of pores, as observed clearly in the idealised case, by splitting pore spaces, which should have been interpreted as throats. It can be argued that the same is valid for the interpretation of Sandstone L. These qualitative observations are quantified with the following data for extracted porosity and pore and throat numbers. Table 1 provides data for the geometry of the analysed images and the volume fraction of the pore space (porosity) as determined by Avizo by counting pore space voxels and as calculated from the irregular pore networks extracted by the two algorithms. As expected, the "real" porosity calculated by Avizo is dependent on the resolution: Higher resolution images (Sandstone H and Carbonate H) show higher porosity because a number of additional features of the pore space are revealed. Evidently, the constructed networks omit fractions of the real pore spaces. For the watershed algorithm (WA), these fractions are relatively small, showing that WA maps the geometry of the real pore space to the geometry of the network with acceptable accuracy on average. For the maximum ball (MB), the omitted pore space is particularly large, confirming the qualitative observation from Fig. 6. The reason for this is that MB reports the inscribed radius of a pore, while WA computes the effective hydraulic radius. In addition, the MB algorithm also reports the volumes of the features as extracted from the tomography and they could be tracked separately from the underlying shapes of the pores and throats in the extracted networks.
Tables 2 and 3 provide data for the pore and throat number statistics, respectively, illustrating how the two algorithms interpret differently the notions of pore and throat. The results for the idealised model confirm the observation from Fig. 6 that both algorithms interpret a large number of regions in the original throats as pores, which split original 1 3 throats and lead to increased numbers of pores and throats in the resulting network. This affects the coordination spectrum and the overall size distributions, which will be presented and discussed in Sect. 4.2. Further, it is noted that MB consistently produces about twice as many pores and 2-3 times more throats than WA for each material and each resolution, suggesting that the potential throat-splitting effect is much stronger with MB. This difference is also due to the fact that the WA employs additional filtering before pore network extraction, which eliminates some of the small features. Such differences in the reported number of pores/throats from different algorithms have been previously observed (Dong et al. 2008).
In line with the porosity results, the number density of resolved features increases significantly on going from low resolution to high with the pore density rising by a factor of 4× for the carbonate (2× higher resolution) and 70×-80× for sandstone (6× higher resolution) and throat density increasing by a factor of 10 for carbonate and 110-180 times for sandstone. This suggests that the sandstone sample possesses substantially more complex pore space across length scales. Details on the geometric and topological properties of the pore networks constructed by the two algorithms are given in Sect. 4.2. These include data typically used in the construction of regular pore network models, such as the pore coordination number spectrum (distribution) required for connectivity representation, the cumulative distribution functions of the pore and throat sizes required for random generation of pore and throat sizes obeying the given distribution (e.g. Jivkov and Xiong 2014), as well as data not presently used in regular network models-the cumulative distribution function of the throat lengths. The latter is used for illustration of the network characterisation capabilities of the two algorithms.

Topological and Geometric Characteristics
Properties of the idealised pore network, as constructed and as interpreted by the two algorithms into PNMs, are presented in Fig. 7.
Regarding the pore coordination number distribution, in Fig. 7a, the PNMs constructed by both algorithms have peaks at coordination number two, which are not present in the original medium. This is clearly due to the introduction of significant number of pores within original throats, creating a strong bias towards pores with coordination of two. As a result, the average coordination numbers determined by the two algorithms are more than two times smaller than the "real" case. In addition, the MB model has recognised a much 1 3 smaller number of isolated pores (coordination number zero) than the WA, apparently due to interpreting isolated pores as chains of pores and throats. Figure 7b indicates relatively good agreement between the prescribed cumulative distribution of pore sizes and that for the extracted models. The WA curve is shifted to the left because this algorithm scales all physically determined pore sizes in order to match porosity. The MB curve starts higher than the prescribed one, reflecting the significant number of small pores introduced by this algorithm within prescribed throats. Figure 7c shows similarly good agreement between the prescribed cumulative distribution of throat sizes and that for the extracted models, with the exception that MB finds throats with larger cross sections because it interprets some of large pores as throats. Figure 7b, c suggests that the statistics for pore and throat sizes could be considered as sufficiently reliable for random generation of pore and throat sizes from this distribution. However, the outcome should be considered with caution as the average pore and throat sizes, depicted in the plots, do not agree well with the prescribed values. Simplified pore network models use the average pore coordination number and the average pore and throat sizes instead of the distributions presented here. In such case, the pore space interpretation by both algorithms will be inaccurate. Further, Fig. 7d shows significant differences between the throat lengths inferred by the PNMs and the prescribed ones. This is the strongest evidence of how the two algorithms (MB more than WA) interpret many throat regions as pores, increasing pore and throat numbers, reducing throat lengths, and as a consequence biasing the results for pore coordination spectrum and size distributions.  Figure 8a shows that the pore network extracted by MB identifies fewer isolated pores than the WA (possibly for the same reasons as with the idealised model) and larger number of pores coordinated by two to six throats. Both models are in agreement for pores with coordination larger than six. While such spectra can be used for future construction of regular pore networks, there is nothing to suggest that they are not biased towards smaller coordination numbers similarly to the result with the idealised model. Hence, the use of these spectra or the resulting average pore coordination number in regular pore network constructions may introduce inaccurate connectivity representation. Figure 8a, b shows that the two models are in good agreement about the cumulative distributions of pore and throat sizes. Due to the possible bias in the coordination spectrum, i.e. misinterpretation of pores and throats, the use of these distributions in pore network construction may also introduce error, albeit smaller than the connectivity inaccuracy. Figure 8d indicates that similarly to the idealised results, MB is consistently producing shorter throats due to splitting pore volumes, interpreted as throats by WA, into pore-throat-pore chains.
The data for sandstone extracted from the lower resolution image, in Fig. 9, show substantial differences between the algorithms. For example, Fig. 9a shows that MB recognises around 5% of the pores to be isolated, whereas WA reports 23%. The explanation for this is the same as in the previous two cases-MB considers as 1 3 pore-throat-chains some pore volumes, which WA considers as throats. At the same time MB network has a larger number of highly coordinated pores than the WA. The possibility here is that for largely coordinated pores, some of the coordinating throats interpreted by MB are considered as part of the pore by WA.
Apart from this, Fig. 9b-d shows that MB produces consistently more and smaller pores and throat radii as well as throat lengths. In general, the results from this lower resolution do not provide a clear choice about which spectrum and size distributions are to be used in regular network constructions. However, a comparison with Fig. 8 illustrates the strong effect of the imaging resolution on the extracted pore network statistics. The interpretation of the more complex pore space revealed with the higher resolution (Sandstone H) by the two algorithms is in better agreement than the interpretation of the lower resolution pore space. In particular, a comparison between Figs. 8a and 9a shows that the higher resolution images provide substantially richer pore space detail yielding smoother coordination number distributions and nearly two times larger average coordination numbers. This makes it difficult to consider a coordination spectrum and distributions spanning the two length scales.
The interpretations of the two algorithms of the carbonate samples are given in Appendix A.  Figures 10 and 11 show relations between pore topological and geometric characteristics in the analysed rocks, which have not been studied in detail. The existence of a relationship between coordination number and averaged pore volume was suggested by a previous study (Sok et al. 2002) for Fontainbleau sandstone imaged at a voxel resolution of 5.7 μm. Specifically, Fig. 10 illustrates the dependence of pore size on coordination number by presenting the average pore radius of all pores with given coordination number. Observed are approximately linear relations for both rocks and at both resolutions. It is important to note that more substantial deviations from linearity occur for a small fraction of very large pores with large coordination numbers. These data points are excluded from these plots. For example, the number of pores in Sandstone L, in Fig. 10a, with a coordination number between 0 and 9 is 99.4% of the total number of pores; hence, the average values calculated for CN > 9 are not reliable because there are simply not enough features contained within the volume at this resolution and field of view. Hence, the linear relations between pore size and coordination number can be assumed as valid for both materials. This provides very useful additional information for future development of regular network models. In the existing construction strategies, the connectivity and the size of any given pore are not linked.

Geometry and Topology Relationships
Finally, Fig. 11 shows relations between pore radii and average radius of the throats by which these pores are coordinated. The earliest attempt to study these relationships has been made by Lindquist et al. (Lindquist et al. 2000) with Fontainbleau sandstone. It is worth noting that the MB model is present with consistently smaller feature sizes which is due to the difference in the pore and throat definition with the WA model. Furthermore, at least 99% of all accounted pores are in the region where the relations are linear as depicted by the fitting lines. For example, the throats in the linear region for sandstone, in Fig. 11a, b, contribute 98.8% to the total number of throats. The deviation from linearity is observed only for a small fraction of pores-the largest and the smallest-the numbers of which are insufficient for deriving reliable statistics for pore-throat

Evaluation of Flow Properties
This section aims to evaluate the performance of the two pore network extraction algorithms in terms of the accuracy of their permeability predictions by comparison with a direct approach.

Direct Evaluation of Permeability
The direct evaluation of the permeability of the segmented 3D images was undertaken using the Avizo's XLab Hydro module. The computational solver in this module has been validated previously against theoretical models and standard glass bead pack models Relations between average radius of coordinating throats and radius of coordinated pore and the effect of image resolution: a WA results for sandstone; b MB results for sandstone; c WA results for carbonate; and d MB results for carbonate 1 3 (Zhang et al. 2011). The computational framework used in that module is the Absolute Permeability Experiment Simulation (APES) which utilises the finite volume method (Harlow and Welch 1965) to solve the Stokes equations for the velocity and pressure field. The Stokes equations are given by: where u is the fluid velocity (m s −1 ), μ is the fluid dynamic viscosity (Pa s) and p is the pressure (Pa).
The flow in the APES simulation is driven by the pressure gradient across two opposite boundaries of the studied geometry by introducing an accommodation zones at the inlet and outlet of the sample to ensure that there is a consistent pressure field. The module also introduces a one-voxel-wide impermeable boundary around the other faces of the geometry in order to prevent loss of fluid. The algorithm then builds the finite volume mesh directly onto the voxelised pore domain of the 3D sample and solves for the velocities and pressures in the system. The process is illustrated in Fig. 12. Once the equations are solved and the volumetric Fig. 12 Schematic illustration of the flow modelling using XLab for Sandstone H in the X direction: a voxel rendering of the pore space structure; b the flow field of a virtual fluid is displayed as streamlines and colour-coded according to the relative velocity; and c the resulting pressure field 1 3 fluxes q total (m 3 s −1 ) through the permeable boundaries are quantified, the APES utilises Darcy's law to evaluate permeability k (m 2 ): where μ is the fluid viscosity (Pa s), L is the distance between the two opposite permeable faces (m) and A is the area of the permeable face (m 2 ). Apart from the parameters described above and the segmented pore space, the module also takes the convergence error criterion (CEC) as an input, on which depends the precision of the outputs of the simulation. We conducted a sensitivity study in order to assess the influence of CEC on the permeability estimation and found that the value of 1 × 10 −7 was a good compromise between precision and computational time. A single simulation in any of the directions was found to take in the order of days, depending on the volume and complexity of the segmented pore space. In this study, the permeability of each segmented sample was evaluated in all directions (X, Y and Z) in order to capture any anisotropy at the length scale of the samples.

Evaluation of Permeability from the PNMs
Even though the outputs from the different network extraction approaches vary, the method (Valvatne et al. 2004;Jivkov et al. 2013) used to evaluate the permeability is similar in order to compare the final results. Mathematically, the output from the pore network extraction algorithms could be treated as 3D graph with the bonds being the pore-throats-pore assemblies and the nodes being the centres of the pores where bonds intersect. Mass transport is allowed to take place along the bonds, with the flow obeying the mass conservation at every node i: where j represents all bonds connected to node i. It should be noted that this is strictly true only for incompressible flow and that the pressure changes in the nodes due to the flow are negligible. The flow rate q b,i−j passing through a bond that connects node i and node j is described as follows: where c b,i−j is the bond conductance, L i−j is the geometric distance between node i and node j and P i/j are the pressures at nodes i/j, respectively. This equation assumes that flow is laminar and incompressible and that the radius of the bond is considerably smaller than the bond length. These assumptions are consistent with the presumptions for mass conservation at the nodes.
Since our model takes into account the contribution of pores to total flow as well as the throats, we describe the conductance of a bond (pore-throat-pore assembly) as a simplification of three cylindrical throats in series using the harmonic mean of the individual conductance: ∑ j q n,i−j = 0, where t indicates the connecting throat. In addition, the hydraulic conductance c is given by Hagen-Poiseuille equation: where k is ½ for a circular throat and G is the dimensionless shape factor which is the ratio of the throat area to its perimeter squared. Combining Eqs. (3) and (4) results in a system of linear equations, with the unknown being the pressures at the nodes. The last step, before solving the system, is to establish a criterion for the boundary nodes. We define the boundary nodes by first adding the pore radii to the respective coordinates of the pore. If the resulting coordinates lay outside of the predefined sample boundaries, they are marked as boundary nodes. Any pore that obeys this condition would be open to the inflow/outflow in the system. The pressures at the nodes residing in the permeable boundaries in the geometry are fixed. By solving the system, the pressures and consequently the flow rates in the whole geometry are obtained. The flow rates through the bonds that are connecting a boundary node to a node in the interior of the system are summed up to calculate the total flux coming in and out of the system. Note that mass conservation yields Darcy's law (Eq. (2)) is then utilised to find the permeability of the network. This procedure has been implemented in MATLAB R2014a with typical computational times for the materials used in this paper being in the order of seconds.
The methodology described above used to calculate the sample permeability is the same for both the WA and the MB algorithms, apart from the fashion in which we calculate the conductance of the individual bonds. This difference arises from the nature of the individual approaches. The WA algorithm, for example, reports the pore and throat radii, coordinates for the pores and the connectivity of the individual pores. The conductance for a bond c b,i−j is again calculated using the harmonic mean of the individual contributors: where R p,i/j is the radius of the pores, L t,i−j and R t,i−j are the length and the radius of the throat and CC is a global coefficient which predetermines the contribution of the pore radius to the area of the cylinder. The calibration for the value of CC is given in Sect. 4.4.3. We find this value by calculating the values of permeability in all principal directions for the idealised pore network (not the voxelised version of it) with values of CC ranging from 0 to 1 with an increment of 0.0001. Then, we take the value of CC for the lowest value of the second norm between the estimates of permeability reported from XLab (P XLab ) and MATLAB (P Mat ). Note that P Xlab = [PermX; PermY; PermZ] and P Mat = [PermX(CC); PermY(CC); PermZ(CC)], and we calculate the second norm as: where k = 1, 2, 3 indicates the direction in which the geometry is evaluated.
q total,inlet = q total, outlet . (8) In contrast to the WA, the MB algorithm has more detailed outputs. The most significant difference is the fact that MB reports shape factors for the individual pores (G p,i/j ) and throats (G t,i-j ) as well as pore lengths (L p,i/j ). We calculate the bond conductance using: The shape factor for a circular cross section is 1/(4π), which is what is used in calculating the conductance of the bonds in the WA model. With the MB, these shape factors are calculated as an average over the length of the reported feature and they reduce as irregularity increases (Valvatne et al. 2004). More information on the definition of the shape factors and pore lengths can be found in Dong et al. (2009).
It should be noted that the approach used to calculate the conductivities with both models and the inherent simplifications does not account for the effects of the flow on the pressure when passing through throats with different cross sections. However, we assume that these influences are negligible at low flow velocities.

Permeability Results
As described in the previous section, prior to calculating permeability values with the extracted irregular pore networks by WA, a calibration of the pore contribution parameter CC in Eq. (8) has been performed using the idealised network. Figure 13 shows the calculated permeability in all directions of this network with changing CC from 0 to 1, as well as the permeability values calculated by XLab on the voxelised representation as described in Sect. 2.1. The value that provides the lowest second norm between the two estimates, in Eq. (9), is CC = 0.63. This value is used for calculating the permeability of all WA extracted networks. It is noteworthy that this value is close to the value used by Dong and Blunt (2009). Table 4 summarises the results of all permeability calculations. For the two rocks, the table also provides the ranges of experimentally determined permeability values. The (10)   and their connectivity. The pore networks of Sandstone L and Carbonate L extracted from MB algorithm tend to predict values smaller than the respective prediction from XLab, whereas this trend is not present with the WA pore networks. This trend is also observed for the other low-resolution sample: Carbonate L. It is interesting to point out that XLab and WA report zero permeability for Sandstone 2 in X direction, whereas MB calculates it as 8.64 × 10 −14 m 2 . The segmented geometry does not have a percolation path in the X direction, which indicates that the MB has introduced links between pores that are not present. If we exclude this sample, the average difference between MB and XLab is 33% and that between WA and XLab is 40%.
In terms of overall performance, both methods provide reasonable predictions, despite the differences in the interpretation of the same pore space. This outcome is expected since the main observation for both algorithms is that they tend to interpret pore volumes, which are potentially throats, by systems of pores and throats, MB more so than WA. While this introduces biases in the pore coordination spectra and the pore and throat size distributions, it does not affect the simple laminar non-reactive flow through the constructed pore networks. However, the potentially incorrect interpretation of some pore volumes may result in inaccurate calculations of more complex reactive flows in the constructed irregular pore networks.

Discussion
In this study, we have compared the performance of two of the leading algorithms for pore network extraction. The key findings are that MB tends to identify a larger number of features (both pores and throats) than WA. In most cases, the pore sizes from MB are much smaller than the ones from WA. Further, MB throats tend to be smaller and shorter than the ones reported by WA. These differences arise from the manner in which both algorithms produce the equivalent networks. Taking the pores reconstruction for example, MB inscribes a sphere with maximum radius inside a cavity, whereas the WA measures the volume of the cavity and reports a sphere with the equivalent radius to match the volume. This explains the discrepancies between the average pore and throat radii as reported from the two algorithms.
The comparison between the extracted statistics and the known statistics for the idealised model shows interesting trends in both models. The spike at CN = 2 and the generally shorter throats indicate that both models interpret the long throats as a series of pores joined by shorter throats. Over-segmentation is a ubiquitous challenge for the watershed approach (Thompson et al. 2006) and strongly suggests the need for pore merging. Regarding isolated pores (CN = 0), the MB and WA report 1 and 169 such pores, respectively, whereas the idealised geometry actually contains 97 such pores. This implies that MB tends to introduce links that are not present in the original segmented pore space. This observation is also supported by the permeability results for Sandstone L in the X direction, which has no percolation path between the permeable boundaries. MB reproduces more accurately the cumulative distribution of pore sizes of the idealised model, but tends to extract smaller features for the rest of the material in comparison with the WA. This could be due to the fact the WA employs additional filtering that removes any small features that might not contribute to the permeability of the network but influence the model construction. In general, both models overestimate the number of pores and throats and underestimate throat lengths and average coordination numbers. The medial axis algorithm (Lindquist et al. 1996) is reported to generate topologically equivalent skeletons of the porous media.

3
However, it suffers from unambiguously identifying pore locations, which can also result in skewing the geometrical and topological statistics of the extracted networks, as observed with the MB and the WA. Clearly, there is a need to recalibrate the coordination spectrum and the distributions of pore and throat sizes if the results were to be considered reliable for subsequent physically realistic regular pore network construction. This could be achieved by a series of synthetic models and is a subject of ongoing work to be reported in future.
The permeability of real materials ranges across many orders of magnitude. For example, experimental measurements of Hollington sandstone's permeability report values in the range 3.0 × 10 −14 m 2 to 3.3 × 10 −12 m 2 (Makse et al. 1996) which is in agreement with the pore network predictions from two models. Similarly, the values for carbonate rock reported range from 7.0 × 10 −15 m 2 to 7.0 × 10 −13 m 2 (Ehrenberg and Nadeau 2005), which is also in agreement with the values calculated in this work. In terms of porosities of the equivalent networks, MB tends to produce networks with significantly lower porosities than the segmented pore space. The reason it is still predicting the permeability of the sample is the use of shape factors, calculated from the voxelised representation, which are not accounted for in the porosity estimation. WA, on the other hand, also underestimates porosity but not as much as MB. WA does not output any shape factors for the reported features, and hence the permeability is calculated based on the feature sizes.
The fact that both models report different topological and geometrical statistics but still manage to predict permeability with reasonable accuracy shows that this property is fairly insensitive to over-segmentation of throats. It can be concluded that the study of simple flow through irregular network constructed directly on images is not affected by the decision of what constitutes a pore and a throat, as long as the connectivity is close to the real pore space. However, the statistical information extracted by the algorithms should be used with caution for constructing synthetic pore network, clearly for the materials studied in this work. Further, even for the image-based networks, the over-segmentation may lead to erroneous simulations for more complex transport processes, for example reactive transport. For example, modelling a sorption process using the MB network might yield faster blocking (unless the individual volumes of the features are tracked separately) due to the generally smaller porosity than the WA. Depending on the length scales of the model, simulating bacterial proliferation or colloidal transport in porous media using both networks could have very different results. A bacteria cell with a size of a couple of microns might not be able to pass through small pores and throats, resulting in parts of the network being physically inaccessible. Furthermore, an experimental study reports bacteria (Delftia acidovorans) preferentially establishing in the pores of the synthetic media used (Zhang et al. 2010). Modelling the process using the extracted networks may result in unrealistic results, given that the models report elevated number of pores.
When interpreting the relationships between coordination number and the average pore radius and pore radii and average throat radii, one should take into account the fact that at different resolutions the tomography scans capture different numbers of features of a certain size. The number of identifiable pores and throats drops off as their sizes approach the resolution limit, as shown by Lowe et al. (2015). This is due to the resolution-limiting factor and the noise filtering applied to the samples and does not necessarily reflect reality. The significance of these relationships for the carbonate and the sandstone at the resolutions imaged is that: the more coordinated a pore is, the larger on average it is; the larger a pore is, the wider on average are the throats coordinating this pore. Whether such relationships hold for other materials or for the same materials at different observation windows is a question for further investigation. Depending on the studied problems, these relations might have to be considered, especially if the pore network model is not based on actual 1 3 tomography scans but on extracted statistics. Building a network that does not obey them might result in different mixing and dispersion coefficients from the ones observed in the real sample. These relationships have been shown to be important for predicting the relative permeability in regular pore networks (Arns et al. 2004). The significance of these relationships for predicting single-phase permeability remains to be investigated.
It is interesting to point out that the methodology used for the construction of the idealised model, namely using Voronoi tessellation and scaling the pore and throat sizes to the sizes of the resulting cells and faces, generated a medium that possesses such linear relations between pore sizes, coordination numbers and throat sizes. This suggests that a combination of dynamically constructed Voronoi tessellations with one of the two algorithms could provide an improved approach for pore network extraction.
In terms of computational speed and resources for evaluating the permeability of segmented pore space, this study showed significant benefits of using PNM over XLab. The computational time needed for the extraction and the calculation of permeability of the equivalent network is in the order of minutes, whereas using XLab computation takes days. Similar observations are made by Varloteaux et al. ) but for modelling reactive transport. Comparing the performance of the WA and the MB in terms of permeability computational time, the WA is faster than the MB due to the fact that it outputs less pores and throats, which ultimately results in less variables that need to be evaluated.

Conclusions
In conclusion, our study outlined the importance of understanding what characterisation approach is used to describe the pore space of a medium. Different approaches interpret the same pore space differently and in addition use alternative methods for calculating permeability. Both methods tend to over-segment longer throats and represent them as pore-throat-pore assemblies. The elevated number of pores could yield unrealistic results, if the modelled process is particularly relevant to the pores, i.e., precipitation or biofilm growth. Additionally, our results suggest that in some cases the MB may introduce throats, whereas the WA could remove existing small throats that are present in the segmented pore space. This effect should be considered in transport modelling directly onto the extracted networks, as it may result in including/excluding certain parts of the network, which would affect the global concentration of solutes. Further, the topological and geometric statistical characteristics of the resulting pore networks depend on the window of observation. The results presented do not suggest repeatability of these characteristics at different length scales; i.e. based on the results, one cannot claim a fractal nature of the pore space for extending potential models beyond these windows of observation. Despite these differences, the networks generated by both algorithms predict permeability values for single-phase flow with acceptable accuracy (33% and 40% on average, for MB and WA, respectively). Additionally, the permeability estimates from the two models are in better agreement for the fineresolution samples, as opposed to the coarse-resolution data sets. The computational effort required for evaluating flow properties of the extracted networks is significantly less than direct fluid mechanics calculations. This implies the use of pore networks will continue to play important role in understanding structure-properties relations in porous media; however, improved approaches for pore space interpretation need to be developed. Our study does not claim that pore network modelling is not a reliable method for analysis of transport, but that the two algorithms considered are not sufficiently advanced to interpret accurately the pore spaces studied here. It is possible that the algorithms provide better results for other 1 3 porous structures with simpler pore space, but we emphasise the necessity for developing improved algorithms, which is a subject of ongoing work. Finally, the linear relationships between pore size, pore coordination number and radii of coordinating throats should be used in future developments of regular pore networks to improve their physical realism.

3
Possibly due to the smaller difference between the two resolutions, two times for carbonate rock, the difference between the results at these resolutions is much smaller than for sandstone. Regarding the coordination spectra, in Figs. 15a and 16a, both algorithms peak at the same numbers, 2 and 1, respectively, with WA interpreting a much larger fraction of pores as isolated than MB. This is consistent with the previous observation of the tendency of MB to treat isolated pores as systems of pores and throats. Nevertheless, the calculated average coordination numbers by the two algorithms are in a relatively good agreement. The tendency of MB to interpret many pore volumes as chains of pores and throats, where WA interprets these as throats, leads to the results shown in Figs. 15a-c and 16a-c. At both resolutions, MB derives larger numbers of smaller pores and throats, as well as shorter throats, than WA. The results at the two resolutions are not very dissimilar, suggesting that the carbonate has a less complex pores system within this window of observation than the sandstone. Unlike the sandstone, the results for carbonate could be reconciled to span this window of observation in terms of coordination and size distributions.