A Novel Regional‐Minima Image Segmentation Method for Fluid Transport Simulations in Unresolved Rock Images

Unresolved digital rock images are often used to avoid high computational costs and limited field of views associated with processing fine‐resolution rock images. However, segmentation of unresolved images using classical methods is suboptimal due to the presence of the partial‐volume effect. Suboptimal segmentations can significantly influence the geometry and effective properties of the reconstructed models. This study reveals that partial‐volume pixels with high pore fractions remain as regional minima in intensity levels in unresolved images. By identifying these regional‐minima pixels, we can effectively extract pore space obscured by the partial‐volume effect. Based on this observation, we propose a novel segmentation method capable of identifying these regional‐minima partial‐volume pixels and converting them to pure pore pixels, thereby binarizing the digital rock images. The method is validated on sandstone and carbonate rock samples. Our method demonstrates a notable improvement in modeled permeability accuracy, surpassing 50% compared to the thresholding method and over 30% compared to the watershed method. Moreover, models segmented by this approach exhibit smaller pore and throat sizes compared to the substantially overestimated results obtained by classical methods. These findings suggest that the regional‐minima segmentation method effectively corrects for the partial‐volume effect and preserves more detailed pore structures. Consequently, it enhances the quality of binarized rock geometries, leading to improved accuracy in fluid‐flow simulations.


Introduction
In recent decades, Digital Rock Physics (DRP) has emerged as a non-destructive approach for investigating the flow and transport characteristics within porous geometries, such as rocks and soils (Andrä et al., 2013a;Farhana Faisal et al., 2019;Goldfarb et al., 2022;R. Li et al., 2022;Wu et al., 2019).In DRP analysis, the quality of digital rock images is of vital importance as it directly influences the characterization of pore structures and, consequently, the accuracy of reconstructed rock models during the image segmentation process.High-resolution images are favorable for resolving fine rock features, such as micro pores and grain interfaces, leading to more accurate rock models.However, utilizing high-resolution images also comes with the drawback of a larger model size, which demands considerable computational time for numerical simulations.Typically, the size of a scanned image does not exceed 4,000 pixels in each direction (Bazaikin et al., 2017).Consequently, higher image resolution results in a narrower field of view, potentially introducing uncertainties in terms of sample representativeness.Alternatively, unresolved images with lower image resolutions are often employed in DRP analysis.But a coarser resolution fails to capture sub-resolution details that finer resolutions can resolve.Moreover, the presence of image artifacts, such as image blur and noise, can cause the rounding and smoothing of pore space (Saxena et al., 2017;Schlüter et al., 2014).These issues, when combined, lead to a partial-volume effect that can significantly impact the precision of reconstructed rock geometries and fluid-flow simulations.
The partial-volume effect typically arises due to low image resolution that is insufficient for resolving delicate pore structures.This effect leads to partial-volume pixels that encompass both solid and pore phases.The primary challenge posed by the partial-volume effect is the blurring of phase boundaries (Schlüter et al., 2014).In essence, these boundaries do not exhibit distinct, abrupt intensity changes but rather manifest as gradual transitions spanning across multiple pixels.The partial-volume effect has been recognized as the major source of uncertainty hindering the success of image segmentation.Even a small fraction of the partial-volume phase can have a considerable impact on flow distribution and the modeled permeability (Gerke et al., 2015).Some computational techniques have been developed to mitigate image artifacts and enhance image quality, and these methods have shown significant success (Andrä et al., 2013b;Buades et al., 2011;Culligan et al., 2006;Kaestner et al., 2008;Schlüter et al., 2014;Sheppard et al., 2004).However, these tools must be used with great care as they unavoidably lead to the loss of some valuable rock structural information from the raw images (Schlüter et al., 2014).Therefore these tools can only be employed to partially mitigate the presence of artifacts.Furthermore, while some artifacts like noise, ring artifacts, and intensity variations can be effectively removed using these tools, addressing the partial-volume effect caused by image blur is a more challenging task (Schlüter et al., 2014).However, reconstructing digital rock models for fluid-flow simulations requires the images to be segmented into pore and solid phases.Thus the partial-volume phase has to be converted to pore and solid phases using a suitable segmentation technique.This requires an accurate understanding of the pore-solid spatial distribution within the partial-volume phase.
Commonly used segmentation methods are categorized into global thresholding and local-adaptive methods.The global thresholding imposes an intensity threshold to divide the image into the pore and solid phases, which considers no presence of the partial-volume phase.Many studies have reported that the thresholding method yields suboptimal results when applied to unresolved images in terms of fluid-flow simulations (Al-Kharusi & Blunt, 2007;Bultreys et al., 2016;Korneev et al., 2018;Scheibe et al., 2015).On the contrary, local-adaptive methods, such as the hysteresis method, watershed-based method, and converging active contour method, have been reported to produce more accurate segmentation results by taking into account local statistics (Iassonov et al., 2009;Schlüter et al., 2014).One of the most widely used local methods is the watershed segmentation method (Misbahuddin, 2020;Vincent & Soille, 1991).This method initiates by flooding initial markers (seeds) that designate distinct regions within the images.Watershed lines are then identified at places where two different regions converge, which will act as boundaries between different regions.The process continues until the entire images are successfully segmented.Though the watershed method tends to yield more accurate segmentation results compared to the thresholding method, it does come with its own set of limitations.Notably, the pore space extracted using this method often tends to be overestimated.Moreover, the watershed method necessitates the placement of seeds and precise edge detection.This introduces a notable user bias and uncertainties in the segmentation outcomes.
Recently, machine-learning or deep-learning-based segmentation methods have gained significant popularity (Ar Rushood et al., 2020;Badrinarayanan et al., 2017;Ronneberger et al., 2015;Sun et al., 2023;Y. D. Wang et al., 2021).With substantial training, this technique demonstrates effectiveness and often yields superior segmentation results compared to traditional segmentation methods.However, a major challenge associated with this technique is the absence of ground truth segmentation results for different digital rock images (X.Li et al., 2023;Manzoor et al., 2023;Reinhardt et al., 2022).The generation of ground truth segmentation results typically relies on conventional segmentation methods.However, generating ground truth results, especially in unresolved rock images with severe partial-volume effect, is a process susceptible to strong user bias.Furthermore, to the best of our knowledge, there is currently no absolute way that can definitely achieve ground truth for unresolved rock images.Hence, the use of deep-learning techniques still relies on the supplementation of conventional segmentation methods.Additionally, achieving effective training for deep-learning tools underscores the importance of identifying an optimal conventional segmentation approach, particularly in the context of unresolved rock images.
In this study, we introduce a novel image segmentation method designed to enhance the reconstruction of rock geometries from unresolved digital rock images for fluid-flow simulations in diverse rock types (Section 2).This method starts from a Gaussian distribution process, which can effectively discriminate between different phases in various rock types while minimizing subjectivity (Section 2.2).It then identifies partial-volume pixels with regional-minimal intensities and subsequently converts them into pure pore and solid pixels through an optimization process based on the principles of maximum-entropy theory (Sections 2.3 and 2.4).In contrast to classical segmentation methods, which often face challenges related to strong user bias and uncertainties, our method minimizes user involvement by offering only a few options for users to determine beforehand.To validate this method, 9 digital rock models of different rock types are reconstructed from unresolved rock images (Section 2.1).We compare the extracted pore-and throat-size distribution using the maximal-ball algorithm and the permeabilities computed by the Lattice-Boltzmann method, using rock models reconstructed by our method and the conventional thresholding and watershed methods (Section 3).The results demonstrate that our method improves the accuracy of fluid-flow simulations in both sandstone and carbonate rock samples, even when utilizing unresolved images at resolutions several times lower than the high-resolution reference.In Section 4.1, we delve into the impacts of the partial-volume effect in unresolved images and reveal that partial-volume pixels with regional-minimal intensities are more likely to represent pore space.Meanwhile, uncertainties and limitations of the proposed method are discussed in Section 4.2.Section 5 provides a summary of the main findings of this study.

Digital Rock Models
The digital sandstone models utilized in this study comprise four paired Bentheimer sandstone models (BS1 ∼BS4) with resolutions of 2, 6, and 18 μm, respectively.Additionally, the Doddington sandstone model (DS) is included at a resolution of 5.3838 μm, and a Mount Simon Sandstone (MS) at a resolution of 5.22 μm.For the carbonate rock models, the Estaillades carbonate (EC) model is employed at a resolution of 3.1 μm (downsampled to 12.4 μm using a bilinear filter in Fiji (Schindelin et al., 2012)).The Indiana limestone (IL) and Middle Eastern carbonate (MEC) models are utilized at resolutions of 2.68 and 10.72 μm, respectively.Specifically, for the IL and MEC models with the resolution of 2.68 μm, we extract a 10,00 3 cubic sub-volume from their original full-core image data sets.Conversely, we use the full-core images for the IL and MEC models at the resolution of 10.72 μm.Table 1 provides a summary of these digital rock models along with their properties.Figure A1 displays the 3D views of these digital rock models.For simplicity, all rock models are treated as pure rock models, and the presence of any additional minerals within them is neglected, except for the MS model.The MS model includes a high proportion of pore-filling clay and feldspar minerals, making it an ideal model for testing the application of our segmentation method in complex porous mediums.
For the Bentheimer sandstone models, an automatic Otsu's thresholding method is applied to segment the 2 μm rock images.Given that the BS images at different resolutions are paired images, the porosities, extracted pore networks, and modeled permeabilities of the 2 μm models serve as the ground truth to validate our segmentation method when applied to the 6 and 18 μm models.Similarly, the segmented IL and MEC models at the resolution of 2.68 μm serve as ground truth for validating the segmentation methods applied to the paired 10.72 μm models.The segmented ground truth images for IL and MEC models are obtained directly from Alqahtani et al. (2021).As the IL and MEC models are full-core, a 10,00 3 sub-volume is extracted from the original data set.Subsequently, a paired sub-volume is extracted from the lower-resolution models after segmenting the rock images using different segmentation methods.We compare extracted pore networks from the sub-volumes at different resolutions to validate the quality of our segmentation methods.However, for permeability comparison, we assess the modeled permeability of the low-resolution full-core models against their laboratory measurements.For DS, MS, and EC models, given the absence of high-resolution paired images to generate ground truth, we validate the quality of the segmentation results by comparing their modeled permeability to the experimental measurements.
For comparisons, all the digital rock models are segmented by different methods, including the thresholding, watershed, and the regional-minima segmentation method proposed in this work.In the thresholding method, an intensity threshold value is selected.Pixels with intensities below it are considered pore while those above the threshold are deemed solid.The threshold value is initially determined using Otsu's method.However, if Otsu's selection leads to a porosity higher than the ground truth or experimentally measured porosity, a lower threshold is manually set.Moreover, for the carbonate rock models, a multiphase Otsu's segmentation method is employed to categorize the rock images into macro pore space, micro pore space, and pure solid phase.In the watershed segmentation method, threshold values are selected to distinguish between pore and solid phases, as described by Alqahtani et al. (2021).The chosen threshold values are conservative to ensure that only certain pore and solid phases are extracted.These extracted pixels serve as seeds for the watershed transformation.By using the image gradient with the Canny method (Canny, 1986), watershed lines are detected, enabling the watershed algorithm to segment the rock images into pore and solid phases.Furthermore, the watershed segmentation method is versatile and can be applied to multiphase systems.In our case, it will be used to segment the carbonate rock models into macro pore, micro pore, and pure solid phases.It is important to note that the thresholding segmentation is executed in Fiji.The watershed segmentation is carried out in Dragonfly software, Version 2022.2 for [Windows].Comet Technologies Canada Inc., Montreal, Canada; software available at https://www.theobjects.com/dragonfly.In addition, the segmented Doddington sandstone model using the watershed method is retrieved from Andrew et al. (2014).

Regional-Minima Image Segmentation Method
In this section, we introduce a novel segmentation workflow that incorporates regional-minima features in rock images.The method begins by distinguishing different phases within the rock images through a Gaussian distribution process.Following this, the method extracts regional-minimal intensity pixels from search windows constructed around each pixel.Ultimately, an optimization process, based on the principles of maximum entropy theory, is employed to determine the optimal segmentation parameters for this method.Further details regarding these steps and the validation of this method are discussed in the subsequent sections.

Phase Separation
The first process of the regional-minima segmentation method is to distinguish among the pure pore phase, partial-volume phase, and pure solid phase through a Gaussian distribution method.Consider a bi-modal digital rock model composed of the pore and another pristine solid phase ∏ i , i = p, s.Ideally, the pure pore and pure solid intensity histograms can be assumed as univariate normal (Goldfarb et al., 2022).The pore and solid histograms are represented by: where f i (x) denotes the probability of pixels having an intensity x; a is the peak height of the Gaussian distribution; μ is the mean intensity; σ is the standard deviation; i = p and i = s represent the pore and solid phase.A sufficiently high image resolution can accurately depict the bi-modal nature of the rock sample.When the image resolution decreases, the partial-volume effect induced by image blur can lead to skewed distributions for the two phases (Schlüter et al., 2014).Since the partial-volume effect is inherently linked to the image resolution, it is anticipated that at a moderately low image resolution, some pixels may still conform to Gaussian distributions.We apply a Gaussian fit algorithm embedded in Matlab to find the two Gaussian distributions for the existing pure pore and solid phase separately.Some typical examples of the Gaussian fit process are illustrated in Figure 1.The figure shows the image intensity histograms along with the Gaussian fits for pore and solid phases for BS1_6 μm, BS1_18 μm, IL, and EC models.It is important to note that in BS1_18 μm, only the Gaussian fit for the solid phase can be identified because the image resolution is not fine enough to exhibit the bi-modal nature of the Bentheimer sandstone.Additionally, in EC, where the micro pore space is a dominant phase situated between the pure pore and solid phases, we employ a Gaussian fit for the micro pore space instead of the pure solid phase.
By comparing the Gaussian distributions and the image intensity histograms, we note the following points: 1.I p : the intersection between the pore Gaussian fit and the image histogram.2. I m : (μ p + 3σ p ), the three standard deviation greater than the mean of the pore phase Gaussian fit, accounting for 99.73% of the entire distribution set.3. I s : the intersection between the solid Gaussian fit and the image histogram.
The selection of three intensity values is common across different rock types.However, for carbonate rock models, an additional value I n is determined as the average of I p and I s .In cases like EC, where the mixture of Water Resources Research 10.1029/2023WR036855 micro pores and solid is prevalent, I n is specifically set to be the intersection between the Gaussian fit for micropores and the image histogram.
With these threshold values, the image histogram is partitioned into distinct regions, as summarized in Table 2.In both rock types, I p is the point where the actual image histogram deviates from the pore phase Gaussian distribution.This value is considered the highest intensity that the pure pore pixels possess.It defines the boundary between the pure pore phase and the macro-mixed phase.Likewise, I s signifies the boundary between the pure solid phase and the mixed phase.It is the point where the image histogram deviates from the Gaussian distribution of the pure solid phase.Ideally, the pure pore and pure solid phases follow the Gaussian distributions.As the image resolution coarsens, pixels situated between the macro pore center and solid grain center are the first to be partialvolume pixels due to image blur.While pixels residing within macro pore center and grain center are less influenced by the partial-volume effect.This is because these central pixels are not in direct contact with their counterpart pixels.Therefore, the central pixels continue to adhere to Gaussian distributions, resulting in the presence of the pure pore regions (S1 and C1) and the pure solid regions (S3 and C4).Pixels in the partial-volume phase encompass both pore and solid phases, resulting in the averaging of their intensities.Thus the partial-volume phase is positioned between the pure pore and solid regions.The separation between the partial-volume phase and the pure phases can therefore be captured by the deviations observed between the actual histogram and the Gaussian distributions.In certain cases, such as models like EC, both the solid phase and micro-mixed phase can be distinctly fitted by Gaussian distributions.I s is set as the lowest intensity value between the two phases, rather than relying on the deviation between the solid Gaussian fit and the image histogram.Naturally, given the coarse image resolution, a mixed phase between the micro-mixed phase and pure solid phase exists.However, as the pore fraction within the micro-mixed phase is minimal, and no pore fraction is present in the pure solid phase, the small pore fraction in this mixed phase can be neglected.Hence, for simplicity, the lowest intensity value between the micro-mixed and pure solid phases is set to I s .
The presence of image artifacts, such as noises, can impact the shape and characteristics of the image histogram.In particular, the pore phase is more susceptible to the influence of image artifacts, due to the relatively lower number of pore pixels compared to solid pixels.In the pore phase, image noises manifest as scattered dots that appear slightly brighter than the surrounding pure pore pixels.This presence of noises causes an early departure from the pore phase Gaussian distribution in the histogram.To solve this problem, we introduce another threshold value, I m which is the three standard deviation greater than the mean of the pore phase Gaussian fit (μ p + 3σ p ), accounting for 99.73% of the entire distribution set.We allow the threshold that separates the pure pore and macro mixed phase to be adjustable within the range of I p ∼ I m .An optimization process is later applied to determine the final I p value.This active region (S1′ and C1′) is therefore called the potential pure pore phase.However, in models where pore phase Gaussian distribution cannot be identified, manual adjustment of I p and the upper boundary I m becomes necessary.Typically, it is assumed that the upper boundary I m will not exceed (μ s 3σ s ) because, beyond this point, the pure solid phase is present.Thus for models where manually setting I m is challenging, such as BS1_18 μm ∼ BS4_18 μm in this work, I p is set manually in a conservative manner and I m is fixed as (μ s 3σ s ), as depicted in Figure 1e.
Between I m and I s is the partial-volume phase.In this region, pixels contain both pore and solid phases.Depending on the source of the pore phase within the partial-volume pixels, this mixed region is subdivided into the macro-mixed phase (S2 and C2) and micro-mixed phase (C3).In the macro-mixed phase, the pore phase fraction contained in each pixel is part of macro pores, therefore these macro-mixed pixels are located at the boundary of macro pores, forming a transition region between the solid and macro pores (see Figures 1c, 1f, 1i, and 1l).While micro-mixed pixels encompass entire pores since their pixel size is larger than the pore size.In sandstone models where macro pores are dominant, we assume that only the macromixed phase is present.While in carbonate rock models where micro pores are abundant, both the macromixed and micro-mixed phases are present.A threshold value, I n , is adopted to separate the macro-and micro-mixed phases in carbonate rock models.In EC where the Gaussian distribution of the micro-mixed phase can be clearly detected, I n is set to be the deviation between the image histogram and the Gaussian fit (see Figure 1k).In this case, the "pure" micro-mixed phase is separated from the macro-mixed phase and the mixed transition of the two phases.In IL and MEC where the Gaussian distribution for the micro-mixed phase is not identifiable, an estimation value of I n is adopted, which is the mean of I p and I s (see Figure 1h).

Regional-Minimal Pixels Identification
After separating the images by the Gaussian distribution process, we develop an algorithm to systematically identify regional intensity minima within the model.Specifically, for a given pixel in the model, e(i, j, k), a square search window of dimensions 3 × 3 is constructed for the 2D scheme (RM2D), or a cubic search window of dimensions 3 × 3 × 3 for the 3D scheme (RM3D), as illustrated in Figure 2a.In RM2D, the range of pixel indices within this search window is defined as (i 1: i + 1, j 1: j + 1, k: k).In RM3D, the pixel indices range is (i 1: i + 1, j 1: j + 1, k 1: k + 1).Within the search window, we proceed to sort all intensity values and extract the n lowest values, where n can vary from 1 to 8 in RM2D and 1 to 26 in RM3D, as shown in Figure 2b.These extracted values are initially considered temporary regional minima, while the remaining pixels are regarded as maxima.It is important to note that any pixels with intensities below the threshold value I p are primarily classified as regional minima and are thus excluded from further consideration when identifying regional minima.Similarly, pixels with intensities above I s are treated as regional maxima and excluded.Furthermore, the n regional minima are unique intensity values, meaning that identical pixels with the same intensity are counted as a single regional minimum.Therefore, the total number of regional-minimum pixels within a search window may be less or more than n.The construction of the search window means that each pixel is shared by 8 search windows in RM2D and 26 search windows in RM3D in addition to the one centered on the pixel itself, as depicted in Figure 2c.A pixel is recognised as a regional minimum only if it retains this characteristic in all the overlapped search windows.In other words, a pixel is not a regional minimum if it is identified as a maximum in any of the overlapped search windows.This process is repeated across all pixels in the data set.As shown in Figure 2d, the original image histogram of BS1_6 μm is then converted into two distributions: regional minima and regional maxima.

Segmentation Parameters Optimization
Following the extraction of regional-minima distributions, we will treat regional-minima pixels as pure pore pixels.To match ground truth porosity, we develop an optimization algorithm to optimize the I p and the upper boundary of the regional-minima distribution.In this optimization process, we introduce two parameters, Tp min and Tp max .Tp min is the adjustable value of I p within the range of I p ∼ I m , and Tp max is the highest intensity of the regional minima.Specifically, pixels categorized as regional minima with intensities below Tp max are designated as pure pore pixels, while all other pixels in the model are classified as solid pixels.
Following the acquisition of the regional minima distributions (n1 ∼ n8 in RM2D and n1 ∼ n26 in RM3D), the first step involves adjusting the values of Tp min and Tp max to preserve the ground truth porosity derived from highresolution models or experimental measurements.Given that Tp min and Tp max are both adjustable parameters, multiple Tp min Tp max pairs can be identified that correspond to the reference porosity.Furthermore, it is crucial The intensity values of pixels in the search window are sorted, and the lowest n values are selected to be potential regional minima, while the rest are treated as regional maxima.(c): For each pixel e, it is shared by 8 search windows in 2D and 26 search windows in 3D centered on its neighboring pixels and one search window centered on itself.A pixel is considered a regional minimum only if it is identified as a regional minimum in all the search windows.(d): An example of extracted regional minima and maxima distributions for the BS1_6 μm model.to determine the optimal regional-minima distribution for image segmentation, which is the value of n.In this work, the information theory is invoked to identify the Tp min Tp max pair and n value that optimize the pore phase distribution.Prior research has developed a maximum entropy function using information theory to determine a threshold that separates original intensity images into pore and solid phases (Kapur et al., 1985;Sahoo et al., 1988).In this step, our objective is to optimize the pore phase distribution through the information theory.
To achieve this, we treat the pore phase as a distinct system and compute its entropy value through: where H p (n) is the computed pore phase entropy value using the regional minima distribution with a factor of n; I min is the minimum intensity values within the model; P p (I ) denotes the probability of pore pixels having an intensity value of I in the pore phase; P p (I min : Tp max ) is the cumulative probability of pore pixels, which is equivalent to the porosity.
In the optimization process, we iteratively increase the value of Tp min within the range of I p ∼ I m .For each Tp min value, we can determine an appropriate Tp max that results in a porosity that matches the reference porosity.Thus this step ends up with variable Tp min Tp max pairs.For each pair of Tp min and Tp max , we calculate the entropy through Equation 2. As a result, for each n value, we identify all potential Tp min -Tp max pairs that achieve the reference porosity.Among these pairs, we choose the one with the highest entropy value as the optimal separation scheme for the given n value.Ultimately, for different n values, we select the one that maximizes the entropy value to be the optimized regional-minima distribution.A typical example (BS1_6 μm) of the optimized n value is illustrated in Figure 3. Figure B1 includes the pore phase entropies for all models used in this work.A common trend in entropy values reveals that the pore phase entropy initially increases and then decreases with the increase of the n value.Therefore, it is possible to identify a specific n value that attains the peak entropy value while maintaining a matched porosity, representing the optimal segmentation scheme.However, in some cases where  It is important to mention that in the proposed method, regional-minima pixels are treated as pure pore pixels.However, in reality, pixels in the partial-volume phases are not pure pore pixels but contain both pore and solid phases.Thus, treating them as pure pore pixels may lead to over-segmentation.This is particularly evident in carbonate rock models, where individual micro pores are dispersed within the solid framework.Oversegmentation not only enlarges the pore size but also bridges pores.Considering that the contribution of micro pore space in carbonate rocks to flow is several orders of magnitude less than macro pore space (Soulaine et al., 2016), we only treat the regional-minima pixels below I n as pure pore pixels in carbonate models, while those above I n are regarded as non-permeable micro pore space.

Pore Network Extraction
In this section, we employ a maximal-ball algorithm to extract a network model that preserves topological equivalence from the digital rock models segmented through various image segmentation approaches.The detailed explanation of the maximal-ball algorithm is referred to (Dong & Blunt, 2009;Dong et al., 2007).The use of the maximal-ball algorithm for constructing network models has been successfully validated in prior studies (Al-Kharusi & Blunt, 2007;Blunt et al., 2013; J. Q.Wang et al., 2015).In essence, this algorithm aims to find the largest inscribed spheres for each pore pixel and subsequently identifies the pores and throats in a binarized rock model.The process commences by choosing a pore pixel within a lattice model, from which the largest inscribed sphere that just contacts the solid phase is determined.This procedure is iteratively applied to all pore pixels within the model, resulting in a collection of overlapping spheres.Spheres that remain uncovered by others are referred to as maximal balls, representing the pores within the network.Throats, on the other hand, are the enclosed smaller spheres that establish connections between these maximal balls.Ultimately, this method allows us to generate a pore network that maintains topological equivalence with the original model, facilitating a comprehensive analysis of the pore and throat size distributions.

Single-Phase Flow Simulation
A critical step in DRP following the imaging process involves the prediction of various rock geophysical properties using different numerical simulation methods.The Lattice-Boltzmann method is a widely used direct approach for fluid-flow simulations (Eshghinejadfard et al., 2016;Golparvar et al., 2018;Maier et al., 1996;Roslin et al., 2019).It is particularly feasible for structured grid systems and can be directly applied to latticebased digital images, treating each pixel as a grid element.This method is therefore well-suited for handling complex rock geometries.In this work, a standard LBM-BGK is used to model the single-phase flow and calculate the absolute permeability (Qian et al., 1992).

Segmented Rock Images
Figure 4 depicts the 2D cross-sections of the segmented results for the BS1 models.Additional segmented results for BS2-BS4 models are provided in Figures D1-D3 in the Appendix.In these figures, the 2 μm segmented image is used as the ground truth for the 6 and 18 μm images.The 6 and 18 μm images are segmented by using the thresholding (TH), watershed (WS), and regional-minima segmentation method in 2D (RM2D), and 3D (RM3D).In addition, Figure D4 in the Appendix displays 2D cross-sections of the segmented results for the Doddington sandstone, EC, IL, and MEC models.Similarly, we use the thresholding, watershed, and the proposed method in 2D and 3D to segment these models.Figure 5 shows the segmentation results generated by different methods of the multi-mineral Mount Simon sandstone model.
When comparing the segmented results generated by different methods with ground truth at a resolution of 2 μm, it is evident that as image resolution decreases, an increasing number of rock structures become unresolved (as exemplified in Figure 4).A resolution of 6 μm can still resolve macro pore space, and the primary challenge to the accuracy of segmentation methods is the ability to extract grain interfaces.Due to their relatively small size compared to macro pores, they are easily influenced by the coarsening of the image resolution.Most of these pore structures are not well captured by the thresholding and watershed methods, leading to the merging of grains and The parenthetical text indicates the segmentation methods, including thresholding (TH), watershed (WS), and regional-minima method in 2D (RM2D), and 3D (RM3D).The BS1_2 μm segmented by Otsu's thresholding method serves as the ground truth of the BS1 model.Red rectangles highlight pore structures that the regional-minima segmentation method can extract, whereas classical methods cannot.The white area denotes the solid phase, while the black represents the pore phase.
Figure 5. 2D cross-sections of the MS model with different segmentation methods, including thresholding (TH), watershed (WS), regional-minima method in 2D (RM2D), and 3D (RM3D).The white area denotes the solid phase, while the black represents the pore phase. Water the formation of a compact solid framework throughout the domain.In contrast, the regional-minima segmentation method successfully captures grain interfaces that are missed by classical methods.Consequently, the proposed method can extract detailed rock structures more accurately, as highlighted in the red rectangles in Figure 4.
Notably, the RM2D and RM3D methods tend to produce over-segmented pore space compared to classical methods.This tendency is particularly evident in 18 μm models, where the proposed method yields more small pore spaces even within the solid framework.This may be attributed to image artifacts, such as image noises being treated as regional-minima pixels and classified as pore pixels in the proposed method.Another possibility is that pixels do not correspond one-to-one in images at different resolutions, even though they are paired images.For instance, a 6 μm pixel represents a cubic volume of dimensions 3 × 3 × 3 in the corresponding 2 μm image.Within the 3 × 3 × 3 volume, rock structures may exhibit significant variations, but the 6 μm pixel cannot capture this complexity.Instead, the 6 μm pixel represents an average of the entire 3 × 3 × 3 cubic volume.As a consequence, the 6 μm image appears different from the 2 μm image.As the image resolution coarsens, the discrepancy becomes more pronounced.Consequently, the 18 μm images exhibit significant differences compared to the 2 μm images.
More interestingly, in the multi-mineral MS model, the proposed method can also extract a lot of detailed structures, such as grain interfaces and small pores within solid grains and clay minerals (as highlighted in red rectangles in Figure 5).However, classical segmentation methods only extract macro pore space, missing almost all small pores within the solid framework and grain interfaces.The pore-filling clay has lower grayscale intensities than quartz and feldspar phases.This indicates that the clay pixels will be naturally recognised as regional-minima pixels if they are surrounded by quartz and feldspar pixels in the proposed segmentation method.Nevertheless, as clay has higher intensity levels than pure pore pixels, these regional-minima clay pixels will not be categorized as pure pore pixels in the optimization process.Therefore, this demonstrates the robust performance of the proposed method in handling the complexities of a multi-mineral porous medium.For carbonate rock models, as exemplified in Figure 6, grain interfaces are less abundant compared to sandstone models due to well-cemented solid grains.However, the regional-minima segmentation method still extracts more small pore space.Furthermore, some small pores near macro pores are observed to be connected to macro pores, forming extensions of the macro pore space (highlighted in rectangles in Figures 6e and 6f)., watershed (WS), regional-minima method in 2D (RM2D), and 3D (RM3D).The white area in panels (c)-(f) represents the pure solid phase, the gray denotes the micro-mixed phase, and the black is pure pore.The green in panels (e) and (f) indicates regional-minima pixels within the micro-mixed phase, identified through the regional-minima segmentation method.

Pore Network Statistics
We utilize the maximal-ball algorithm to extract pore networks from all rock models following the image segmentation process.The pore network statistics of these models are summarized in Tables C1 and C2, including the number of extracted pores and throats, average coordination number, as well as their average radius.Figure C1 illustrates the relationship between pore volume fraction and pore/throat radius.
As the image resolution decreases, it is evident that the total number of extracted pores and throats decreases.Furthermore, the average pore and throat radius increases as image resolution coarsens, which aligns with findings reported in Alyafei et al. (2015).This indicates that macro pores are smoothed out, and many small pores become unresolved as the image resolution coarsens (Gooya et al., 2016).In particular, the pore and throat size distribution curves of conventional methods are skewed toward larger size ranges, with limited data pertaining to smaller sizes.In contrast, the regional-minima segmentation method generates more pores and throats with average sizes smaller than those in classical methods.Notably, in carbonate rock models, we find that the coordination number generated using models segmented by the proposed method is consistently higher than with classical methods.This indicates that pores have more connections to other throats.

Absolute Permeability
For BS1-BS4 models, we compare the computed absolute permeabilities along the X-, Y-, and Z-axes of the 6 and 18 μm models to the 2 μm ground truth models.The calculated permeabilities along three axes are listed in Tables E1-E4.The permeability error is calculated by: where E i denotes the permeability error along the i axis, K i is the permeability of the low-resolution models, and K GT,i is the permeability of the 2 μm models.To quantitatively assess the improvement achieved by the proposed method, Table E5 provides a summary of the average errors observed with various segmentation methods across the four Bentheimer Sandstone models.The average error across BS1-BS4 models, E AVG , is calculated by: For MS, DS, EC, IL, and MEC models without high-resolution models as ground truth results, we compare their permeabilities along the Z-axis to their experimental results as this direction aligns with their experimental measurements.The permeabilities of EC, IL, and MEC models are calculated considering the micro pore space as a no-flow phase.The results and errors are listed in Table E6.
Overall, the proposed method in either 2D or 3D scheme shows more accurate permeabilities than thresholding and watershed segmentation methods.Relative to the thresholding method, the proposed approach improves permeability accuracy by over 60% for all four Bentheimer Sandstone models at a resolution of 6 μm, and by over 50% at 18 μm.Compared to the watershed method, the proposed technique exhibits improvements of over 40% for the Bentheimer Sandstone models at 6 μm resolution and over 25% at 18 μm.Notably, for 18 μm models, determining seeds for the watershed transform is hindered by the extremely coarse image resolution.As a result, the watershed segmentation method tends to generate lower porosities than other methods.Certainly, different choices of seeds may yield better results, but this demonstrates the strong user bias that exists in the watershed method.
The proposed method also provides more accurate permeabilities for other sandstone models.In the MS model, the modeled permeabilities are compared to the reported permeabilities (ranging from 10 to 50 mD) for samples with porosities ranging from 8% to 12% (Frailey et al., 2011).Given the wide range of reported permeabilities, all segmentation methods exhibit good agreements.In the DS model, we observe that the thresholding and watershed methods significantly overestimate the permeability, while our method aligns perfectly with the measured permeability.
In carbonate rock models, we also find consistent results.In the EC model, the thresholding and watershed methods yield much lower permeabilities than the measurement.This discrepancy arises because the EC model

Water Resources Research
10.1029/2023WR036855 contains a considerable amount of micro pores.Some of the micro pores that are connected to macro pores can provide additional flow paths.However, classical methods cannot extract these micro pores, resulting in much lower permeabilities.In contrast, the proposed method extracts many small pores around macro pore space (see Figure 6).These small pores may connect to each other and macro pores, forming additional flow paths other than macro pores.Thus, the proposed method demonstrates more accurate permeabilities in the EC model.For the IL and MEC models, all segmentation methods demonstrate good agreement with the experimental results, indicating that micro pore space contributes very little to fluid flow in these models.However, we observe that the proposed method in the 2D case yields less satisfactory results compared to other segmentation methods, particularly in the MEC model.This discrepancy can be attributed to the regional-minima pixel extraction process in RM2D, where we employ a 3 × 3 search window (as described in Section 2.2), with the n index ranging from 1 to 8. Since the range of the n index is relatively coarse, the optimized case often fails to be the most optimal one.In contrast, in the RM3D case, a search window of dimensions 3 × 3 × 3 allows the n index to vary from 1 to 26.Consequently, we can finely tune the n index and identify the optimal case during the optimization process.Therefore, RM3D typically yields more accurate permeabilities compared to RM2D.

Regional-Minima Features
In this section, we delve into the specific influence of image resolution on the rock images themselves and elucidate why the proposed method generates improved segmentation results.Taking the BS1 images at different resolutions as an example, we select three straight lines at identical locations on these normalized images, denoted as A 0 A 1 , B 0 B 1 , and C 0 C 1 , as illustrated in Figures 7a-7c.The images are normalized by dividing the intensity of each pixel by the mean intensity value of the entire image.In Figures 7d-7f, we plot the normalized intensity levels along these lines.These lines intersect both the solid and pore phases and cover phase boundaries, allowing us to observe the intensity transition from the solid to the pore phase.Specifically, line A 0 A 1 captures a macro pore, line B 0 B 1 crosses a grain interface, and line C 0 C 1 encompasses two larger pore gaps.Figure 7d demonstrates that higher image resolution results in a steeper transition between the solid and pore phases, whereas lower resolutions exhibit broader and smoother solid-pore transition bands.This phenomenon highlights that decreasing image resolution tends to smooth out the solid-pore interfaces, resulting in a significantly blurry interface between the solid and pore phases.The smoothing of solid-pore interfaces has been identified as a dominant issue affecting the accuracy of fluid-flow simulations in sandstone rocks as image resolution reduces (Gooya et al., 2016;Schlüter et al., 2014).The presence of these blurry transitions makes it increasingly challenging to precisely locate the boundaries of macro pores at lower resolutions.
Additionally, as the image resolution coarsens, the intensity contrast between peaks and valleys surrounding narrow pore structures becomes flattened and diminished, as demonstrated in Figures 7e and 7f.This phenomenon results in a less distinct separation between the solid and pore phases.This flattening effect is responsible for the unresolved pore space in the constructed models.As depicted in Figure 7e, an intensity valley is clearly presented for the 2 μm image.However, as the resolution decreases, the 6 μm image presents a similar valley that is less distinguishable from the solid phase.In the 18 μm image, the valley disappears, and the intensity in this region becomes comparable to the surrounding pure solid pixels.This indicates that distinguishing small pore structures, such as scattered intra-granular pores and grain interfaces, from the solid phase becomes more difficult as image resolution decreases.The reduced contrast between the pore and solid phases in fine pore structures is primarily due to the presence of the blurry transition.These blurry transitions can have broader widths than the narrow pore structures, causing the pixels covering the narrow pore structures to be classified as transition pixels with significantly higher intensity values than the pure pore pixels.Extracting these structures using the thresholding method would require an extremely high threshold value.Moreover, the wider transition makes it harder to determine the separation between pore and solid phases in the watershed transform.As a result, many of these narrow pore spaces are lost when using the thresholding and watershed methods.However, despite the decrease in resolution causing a flattening of the intensity valley, we observe that the intensity levels of the pixels within the valleys remain lower than those of the surrounding pure solid pixels.Therefore, by incorporating the regional-minima characteristics in the segmentation process, we can successfully extract these narrow pore structures in unresolved images.As a result, the proposed segmentation method can yield more reliable segmentation outcomes and accurately distinguish between the pore and solid phases.

Uncertainties and Limitations
Like any other image processing method, the proposed method is subject to some common uncertainties.
Particularly, image artifacts, especially image noise, can significantly influence our method.Image noise may appear darker than the surrounding high-intensity solid phase or brighter than the surrounding low-intensity pore phase.These noise artifacts may be treated as regional minima and maxima pixels, thereby affecting the segmentation result.Therefore, a proper denoising process is essential before applying our segmentation method to ensure accurate outcomes.
Despite the common uncertainties, there are inherent uncertainties in the proposed method.The initial Gaussian fit process introduces some uncertainties to the segmentation results, as the quality of the fit may vary.While pure phases in the image histogram typically align well with Gaussian distributions, the selection of intensity cutoffs to separate the image into different phases introduces certain user bias.Additionally, in rock models where pure pore phases are not clearly present, manual selection of threshold values is necessary, as seen in the 18 μm BS models used in this work.Similarly, for carbonate rock models like IL and MEC models, the separation of macro-mixed and micro-mixed phases relies on estimation.However, despite these uncertainties, the subsequent steps in the proposed method are conducted without additional manual settings.Therefore, the method is capable of achieving good results while minimizing user bias.
The proposed method has several limitations.First, it requires a sufficiently high image resolution to accurately cover pure pore and solid phases in the images.If the resolution is too low, manual adjustment of intensities may be necessary to distinguish between different phases, introducing additional uncertainties.Second, reference porosity is essential for optimizing segmentation parameters in the proposed method.Without this reference, the segmentation process may be less accurate.Lastly, the method currently only produces binarized models and may not be suitable for multi-mineral segmentation.Although it has been tested on a four-phase Mount Simon Sandstone model, it currently can only extract the pore phase rather than directly segmenting the model into various phases.

Water Resources Research
10.1029/2023WR036855 LI ET AL.

Conclusions
In this study, we proposed a novel segmentation method to overcome the influence of the partial-volume effect in unresolved rock images.This method is based on identifying indistinguishable pore space using regional minima and parameter optimization with maximum entropy theory.
The pore system statistics, such as pore body, pore throat and coordination number, as well as permeability in various rock samples, are compared systematically using the segmented rock images from the proposed approach against that from conventional segmentation methods.The results indicate that the neglected pore space (using the conventional segmentation methods) in the unresolved rock image plays a significant role in the DRP simulation of permeability.Compared with other conventional segmentation methods, the newly proposed method overcomes the influence of the partial-volume effect and extracts the neglected pore space (normally the micropore space) successfully, which results in a more realistic pore system and permeability calculation in unresolved rock images.The proposed method demonstrates a significant improvement in permeability prediction in coarse Bentheimer sandstone images (compared with permeability simulated in higher-resolution models, the proposed method shows 50% improvement against the thresholding method and over 30% improvement against the watershed method in permeability accuracy).The improvement is also found consistent across various rock models tested in this study comparing the simulated permeabilities to their laboratory measurements, including Mount Simon sandstone, Doddington sandstone, EC, IL, and MEC.
One natural extension of the current work is to perform two-phase flow simulations, such as pore network modeling of capillary pressure and relative permeability curves, which may better distinguish the differences between our newly proposed segmentation method and other existing approaches.Unresolved rock images segmented with this novel regional minima-based segmentation method make it possible to develop more accurate digital rock samples from unresolved rock images.This will help us to generate more realistic pore network models and potentially more accurate DRP simulations in large unresolved rock images.[RM2D] and [RM3D] denote the regional-minima segmentation method in 2D and 3D cases, respectively.Hp(n) denotes the pore phase entropy, ϕ is the porosity, and ϕ ref indicates the reference porosity used to calibrate the optimization.The n index corresponds to the number of regional minima chosen from search windows during the regional-minima extraction process.The dotted green line represents the final optimized n index that matches the reference porosity and achieves the highest pore phaser entropy.Note.The parenthetical text represents segmentation methods, including thresholding (TH), watershed (WS), and regional-minima method in 2D (RM2D), and 3D (RM3D). Water

Appendix E: Modeled Permeability
Modeled permeabilities are contained in this section (Tables E1-E6).Note.Permeabilities of 2 μm model segmented by Otsu's segmentation method are taken as ground truth (K GT ) for the validation of 6 and 18 μm models segmented by thresholding (TH), watershed (WS), regional-minima segmentation method in 2D and 3D (RM2D and RM3D).ϕ denotes porosity.K X , K Y , and K Z represent permeabilities along the X-, Y-, and Z-axis.
Corresponding errors are computed by (|K K GT |)/K GT .Note.Permeabilities of 2 μm model segmented by Otsu's segmentation method are taken as ground truth (K GT ) for the validation of 6 and 18 μm models segmented by thresholding (TH), watershed (WS), regional-minima segmentation method in 2D and 3D (RM2D and RM3D).ϕ denotes porosity.K X , K Y , and K Z represent permeabilities along the X-, Y-, and Z-axis.Note.Permeabilities of 2 μm model segmented by Otsu's segmentation method are taken as ground truth (K GT ) for the validation of 6 and 18 μm models segmented thresholding watershed (WS), regional-minima segmentation method in 2D and 3D (RM2D and RM3D).ϕ denotes porosity.K X , K Y , and K Z represent permeabilities along the X-, Y-, and Z-axis.
Corresponding errors are computed by (|K K GT |)/K GT .Note.Permeabilities of 2 μm model segmented by Otsu's segmentation method are taken as ground truth (K GT ) for the validation of 6 and 18 μm models segmented by thresholding (TH), watershed (WS), regional-minima segmentation method in 2D and 3D (RM2D and RM3D).ϕ denotes porosity.K X , K Y , and K Z represent permeabilities along the X-, Y-, and Z-axis.
Corresponding errors are computed by (|K K GT |)/K GT .

Segmentation methods
Average error along axes (%)   Note.The modeled permeabilities are compared to their experimental measurements (Exp).The models are reconstructed by different segmentation methods, including thresholding (TH), watershed (WS), and regional-minima segmentation method in 2D and 3D (RM2D and RM3D).ϕ is porosity.Permeability errors are computed by (|K K Exp |)/K Exp .The experimental permeabilities of the MS are reported from 10 to 50 mD for porosities of 8%-12% (Frailey et al., 2011).The porosities of EC, IL, and MEC models consider only the macro porosity.

Figure 2 .
Figure 2. The workflow of extracting regional intensity minima.(a): A 3 × 3 search window in 2D scheme (RM2D) or a 3 × 3 × 3 search window in 3D scheme (RM3D) are constructed based on the central pixel e. (b):The intensity values of pixels in the search window are sorted, and the lowest n values are selected to be potential regional minima, while the rest are treated as regional maxima.(c): For each pixel e, it is shared by 8 search windows in 2D and 26 search windows in 3D centered on its neighboring pixels and one search window centered on itself.A pixel is considered a regional minimum only if it is identified as a regional minimum in all the search windows.(d): An example of extracted regional minima and maxima distributions for the BS1_6 μm model.

Figure 3 .
Figure 3. (a) and (c) illustrate the pore phase entropies for different n values for the BS1_6 μm model in RM2D and RM3D scheme.The final n value is chosen as the one with maximum entropy while its porosity (ϕ) matches the reference porosity (ϕ ref ).Therefore, n = 6 in RM2D scheme, and n = 16 in RM3D scheme.(b) and (d) are the pore and solid histograms after determining the n value.

Figure 4 .
Figure 4. 2D slices of segmented BS1 model at resolutions of 2, 6, and 18 μm.The parenthetical text indicates the segmentation methods, including thresholding (TH), watershed (WS), and regional-minima method in 2D (RM2D), and 3D (RM3D).The BS1_2 μm segmented by Otsu's thresholding method serves as the ground truth of the BS1 model.Red rectangles highlight pore structures that the regional-minima segmentation method can extract, whereas classical methods cannot.The white area denotes the solid phase, while the black represents the pore phase.

Figure 6 .
Figure6.2D cross-sections of the Estaillades carbonate model segmented using different methods, including Otsu's thresholding (TH), watershed (WS), regional-minima method in 2D (RM2D), and 3D (RM3D).The white area in panels (c)-(f) represents the pure solid phase, the gray denotes the micro-mixed phase, and the black is pure pore.The green in panels (e) and (f) indicates regional-minima pixels within the micro-mixed phase, identified through the regional-minima segmentation method.

Figure 7 .
Figure 7.Comparison of intensity values along three straight lines selected from the BS1 models at different resolutions.The images are normalized by dividing the intensity of each pixel by the mean intensity value of the entire image.Panels (a), (b), and (c) show the normalized images and the three lines on the images at resolutions of 2, 6, and 18 μm, respectively.Panels (d), (e), and (f) depict the normalized intensity values along these lines.Line A 0 A 1 crosses a macro pore, line B 0 B 1 captures a grain interface, and line C 0 C 1 contains two larger pore gaps.The intensity valleys in Panels (d), (e), and (f) represent the pore structure.

Figure B1 .
Figure B1.Entropy optimization of the digital rock models in this work.[RM2D] and [RM3D] denote the regional-minima segmentation method in 2D and 3D cases, respectively.Hp(n) denotes the pore phase entropy, ϕ is the porosity, and ϕ ref indicates the reference porosity used to calibrate the optimization.The n index corresponds to the number of regional minima chosen from search windows during the regional-minima extraction process.The dotted green line represents the final optimized n index that matches the reference porosity and achieves the highest pore phaser entropy.
Corresponding errors are computed by (|K K GT |)/K GT .

Table 1
Digital Rock Models Used in This Work and Their Properties

Table 2
Phases Divided According to Intensity Values Chosen Through the Gaussian Distribution Process Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR036855by University Of Aberdeen The Uni, Wiley Online Library on [26/06/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR036855byUniversityOf Aberdeen The Uni, Wiley Online Library on [26/06/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)onWileyOnline Library for rules of use; OA articles are governed by the applicable Creative Commons Licensethe highest entropy values do not correspond to matched porosity, we opt for the case with matched porosity and slightly lower entropy, as illustrated in FiguresB1o, B1p, B1u, and B1v.

Table E3
Permeabilities of the BS3 Model Computed by Lattice-Boltzmann Method

Table E4
Permeabilities of the BS4 Model Computed by Lattice-Boltzmann Method

Table E6
Permeabilities (Along Z-Axis) of the MS, DS, EC, IL, and Middle Eastern Carbonate (MEC) Models Computed by Lattice-Boltzmann MethodNote.Average error along axes are calculated by (E X + E Y + E Z )/3 for each model, where E X , E Y and E Z are permeability errors along the X-, Y-, and Z-axis.E AVG denotes the average of mean axis errors across the four models.Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR036855by University Of Aberdeen The Uni, Wiley Online Library on [26/06/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License