Digital image processing of saturation for two-phase flow in planar porous media model

Zhi DOU*, Zhi-fang ZHOU, Si WANG, Yong HUANG

School of Earth Sciences and Engineering, Hohai University, Nanjing 210098, P. R. China

Digital image processing of saturation for two-phase flow in planar porous media model

Zhi DOU*, Zhi-fang ZHOU, Si WANG, Yong HUANG

School of Earth Sciences and Engineering, Hohai University, Nanjing 210098, P. R. China

In this paper, the accuracy of estimating stained non-wetting phase saturation using digital image processing is examined, and a novel post-processing approach for calculating threshold is presented. In order to remove the effect of the background noise of images and to enhance the high-frequency component of the original image, image smoothing and image sharpening methods are introduced. Depending on the correct threshold, the image binarization processing is particularly useful for estimating stained non-wetting phase saturation. Calculated saturation data are compared with the measured saturation data during the two-phase flow experiment in an artificial steel planar porous media model. The results show that the calculated saturation data agree with the measured ones. With the help of an artificial steel planar porous media model, digital image processing is an accurate and simple method for obtaining the stained non-wetting phase saturation.

digital image processing; saturation; two-phase flow; planar porous media model

1 Introduction

Non-aqueous phase liquids (NAPLs) have been widely manufactured and used in various chemical process industries. These NAPLs include hazardous organic compounds, such as solvents, chlorinated and fluorinated compounds, and transformer oils. Accidental spills and leakage of these chemicals from various sources (e.g., underground storage tanks) into hydro-environmental systems pose a serious long-term threat to the environment (Delshad et al. 2003; Das et al. 2006; Tsakiroglou et al. 2007; Tzovolou et al. 2009).

Non-wetting phase saturation under different relative permeabilities is usually measured from flow experiments performed on porous media models (Ataie-Ashtiani et al. 2001). However, this method for measuring saturation makes the procedure of the steady-state two-phase flow experiment complex and time-consuming. Using digital image processing, the saturation can be calculated accurately, which helps to improve the accuracy of the results and avoid unnecessary experimental procedures. The correctness of saturation estimation is criticalto the creditability of relative permeability and saturation curves. Therefore, it is necessary to simply and accurately obtain the correct saturation during two-phase flow experiments.

In two-phase flow experiments, the stained non-wetting phase or wetting phase infiltrates the porous media, and the dye patterns are photographed. Then, the saturations of the non-wetting phase and wetting phase are simply obtained through image analysis of the dye photographs. This method has been proven to be very useful for estimating the non-wetting phase saturation during two-phase flow experiments in planar pore network media models (Nguyen et al. 2006). The planar pore network media model is made of a fixed media pattern, which is different from the artificial steel planar porous media model filled with standard sands with random patterns. Even though they have the same porosity, the patterns of standard sands in the artificial steel planar porous media model are different. For immiscible two-phase flow in the artificial steel planar porous media, it is very important to determine an approach to obtain the positions of the standard sand so that the non-wetting phase and wetting phase saturations can be calculated using digital image processing. Traditionally, digital image processing of dye photographs only involves separation of stained and non-stained phases. There seem to have been no reports on the calculation of saturation with digital image processing when the image information includes the non-wetting phase, wetting phase, and standard sand.

Usually, the kind of background noise has inevitable effects on the RGB color space of images (Weijer and Boomgaard 2005; Persson 2005). When the image information is distorted by transmission through a random medium, such as the air, or an imperfect communication system, a filtering method is needed to remove the background noise. The filtering method is termed an adaptive filtering technique, which must adapt itself to the distortion presented in the system so that the effects of the non-useful information can be reduced. Image binarization techniques convert an RGB color space into a bi-level space, which is the first step in most digital image processing techniques. Many binarization techniques used in processing tasks are aimed at simplifying and unifying the image data at hand (Gatos et al. 2006; Badekas and Papamarkos 2007). The simplification can enhance the effect of oncoming image processing. The quality of the binarization is critical for subsequent steps, and the correct threshold is a crucial factor for image binarization processing. However, despite the introduction of various threshold algorithms, the quality of existing methods still limits the performance of character recognition systems. These existing methods determine the threshold value based on the local properties of an image, while the common properties of the image are not found and considered.

This study used a procedure whereby the artificial steel planar porous media model was first used to measure the saturation of the non-wetting phase during the steady-state two-phase flow experiment. Second, the image smoothing processing and the image sharpening processing were used to remove the effect of the image background noise and to enhance the high-frequency component of the original image. The red band of RGB color spaces of theimage was changed to bi-level (white and black) spaces. Third, considering the common properties of the image, the threshold was calculated with post-processing. Finally, with the help of the calculated threshold, the non-wetting saturation was calculated accurately.

2 Methods and materials

2.1 Standard sand and dye materials

Standard sand with a particle size ranging from 0.5 mm to 1 mm was chosen as porous media. The standard sand was made of quartz sand with different particle sizes. The properties of the standard sand in the artificial steel planar porous media model are displayed in Table 1. Purified water was chosen as the wetting phase, and 93# gasoline was chosen as the non-wetting phase. In order to enhance the color contrast intensity of images, the non-wetting phase was stained by Sudan III.

Table 1 Properties of standard sand in artificial steel planar porous media model

In the experiment, the maximum influx rate with a syringe pump was 0.002 3 cm/s. Under this maximum influx rate, the Darcy law can be satisfied.

2.2 Experimental model

The artificial steel planar porous media model was filled with the standard sand with a particle size ranging from 0.5 mm to 1 mm, which was identical to that used in the steady-state two-phase flow experiment. The artificial steel planar porous media were installed above a two-dimensional light transmission cell. The length, width, and height of the artificial steel planar porous media model were 10 cm, 10 cm, and 1 cm, respectively (Fig. 1).

Fig. 1 Experimental device

The details of the steady-state two-phase flow experimental procedure can be found in the literature (Avraam et al. 1994; Avraam and Payatakes 1995). Initially, the porous media were fully occupied by the wetting phase. The non-wetting phase and the wetting phase were injected into the porous media at the same time and at the same influx rate using two syringepumps. Once the hydraulic heads at the inlet and outlet were steady, the non-wetting phase saturation was measured and the image was photographed. Then, the influx rate of the non-wetting phase was increased with the first syringe pump, and the influx rate of the wetting phase was reduced with the second syringe pump accordingly. When the hydraulic heads at the inlet and outlet were steady at the next influx rate, the next non-wetting phase saturation was measured and the image was photographed. Until the wetting phase saturation in the model did not decrease any more, the experiment continued. Before the non-wetting phase was injected, the first image was photographed. The measured saturation of the non-wetting phase is given by Eq. (1).

whereVais the total interior volume of the artificial steel planar porous media model, andVintandVoutare the influx and outflux volumes of the non-wetting phase, respectively.

The intensity of light has an effect on the images of color spaces. Thus, the two-phase flow experiment was performed in a dark room, and a fluorescent lamp was used. In order to make the total number of pixels in every image constant, the camera (a Canon 400D with a focal length of 20 mm) was held at the same location. During the two-phase flow experiment, 13 images were photographed, of which one image was photographed before the non-wetting phase was injected into the artificial steel planar porous media model, and was designated the standard image. The image processing steps with the help of MATLAB are shown in Fig. 2. MATLAB is a high-level technical computinglanguage and interactive environment for algorithm development, data visualization, data analysis, and numerical computation. In this study, the image smoothing processing, image sharpening processing, and image binarization processing were realized with MATLAB 7.8.

Fig. 2 Flow chart of image processing and analysis

3 Image processing and analysis

3.1 Pre-processing of image

In order to remove the background noise, we needed a filtering method that could deconvolve or equalize the distorted signal to produce the original, undistorted signal. Gaussian filtering is a useful tool for image smoothing. The Gaussian filtering function isexpressed as

whereD(u,v) is a radial basis function,uandvare two-dimensional coordinates similar toxandyin thex-ycoordinate system, andσis the scale parameter of Gaussian filtering function.

It can be seen from Fig. 3 that the background noise in circles in Fig. 3(a) is removed with Gaussian filtering. However, the edges of the non-wetting phase and the wetting phase are illegible, which inhibits the oncoming image processing. Therefore, image sharpening is needed. Image sharpening is a method that makes the digital image much more clear by compensating for outlines and sharpening edges of images. The aim of this method is to enhance the high-frequency component of the original image. The normal sharpening method performs high-frequency enhancement for the whole image. In this study, the sharpening processing was based on the Laplace operator:

wherefis the color gradation of the original gray image, andx, andyare the two-dimensional coordinates in thex-ycoordinate system. The Laplace operator is easy to use in MATLAB. In order to improve the results of sharpening processing, the enhanced gray image is obtained by

whereg(x,y) is the color gradation of the enhanced gray image, andh(x,y) is the color gradation of the gray image after pre-processing. As for different images,g(x,y)∈[?255,255], which is necessary to scale the value into the range of [0,255].

Fig. 3 Effect of image smoothing processing on background noise

3.2 Image binarization processing and selection of threshold

The images were processed with image pre-processing. Let us make the followingassumptions regarding the pixels in images: (1) there are three parts of the pixels in every image: one part containing the information of the stained non-wetting phase, one part containing the information of the wetting phase, and the other part containing the information of the porous media and standard sand; (2) the pixels containing the information of the porous media are unchanged during the two-phase flow experiment; and (3) all images have a common threshold. Under these assumptions, the threshold of every image,Ti, was introduced, and the image was processed with binarization under the introduced threshold,Ti. Since the camera was located at an unchanged height, every image had an equal number of pixels. The calculated saturation of the stained non-wetting phase in the image is given by

whereNnw(Ti) andNw(Ti) are the numbers of pixels of the non-wetting phase and wetting phase, respectively, under the same threshold,Ti. If we know the value of the thresholdTi, the non-wetting phase saturation can be calculated with Eq. (5). Considering the assumptions above, if a measuredis recorded, the threshold of all images can be calculated with Eq. (5). In two-phase flow experiments, the measuredat any steady-state moment is obtained by Eq. (1). We substituted Eq. (1) into Eq. (5), and adjusted the threshold until Eqs. (1) and (5) had the same value. The threshold was calculated with the post-processing approach.

4 Results and discussion

Fig. 4 shows the images after binarization processing. Table 2 gives the calculated and measured saturation values. In Fig. 4, the black area is the non-wetting phase, and the white area is the wetting phase and the porous media. During the two-phase flow experiment, the non-wetting phase displaced the wetting phase. The non-wetting phase saturation increased with the influx rate. Obviously, as seen from Figs. 4(a) through (e), when the wetting phase was injected at a small influx rate, the wetting phase occupied the middle of the model, due to the small porosity in the middle. The parts of the porous media and the wetting phase that were not displaced by the part of the non-wetting phase were transformed into the white area. In the white area, the part of the porous media was unchanged. Therefore, when the test area was occupied by the non-wetting phase, the wetting phase, and the porous media, there were two kinds of color, black and white, in the images.

It can be seen from Table 2 that the calculated saturation values agree with the measured ones, and the relative errors are not more than 6%, which is satisfactory for making the saturation curve during the two-phase flow experiment. This indicates that the threshold for binarization is correct and the calculated results are accurate. But we find that most of the calculated results are a little larger than the measured results in Table 2. This discrepancy between the measured and calculated saturation values resulted from the construction of the artificial steel planar porous media model. The height of the model was only 1 cm, the image information did not contain all channels, and a small amount of the channels which had nomain effects on the calculated results were hidden at the bottom of the model. From this point of view, the saturation calculated with digital image processing was allowed to substitute for the measured saturation, which was very helpful to the study of the two-phase flow.

Fig. 4 Images after binarization processing

Table 2 Measured and calculated saturation values

In the steady-state two-phase flow experiment, when the non-wetting phase was injected at a slow influx rate, which means that the influx rate of the wetting phase was faster than that of the non-wetting phase, the non-wetting phase occupied a small part of the standard sand. When the wetting phase saturation in the model could not decrease any more with the increasing influx rate of the non-wetting phase, the wetting phase saturation was irreducible (Fig. 4(l)). Table 2 shows the non-wetting phase saturation with different influx rates, calculated with digital image processing. Usually, for the steady-state two-phase flow experiment, the saturation was obtained at different influx rates. Digital image processing in this study was used to obtain the saturation during the experiment.

5 Conclusions

The present study focused on digital image processing of the saturation of the non-wetting phase in a planar porous media model during the two-phase flow experiment. The effect of background noise on the image information was also considered. Meanwhile, thebackground noise was removed with Gaussian filtering. The negative effect of the image smoothing was removed through the image sharpening which helped to improve the accuracy of the calculated results.

With the help of the post-processing approach for calculating the threshold, effective calculated results were obtained with digital image processing. The results clearly show that it is possible to obtain the saturation with digital image processing in an artificial steel planar porous media model, and the calculated results are satisfactory.

Although the calculated results are obtained during the steady-state moment, from the comparison of the measured and calculated saturations, it is confirmed that digital image processing is also useful for the unsteady-state moment at any time.

Ataie-Ashtiani, B., Hassanizadeh, S. M., Oostrom, M., Celia, M. A., and White, M. D. 2001. Effective parameters for two-phase flow in a porous medium with periodic heterogeneities.Journal of Contaminant Hydrology, 49(1-2), 87-109.

Avraam, D. G., Kolonis, G. B., Roumeliotis, T. C., Constantinides, G. N., and Payatakes, A. C. 1994. Steady-state two-phase flow through planar and nonplanar model porous media.Transport in Porous Media, 16(1), 75-101.

Avraam, D. G., and Payatakes, A. C. 1995. Generalized relative permeability coefficients during steady-state two-phase flow in porous media and correlation with the flow mechanisms.Transport in Porous Media, 20(1-2), 135-168.

Badekas, E., and Papamarkos, N. 2007. Optimal combination of document binarization techniques using a self-organizing map neural network.Engineering Applications of Artificial Intelligence, 20(1), 11-24. [doi:10.1016/j.engappai.2006.04.003]

Das, D. B., Mirzaei, M., and Widdows, N. 2006. Non-uniqueness in capillary pressure-saturation-relative permeability relationships for two-phase flow in porous media: Interplay between intensity and distribution of random micro-heterogeneities.Chemical Engineering Science, 61(20), 6786-6803. [doi:10.1016/j.ces.2006.07.028]

Delshad, M., Lenhard, R. J., Oostrom, M., and Pope, G. A. 2003. A mixed-wet hysteretic relative permeability and capillary pressure model for reservoir simulations.SPE Reservoir Evaluation and Engineering, 6(5), 328-334. [doi:10.2118/86916-PA]

Gatos, B., Pratikakis, I., and Perantonis, S. J. 2006. Adaptive degraded document image binarization.Pattern Recognition, 39(3), 317-327. [doi:10.1016/j.patcog.2005.09.010]

Nguyen, V. H., Sheppard, A. P., Knackstedt, M. A., and Pinczewski, W. V. 2006. The effect of displacement rate on imbibition relative permeability and residual saturation.Journal of Petroleum Science and Engineering, 52(1-4), 54-70. [doi:10.1016/j.petrol.2006.03.020]

Persson, M. 2005. Accurate dye tracer concentration estimations using image analysis.Soil Science Society of America, 69(4), 967-975. [doi:10.2136/sssaj2004.0186]

Tsakiroglou, C. D., Avraam, D. G., and Payatakes, A. C. 2007. Transient and steady-state relative permeabilities from two-phase flow experiments in planar pore networks.Advances in Water Resources, 30(9), 1981-1992. [doi:10.1016/j.advwatres.2007.04.002]

Tzovolou, D. N., Benoit, Y., Haeseler, F., Klint, K. E., and Tsakiroglou, C. D. 2009. Spatial distribution of jet fuel in the vadoze zone of a heterogeneous and fractured soil.Science of the Total Environment, 407(8), 3044-3054. [doi:10.1016/j.scitotenv.2009.01.020]

Weijer, J. V. D., and Boomgaard, R. V. D. 2005. Least squares and robust estimation of local image structure.International Journal of Computer Vision, 64(2-3), 143-155. [doi:10.1007/s11263-005-1840-0]

(Edited by Yan LEI)

This work was supported by the National Natural Science Foundation of China (Grant No. 51079043), the Special Fund for Public Welfare Industry of Ministry of Water Resources of China (Grants No. 200901064 and 201001020), and the Research Innovation Program for College Graduates of Jiangsu Province (Grant No. CXZZ11_0450).

*Corresponding author (e-mail:dz.uriah@gmail.com)

Received Apr. 17, 2011; accepted Jul. 4, 2011