Quantifying fracture geometry with X-ray tomography: Technique of Iterative Local Thresholding (TILT) for 3D image segmentation

This paper presents a new method—the Technique of Iterative Local Thresholding (TILT)—for processing 3D X-ray computed tomography (xCT) images for visualization and quantification of rock fractures. The TILT method includes the following advancements. First, custom masks are generated by a fracture-dilation procedure, which significantly amplifies the fracture signal on the intensity histogram used for local thresholding. Second, TILT is particularly well suited for fracture characterization in granular rocks because the multi-scale Hessian fracture (MHF) filter has been incorporated to distinguish fractures from pores in the rock matrix. Third, TILT wraps the thresholding and fracture isolation steps in an optimized iterative routine for binary segmentation, minimizing human intervention and enabling automated processing of large 3D datasets. As an illustrative example, we applied TILT to 3D xCT images of reacted and unreacted fractured limestone cores. Other segmentation methods were also applied to provide insights regarding variability in image processing. The results show that TILT significantly enhanced separability of grayscale intensities, outperformed the other methods in automation, and was successful in isolating fractures from the porous rock matrix. Because the other methods are more likely to misclassify fracture edges as void and/or have limited capacity in distinguishing fractures from pores, those methods estimated larger fracture volumes (up to 80 %), surface areas (up to 60 %), and roughness (up to a factor of 2). These differences in fracture geometry would lead to significant disparities in hydraulic permeability predictions, as determined by 2D flow simulations.

Fracture apertures are commonly small compared to the two planar fracture dimensions, so in order to quantify fracture apertures in experimental cores, fine imaging resolution is needed. The latest developments in xCT instruments have enabled unprecedented resolution of fractures within geologic media under a variety of in situ conditions [2,9,52]. For example, laboratory and synchrotron-based xCT instruments that have adopted two-stage magnification can achieve micron to submicron voxel sizes at source-sampledetector separations that accommodate relatively small rock cores (e.g., 4-mm-diameter cores) within vessels that allow flow, temperature, and pressure control [52,54]. In contrast, advanced conventional xCT scanner designs, which rely solely on geometric magnification, provide the X-ray flux and source-detector spacing needed to image larger fractured rock cores (e.g., 2 cm dia.) within larger and stronger vessels [55]. The tradeoffs between voxel size, field of view, sample size, and X-ray flux [2,32,33] suggest that both types of instruments will continue to play important roles in core-flow research for the foreseeable future, and regardless of xCT design, larger datasets will continue to accompany advances in instrumentation.
Once high-resolution reconstructed grayscale images have been generated, the images are segmented, in which voxels are classified into void space and mineral phases. Segmentation techniques for images of porous media are well established and generalized [2,9,[56][57][58][59]. Often used are histogram thresholding algorithms (e.g., [18,23,60,61]), which rely on identification of a single grayscale value to partition intensities into different categories. These "intensity-based" methods are typically computationally efficient, automated, and repeatable, and can be very reliable when different features result in distinctive peaks on the grayscale histograms. Another group of segmentation algorithms, which includes edge detection and region-grow, uses spatial information about local variation in grayscale intensity. They are often regarded to be more reliable than thresholding [52], but are computationally intensive and time consuming partly due to the need for human supervision [56]. More sophisticated methods have also been developed, by utilizing additional criteria or by coupling more than one basic segmentation scheme, such as the active contour method [56,62] and the indicator kriging method [63].
For fractured media, the algorithms commonly used for segmentation of porous media are generally applicable [11,14,19,21,52,64]. For example, Karpyn et al. [14] adopted a thresholding-based approach to differentiate water and oil phase from rock matrix, and Deng et al. [21] devised a ternary thresholding method to identify the fracture, reaction-altered rock matrix, and unaltered rock matrix, while Gouze et al. [11] and Noiriel et al. [52] applied the region growing method. However, there are challenges to direct application of those methods to fracture segmentation. In particular, fracture apertures may be smaller than or comparable to the achievable voxel sizes, especially for fractures under confining pressure and around fracture asperities. The voxels of fine fractures are affected by boundary blurring of both fracture walls and therefore are subject to misclassification due to the partial volume effect [18,29,57,65], making applications of thresholding methods difficult. Several methods have been developed to capture fractures that have fine apertures (e.g., [18,66]). One example is the multi-scale Hessian fracture (MHF) filter developed by Voorn et al. [66], which was inspired by blood vessel detection techniques [67]. It is, however, prone to overestimating the apertures of fine fractures when the aperture sizes vary over a wide range, which is common in natural fractures. A different type of challenge in fracture image segmentation arises when the fracture is in granular rock because of the need to distinguish fractures from the void space in the porous matrix. This requires additional post-processing of the binary segmented images, which can be done by application of morphological filters (e.g., dilation and erosion) and topological operators (e.g., the 3D connectivity filter) [6]. For example, Landry and Karpyn [64] used a series of erosion and dilation steps to isolate fractures from porous media.
A final type of challenge is that the 3D image datasets are often enormous because of the improvements in image resolution that have resulted from recent advancements in xCT scanning technology (as described above). Because fractures have large planar dimensions relative to fracture apertures, and fractures typically occur within rock matrices with high spatial variability of mineral composition and rock texture, it is often necessary to image the entire specimen (plus its flow vessel). The resulting large datasets limit the application of algorithms that are computationally intensive and that require human intervention.
In this manuscript, a novel segmentation routine, the Technique of Iterative Local Thresholding (TILT), is introduced. TILT was devised for xCT image segmentation for fractured rock with particular capabilities to address the challenges of (i) segmenting very fine and highly variable apertures, (ii) distinguishing fracture from pores, and (iii) managing large 3D image datasets. TILT combines a novel method for dilation-based masking for local thresholding with a morphological filter for fracture isolation, in an automated and optimized iterative algorithm. We welcome and encourage the scientific community to use TILT and have made it available online. 1 As a testbed for developing TILT, we used two 3D image datasets that exemplify all three challenges listed above. The datasets are from our recently published investigation of fracture geometry alterations due to mineral dissolution caused by reactions with flowing CO 2 -acidified brine [55].
Ultimately, the goal is to accurately and precisely characterize the geometrical characteristics of fractures to enable sound inferences from laboratory experiments. To examine this, geometric quantities were determined from the binary segmented images produced by TILT, including aperture distribution, total void volume, surface area, and roughness. For comparison, several other common segmentation methods were also applied to the same 3D images to gauge the relative uncertainties in fracture characterization that can be introduced during image processing. Finally, fracture flow simulations were conducted on the fracture geometries using a 2D numerical model. The simulation results provide a quantitative understanding of how differences in image processing methods can lead to differences in inferences about important physical processes such as hydrodynamic processes.

TILT
The complete workflow of the TILT routine is shown in Fig. 1. The automated algorithm takes as input the grayscale 3D image and progresses through two stages-initialization and segmentation optimization-to create a binary image of the isolated fracture.
In the initialization stage, the grayscale 3D image undergoes segmentation and post-processing for the purpose of generating an initial isolated fracture to start the segmentation optimization stage. This coarse segmentation is done using a "global threshold" value that is determined based on the histogram of the entire 3D grayscale image. The thresholding is performed using Otsu's method, a well-established Fig. 1 Workflow of the Technique of Iterative Local Thresholding (TILT). The blue box is the initialization stage, with the dotted arrows illustrating the two steps of initialization: coarse segmentation using global threshold and fracture isolation. The red box highlights the iterative fracture segmentation stage. The red solid arrows illustrate three procedures in a single iteration, and the red dashed arrows show the conditions used to determine continuation or termination of the iterations. In steps where more than one option is available, the options in bold are used in the application for the testbed datasets. The gray box and arrows represent the step for the generation of the fracture template. The detailed procedures of the MHF filter are provided in the workflow in the supplementary materials thresholding method which minimizes intra-class variations and maximizes inter-class variations, and has been demonstrated to be a computationally efficient and robust algorithm [59,68,69].
The resulting coarsely segmented image then undergoes post-processing to distinguish the fracture from pores. Postprocessing for fracture isolation can be realized using two approaches, depending on the properties of the specimen being analyzed. In a porous rock matrix, a fracture template is used, while for a rock matrix that is less porous, the 3D connectivity filter can be applied directly. If the fracture template is used, void objects that intersect with the template are classified as fracture, and the remaining void objects are characterized as pores.
The fracture template is generated by employing the MHF filter, which uses geometrical information in the grayscale 3D image and assigns a "sheetness value" to each voxel enabling differentiation of the pore void space and the fracture void space. (Details about MHF are in the next section and in the Supplementary Material(online resource).) The MHF filter preserves void space that has a planar shape, laying the ground for the fracture template. The 3D connectivity filter is then applied to clean off small and unconnected void spaces. In this way, a fracture template is generated without human intervention. The template is used for fracture isolation in the initialization stage and is passed to the optimized segmentation stage. Therefore, it needs to be generated only once.
Application of TILT to the testbed datasets used the fracture template generated by the MHF filter for the postprocessing because of the porous nature of the Indiana limestone. The MHF filter is computationally costly. The computational time attributable to template generation averaged out to approximately 2 h per 3D image. For the two testbed datasets, i.e., before-and after-reaction xCT scans of the same fractured core, the same fracture template was used to locate the fracture.
In the segmentation optimization stage, a single iteration in the algorithm consists of three procedures: mask generation, image segmentation, and fracture isolation. TILT generates custom masks to delineate a tight domain for local thresholding. It is an alternative to the common approach of a moving window of regular geometry, which can introduce artifacts such as sharp corners. In TILT, the masks are created by morphological dilation of the fractures in binary 2D image slices. (For the testbed application, we used a disk radius of ten times the voxel size.) Local thresholding in each iteration is accomplished using Otsu's method, which was described above. Thresholding can be conducted in either 2D or 3D. In 2D, the thresholds used for segmentation of each 2D image are calculated based on the intensity histograms of each 2D mask, while for 3D processing, the threshold calculation uses the intensity histogram of the 3D imagethe stacked fracture masks. 2D image processing offers computational advantages, but it may introduce bias in segmentation by not including the surrounding pixels in neighboring slices [59,61]. On the other hand, selecting a single threshold for a large stack of images may be problematic because of variations in brightness and contrast over the whole domain. Comparing the 2D and 3D segmentation results provides a way to gauge the magnitude of such biases. Although 3D segmentation could be done using a 3D mask (generated from dilation in 3D), in this analysis, we chose to conduct 3D segmentation using the stack of 2D masks because this approach is more efficient and the stacked masks still preserve all critical information from surrounding voxels at the fracture boundaries.
Post-processing for fracture isolation is then conducted using either the fracture template approach or the 3D connectivity operator. In the application to the testbed dataset, the fracture template generated by the MHF filter was used, as it was used in the initialization stage. The fracture void space identified at the end of an iteration is used as the basis for new mask creation for the next iteration.
The iterative algorithm generates an optimized local threshold, and this allows for practical application to large datasets. Furthermore, this prevents biases and inefficiencies introduced by human selection of a pre-defined threshold value. The iterations increment the threshold and continue until the volume difference of the fracture between two consecutive iterations is smaller than a threshold specified by the user. For the testbed datasets studied, the TILT routine usually required three iterations to converge on a 1 % volume difference.

Other segmentation methods
Here, we present the comparative segmentation methods. The simplest method isolates the fracture region by cropping the images with bounding boxes defined by visual examination of the fracture undulations. The second method that was studied uses contrast enhancement. Here, the contrast-limited adaptive histogram equalization (CLAHE) algorithm [70] was used with an 8 × 8 × 8 kernel. Interpolations were conducted to avoid artifacts at the boundaries. In these two methods, the thresholds were calculated from Otsu's method, as in TILT. For both methods, segmentations were performed in either slice-by-slice fashion, i.e., 2D, or once for the entire 3D dataset.
Two spatial information-based segmentation methods were also compared. One is the active contour (2D)/active surface (3D) method [71,72]. For this method, seeds must be specified, which grow and deform to fill the target object.
This process is user controlled by specifying geometrical factors, which ultimately affect the smoothness of the object's surface [56,72]. The user also has the choice of either image intensity or intensity gradient as the basis for object growth. In this study, 3D application of this method was realized in ITK-SNAP [72]. Because of extremely high time and computational costs [71], the application of intensity gradient in this work was only for a sub-section of the scans. Application of the intensity gradient-based active surface method on 100 images of the fractured core took more than 6 h on a workstation with a six-core processor and 48 GB of memory to produce reasonable results. Segmenting the entire stack of images (˜1800) using the intensity-based active surface method was completed within 2 h on the same workstation, although a threshold needed to be pre-defined. By comparison, processing the entire core using TILT, the iterations needed to generate the optimized threshold and isolated fracture volume were achieved within 3 h.
The other spatial information-based method used was the MHF filter. In this context, the goal was to explore its capability as a stand-alone method, not when used to generate the fracture template that is used in the TILT iterations for fracture isolation. The segmentation routine using the MHF filter first calculates the Hessian matrix of the grayscale intensity and determines the eigenvectors and eigenvalues, which provide information of local principal direction and magnitude of curvatures, respectively. Voorn et al. [66] devised a sheetness parameter, A s , for fracture identification by using the eigenvalues to indicate a planar object. The subscript s refers to the Gaussian scale used to derive the Hessian matrix. For a fracture of aperture a, a Gaussian scale of a 2 captures the aperture most efficiently. In the MHF filter, multiple Gaussian scales that encompass the range of fracture apertures should be applied. At each voxel, A s is calculated for a range of s values, normalized to [0,1], and the value of A s that is preserved is the maximum value. Segmentation of the A s matrix then generates a binary image. In the current study, a scheme with the same concept as Otsu's method was used to segment the A s matrix to generate a binary image of the fracture (Supplementary Material(online resource).
One of the limitations of the MHF filter is that the range of Gaussian scales used has a strong impact on the segmentation results; small scales introduce more noise, while large scales may broaden the features unrealistically. In the testbed dataset for this study, with large variations in fracture apertures, enlarging the range of Gaussian scales to capture large apertures not only increases computational burden but also results in unrealistic broadening of narrow apertures. For instance, processing the entire core of approximately 1800 images took approximately 18 h, with Gaussian scales ranging from 1 to 12, and approximately 24 h with a range of 1 to 16. Therefore, this technique was applied for only the 3D image of the initial experimental scan.
For fracture isolation, post-processing of the binary images produced by the MHF filter was completed using the 3D connectivity operator, while the same post-processing procedure as in TILT was applied on results of other segmentation methods.

Data-acquisition, reconstruction, and pre-processing
The xCT images used as a testbed for this study are reconstructed from xCT scans of specimens of Indiana limestone, which is very porous and comprised of carbonates that are readily soluble in acidic solutions. The rock sample was purchased from Kocurek Industries (product ID B-101a). It was from the Bedford formation, with porosity of˜14 % and brine permeability of˜3 mD. The experimental cores were 5 cm long and 2.5 cm in diameter and were fractured for the experiments using the modified Brazilian method. During an experiment, the fractured core was placed in a high-pressure carbon fiber core holder, and was exposed to CO 2 -acidified brine with pH of 3.3. The fractured core was scanned in this vessel at different stages of the flowthrough experiment. Details about core preparation and experimental conditions are in Deng et al. [55].
The imaging was done using a conventional X-ray CT scanner (North Star Imaging M-5000) at the National Energy Technology Laboratory (NETL) in Morgantown, WV. The Supplementary Material (online resource) describes data acquisition of the attenuated X-ray beam, specimen rotation, and reconstruction of the 3D image. Each 3D image consisted of approximately 1800 image slices, with a total of more than one billion voxels, with voxel sizes of approximately 28 µm. Two 3D images were used as a testbed for this study; these are the initial and final scans from one of the core-flow experiments, the "low P CO2 " experiment, representing before-and after-reaction geometries. Figure 2 shows cross sections from the two reconstructed grayscale 3D xCT images. The attenuation contrast between the rock and the water is sufficient to visually distinguish the fracture from the rock matrix. It is also visually evident that before the experiment the fracture was relatively fine, and after the reactive flow experiment, carbonate dissolution had enlarged the fracture. However, the enlargement is not uniform; there are preferential flow pathways separated by regions where the fracture aperture remains very small, as discussed by Deng et al. [55]. The close-up images reveal the large variation in fracture aperture and the need for a segmentation algorithm that can characterize both fine and enlarged apertures.  Figure 2b shows the histogram of the grayscale intensities from the entire reconstructed 3D xCT image from the scan of the unreacted core. The distribution appears to be uni-modal because the void space accounts for such a small portion of the total volume that the intensities corresponding to mineral matter dominate the histogram. This illustrates that it is not feasible to segment the 3D image simply by partitioning the histogram to determine a global threshold.
Grayscale images with ring artifacts required preprocessing, for which we adopted the combined wavelet-Fourier filtering method [73]. Details of the destripe and denoise procedures can be found in the Supplementary Material (online resource).

Geometrical characterizations and permeability simulation
Fracture vertical aperture, which is the measure of the distance between two fracture walls, was calculated at each location directly from the binary 3D image by counting voxels, column by column. As a result, 2D aperture maps and aperture histograms were generated. Fracture volume was calculated by summing the voxels that are classified as fracture. To quantify fracture surface area, we employed the iso-surface method [23]. The iso-surfaces of the fractures were generated from the binary 3D image using an open-source 3D mesh generator iso2mesh with minimum smoothing to preserve the roughness [74], and the sum of the areas of the surface meshes provides an estimate of the fracture surface area. Surface roughness was quantified by the surface Z 2 parameter, which has been used and documented in details in previous studies [21,75]. The surface Z 2 parameter is the root mean square of the elemental surface slope, and has a larger value for a rougher surface for which the local surface slope varies substantially.
Flow modeling was performed on the 2D aperture maps to estimate fracture permeability. The 2D steady-state local cubic law model [21,76] was selected for this purpose, and the same numerical schemes as in Deng et al. [55] were used. In the application, there was no smoothing; i.e., the resolution of the flow simulation model was the same as the resolution of the image.

Results and discussion
This section presents the results of identification and characterization of the fractures of the two 3D images described in Section 2.1. Throughout, the segmentation results from TILT are presented along with the results from the other segmentation methods. In Section 3.1, we examine the separability of rock and void peaks on the grayscale intensity histogram for TILT and for the other methods. This separability is key to minimizing error in grayscale thresholding. In Section 3.2, we compare the different methods through visual inspection of the fracture in binary segmented images. This illustrates differences in the methods' abilities to characterize fine features, aperture variation, and surface roughness. In the remaining sections, we provide quantitative comparisons of geometrical properties. In Section 3.3, we examine fracture geometries by evaluating the distributions of the fracture aperture and the total void volume of the fracture. Aperture distributions are relevant because these distributions provide insights about permeability evolution [20], channelization [55], geomechanical strain [77], and fluid imbibition and drainage [78] in fractured core experiments. As to the fracture void volume, it provides insights about mineral dissolution in reactive flow experiments [40]. In Section 3.4, we examine how segmentation affects estimates of fracture surface area and surface roughness. Fracture surface area is a key input parameter for reactive transport modeling [79]. Fracture surface roughness is important in flow modeling because it leads to an effective hydraulic aperture that is smaller than the actual mechanical aperture, thus reducing fracture permeability [21]. Finally, in Section 3.5, we examine the extent to which differences in fracture geometry lead to differences in estimation of fracture permeability, a key hydrodynamic variable. This gives us a sense of the importance of the quantitative differences observed in the fracture geometry estimates.
It must be emphasized that application of TILT to these datasets is for purposes of illustration of the TILT routine and its performance in comparison to other segmentation methods. It is not an attempt at benchmarking.

Grayscale intensity histograms
The separability of features on the grayscale intensity histograms provides valuable insights for evaluating the performance of different segmentation methods. Figure 3a and b are the histograms of the grayscale intensities of the entire unprocessed 3D datasets from the initial and final scans of the experiment. Neither histogram shows distinctive peaks corresponding to rock matrix and void spaces and therefore presents a segmentation challenge. For the scan after reaction, in which the enlargement of the fracture has increased the proportion of void space, a very small peak of smaller intensities is present, but the signal is weak. Figure 3c-h of are histograms resulting from applications of the two comparative thresholding methods and TILT in 3D. Applying a rectangular bounding box to crop the images amplifies the signal from void space, but the improvement is not substantial. Similarly, contrast enhancement using the CLAHE algorithm offers only limited improvement. In contrast, TILT enhances the void signal substantially and generates intensity histograms with a clear bimodal distribution for subsequent thresholding. Even for the scan before reaction, in which the fracture is very fine, TILT successfully improves the separability of intensities of rock matrix and void spaces.
Thresholding algorithms (like Otsu's method) are most accurate when the intra-class variances are roughly equal. When the two classes have imbalanced variances in their intensity distributions, the calculated threshold is biased such that intensities belonging to the class with larger variance are misclassified as the class with the smaller variance [80]. Therefore, the type of error introduced in thresholding histograms such as those in Fig. 3c-f is to misclassify rock as void, i.e., to overestimate the void volume. In particular, the voxels that will be misclassified are those with intermediate grayscales such as those at the fracture surface and in sub-voxel apertures, which are subject to boundary blurring. Because TILT amplifies the signal from the fracture voxels and reduces the differences in the intra-class variances of void space and rock, it alleviates this bias. Figure 4 shows a fracture cross section of the 3D grayscale xCT image before reaction and binary images of the isolated fracture resulting from different segmentation methods. By visual inspection, the isolated fractures determined by the 2D and 3D TILT routines match the grayscale image well, as do several of the other methods. Unlike the segmentation results from the stand-alone MHF filter, TILT preserves aperture variations. For each of the thresholding methods, the isolated fractures from the 2D and 3D approaches are in good agreement with each other. The active surface method produces smoother fracture-rock surfaces, largely because of the user-specified geometrical regulations. The result from the edge-based active surface method reveals evident distortions, which are attributable to the complex geometry of the fractured porous media. With this method, edge blurring and edge preservation are difficult to balance. Even with careful initiation of the seeds and large time commitment to allow their slow growth, such distortions are inevitable. Therefore, results of the edge-based active surface method are not discussed further.

Fracture aperture distributions and volume
The distributions of fracture apertures resulting from the different segmentation methods are shown in Fig. 5 for the initial and final scans. While all the methods are effective Fig. 3 Grayscale intensity histograms of the 3D (a), (b) original xCT images, (c), (d) rectangular-cropped images, (e), (f) contrast-enhanced images, and (g), (h) images in the masks used for segmentation at the last iteration of TILT with 3D thresholding, for the initial and final scans of the experiment. The red dashed lines correspond to the threshold values determined from Otsu?s method in capturing the difference between the apertures before and after the reaction, it is the differences between the distributions resulting from the different methods that are discussed here. The distributions of apertures resulting from application of TILT (either 2D or 3D) peak at smaller values, and aperture distributions determined by thresholding images cropped by rectangular boxes and the CLAHE algorithm peak at larger values. As described above, when there is imbalance in the intra-class variances, thresholding is prone to misclassify rock as void, especially for the voxels at the fracture surface and in sub-voxel apertures (as was discussed in Section 3.1). Therefore, the aperture distributions produced by TILT, which peak at smaller aperture values, are expected to be more accurate. Now we focus on identification of zero apertures, which indicate fracture contact areas. These areas are critically important in geomechanical modeling because contact areas support normal stress and constrain fracture slip [77]. It is observed that TILT consistently predicts a higher percentage of zero apertures than do the other two thresholding Fig. 4 a Cross section illustrating the isolated fracture before reaction of the grayscale xCT image and binary images determined using b multi-scale Hessian fracture (MHF) filtering, c edge-based and d intensity-based active surface method, e 2D and f 3D thresholding with rectangular cropping, g 2D and h 3D thresholding with contrastenhanced images, and i 2D and j 3D TILT, from the initial scan of the experiment methods. The reason is, as explained above, that TILT is not inclined to misclassify intermediate grayscale values as void as do the other thresholding methods. Figure 5e compares the aperture distribution generated by TILT for the initial experimental scan with those generated by the MHF filter and the intensity-based active surface method. Both methods are highly sensitive to user-specified parameters. Fracture apertures determined by the MHF filter are largely affected by the ranges of Gaussian scales used (Supplementary Material(online resource)). Because of the sensitivity of the fracture aperture to the value of s, it is not unexpected that most of the resulting aperture values fall between zero and twice the largest value of s. It implies that though the MHF filter provides a useful tool in Fig. 5 Aperture distributions of a, c 2D and b, d 3D thresholding methods for the initial and final scans, and e the MHF filter, the TILT routine, and intensity-based active surface area method (ITKSNAP-int) for before reaction scan (x-axis was truncated at 60) fracture detection, and it has been used to distinguish fracture and pores in the TILT routine and post-processing of other segmentation schemes, it is not ideal for quantification of apertures with substantial spatial variations. Large apertures may be underestimated if the range of s values is not large enough, while narrow apertures can be overestimated when a large range is used, which is evident in Fig. 4b.
Results generated from the intensity-based active surface method are, inevitably, affected by the user-specified thresholding values. In this work, the threshold values used were the same as the values generated by the TILT routine. As a result, the aperture distributions of the two techniques show good agreement for apertures above 4 voxels (Fig. 5e). There are still some differences in small apertures, because the user-specified parameters on geometrical regulations tend to reduce curvatures that may be caused by fine fracture apertures. Figure 6a shows the estimates of fracture volume resulting from the different segmentation methods. The trends in fracture volume are identical to the trends in fracture apertures. The volumes estimated by the rectangular-cropped and contrast-enhanced schemes are consistently larger than those resulting from the 3D TILT routine by +77 and +63 % for the initial scan, and by +16 and +36 % for the final scan. The large range of volumes for the initial scan derives from the large proportion of fine apertures, which are more susceptible to overestimation, as discussed in Section 3.1. After the reaction, in which the fracture had enlarged apertures, the difference between the rectangularcropped method and the TILT routine is smaller, as the masks used for TILT grow with the fractures and approach the size of the bounding boxes from the rectangular-cropped scheme.
Fracture volumes determined by the MHF filter and intensity-based active surface method are more comparable with the results of the TILT routine, showing a difference of −18 and −8 % respectively. Comparing the fracture volumes generated by 2D and 3D processing, the differences are smaller than those among different segmentation schemes, within 10 % for the initial scan that has narrower apertures and 3 % for the final scan.

Surface area and roughness parameters
Shown in Fig. 6b and c are the surface areas and roughness parameters resulting from the different methods. For both, the rectangular-cropped and CLAHE methods generate estimates that are larger than those from TILT. As was discussed in Section 3.1, the rectangular-cropped and CLAHE methods are prone to overestimate void space. This not only leads to larger fracture aperture estimates as shown in Section 3.3 but also increased connectivity between fracture and pores, thereby larger surface area and roughness estimations.
For both surface area and roughness, the range of values is larger for the initial scan than that for the final scan of the experiment. For instance, the surface area determined from the segmented 3D rectangular-cropped images is 5.5 × 10 −3 m 2 , which is more than 50 % higher than that of the 3D TILT routine, while for the surface Z 2 parameter, the 3D CLAHE scheme estimate is larger than the 3D TILT routine estimate by a factor of 2. This implies that there is greater uncertainty in quantifying fracture surface area and roughness for fractures with significant amounts of fine apertures. It is because, as discussed in Section 3.1, for these fractures, the signal from the void space is so weak that the differences in variances of the classes are significant, and the thresholding bias is stronger. Such uncertainty may cloud subsequent interpretations, especially in the case of surface area, which shows opposite trends among different methods. In addition, the uncertainty may be carried into reactive transport modeling, where reaction rate is proportionally related to surface area. Fracture surface areas and roughness parameters derived from the MHF filter are smaller than those of the TILT routine. It is unsurprising given the smooth fracture-rock interface and lack of aperture variations for fracture geometries generated by the MHF filter (Fig. 4). Even though the intensity-based active surface method used the same threshold values as the TILT routine, unlike the aperture distributions and volume estimates, the surface area and roughness parameters are quite different. It is primarily because such parameters are sensitive to local details that are subject to modifications by the geometrical regulations of the active surface method.

Fracture permeability
In this section, we examine the extent to which the differences in geometry parameters derived from the different segmentation methods would affect estimates of fracture permeability. Figure 6d shows the estimated fracture permeabilities for the initial and final scans. The permeabilities resulting from the rectangular-cropped and CLAHE schemes are higher than those from TILT, by 155 and 122 % for the initial scan and 29 and 66 % higher for the final scan in the case of 3D processing. The explanation is not straightforward because the rectangular-cropped and CLAHE schemes generate higher estimates for both fracture volume and roughness. Larger fracture volumes translate to larger fracture permeability, but higher roughness corresponds to lower fracture permeability [21]. The two effects counteract each other.
The estimated fracture permeability determined by the intensity-based active surface method is in agreement with those of the TILT routine, in contrast with the MHF filter. For instance, fracture volume of the initial scan estimated from the MHF filter is 18 % lower than that of the TILT routine, while the estimated permeability is 36 % higher. The difference can be explained by the overestimation of narrow apertures and underestimated roughness.
The estimated fracture permeabilities span a range of a factor of 3 from the different segmentation methods. This range is significant, as it is comparable to the range caused by adopting different modeling approaches for the estimation of hydrodynamic properties [21]. The differences in permeability derived from the 2D and 3D approaches are more significant than those in fracture volume. For example, for the initial scan, the volume difference between 2D and 3D approaches of TILT is˜10 %, while the percentage difference in permeability is˜44 %. The findings indicate that the fine features treated differently in the 2D and 3D approaches could disproportionally affect the flow field in the fractures. Therefore, when fracture geometry is to be used for permeability determination, 3D segmentation is recommended.
Among all the segmentation methods used, fine fracture apertures of sub-voxel scale were not treated specially, and therefore are subject to misclassifications as discussed in Section 3.1. In order to understand the importance of quantifying sub-voxel apertures in affecting flow, we conducted additional simulations. Because of the lack of ground truth, there is no way to differentiate whether the contact points in the segmented images are real or results of misclassification. As a worst case scenario, we assumed that the contact areas in the current segmentation results are all mischaracterized fine apertures, and assigned an aperture value in these areas. The values tested ranged from 0.1 to 30 μm, which is the size of one voxel. An aperture of 0.1 μm gave an effective permeability of 0.8 mD, which is on the same order of magnitude as the permeability of the rock matrix. Simulations using modified aperture maps gave hydrodynamic properties in great agreement with one another, with indiscernible differences (Supplementary Material (online resource)). It implies that for fracture geometries in this study, fine apertures play a very limited role in controlling fluid flow.

Conclusions
In this study, a new automated routine has been developed for segmentation of 3D xCT images of fractures. The method represents an advancement in three regards. The use of dilation for custom-mask generation significantly enhances separability of the grayscale intensities of void space and rock matrix, and therefore reduces potential thresholding bias and enhances the confidence in the classifications of the intermediate grayscale intensities. Second, the combination of adaptive local thresholding with post-processing algorithms for fracture isolation enables characterization of highly variable fractures in rocks that have a porous matrix. Finally, optimized iteration provides an automated platform that minimizes human intervention and enables processing of large image datasets.
TILT outperformed the other segmentation methods in terms of classification accuracy, and therefore produced fracture geometry estimates of higher confidence. The rectangular-cropped and CLAHE methods result in larger fracture apertures and higher fracture-pore connectivity. As a result, the determined fracture volume, surface area, and roughness parameters can be 80 % higher than those from the TILT routine, which led to estimated fracture permeabilities as much as 150 % larger. The MHF filter provides an effective tool in differentiating fractures from pores, but as a stand-alone segmentation method, its quantification capacity is poor in the case of fractures with varying apertures. The results are sensitive to a range of user-specified parameters, which is also true in the case of the active surface method. In our study, the 2D and 3D approaches did not differ significantly in volume and aperture distributions (within 10 %) but showed evident differences in estimated fracture permeabilities, which can be as high as 44 %. This comparison highlights the importance of fine features in permeability estimation, and therefore supports the option of 3D segmentation as it preserves information in all three dimensions. In comparison, sub-voxel scale apertures, which are very challenging in image processing, have very limited impacts on the flow field in fractures. How much the estimated values vary among different segmentation schemes largely depends on the parameters of interest and the geometries of the fractures. For fractures with larger regions of fine apertures, the geometrical quantifications have larger uncertainties. Geometrical parameters such as surface area and roughness are in general more sensitive to segmentation uncertainties, and so is fracture permeability estimation.

Supplementary Material (online resource)
The supplementary material provides details on four subjects. The first section documents data acquisition and reconstruction in detail, and includes a graphical presentation of the xCT system used to collect the image data for this study. In the second section, we present the pre-processing procedures for ring artifacts reduction and noise removal. Third, we demonstrate the workflow of the MHF filter with examples, and discuss the sensitivity of the results of the MHF filter to the range of Gaussian scales. Lastly, we present the permeability estimates for different sub-voxel aperture values.