Multiscale characterization of pore spaces using multifractals analysis of scanning electronic microscopy images of carbonates

Pore spaces heterogeneity in carbonates rocks has long been identified as an important factor impacting reservoir productivity. In this paper, we study the heterogeneity of carbonate rocks pore spaces based on the image analysis of scanning electron microscopy (SEM) data acquired at various magnifications. Sixty images of twelve carbonate samples from a reservoir in the Middle East were analyzed. First, pore spaces were extracted from SEM images using a segmentation technique based on watershed algorithm. Pores geometries revealed a multifractal behavior at various magnifications from 800x to 12 000x. In addition, the singularity spectrum provided quantitative values that describe the degree of heterogeneity in the carbonates samples. Moreover, for the majority of the analyzed samples, we found low variations (around 5 %) in the multifractal dimensions for magnifications between 1700x and 12 000x. Finally, these results demonstrate that multifractal analysis could be an appropriate tool for characterizing quantitatively the heterogeneity of carbonate pore spaces geometries. However, our findings show that magnification has an impact on multifractal dimensions, revealing the limit of applicability of multifractal descriptions for these natural structures.


Introduction
Studying carbonate reservoirs is crucial for the oil industry as they contain around half of the hydrocarbon reserves in the world (Moore, 2001). Unlike sandstones which have mainly homogeneous intergranular pores, the diagenesis of carbonate rocks combined with other geological processes create heterogeneous and complex pore spaces with sizes ranging from micro-meter to centimeters (Zhang et al., 2004). In Correspondence to: M. S. Jouini (mjouini@pi.ac.ae ) addition, these heterogeneities can be observed over several scales ranging from micro-meter to hundred of meters. Moreover, the geometry and the heterogeneity of carbonate pores can greatly impact the rock physics properties. For example, recent experiments reveal that pore spaces heterogeneity can cause variations in elastic properties values with as much as 40 % change for the same total volume of pores (Xu and Payne, 2009). Furthermore, pore spaces geometries and heterogeneities can highly impact the transport properties (Payne et al., 2010;Weger et al., 2009). These observations indicate that the relationship between pore spaces structure and rock properties can be highly complex, which makes their modeling a challenging issue. Nowadays, image acquisition and analysis techniques represent powerful nondestructive tools that make it possible to quantify the total porosity volume, morphologies and size distribution of pore spaces in porous media. This digital description can be used as input for mathematical models analysing pore distributions from images in order to provide some quantitative descriptors for heterogeneity. This characterization needs a good specimen preparation and acquisition technique to produce representative images. Scanning electron microscopy (SEM) has been widely used as imaging acquisition technique providing high resolution images of rock samples and suitable distinction between solid particles and pores (Bogner et al., 2007;Joos et al., 2011;Nadeau and Hurst, 1991). These images consist on grey level pixels, where pore spaces are represented by low levels whereas grains are represented by high level ones. In addition, the measurement of geometrical parameters requires a robust and reproducible segmentation algorithm in order to separate the void phase (pore spaces) from the solid phase (grains).
The aim of this paper is to study pore spaces distribution in SEM images of carbonates rocks using a multiscale characterization. To accomplish this, we first provide an image analysis procedure to automatically detect pore spaces (Fig. 1). Then we use the concept of fractality as a  quantitative description method of the heterogeneity of pore spaces geometry. Finally, as SEM images were acquired at several magnifications, we analyse carbonate pores spaces at different scales to assess the magnification impact on the sample porosity and the heterogeneity descriptors.
In order to solve the segmentation problem of grey level images, authors have used either statistical or spatial methods (Electron, 2004;Ramlau and Ring, 2007). The histogrambased technique (Sahoo et al., 1988;Glasbey, 1993) is one of the most frequently used approaches in statistical meth-ods and proposes the finding of suitable thresholds in order to segment grey level images, depending on histogram shapes. A simple way to obtain a binary image is to choose a suitable threshold for the grey level distribution. Heilbronner and Keulen (2006) adjusted the threshold using a subjective criterion based on a visual comparison between the input SEM image of rock sample and the segmented image result. Recently, Joos et al. (2011) proposed the use of a single value threshold computed automatically using the Otsu's algorithm (Otsu, 1979) in order to segment SEM images where histogram clearly showed a bi-modal behavior related to pore spaces phase and the solid one. Vogel (1996) suggested a more sophisticated segmentation technique based on bi-level thresholding to segment SEM images of soil where there was a partial merge between the pore spaces and grains modes in the histogram. Tsai (1995) proposed computing the maximum curvature to detect the suitable threshold in order to solve the highly complex segmentation problem of unimodal distributed histograms in a porous media. Other authors used spatial methods based on the analysis of local geometrical information in order to segment grey level images. Videla et al. (2007), Malcolm et al. (2007), and Jorgensen et al. (2010) proposed the use of a spatial method based on the watershed technique in order to segment micro tomography grey level images.
Since there is no general solution to segment grey level images, a suitable answer to this problem is to combine statistical and spatial image segmentation techniques to take advantage of both methods (Schulter et al., 2010). The proposed methodology in this paper is based on a combination of statistical and spatial methods using mainly watershed algorithm, bi-level segmentation, morphological operations and region growing techniques (Pratt, 2007;Dougherty, 1992;Henk, 1998;Soille, 1999). The main output is a binary image (black and white) where pixels are classified as pore spaces or grains. It is then possible to compute some statistics related to the distribution of pore spaces surfaces and their aspect ratio. In addition, pore spaces of carbonates rock in general have complex shapes. Thus, heterogeneity quantification can be an interesting characteristic to compute over the several magnifications SEM acquisitions.
Recently, the fractal approach was introduced to characterize complex and heterogeneous geometries. This technique is a powerful tool, proposed by Mandelbrot (1983), to describe quantitatively properties of complex morphologies. Fractal analysis is widely recognized to describe and analyze the scale variability of geometrical distributions over a wide range of scales by computing a single exponent called fractal dimension (Schertzer and Lovejoy, 1994). The technique is independent of the image source and the real scale of the features. Moreover, it has been used in many applications of sciences for describing complexity and self-similarity in nature as proposed by Evertsz and Mandelbrot (1992). Multifractal is generalization of a fractal approach applied when a single fractal dimension is not able to describe system geometry and instead computes a whole spectrum of exponents (called singularity spectrum). Several studies of fractals and multifractals related to the analysis of pore spaces distributions in porous media and sedimentary rocks can be found in earth sciences (Block et al., 1991;Grout et al., 1998;Wong and Howard, 1986;Flavio et al., 1998;Krohn, 1988;Hansen, 1988;Piñuela et al., 2007;Tarquis et al., 2003Tarquis et al., , 2006. Muller and McCauley (1992) compared multifractal spectrum of segmented SEM images with a shape factor related to pores, allowing them to distinguish soil groups associated to different soil structures. Muller (1996) demonstrated that multifractal analysis can be useful to characterize different geological chalk environments. In particular, he observed that the multifractal spectrum of pore space is correlated with air permeability values measured from the corresponding core samples. Dathe and Thullner (2005) analyzed pore spaces in porous media using, respectively, monofractals and multifractal approaches where the authors tried to establish a relationship between fractal dimensions of pore spaces and solid phase. In this study, we use the multifractality concept on SEM images of carbonate rocks in order to achieve two main goals: the first is to evaluate pore spaces heterogeneities; the second is to assess the effect of magnification on multifractals descriptors and the estimated porosities.

Image processing and segmentation
In this study, we used 60 SEM images from 12 samples of a carbonate reservoir in the Middle East taken at several depths and locations. Five (5) SEM images of each sample were acquired at different magnifications by focusing at each step on either the central part or a random location of the sample. The size of images was 1024 × 1024 pixels; magnifications used for acquisitions ranged between 800x and 12 000x. Pore spaces were detected using the segmentation workflow (explained presently) and porosities were estimated at every magnification scale. We developed the image segmentation algorithm and the multifractal analysis codes using Matlab scripts. In order to segment grey level of SEM images, we used an image processing technique based on applying the watershed algorithm as a main procedure combined at different steps with statistical and other spatial methods as described by the whole workflow in Fig. 2. Watershed transformation has been widely used in image segmentation problems (Rambabu and Chakrabarti, 2007;Bieniek and Moga, 2000;Vincent and Soille, 1991). In this method, a grey-level image may be considered as a topographic relief where pixel value represents a height. Intuitively, watershed of a relief corresponds to the limits of the adjacent catchment basins of water drops. The main idea of this technique is to find local minimums which will be sources of a continuous water rising flood. Then in order to prevent merging water coming from different sources, a barrier is built at each point of contact. The process ends when the maximum grey level is reached and the union of all those barriers composes the watersheds (Fig. 3). Unfortunately, this transformation creates an over segmentation of the image in most cases. In order to prevent this problem, we used an efficient strategy based on using controlled markers (Meyer and Beucher, 1990). These markers aim to help the segmentation process by detecting some zones in pore spaces and in grains which will represent sources in addition to local minima of the gradient of the image used in the watershed algorithm itself. Markers corresponding to pore spaces were detected using a bi-level segmentation and a region growing technique. Markers related to grains were obtained using morphological operations (Fig. 4). The bi-level segmentation technique combined with a region growing is a useful segmentation tool for grey level data with histogram partial merge between two main phases. This technique aims in a first step to find two threshold levels G min and G max . Hence, in our images segmentation of all grey levels below the lowest threshold G min were detected as pores and all grey levels above G max were detected as grains. These thresholds were computed using gradient masks. Sobel and the Laplace filters were used to compute histograms of grey levels corresponding to detected edges and deduce from them these two thresholds. The second step consisted of using a growing region algorithm to segment the fuzzy interval zone (grey levels between G min and G max ) using connectivity between pixels belonging to pore spaces. Furthermore, morphological operations were applied on the 60 SEM images (from 12 carbonate samples) in order to detect grains markers used as sources for the watershed algorithm. The detected markers represent connected blobs of pixels inside each of the grains zone (Fig. 4c). A series of erosion, opening by reconstruction and closing by reconstruction were used to create flat maxima inside each object. Finally, pixels connectivity was used to extract these markers. The final segmentation result of the workflow is a binary image where black and white pixels denote respectively grains and pore spaces (Fig. 4e). This result depends on the quality of the detected markers of pore spaces and grains. The next step was to give labels to pore spaces by using connectivity algorithms which allows extracting and evaluating surfaces and aspect ratio of every pore space. Consequently, the final result of this images processing is given by the total pore space, porosity and the pore aspect ratio distribution.

Multifractal model
As part of our images analysis, we also studied pore spaces morphologies and heterogeneities based on multifractals. Unlike mathematical definition of ideal fractals, pore spaces structures reveal fractal properties only in a statistical sense, which means we need a statistically representative number of samples to be analyzed. In order to compute a fractal dimension, it is necessary to define a measure in the digital images which is closely associated with pore spaces geometries in images. A widely used measure to study spatial structures is based on the box-counting method (Posadas et al., 2001). The principle is to cover the binary image by a regular square grid partitioning the space into boxes of size ε. This process is repeated for different box sizes, which is equivalent to analyzing the studied geometry or structure at different scales. The main equation for fractal theory establishes a relation- ship between the number of boxes N(ε) needed to cover the analyzed structure and their size ε: where D 0 is the fractal dimension: The generalization of fractal approach is the multifractal analysis, and is applied when a single fractal dimension is not able to describe system geometry. These objects can be more completely characterized by a spectrum of fractal dimensions. Multifractals are measured using a probability distribution P i in each ith box as: where ε is the box size and α i is the Lipschitz-Hölder exponent characterizing the singularity strength in the ith box. The factor α i allows quantifying the distribution complexity in spatial location (Posadas et al., 2005). The number of boxes N(α), where the probability P i has singularity strengths between α and α+dα, can be related to the box size ε as: where f(α) is the singularity spectrum of boxes characterized with singularity α. For mono-fractal pore distributions, α remains constant and the multifractal spectrum is composed of a single point. For multifractal pore distributions, the spectrum can be represented by a curve with a wide range of values for α. This interval increases with the increase of the distribution heterogeneity (Posadas et al., 2005;Xie et al., 2010). Multifractals generalized dimensions D q of the q −th order (Hentschel and Procaccia, 1983) are defined as: where µ(q,ε) is the partition function (Chhabra and Jensen, 1989): Moreover, this partition function scales with box size ε as: where τ (q) is the q th mass exponent and f(α) the singularity spectrum (Halsey et al., 1986), are related as: D 0 , D 1 and D 2 can be computed using the following equations (Posadas et al., 2005): where C(ε) represents the correlation sum. A first analogy with physical properties appears as D 1 and corresponds to the information dimension (Shannon and Weaver, 1949). D 0 is called capacity dimension and provides an average related to the information of the analyzed structure distribution (Voss, 1988). D 2 is the correlation dimension. Equality for these three dimensions occurs only when the distribution is ideally monofractal (Korvin, 1992).

Image analysis
Porosity was computed as a percentage of the segmented pore area to the total image size for the analyzed samples. The segmentation process provides a binary image where black pixels denote the detected grains and the white pixels denote pore spaces. In addition, as no ground-truth information was available, the resulted detected pores were analyzed visually in order to validate the accuracy of the segmentation process. Figure 5 provides an example of one of the SEM images (sample S 1 ) and the associated segmentation result. The analyzed zone is a square of 45 µm × 45 µm, composed mainly of calcite and revealing a homogeneous distribution of grains and pores. Detected pore spaces in images are in agreement with the expected result based on visual observation. Moreover, these pores are regular and homogeneous and are mainly micro-pores and meso-pores, with a majority Table 1. Estimation of capacity dimension D 0 , information dimension D 1 and correlation dimension D 2 for samples S 1 to S 12 . R 2 is the coefficient of determination computed for the linear interpolation in the loglog plot.   of pore sizes ranging between 0.1 to 1 µm 2 (Fig. 6a). Furthermore, we provide the aspect ratio distribution for the same sample in Fig. 6b. This statistic refers to the ratio between the minor and major axes of an ellipsoidal pore. The aspect ratio distribution for sample S 1 shows a maximum close to 1. However, this value is not representative of the whole sample but is generated by the large number of detected micro-pores with sizes close to the image resolution limit.

Determination of multifractals parameters
The choice of suitable box sizes and range of moments orders are crucial for the multifractal analysis (Saucier and Muller, 1999). In order to find these ranges, we have used two criteria based on assessing the linear behaviors of the mass exponent τ (q) as a function of q and the partition function µ(q,ε) as a function of ε (Dathe et al., 2006). First, we assessed the variation of the mass exponent τ (q) for several ranges of moments q. We observed that τ (q) behaves nonlinearly when considering an interval of moments wider than q [−5, +5].
Thus, we chose to restrict our study to this interval where the mass exponent τ (q) versus q is showing a linear relationship with a coefficient of determination R 2 above 0.94 (Posadas et al., 2003). Moreover, we found that a moment step equal to 0.5 allows showing accurately the multifractal behavior in the range q [−5, +5]. In addition, we applied the boxcounting method to the analyzed images (1024 × 1024 pixels) using box sizes ranging from 1 to 1024 and assessed the partition function µ(q,ε) as function of box sizes. We found a linear behavior relationship with coefficient of determination R 2 above 0.9 for box sizes range ε [8, 512] and moments range q [−5, +5]. These ranges were suitable for the majority of analyzed images. Values of multifractal coefficients were determined first by computing slopes of logarithmic values of the partition function over logarithmic values box size as presented in Eq. (5), and second by using the relation presented in Eq. (8). Figure 7 shows the linear behavior of the partition function µ(q,ε) for q = 0 and q = 1 of samples S 1 and S 2 . Table 1 provides the capacity dimension D 0 , the information dimension D 1 , the correlation dimension D 2 ,  image porosities, and the respective coefficients of determination R 2 of the partition function for the twelve analyzed samples for the magnification 12 000x. Furthermore, in order to validate the choices made for the moments and the box sizes ranges, we verified that the singularity spectrum f(α) is a concave parabola for each analyzed image showing a multifractal behavior. In addition, we also verified that the spectrum touches the internal bisector (f(α) = α) of the axis (Riedi, 1997). The point of intersection between the tangent and the graph of f(α) corresponds to f(α(1)) = α(1) = D 1 (Evertsz and Mandelbrot, 1992). Moreover, we have investigated the possibility of existence of a relationship between the variation of the total porosity and multifractal dimensions. Figure 8 reveals a direct linear fitting relationship with a coefficient of determination R 2 equal to 0.67 between porosity and the capacity dimension D 0 . This observation suggests that high porosities correspond to high capacity dimensions. This is an expected result, as for high porosity images the pore spaces structure is more compact, which increases the capacity dimension. In addition, the same increasing behavior is observed in Table 1 for the information dimension D 1 and the correlation dimension D 2 . Neverthe-less, the generalized dimensions D q for q ≥ 2 do not show a clear tendency that can be correlated to the total porosity. In relatively homogeneous samples, the generalized dimension D q showed a slow rate decrease leading to a convergence toward a constant value for moments q≥0, as illustrated in Fig. 9 for sample S 1 . In this example, generalized dimensions D q are constant for q≥2. This behavior was also observed for homogeneous samples S 4 , S 10 and S 11 . Moreover, the generalized dimensions D q for the heterogeneous samples S 2 , S 3 , S 5 , S 6 , S 7 , S 8 , S 9 and S 12 revealed a continuous decrease for 0≤q≤5. These two observations corroborate that generalized dimensions D q can provide an indicator of data heterogeneity degree, since their values would converge to a constant faster in images where pore distributions were less heterogeneous. Furthermore, Posadas et al. (2005) have shown that homogeneous multifractal distributions have a narrow concave f(α)-spectra, whereas the opposite is true for heterogeneous structures. Thus, the wider is the magnitude of change of f(α) around the value f(α(0)), the higher is the heterogeneity of analyzed structures. A visual interpretation of our segmented images qualitatively reveals more heterogeneity in samples with wider f(α)-spectra, which corroborates previous studies. Figures 10 and 11a provides the segmented images of samples S 1 , S 2 , and S 3 . Figure 11b and c shows respectively their f(α)-spectra and their mass exponents τ (q) with coefficients of determinations, respectively, 0.998, 0.995 and 0.999. For instance, sample S 1 f(α)-spectra showed a narrower shape compared to samples S 2 and S 3 , with a computed interval of α [1.48, 1.87] and α = α maxα min = 0.39. Whereas, the singularity spectra revealed higher magnitudes of alfa changes of α [1.64, 2.17] and α = 0.53 for sample S 2 , and α [1.72, 2.24] and α= 0.51 for sample S 3 .

Scale effect
The magnification at which images are investigated in carbonate rocks may affect the physical interpretation, in terms of structures and distributions of the visualized pore spaces. Patterns that can be distinguished at a certain image size and magnification may be not present at others. Zhi-bin Liu et al. (2005) studied the magnification effects on the interpretation of SEM images of expansive soils and provided a range of validity for fractal dimensions variations depending on image magnification. Therefore, we studied the magnification effect on SEM images in order to evaluate the variations in multifractal dimensions and the total porosity.

Scale effect on multifractal dimensions
Dimensions were computed for the 12 samples at magnifications ranging from 800x to 12 000x. We found that, for most of analyzed samples, the variation of multifractal dimensions was relatively low, around 5 % in intermediate magnifications (1700x to 12 000x). Figure 12 shows    Figure 13. Results of multifractal analysis for sample S1. Plots of multifractal dimension Dq 2 versus moment q, qє [-5, 5] the sample S 1 acquired at five magnifications, and Fig. 13 provides its associated multifractal dimensions. The average variation of multifractal dimensions between magnifications M 1 (1700x) and M 4 (12 000x) was equal to 5.76 % for q≥0 for sample S 1 . Table 2 recapitulates these results for all samples and provides the variation of alpha ( α) for all images having multifractal behavior. In Fig. 14, we provide the multifractal dimensions results for all the other analyzed samples at these five magnifications. For the majority of samples, the lowest magnification M 5 (800x) showed significant variations in multifractal dimensions between 2.3 % and 13.95 % for the worst case when compared to the other scales. This magnification corresponds to the highest scale of visualization and represents the lowest magnification from which detected patterns were most likely related to macro-pores and megapores with an equivalent pore radius ranging, respectively, from 2 to 10 µm and from 10 to 100 µm. Furthermore, we found that no clear multifractal behavior could be detected for three samples at the highest magnification S 3 , S 3 and S 8 as their respective spectra f(α) failed to show a concave shape (Riedi, 1997). The main reason is that the detected patterns at this magnification showed very irregular pores, which are unrepresentative of the general pore spaces distribution observed for the other magnifications.
These findings show that magnification has an impact on multifractal dimensions, revealing the limit of applicability of multifractal descriptions for natural structures, which seem to have multifractal behavior for a limited range of observations scales. Moreover, it corroborates the study of Zhibin Liu et al. (2005). In our case, the ratio between magnifications providing significant estimations for multifractal dimensions was 0.14 (ratio between magnifications 1700x and 12 000x). This result indicates that future studies should be orientated to investigate self similarity for higher scales and its limits.

Scale effect on porosity
We also studied the scale effect on porosity estimation from the segmented SEM images. We found that porosity had an increasing behavior with image magnifications for nine samples and a decreasing behavior for the other three. Figure 15 shows the porosity estimations for samples S 1 , S 2 and S 3 for the five magnifications. First, in terms of image analysis, the decreasing behavior was expected for samples composed mainly of micro-pores and meso-pores, with an equivalent pore radius ranging, respectively, from 0.1 to 0.5 µm and from 0.5 to 2 µm, like samples S 1 and S 2 . Indeed, for lower magnifications, micro-pores and meso-pores were less visible and therefore less detectable by the segmentation algorithm. Second, the increasing behavior of the other samples was related to the detection of macro-pores and megapores. Sample S 3 , for example, revealed at the lowest magnifications some pores with an equivalent radius larger than 12 µm. Finally, we investigated the possibility of a correlation between the variation of multifractal dimensions for the five magnifications and the variation of porosities for these magnifications, but unfortunately we did not find any clear correlation between them.

Conclusions
The image analysis techniques presented in this paper allow studying pores morphologies using a grey level image as input data without any a priori required for images source. The proposed study showed how useful the multifractal analysis could be when investigating spatial features in highly complex geometry images. This paper indicates that the pore space of the analyzed carbonates samples is multifractal and their heterogeneity can be quantified through the variation of the singularity α. In addition, multifractality behavior was detected over a limited range of magnifications (magnifications between 1700x and 12 000x). Thus, we could provide a range of scales for the multifractal modeling validity. We also focused on the scale effect on the porosity estimation from the segmented SEM images. We found that porosity showed an increasing behavior with image magnification, increasing for nine samples and a decreasing for the other four samples. However, no clear correlation could be found between the variations of multifractal coefficients and the variations of porosities due to the magnifications effect. Finally, future acquisitions of thin sections and X-ray computed tomography scans for the same samples will be conducted in order to study the behavior of detected carbonates pore spaces at higher scales and to compare them to our actual results.