Characterization of Pore Geometry of Indiana Limestone in Relation to Mechanical Compaction

Characterization of Pore Geometry of Indiana (USA) Limestone in Relation to Mechanical Compaction — We studied the pore structure in intact and inelastically compacted Indiana limestone using X-ray microtomography imaging. Guided by detailed microstructural observations and using Otsu’s global thresholding method, the 3D images acquired at voxel side length of 4 μm were segmented into three domains: solid grains, macropores and an intermediate zone dominated by microporosity. The macropores were individually identified by morphological processing and their shape quantified by their sphericity and equivalent diameter. Our new data revealed a significant reduction of the number of macropores in hydrostatically and triaxially compressed samples with respect to the intact material, in agreement with previous microstructural analysis. The intermediate (microporosity) domains remained interconnected in compacted samples. Our data suggest that the inelastic compaction in Indiana limestone is manifested as not only a decrease in the volume fraction of the microporosity backbone but also a corresponding decrease in its thickness. Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 67 (2012), No. 5, pp. 753-775 Copyright © 2012, IFP Energies nouvelles DOI: 10.2516/ogst/2012051 Pore2Field – Flows and Mechanics in Natural Porous Media from Pore to Field Scale Pore2Field – Physique des écoulements en milieux poreux naturels : de l'échelle du pore à l'échelle du réservoir IFP Energies nouvelles International Conference Rencontres Scientifiques d’IFP Energies nouvelles ogst120006_Yuntao 3/12/12 15:40 Page 753


INTRODUCTION
The mechanical and transport behaviors of a rock are sensitively dependent on its pore geometry. Hence it is of fundamental importance in rock physics to characterize the pore space and its geometric attributes, including pore size, shape, tortuosity and connectivity. Traditionally the pore structure is characterized using optical and Scanning Electron Microscopes (SEM). However recent advances in 3-dimensional imaging techniques such as X-ray microtomography (microCT) and Laser Scanning Confocal Microscopy (LSCM) have provided enhanced perspective on pore geometry complexity. It has contributed useful insights into the preexisting pore space and how it influences rock physical properties, as well as damage evolution and its relation to the micromechanics of failure.
Conventionally the void space of a rock is considered to be made up of two types of cavities: relatively elongated microcracks and equant pores (Brace et al., 1972;Paterson and Wong, 2005). In a compact crystalline rock, many of the preexisting cracks and pores have relatively narrow apertures down to the submicron range, which are not resolvable by X-ray microCT unless one resorts to a sample that includes so few grains that cannot provide a representative elemental volume. However, LSCM (Fredrich et al., 1995) has been used successfully to visualize these relatively fine features in 3D and, even though the technique can be used only to probe a depth that is typically less than 1 mm or so, it has proved to be viable for characterizing the connectivity and modeling permeability (Fredrich, 1999).
In a clastic rock, the void space is dominated by relatively equant pores connected by throats that are sufficiently large for direct imaging using microCT (Spanne et al., 1994;Auzerais et al., 1996) and the data (mostly for sandstones) have provided useful insights into the 3D geometric complexity (Coker et al., 1996;Lindquist et al., 2000) and its control over permeability (Arns et al., 2004;Fredrich et al., 2006;White et al., 2006;Zhan et al., 2010;Sun et al., 2011). In contrast, the use of X-ray microCT for 3D imaging of the pore space in carbonate rocks is not as straightforward (Okabe and Blunt, 2004;Bauer et al., 2011). Carbonate rocks are generally recognized to have pore geometry that is significantly more complex than other sedimentary rocks such as siliciclastics (Choquette and Pray, 1970;Lucia, 1995). One of the reasons for the geometric complexity is that depositional environment and diagenesis exert significant genetic influence over the development of texture and fabric of a carbonate rock (Folk, 1980), which can in turn modify both the size and connectivity of the pore space in a relatively rapid and drastic manner. The pore size in a carbonate rock may span over a very broad range, with a distribution that is often bimodal, including a significant subset of microporosity with submicron features that can only be resolved under a tool such as the SEM (Pittman, 1971;Anselmetti et al., 1998).
Recent studies have demonstrated that the macropores and micropores exert distinctly different influences on the elastic (Baechle et al., 2008;Knackstedt et al., 2009) and inelastic (Zhu et al., 2010;Vajdova et al., 2010) behaviors of limestone and accordingly realistic modeling of the micromechanics would require a dual porosity model. Although most of the micropores can be readily resolved under the SEM, the quantitative characterization of a statistically representative sample of these fine features can be tedious. Furthermore, the 3D characteristics can only be inferred from multiple 2D sections of the sample. Guided by parallel microstructural observations using optical microscopy and SEM, a first objective of this study is to investigate to what extent X-ray microCT can be used to characterize the macropores and indirectly infer the geometry and spatial distribution of the microporosity. Instead of following the conventional approach to binarize the CT image into void and solid, we segment it into three different domains: other than the solid and void (made up of macropores), there is an intermediate domain made up of voxels containing solid material embedded with micropores.
To investigate the micromechanics of inelastic and failure behaviors, it is important to characterize the evolution of damage. As reviewed by Dresen and Guéguen (2004) and Paterson and Wong (2005), most previous studies were performed by SEM since a primary mechanism is stressinduced microcracking, which is difficult to resolve with microCT. There have also been limited studies using LSCM to probe the 3D characteristics of damage evolution (Fredrich et al., 1995;Wong et al., 2004). In carbonate rocks, the damage may derive not only from cataclastic processes but also crystal plasticity processes such as twinning and dislocation slip which can be activated by stress under room temperature (Fredrich et al., 1989;Baud et al., 2000;Vajdova et al., 2010). For inelastic compaction, however, recent studies of Zhu et al. (2010) and Vajdova et al. (2010Vajdova et al. ( , 2012 of several porous limestones have indicated that it is dominated by the collapse of macropores, particularly in the initial yield phase. Since microCT is quite effective in identifying the macropores and characterizing their geometric characteristics, a second objective of the present study is to characterize and contrast these characteristics in undeformed and deformed samples of Indiana limestone, so as to gain insights into the damage evolution and micromechanics of inelastic compaction.

EXPERIMENTAL PROCEDURE
Indiana limestone has been investigated in a number of rock physics and experimental deformation studies. In particular, brittle faulting and its micromechanics were investigated by Robinson (1959), Wawersik and Fairhurst (1970), Zheng (1989) and Myer et al. (1992). Brittle-ductile transition and compactive yield behavior in Indiana limestone were studied recently by Vajdova et al. (2004). Its pore geometry was investigated by Churcher et al. (1991). Recently Vajdova et al. (2012) conducted a relatively detailed characterization of the pore geometry of an undeformed sample as well as damage in deformed samples using optical and SEM, which provides useful guidance on the segmentation of our microCT image.

Sample Preparation and Mechanical Deformation
Cylindrical samples of the limestone were cored perpendicular to the sedimentary bedding. Indiana limestone cores have been ground to diameter of 18.4 mm and length of 38.1 mm. For mechanical testing, a sample was first dried in vacuum at 80ºC for 48 hours and then jacketed in a thin copper foil and polyolefine (heat-shrinkable) tubings to separate the rock from the confining medium (kerosene). Electric resistance strain gages (TML type PFL10-11) were attached in orthogonal directions to the copper jacket to measure axial and transverse strains. Details of the experimental methodology used in Stony Brook laboratory were described in details by Vajdova et al. (2004Vajdova et al. ( , 2012. The jacketed samples were deformed in the conventional triaxial configuration at room temperature. A sample (IH) was hydrostatically compressed to a maximum pressure of 175 MPa. Two other samples (IC1 and IC2) were triaxially compressed at a confining pressure of 20 MPa to axial strains of 0.5% and 1.5%, respectively. The axial displacement was servo-controlled at a fixed rate (corresponding to a nominal strain rate of 1.3 × 10 -5 s -1 ). The volumetric strain was calculated using the relation ε v = ε ⏐⏐ + 2ε ⊥ , where ε ⏐⏐ and ε ⊥ are the axial and transverse strains measured by the gages, respectively. The samples stressed to different stages of deformation were then unloaded. The samples and their deformation histories are summarized in Table 1.
Upon retrieval from the apparatus, the samples were encapsulated with epoxy and two thin-sections were then cut along the central vertical plane of each sample. One thin-section was used primarily for optical microscopy and the other was coated with gold for SEM observations. For microCT imaging, a cylindrical subvolume (of diameter 4 mm and length 7 mm) was also extracted from one half of the sample by coring perpendicular to the cut surface (Fig. 1). For reference, thin sections and subvolume for microCT imaging were also obtained from an undeformed sample (I0-CT) of Indiana limestone.

X-Ray Microcomputed Tomography
and Microstructural Analysis X-ray microCT imaging was performed at the High-Resolution CT Facility at the University of Texas at Austin (UTCT), following procedures detailed by Ketcham and Carlson (2001). The system employs a 200 kV microfocal X-ray source capable of a < 10 μm focal spot size. Scanning specifications were varied to achieve the best imaging for each sample. X-rays were set to 60 kV and 10 W without any filter and a 4X objective was employed. Acquisition times ranged from 50 to 60 s per view and from 476 to 481 views were gathered over an angular range from 190°to 192°. Reconstructed slice images were of size about 1 000 × 1 000 voxels. Each sample spanned 934 image slices and the voxel side length was 4 μm (corresponding to a voxel size of 64 μm 3 ). Microstructure of the intact and deformed samples was also studied under optical microscope and SEM on thin-sections. For SEM observations, the gold-coated thin-sections were studied using a LEO 1550 microscope with a voltage up to 10 KV. All SEM micrographs presented here were acquired in the backscattered electron mode. Optical microscopy was performed using a Nikon polarizing microscope. We also used an Epson Perfection™ V700 photo scanner to characterize the macropore size statistics. The thin sections were scanned at a resolution of 1 800 dpi or higher. Our experience is that the scanner can be used to resolve the macroporosity as effectively as an optical microscope (Zhu et al., 2010;Vajdova et al., 2010), with the advantage that it can cover the whole area of the thin section, thus circumventing the need to assemble a mosaic of numerous optical micrographs. The macropores were identified using a brightness thresholding approach (Russ, 1990) and the binarized image was then analyzed using ImageJ, a public domain image processing program developed at the National Institutes of Health. The area of each individual pore was measured and the equivalent diameter was evaluated.

GEOMETRY OF PREEXISTING PORE SPACE
Indiana limestone used in this study is a gray, fossiliferous limestone that is quarried in the Bedford-Bloomington area in Indiana. It was formed during Mississippian period and belongs to the Salem limestone formation (Patton and Carr, 1982). Its average modal composition is: calcite 97.1%, magnesite 1.2%, silica 0.8%, alumina 0.7%, iron oxide 0.1% and undetermined material 0.1% (Indiana Limestone Handbook, ILI, 2007). The allochems (that include fossils, ooids and some peloids) constitute ~65% of the bulk rock volume. Many of the ooids contain fossils and fossil fragments which are coated with concentric layers of calcite.
In terms of rock fabric, Indiana limestone can be classified as a calcite-cemented grainstone. Due to the material heterogeneity associated with three very different solid constituents (allochems, micrite, sparite), the grain sizes in our Indiana limestone samples range from < 5 μm for the pore-lining micritic material to > 300 μm for the allochems. The porosity of Indiana limestone is quite variable (ILI, 2007) and a range between 12% and 20% have been reported in the literature. We calculated the total porosity from the density of a vacuum dried sample, assuming a solid composition of 100% calcite with density of 2 710 kg/m 3 . The initial porosities so determined for our Indiana limestone samples range from 16.5% to 17.6% (Tab. 1).  Sample preparation for μCT and microstructural analysis: 1) intact and deformed samples of 38.1 mm length and 18.4 diameter were cut in two halfs. One half was impregnated with high viscosity epoxy and 2) a petrographic thin section was prepared; 3) a small sample of 7 mm length and 4 mm diameter for μCT imaging was cored in the impregnated half of the sample. Figure 2 illustrates the microstructure of intact Indiana limestone (I0-CT) observed under an optical microscope. The allochems are commonly coated with micrite cement around their rims and the interparticle porosity is made up of relatively large pores (areas with lightest color). Some of these pores are partially filled with sparry cement. In addition to macroporosity, micropores occur especially within allochems. High microporosity is also concentrated along boundaries of allochems. Extensive SEM observations of the micropores were presented by Vajdova et al. (2012). We also scanned the thin-section (covering an area of 38 × 18 mm 2 ) at a resolution of 1 800 dpi and after the scanned image was binarized, it was analyzed to determine the area of each individual pore, from which the equivalent diameter was evaluated. Given the thickness of the thin section and limited resolution of the scanned image, following Anselmetti et al. (1998), we arbitrarily used a lower limit of 33 μm on pore size, so that a macropore is defined to be a feature with equivalent diameter of 33 μm or larger. From geometric probability theory, the 2D measurement of areal porosity (sum of the pore areas normalized by the total area) can be used to estimate the macroporosity of the sample if it is isotropic (Underwood, 1970). For an undeformed sample (sample I0-TS in Vajdova et al., 2012), the macroporosity was so determined to be about one-third of the total porosity of 16.5%, with the implication that the microporosity represents a very significant fraction. Detailed description of the microstructural data from 2D measurements was presented by Vajdova et al. (2012).

Segmentation of MicroCT Image into Three Domains
Certain first-order geometric characteristics of the pore space can be inferred from the global histogram of the X-ray attenuation coefficient. Figure 3a shows the microCT data for a Fontainebleau sandstone with porosity of 15%, comparable to that of our Indiana limestone sample I0-CT. The 8-bit data were acquired by Lindquist et al. (2000) from a cylindrical sample of diameter 4.53 mm and length 2.91 mm, with a voxel dimension of 5.7 μm. The X-ray attenuation coefficients were digitized and mapped onto a gray level in the range [0, 255] and the relative frequency corresponds to the number of voxels normalized by the total number of voxels in the sample (or equivalently, the voxel volume normalized by the total sample volume). One can readily recognize in the histogram two distinct modes, with relatively symmetric distributions on either sides of the peaks. The two modes are separated by a broad range of gray level, with a well-defined intermodal valley in between. Using any one of the standard protocols in image analysis, one can select an appropriate threshold that falls on the intermodal valley, which would then provide a criterion for binarizing the CT-image into void and solid domains, that corresponds to voxels with gray levels below or above this threshold. For comparison, we show in Figure 3b the global histogram for our 16-bit CT image of Indiana limestone sample I0-CT. Although two peaks can be recognized, unlike the sandstone histogram (Fig. 3a) the distributions on either sides of the peaks are not symmetric. Instead of a well-defined intermodal valley, there is a long tail to the left of the second peak that persists to relatively low gray levels to the right of the first peak. This suggests that there is an intermediate domain with gray levels between those for the void and solid domains, with the implication that voxels in this domain are embedded with micropores that are smaller than the dimension of 4 μm. It is likely that delineation of this intermediate domain would provide useful information on the spatial distribution of microporosity in the limestone. Before one can segment the CT image into three distinct domains, it is necessary to first select two thresholds that lie between the two modes in the global histogram. In this study, we determined the thresholds using Otsu's (1979) method that maximizes the between-domain variance. If the two global thresholds so selected are denoted by I * void and I * solid (with I * solid > I * void ), then the following criterion on the gray level I of a voxel can be used to identify to which of the three domains it belongs: As illustrated in Figure 4a, the global histogram for I0-CT has two modes corresponding to gray levels of I p void = 16 125 and I p solid = 27 848 and the two thresholds determined using Otsu's method in the intermodal region were I * void = 19 472 and I * solid = 27 287. Since Indiana limestone is basically monominerallic (with a modal composition of 97% calcite) and the CT attenuation coefficient is primarily dependent on the density, the gray level of a voxel is expected to be linearly related to its local porosity ϕ. Accordingly we adopted the following relation for inferring the local porosity from the CT measurement: (2) The local porosity so inferred is plotted in Figure 4b as a function of the voxel volume normalized to the sample volume. The distribution follows an approximately bilinear trend, with significant proportion of the voxels inferred to have local porosities < 20%. In Equation (2), we selected I p void and I p solid as the limits and this choice would have been unambiguous if the two peaks corresponding to void and solid are given by Dirac delta functions. However, the histogram for the CT data invariably deviates from such an idealized distribution and hence there are other choices one can make and alternative expressions which one can propose for inferring the porosity and it is a moot question which will give the best approximation. It should be noted that if we perform an integration over an appropriate range of grey level, our approach is consistent with that of Bauer et al. (2011) who did not explicitly consider the local porosity distribution, but instead identified I in Equation (2) with the mean grey level of the image and used the expression for inferring the global porosity. According to our segmentation analysis, the microCT image can be separated into three distinct domains (Fig. 4a): voids with dimensions larger than the voxel resolution (corresponding to the macropores) that occupy 4.96% of the bulk sample volume, solid grains that occupy 47.84% of the bulk volume and an intermediate zone Distribution of X-ray attenuation (relative frequency) for a) an intact sample of Fontainebleau sandstone (Lindquist et al., 2000) and b) our 16-bit CT image of intact Indiana limestone (sample I0-CT).
embedded with micropores) that occupies 47.16% of the bulk volume. From Equation (2), we would infer that the intermediate zone includes numerous voids that occupy 10.79% of the bulk volume. In other words, from the CT image the limestone sample is inferred to have a total porosity of 15.75%, that is partitioned between a macroporosity of 4.96% and microporosity of 10.79%. It should be noted that these values are in good agreement with other independent measurements (Tab. 1). From density measurement of the dry sample, the total porosity of I0-CT was determined to be 18.1%. From 2D measurements on scanned image of the thin-section (of another sample I0-TS of total porosity 16.5%), the macroporosity was inferred to be 5.2% (Zhu et al., 2010;Vajdova et al., 2012). The difference between the two would give a microporosity of 11.3% for I0-TS. It should be emphasized that, for this stage of our analysis, the segmentation into three phases was conducted on the microCT image of Indiana limestone that had not been subjected to any filtering. However, we have unpublished CT data for Majella limestone that has a histogram with an "elbow" in the intermediate zone which is not as well defined and for this limestone with a porosity of 31%, we found it necessary to apply appropriate filtering to accentuate the contrast among the solid, void and intermediate zone before we could segment the image into three phases. Indeed, Bauer et al. (2011) concluded that such pre-segmentation filtering was absolutely necessary in their recent study of Estaillades and Lavoux limestones (with porosities of 24.7% and 28.7%, respectively). This study was kindly pointed out to us by one of the reviewers, since we were not aware of it during our own investigation. Most probably in Estaillades, Lavoux and Majella limestone which are significantly more porous than Indiana limestone, the intermediate zone covers a larger volume fraction and the transition is more gradual, which would necessitate the use of filtering to enhance the contrast before one can perform the segmentation.
Understandably our preference is to minimize any unnecessary filtering unless it is absolutely essential and the quality of our data is such that we do not see any justification or advantage for filtering the data during this first stage of segmentation. Nevertheless, in a later stage during our morphological analysis when we attempted to isolate the macropores and analyze the connectivity of the intermediate zone, we did find it necessary to first apply a median filter to the CT-image. Accordingly it is useful to perform an analysis of how our segmentation results and porosity estimates are sensitive to filtering. This analysis is detailed in Appendix.
In Figure A1 in Appendix, we present the histograms for the data in Figure 4a after they have been subjected to median filters applied to subvolumes of 5 3 and 11 4 neighbor voxels, respectively. It can be seen that since the histograms are qualitatively similar, not surprisingly the thresholds determined for the unfiltered (Fig. 4a) and two filtered (Fig. A1) images almost coincided with one another. However, there is an overall trend for the right peak (corresponding to solid) to attain higher value with increased filtering. As a result, solid volume fraction is predicted to be higher in a segmented image after filtering, whereas total porosity is lower. For example, for the filtered image in Figure A1a  to Estaillades and Lavoux limestones) the solid volume and total porosity have values of 54.15% and 14.08%, about 6% higher and 2% lower than the corresponding values in the unfiltered image (Fig. A1).
A 3D visualization of the gray level distribution is presented in Figure 5a. Whereas many of the macropores were observed to be isolated features (Fig. 5b), most of the solid and intermediate domains seem to be interconnected. In Figure 6,   we present five slices of a subvolume, which show many features qualitatively similar to the optical micrographs in Figure 2. Many of the solid grains and macropores seem to be surrounded by layers of the intermediate domain. It should be noted that since these images have not been filtered, there are many isolated, relatively small blobs (especially in the intermediate domain) which may correspond to random noise.

Geometric Characterization of the Macropores and Microporosity Domain
To characterize geometric attributes of the pore space from the X-ray tomography images, it is necessary to perform morphological processing on the microCT image. As a first step, we applied a median filter, which is considered to be an effective method to remove random noise from an image (Russ, 1990). The basic idea of median filtering is to analyze the gray level distribution in a neighborhood that surrounds a voxel and then replace the gray level of the voxel by the median value of the distribution. In this study the median filter was applied to a neighborhood corresponding to a cubical array of 11 × 11 × 11 voxels. Such a procedure can potentially reduce the random noise while preserving the edges. The surface roughness is reduced while the volume of an isolated blob would remain basically unchanged. However, it may also shrink or even remove some of the relatively small blobs. The impact of the size of the median filter on the results is discussed in Appendix. We next used the two global thresholds (obtained using Otsu's method as described above) to segment the filtered image into three distinct domains. These thresholds were determined for the unfiltered image but as discussed above, they almost coincide with thresholds obtained for the filtered images). Using morphological concepts of connectivity on the segmented image, we identified individual clusters of interconnected voxels that correspond to blobs of solid, macroporosity or microporosity. The last step was to map out the surfaces of these blobs using a "marching cubes" algorithm (Lorensen and Cline, 1987). The algorithm considers an imaginary binary cube created from eight adjacent vertices. Each of these vertices is assigned a binary value of one or zero according to the gray level and domain type. It can be shown that there are fourteen possible configurations for polygonal isosurfaces that pass through the imaginary cube (including the configuration with no surface). For a given combinaison of eight vertex values, a unique polygonal isosurface can be determined and if this process is repeated for all the voxels, the union of these polygons that have been mapped out for all the marching cubes would correspond to the totality of surfaces that enclose the domain.
According to our analysis, a significant proportion of the macroporosity is associated with isolated voids. We present here six macropores mapped out using the marching cubes algorithm. Figure 7a,b shows two relatively small pores with shapes that are approximately spherical and ellipsoidal, respectively. Figure 7c shows two quasi-ellipsoidal blobs connected by a channel with a diameter on the order of 10 μm. The polygons covering the surfaces are also shown in these figures. Figure 7d-f shows pores with increasing geometric complexity, with multiple blobs connected to one another by elongated channels. The surfaces of each of these three complex pores are covered by a relatively dense mesh of polygons, which is not shown in the figure because individual polygons are not readily resolvable.
For an isolated macropore with volume V and surface area S, we used two parameters to characterize its geometry: The equivalent diameter D is the diameter of a sphere with volume V identical to that of the macropore. The sphericity ψ (Wadell, 1935) is the ratio of the surface area of a sphere (with the same volume as the macropore) to the surface area S of the macropore. It is a measure of roundness, such that ψ = 1 for a sphere and a ψ value significantly less than unity is associated with a highly nonspherical geometry or relatively rough surface. This is illustrated by the six macropores in Figure 7 with sphericity ranging from 0.21 to 0.93.
From our morphological analysis, 872 isolated macropores that contribute to a porosity of 3.35% were identified and their equivalent diameter distribution is shown in Figure 8a. Most of these isolated pores have equivalent diameters less than 200 μm and there seems to be an overall trend for a more significant proportion of the larger pores to deviate from a spherical geometry. For comparison, we show in Figure 8b the data of Vajdova et al. (2012), which were obtained on 2D scanned image of the thin-section of I0-TS. They did not characterize the pore shape and the equivalent diameter was calculated for a circle with area identical to the projected area of the pore observed on the 2D section. Both sets of data fall on similar ranges of equivalent diameter and they follow qualitatively similar trends with number of pores decreasing with increasing diameters. Using stereological concepts (Underwood, 1970) and assuming isotropy, the 2D and 3D measurements (of number per unit area n A and per unit volume n V ) of a system of spherical objects are related by: where D denotes the mean diameter of the spherical objects.
Although many of the macropores deviate from spherical geometry, the 2D and 3D measurements of pore size statistics are in approximate agreement with the stereological relation (4). For example, if we consider equivalent diameters in the range of 50-100 μm, the number of macropores per unit area is n A~ 2.5 mm -2 (Fig. 8b) and using a mean diameter of D -= 75 μm, we would expect n V~ 33 mm -3 . In comparison, a value of n V = 11 mm -3 was obtained from our microCT data analysis (Fig. 8a).
For the intermediate domain, our morphological analysis could identify very few isolated clusters of microporosity that altogether occupy < 1% of the voxel volume in this domain. In other words, the microporosity domain in the limestone sample is basically made up of one interconnected backbone. We will denote the voxel volume occupied by this microporosity backbone normalized by the bulk sample volume by v bb and the ratio between the backbone volume and surface area by ρ. In sample I0-CT, these two geometric parameters for its microporosity backbone are: v bb = 41.09% and ρ = 26.3 μm.

Mechanical Data
In presenting the mechanical data, we have adopted the convention that compressive stresses and compactive strains (i.e. shortening and porosity decrease) are positive. The maximum and minimum principal stresses will be denoted by σ 1 and σ 3 , respectively. In Figure 9a, b, we present the stress-strain curves for the three deformed samples investigated in this study.
Sample IH was hydrostatically compacted to a maximum confining pressure of 175 MPa (Fig. 9a). It has a total porosity of 14.6%, lower than those of the other three samples investigated in this study. As noted by Vajdova et al. (2004), the hydrostat has an inflection point at P* ~ 60 MPa, which marks the onset of inelastic compaction and pore collapse (Vajdova et al., 2012). Hence the mechanical data suggest that the sample would have accumulated as much as 2% of permanent volumetric strain.
Samples IC1 and IC2 were triaxially compressed at confining pressure of 20 MPa. While the initial porosities of these two samples were comparable (17.6% and 17.3%), they were somewhat lower than the total porosity of I0-CT (18.1%). Guided by the data of Vajdova et al. (2004) on the critical stress C* for the onset of shear-enhanced, sample IC1 was deformed to just beyond C* before it was unloaded. In contrast, sample IC2 was deformed well beyond C* and therefore expected to have sustained significant shear-enhanced compaction and pore collapse.

Microstructural Observations
The inelastic damage induced by hydrostatic and triaxial compression has been characterized from optical microscopy and SEM observations. Damage in the inelastically compacted samples is primarily associated with collapse of macropores. The macropores were divided into three groups according to their sphericity: sphericity > 0.78 (blue), 0.59 < sphericity ≤ 0.78 (green), 0 < sphericity ≤ 0.59 (brown); b) Size distribution of pores resolved from a scanned image at a resolution of 2 400 dpi (Vajdova et al., 2012). The number of pores per unit area n A is plotted versus equivalent diameter.
In the sample IH which was hydrostatically compressed to an additional 100 MPa beyond the critical pressure P* (Fig. 9a), macropores at different stages of collapse were observed (Fig. 10a). Intragranular cracking at an early stage can be observed in allochems, such as the one marked "A" in Figure 10b with numerous microcracks near the segment of the allochem boundary along which the micritic cement had delaminated. In this hydrostatically compacted sample, cataclastic pore collapse seems to involve largely the cement, with most of the allochems remaining relatively intact. Figure 10c illustrates the later stage of collapse of a macropore within cement surrounded by five allochems. Intense microcracking had developed in the cement surrounding the macropore but the damage mode and intensity were heterogeneously distributed around the pore, possibly due to the mechanical contrast between the micrite and sparry cement. Numerous cracks had emanated from the micropores and coalesced. Large and small spalled fragments have imploded into the cavity. Unlike a micritic limestone such as Tavel limestone which has many quasi-spherical macropores (Vajdova et al., 2010), the geometry of a macropore in an allochemical limestone can be quite irregular, adding to the heterogeneity of the spatial distribution of cataclastic damage in the periphery of a collapsed pore (Fig. 10b,c).
In sample IC1 triaxially compressed at 20 MPa confining pressure to 0.5% of axial strain, no significant damage was visible, most probably because this sample was not deformed to a high enough amount of plastic strain. In the sample IC2 triaxially compressed at the same pressure but to a higher amount of strain, the development of shear-enhanced compaction was manifested by pervasive collapse of macropores (Fig. 11a). Cataclastic pore collapse is mostly confined to within the cement and the allochems are relatively intact in this sample (Fig. 11b). Vajdova et al. (2012) showed that cracking within the allochems is typically occurring after accumulation of a significantly larger amount of plastic strain. The collapse of macropores in sample IC2 is accompanied by cataclastic damage in their peripheries analogous to those observed in the hydrostatically compacted sample IH (Fig. 11c).

MicroCT Imaging of Pore Collapse
The global histograms of the three deformed samples (Fig. 12) are qualitatively similar to the undeformed sample (Fig. 3b). Again even though two peaks can be recognized, unlike the sandstone histogram (Fig. 3a) the distributions on either sides of the peaks are not symmetric and there is a long tail to the left of the second peak that persists to relatively low gray levels to the right of the first peak. Accordingly, we follow the same protocol as before to segment each of the CT image into three distinct domains, by determining two thresholds I * void and I * solid using Otsu's (1979) method. Again, it should be emphasized that, for this stage of our analysis, the quality of the data is such that segmentation into three phases was conducted on the microCT images of the deformed samples which had not been subjected to any filtering.
With reference to the histogram in Figure 12a, sample IH has two modes corresponding to gray levels of I p void = 12 831 and I p solid = 28 422 and the two thresholds determined using Backcattered SEM images of sample IH deformed hydrostatically beyond P*. a) Pervasive collapse of macropores; b) Intragranular cracking at an early stage can be observed in allochems, such as the one marked "A"; c) Collapsed macropore within a cemented region that is surrounded by five allochems.

Figure 11
Backcattered SEM images of sample IC2 triaxially deformed at 20 MPa of confining pressure. Direction of σ 1 is vertical. a) Pervasive collapse of macropores; b) Early stage of the collapse of a relatively small macropore; c) High magnification view of a collapsed macropore with cataclastic damage in its peripheries.
Otsu's method in the intermodal region were I * void = 17 771 and I * solid = 27 826. Similarly sample IC1 (Fig. 12b) has two modes at gray levels of I p void = 13 297 and I p solid = 28 440 and the two thresholds determined using Otsu's method were I * void = 18 126 and I * solid = 27 520. Sample IC2 (Fig. 12c) has two modes I p void = 13 034 and I p solid = 29 045 and the two thresholds were I * void = 17 762 and I * solid = 28 294. According to our segmentation analysis, the microCT image of each deformed sample can be separated into three distinct domains and the partitioning of voxel volume and porosity among the three domains is summarized in Table 1. Five slices from a sub-volume of samples of Indiana limestone a) IH deformed hydrostatically beyond P*, and IC1 b) and IC2 c) triaxially deformed beyond C* at 20 MPa of confining pressure. The partitioning between solid grains, macroporosity and the intermediate zone dominated by microporosity is shown.
Using Equation (2), the local porosity in the intermediate zone was calculated and plotted in Figure 13a. The inferred porosity distributions for the three deformed samples almost coincide with one another. For comparison, we also include the distribution for I0-CT, which seems slightly higher. This discrepancy between the undeformed and deformed samples is possibly related to the backbone geometry of the microporosity domain in these samples.
For each of the three deformed samples, after appropriate filtering as discussed above, our morphological analysis could identify only very few isolated clusters of microporosity that altogether occupy 1% or less of the voxel volume in its intermediate domain, which indicates that the domain is dominated by one interconnected backbone that percolates the sample. The voxel volume fractions v bb of IH, IC1 and IC2 are respectively 34.07%, 35.44% and 31.20%, which are significantly lower than the value of 41.09% for the undeformed sample (Fig. 13b). This decrease in the volume fraction of the microporosity domain is possibly related to inelastic compaction sustained in the samples.
Furthermore, the volume-to-area ratio ρ in the deformed samples IH, IC1 and IC2 has values of 21.2 μm, 21.9 μm and 21.6 μm, respectively (Fig. 13b). These are also significantly lower than the value of 26.3 μm in the undeformed sample I0-CT. If the microporosity backbone can be approximated as a thin-walled shell that percolates the sample, then the mean thickness of the shell is given by 2ρ. Our data therefore suggest that the inelastic compaction in Indiana limestone is manifested by not only a decrease in the volume fraction of the microporosity backbone, but also a corresponding decrease in its thickness from ~52 μm to ~42 μm.
In Figure 14, we present slices of segmented subvolumes of the three deformed samples. The spatial distribution of the three domains seems qualitatively similar to that in the undeformed sample (Fig. 6). However, the macropores appear to be somewhat smaller and to an extent this is supported by the macropore size distributions of the deformed samples. Among the three, IH (Fig. 9a) and IC2 (Fig. 9b) were not unloaded until well after inelastic yielding and therefore expected to have significant damage and pore collapse, which should be reflected by corresponding decreases in the number of macropores observed in the microCT images (Fig. 15a, c). In contrast, IC1 was unloaded soon after the onset of shear-enhanced compaction and expected to have sustained relatively little damage and accordingly one would expect its pore size distribution (Fig. 15a) to be comparable to that in the undeformed sample I0-CT. Yet a comparison of Figure 8a with Figure 15a indicates a fundamental difference in the pore size distributions of I0-CT and IC1, which we interpret to be likely due to variability between samples and not mechanical deformation. We identified significantly more macropores in IC1 (1 146 versus 872 in I0-CT), many of which have relatively high sphericity and small equivalent diameter. If we defined a complex macropore as one with Distribution of equivalent diameter of the macropores identified by morphological analysis in the deformed samples of Indiana limestone IH a); IC1 b) and IC2 c). The number of pores per unit volume n V is plotted versus equivalent diameter D. The macropores were divided into three groups according to their sphericity: sphericity > 0.78 (blue), 0.59 < sphericity ≤ 0.78 (green), 0 < sphericity ≤ 0.59 (brown).
Using the marching cubes algorithm, the surfaces associated with the three domains were mapped out. Our morphological analysis of the microCT images shows that both the solid and intermediate domains are basically interconnected. In particular, the intermediate domain has a percolative backbone that occupies > 99% of the voxel volume, with a characteristic thickness (inferred from twice the volume-to-area ratio) on the order of 53 μm. As for the void domain, a significant proportion of its void space can be described as isolated macropores (Fig. 7). While the pore size statistics show qualitatively similar trend as the 2D data (Fig. 8), our new data underscore some of the 3D complexities of the macropores. We observed a number of complex macropores that extend over a length of > 1 mm (Fig. 7c).
Our observations point to two important questions related to the mechanical and transport properties of porous carbonate rock. First, what are the elastic and inelastic responses of these macropores to an applied stress field? Typically the cavities in a rock are modeled as one of two end-members: equant pores or elongated cracks. It seems that a complex pore in Indiana limestone is more akin to a cluster of equant pores connected by cracks and a deeper understanding of its deformation behavior would require systematic numerical simulation on the realistic geometry acquired from microCT. Second, how are the macropores and micropores connected among themselves to provide a percolative conduit for fluid transport? Since our observations indicate that the macropores are not interconnected, they by themselves cannot provide a percolative path for a fluid, which must therefore traverse parts of the interconnected microporosity domain to flow through. Mercury injection data on Indiana limestone (Churcher et al., 1991) indicate that the micropores are connected by throats with constricted apertures that are submicron. To quantitatively simulate the transport processes, more elaborate and systematic studies of these issues are warranted beyond the morphological analysis we pursued here. Our focus here is inelastic compaction and damage development and therefore, these issues are beyond the scope of this study. Nevertheless, there are recent investigations notably by Okabe and Blunt (2004), Al Kharusi and Blunt (2007) and Bauer et al. (2011) which have advanced the understanding on this front.
This study represents one of the first attempts to characterize the damage accumulated in an inelastically compacted carbonate rock using microCT imaging. Previously microstructural observations on thin-sections have shown that cataclastic collapse of the macropores is the dominant micromechanical process (Zhu et al., 2010;Vajdova et al., 2010Vajdova et al., , 2012. For three hydrostatically and triaxially compacted Indiana limestone samples, we have demonstrated that the CT image can also be segmented into three distinct domains.
Additional insights were gained from our 3D data. In all three compacted samples, the microporosity domains remain sphericity ψ < 0.3, similar to the example in Figure 7f , then our analysis shows that IC1 has only two such pores that contribute to 0.15% porosity. In contrast, I0-CT has six relatively large complex pores that take up as much as 1.66% porosity.
Even if we take into account sample variability, since the numbers of macropores identified in IH (359) and IC2 (787) are lower than either I0-CT or IC1, we consider it likely that these decreases are related to inelastic damage. More pronounced decreases were inferred in pores with equivalent diameters of 100 μm or less (Fig. 15a,c), which are probably related to the collapse of this subset of macropores of relatively high sphericity, as documented in SEM observations (Fig. 10).
However, the scenario for the complex macropores may be quite complicated. Our analysis shows that sample IH has three relatively large pores (with ψ < 0.3) that take up as much as 2.63% porosity and IC2 also has three large complex pores that contribute 0.96% porosity. It is possible that these features are related to mechanical damage but sample variability cannot be ruled out either. This issue cannot be resolved in the absence of a more comprehensive investigation into samples deformed to different stages of deformation.

CONCLUSIONS
Given the geometric complexity of its pore space, standard techniques for binarization of microCT data of a porous siliciclastic rock are not directly applicable to a carbonate rock. Hence most investigations of its pore structure have been done on 2D thin-sections, synthesizing observations on different scales using the optical microscope and SEM. Additional information on geometry of the pores and throats may also be gained from mercury injection (Churcher et al., 1991) andNMR imaging (Gleeson andWoessner, 1991;Kleinberg, 1999). In this study, we developed a methodology that can extract 3D information from microCT data on the partitioning of porosity and pore size statistics in Indiana limestone.
The microCT image of Indiana limestone can be partitioned into three distinct domains by identifying two thresholds using Otsu's (1979) method. There is an intermediate domain with attenuation coefficients that fall between those for the solid and void domains. Guided by microstructural observations and using a dual porosity framework, this intermediate domain was interpreted to be made up of voxels of solid embedded with micropores below the microCT resolution (Fig. 6). Using a linear interpolation scheme, the local porosity in this domain can be inferred and in agreement with previous microstructural observations on carbonate rocks (Baechle et al., 2008;Zhu et al., 2010). The cumulative microporosity in this intermediate domain was estimated to be more than half of the total porosity in Indiana limestone (Tab. 1).
interconnected. Our data suggest that the inelastic compaction in Indiana limestone is manifested by not only a decrease in the volume fraction of the microporosity backbone but also a corresponding decrease in its thickness bỹ 10 μm (Fig. 13b). Even if we take into account sample variability and uncertainty related to the morphological analysis, the data suggest an overall decrease in the number of macropores in the compacted samples that are related to inelastic damage. To define a more definitive trend, a systematic investigation of samples deformed to different stages of damage is warranted. Although, we have focused on Indiana limestone in this study, it is likely that the methodology developed here can be applied to other porous limestones, which we intend to pursue in future work.
In the future, such analysis could potentially be pushed further using local microtomography, where only a certain part of a sample is imaged before and after deformation. However due to the intrinsic difficulty in performing loading/ unloading cycles on a sample without damaging it, we believe that such technique would only be useful to better understand mechanical compaction if a sample would be imaged in situ. To our knowledge, no existing in situ imaging system can yet reach the stress levels at which rocks such as Indiana limestone typically deform inelastically.

APPENDIX: IMPACT OF THE MEDIAN FILTER WINDOW SIZE ON THE THRESOLDS AND PORE SIZE DISTRIBUTION
To eliminate random noise, a standard procedure in morphological analysis, used in this study, is to first apply an appropriate filter to the microCT image. Accordingly the macropore size statistics presented on intact and deformed Indiana limestone may vary with the particular filtering that was applied to the CT data. The results presented in this study were all obtained using a median filter on an 11 × 11 × 11 voxel volume. In this section, we examine the impact of the median filter size on the CT data and the pore size statistics. In Figure A1, we present the histograms for the data in Figure 4a (intact sample I0-CT) after they have been subjected to median filters of applied to subvolumes of 5 3 and 11 3 neighbor voxels, respectively. It can be seen that since the histograms are qualitatively similar, not surprisingly the thresholds determined for the unfiltered (Fig. 4a) and two filtered (Fig. A1a) images almost coincided with one another. However, there is an overall trend for the right peak (corresponding to solid) to attain higher value with increased filtering (Fig. A1b). The same trends were also observed in the deformed samples (Fig. A2). As a result, porosity is predicted to be lower in a segmented image after filtering (Fig. A3a), whereas the solid volume fraction higher (Fig. A3b). For example, for the filtered image in Figure A1a subjected to a 53 median filter the solid volume and total porosity have values of 54.15% and 14.08%, about 6% higher and 2% lower than the corresponding values in the unfiltered image (Fig. A1a). Our analysis shows that in fact both the macro and microporosities decrease in the filtered images (Fig. A3b, c). The sensitivity of pore size distribution to this type of median filtering is illustrated in Figure A4, which compares the size distributions obtained using median filters on five different voxel volumes ranging from 3 × 3 × 3 to 11 × 11 × 11. Significant differences were observed for equivalent diameters < 60 μm (comparable to the dimension of 11 voxels). Hence, it is possible that the results presented earlier have underestimated the number of these smaller macropores, but then some of these smaller features might be associated with noises in the microCT image. This is an issue that warrants more thorough investigation, but as long as identical filtering was applied to all the samples, comparison of pore sizes in undeformed and deformed samples should still provide useful insights into damage accumulation. Sensitivity of the gray level distribution to the type of median filtering for the intact sample I0-CT. a) Gray level distribution for the raw data (blue) and after median filters of 5 × 5 × 5 (green) and 11 × 11 × 11 (red); b) detail of the first peak of the distributions (macropores); c) detail of the second peak of the distributions (solid phase).  Figure A3 Influence of the median filters on the inferred values of a) total porosity, b) solid grain fraction, c) macroporosity and d) microporosity. Results on the intact sample I0-CT and the deformed samples IH, IC1 and IC2 are presented for median filters between 3 × 3 × 3 and 11 × 11 × 11. Sensitivity of the gray level distributions to the type of median filtering for samples a) IH, b) IC1 and c) IC2. The gray level distributions are presented for the raw data (blue) and after application of median filters of 5 × 5 × 5 (green) and 11 × 11 × 11 (red). Median filter 9 x 9 x 9

IH
Median filter 11 x 11 x 11

Figure A4
Sensitivity of pore size distribution to the type of median filtering for the intact sample I0-CT. The number of pores per unit volume n V is plotted versus equivalent diameter D for median filters between 3 × 3 × 3 and 11 × 11 × 11.