X-ray CT analysis of pore structure in sand

The development of microfocused X-ray computed tomography (CT) devices enables digital imaging analysis at the pore scale. The applications of these devices are diverse in soil mechanics, geotechnical and geoenvironmental engineering, petroleum engineering, and agricultural engineering. In particular, the imaging of the pore space in porous media has contributed to numerical simulations for single-phase and multiphase flows or contaminant transport through the pore structure as three-dimensional image data. These obtained results are affected by the pore diameter; therefore, it is necessary to verify the image preprocessing for the image analysis and to validate the pore diameters obtained from the CT image data. Moreover, it is meaningful to produce the physical parameters in a representative element volume (REV) and significant to define the dimension of the REV. This paper describes the underlying method of image processing and analysis and discusses the physical properties of Toyoura sand for the verification of the image analysis based on the definition of the REV. On the basis of the obtained verification results, a pore-diameter analysis can be conducted and validated by a comparison with the experimental work and image analysis. The pore diameter is deduced from Young–Laplace’s law and a water retention test for the drainage process. The results from previous study and perforated-pore diameter originally proposed in this study, called the voxel-percolation method (VPM), are compared in this paper. In addition, the limitations of the REV, the definition of the pore diameter, and the effectiveness of the VPM for an assessment of the pore diameter are discussed.


Introduction
The estimation of pore dimensions and pore networks in soil is one of the most important studies to evaluate the mechanical and hydrodynamic properties for soil science (Topp and Miller, 1966), soil mechanics (Carman, 1939;Mualem, 1976;Dullien, 1992;Chanpus, 2004), geotechnical and geoenvironmental engineering, and petroleum engineering (Chatzis et al., 1983;Culligan et al., 2006;Gharbi and Blunt, 2012).Carman (1939) considered that water did not seep into straight channels but around irregularly shaped solid particles.Topp and Miller (1966) reported on the waterpressure changes, which shifted the main-branch wetting and drying curves, as well as the families of "rewet" and "redry" scanning loops in the water retention curve (WRC) of glassbead media.Mualem (1976) proposed a simple model that predicted the unsaturated hydraulic conductivity curves by using WRCs.Chatzis et al. (1983) evaluated the effects of the particle size, the particle-size distribution, the macroscopic particle size, the macroscopic and microscopic heterogeneities, the microscopic dimensions such as the ratio of the pore body to the pore-throat size, and the pore-to-pore coordination number on the residual oil under water-wet conditions.In fact, it is difficult to define the pore dimensions in grains because the pores are surrounded by grains and are thus not isolated.The shape of a soil particle is not spherical but rather a complicated shape; thus, the pore dimensions are only defined on the basis of some assumptions.In the past decade, microfocused X-ray computed tomography (CT) scanners have been used to study the pore structure and residual fluid in porous media.Culligan et al. (2006)  an X-ray CT scanner to visualize an oil-water-glass-bead system.Gharbi and Blunt (2012) reported that the efficiency of water flooding as an oil recovery process in carbonates should depend on the wettability and connectivity of the pore space.Figure 1 shows an X-ray CT image of a grain sample in two dimensions.The X-ray CT shows the spatial distribution of density, which enables the soil particles and pores to be distinguished.Locally, the longest and shortest length of the pore, as shown in Fig. 1a, can be measured by using software for image analysis.However, it is partial property of the pore, and the required information is at least a property in representative volume.It should be required not to measure individual pore size by using software, but to estimate them by using a systematic method.Moreover, complicated pores have an aspect ratio, as shown in Fig. 1; thus, a discussion of connectivity of pores will be required for the study of the hydrodynamic issue in soils.The aim of this paper is to propose an evaluation method for the pore size.We now review the current techniques for pore analysis.The most popular methods for measuring pores in soil are the indirect methods of a mercury intrusion porosimetry (MIP) or the air intrusion method (AIM).These methods are based on the concept that the pore structure is assumed to be a straight tube.Through the recent developments of the scanning electron microscope (SEM) and nondestructive testing methods such as CT and magnetic resonance imaging (MRI), the pore structure in soil can be directly measured (Tomotsune et al., 2015).The application of CT to soil mechanics with respect to the localization of the soil under triaxial compression was reported by Desrues et al. (1996).Otani et al. (2000) evaluated the crack width of undisturbed clay using an industrial X-ray CT scanner.Likewise, the applications of CT to geotechnical engineering have been diverse over the past two decades.In particular, advanced CT has been developed to scan at the microscale.For example, Altman et al. (2005) evaluated synchrotron-source X-ray CT for the visualization of flow phenomena in rock samples.Wildenschild et al. (2002;2005a, b) focused on the flow process in the pores between sand particles using a microfocused X-ray CT scanner and successfully visualized in granular materials.Higo et al. (2011) studied the deformation of unsaturated sand using X-ray micro-CT scanner and discussed the pore size and the meniscus in a pore.As for a discussion of a capillary in soil, Andrew et al. (2014Andrew et al. ( , 2015) ) developed a system that conducts a CO 2 flow test and acquires images of rock grains and then discussed the snap-off phenomena in the pore structure.Taylor et al. (2015) have a similar motivation as our research; they visualized the pores between sand particles and originally developed a void segmentation technique.Yang et al. (2015) focused on the issue of geologic carbon storage (GCS) and presented the direct extraction of the pore-size distribution in a porous structure without an approximation or a complex calculation.Further, Yang et al. (2015) tried to predict the WRC on the basis of a comparison of measurements obtained by MIP.Andrew et al. (2014Andrew et al. ( , 2015) ) and Taylor et al. (2015) strictly evaluated the pore size; however, a study on the relationship between the capillary phenomena and the pore size and its distribution in three dimensions is still under discussion.Most likely, it is significant but difficult to define the pore size in three dimensions.To evaluate a structure consisting of a large number of pores, a suitable image analysis method is required.
This paper discusses an evaluation method for pore structure of sand from microfocused CT scan data.In this paper, the authors distinguish a pore from the pore structure.In the first part of this paper, the authors propose the application of a mathematical morphology method for estimating the pores between sand particles.By showing the analysis results of simple subjects, the usefulness of proposed method will be validated.Next, the importance of the selection of a representative element volume (REV) for estimating the grain-size distribution and the averaged pore index, such as the porosity and specific surface, is discussed.The authors show there is an optimum REV in this paper.On the basis of the above fundamental examination of treating CT data, the authors propose a voxel-percolation method (VPM) to evaluate the pore structure of sand.The estimated results are compared with a WRC test, and the effectiveness of proposed method is described.It is concluded that the resolution required for evaluating the pore structure is almost equivalent to that for a pore.The final objective of this research is to develop a general method for soils.As the primary aim of this research, the evaluation of sand is treated in this study because it is natural material.

X-ray CT scan
Table 1 lists the specifications of the microfocused X-ray CT scanner (TOSHIBA TOSCANNER 32300 FPD) installed at the X-Earth Center at Kumamoto University in 2010.The generated X-ray beam is a polychromatic beam with a wide range of frequencies; hence, corrections are made for the beam-hardening effect.Radioscopic images alone are not Table 2 lists the scan conditions selected for this study.An X-ray tube voltage of 60 kV and a current of 200 µA were chosen, and the focus-to-center distance (FCD) was defined as 24.7 mm; hence, the dimensions of one voxel are 5 × 5 × 5 µm 3 in this study.In general, the scan conditions depend on the target of study, and the voltage, current, and FCD should be variable.The scan conditions selected in this study are not absolute conditions, and each user has to find the best conditions to observe the target in each study.The soil tested was Toyoura sand with a dry density of 1.57 t m −3 and specific density of 2.65 t m −3 (i.e., the porosity was 0.41).Toyoura sand has been widely used as benchmark sand for civil engineering in Japan, and it is a quartz-dominated sand (Miura et al., 1989;Ojuri and Fijabi, 2012).The range of grain sizes is 75-850 µm, and the uniform coefficient, defined as D 60 /D 10 , is 1.5; hence, Toyoura sand is categorized as a classified soil with low range of grain sizes.The sample was well-packed into an acrylic mold with a diameter of 10 mm and a thickness of 1 mm.If the diameter of the specimen is greater than 10 mm, e.g., 20-30 mm, the center of the specimen is too far from the X-ray tube, and high-resolution CT images cannot be obtained because of the loss of focus between the X-ray tube and the center of the specimen.In this study, the authors did not use any special region of interest (ROI) and only used Image J to extract the cubic area of interest from the center of the image.There may be some methods for processing an image of a pore in sand from X-ray CT data.In this paper, the applicability of the mathematical morphology is discussed for the definition of the pore structure.The mathematical morphology describes the complicated shape by using a spherical element whose diameter is known (Soille, 2003).The basic operations are available in many image analysis software packages (Luis et al., 2005).
The basic process for image analysis in this study consists of image segmentation to create a binary image from an original image (Sect.3.1.1)and the determination of the pore diameter using a granulometric method based on the mathematical morphology (Sect.3.1.2).

Image segmentation
The most popular method to evaluate the porosity of porous materials from CT data is based on a statistical assumption.The parameters of a distribution function based on the CT data can be determined by an optimization technique (e.g., Kato et al., 2013;Mukunoki et al., 2014).The accuracy of these methods depends on the selection of distribution functions based on the CT data.The required number and type of functions are still under research, and there may be a number of solutions to these issues.As the first step of image analysis in this study, the pore and grain are identified from CT data and Fig. S1 in the Supplement shows a histogram of a CT image of Toyoura sand.The image segmentation method developed by Otsu (1979) was used because the Toyoura sand tested exhibited two distinct peaks, as shown in Fig. S1.The two peaks in the CT values indicate the two phases of the particles and pore space.During this process, the image is analyzed so that the grains are interpreted as black areas and the pores are white areas.By determining the number of voxels for each region, the averaged index associated with the pore, such as the void ratio, can be evaluated.However, by this process, the distribution of the pores according to size cannot be evaluated.

Granulometric method
In the second step of the image analysis, for each pore identified from the binary data, its shape is evaluated by using the mathematical morphology.On the basis of this concept, a unit element B belongs to a pore X.Then, the unit element in X, B x , can be written as follows: (1) It is likely that the shape of the unit element is a square or circle with a symmetric shape.In this study, the target subject is a pore, and we are interested in its size.By obtaining an image of a pore from an X-ray CT scanner, the dimensions of the pore can be determined by using the unit element because the dimensions of the unit element can be modified.In this study, the unit element of a sphere is used; therefore, the element is called a spherical element in this paper.Figure S2 illustrates the spherical element.In this study, 13 spherical element sizes are used.The first three spherical elements for the granulometric method are shown in Fig. S2.The maximum radius r of the spherical element in this study is r = 13.At the beginning of the granulometric image analysis (GIA) process, the maximum size of a spherical element is not known; therefore, only the appropriate number of radii of the spherical element is given.The GIA continues that the radius of the spherical element is increased until there are no pores that fit that element.The smallest element is one voxel (Fig. S2a).Elements with a small radius do not appear exactly spherical in shape.However, an element becomes more spherical as the number of voxels is increased for a specific diameter of that structural element, as shown in Fig. 5.The voxel number corresponding to the diameter of a spherical element, D, can be defined as follows: where r is the number of voxels from the center of the spherical element.As shown in Fig. S2, when r is 1, which is the minimum number, d is 3; therefore, the diameter of a spherical element is 3.The diameter increased on the basis of the center inclusion concept; thus, the diameter of the spherical element should be always an odd number.As r increases further, the shape of the spherical structure becomes more spherical.It is proposed that the diameter of a spherical element is the pore size at the corresponding location in this study.For the generation of a spherical element with different dimensions, an image tool kit called ITK (Luis et al., 2005) was used.
The granulometric method can recreate the complicated pore space by overlapping spherical elements having several different diameters.Initially, the smallest spherical element, i.e., one voxel, should be applied to the analysis area, which is a pore space, and it should occupy the entire pore space.Then, the next larger spherical element with a diameter of three voxels, as shown in Fig. S2b, is applied to the same area of pore space, but the next larger spherical element cannot occupy the entire pore space, as shown in Fig. S2c.Likewise, the entire pore space is scanned by each spherical element.For a more complicated pore shape, more spherical elements having different diameters will be required.In the granulometric method, a spherical element may be partially overlapped; hence, the distribution of the pore diameter indicates that the nonoverlapping part of a spherical element is evaluated.In addition, the sum of the ratios of overlap of spherical elements at each step can be used to evaluate the pore volume and degree of saturation by the VPM, as explained in the following section.

Verification of area measurement using granulometric method
It is noted that a spherical element is referred to as a circular element in this section because the target subject is drawn in two dimensions.At the end of the image analysis using the GIA, the number of voxels in each circular element can be counted.The accuracy of the area measurement by the GIA is verified by comparing the number of counted voxels with the calculated areas of a circle and square.Figure 2a and b show a black circle and a black square rotated 45 • in an area with dimensions of 200 × 200 voxels (i.e., the total number of voxels is 40 000), and Fig. 2c shows the circle and square drawn in an area with dimensions of 400 × 200 voxels (i.e., the total number of voxels is 80 000).The space other than the circle and square is defined as the pore space (see Fig. 2a and b).The color legend indicates the diameter of the circular element: white corresponds to the largest diameter of the circular element.As the diameter is reduced, the color becomes dark.The number of voxels for the smallest circular element is equal to the number of pore structures in total.Table 3 summarizes the area of the pore space obtained from the GIA and the areas calculated for the solid and pore spaces by the following set equation: where S is the number of voxels as the total area counted by GIA, and B(r i ) is the area of the circle/spherical element with a radius from i to n.When i is equal to n, S is equal to the target area/volume.Even if the shape of the pore space is complicated, as in Fig. 2c, the GIA can estimate the total number of voxels (see Table 3).
If the diameter of a spherical element is equivalent to the pore diameter, its pore geometry must be a set of parallel lines or a rectangular shape.The area enclosed in dotted lines, as shown in Fig. 2d and e, was analyzed by the GIA, and both areas were found to be 19 801 voxels.The two artificial images in Fig. 2d and e can provide an interesting discussion.Certainly, the GIA estimates the width in the area of the diamond; despite the fact that the definition of width is vague, GIA produced Fig. 2d.Even if the target image is rotated by 45 • , the same results should be obtained.The results in Fig. 2d and e indicate that the GIA provides the same results, even though the shape of the pore space changed; hence, the GIA can scan the complex pore structure to investigate the pore size.The point is that the spherical element evaluates the space at the four corners, as shown in Fig. 2e.This result indicates the important feature of this image analysis; in short, the spherical element finds small pore spaces and evaluates the pore size by the part remaining after the overlapping process.This behavior, whereby the spherical element finds the small space, is similar to the capillary behavior in porous materials such as sand.Note that the labelled number is given to each circular element and its number is corresponding to the radius of circular element.Besides the location of each circular element is also recoded; hence, the spatial distribution of the overlapping circular element can be visualized as shown in Fig. 2. The GIA is very useful technique because the pore size can be evaluated without measuring the individual pores.However, it should be recognized that the diameter of the circular element is not strictly equal to the pore diame- ter.The authors suggest a diameter for the spherical elements fitting the pore space as the pore size in sand; this will be discussed later in Sect. 4.

Selection of appropriate representative element volume (REV) for analysis of sand
The REV of CT image should be discussed to estimate the geometrical properties of sand by an image analysis.In this study, the porosity and specific surface area of Toyoura sand were evaluated, and the REV was assessed.The detailed steps can be found in Fujiki et al. (2014); hence, the concept is only introduced.A subsampling region is defined; the porosity (n) and specific surface area (S sp ) are calculated from the CT image using the following equations: where V pore is the volume of the pores, V t is the total volume of the specimen, and S ps is the mean surface area of the grains obtained from the CT image.In addition, the means and standard deviations of the porosity and specific surface area are calculated.The above described process is continuously repeated owing to the enlargement in the calculation region until the relative standard deviation (RSD), which is defined as the standard deviation divided by the mean, is less than 1 %.
The dimension of the CT image is 1024 × 1024; however, if the extraction of the cubic area from this CT image is required, the maximum voxel number to be used is 700.The image can be obtained from binary images following the image segmentation process.Figure 3 shows three-dimensional images; the porosity and specific surface area were evaluated to validate the representative volume.The measured porosity of the specimen was 0.41. Figure 4a and b show the analyzed porosity and specific surface area with along with their  RSDs for various subsample sizes.In this analysis, the domain size of the subsample size and the location for the initial calculation of the porosity and specific surface were randomly changed 20 times; then, the RSD for each voxel size was calculated.The porosity and specific surface analyzed from the CT images converged owing to the increase in the voxel size; moreover, a tendency to decrease is observed for both parameters.The behavior observed for porosity results in close measured value.The difference between the porosity and the specific surface is not significant for the tested materials but, generally, the REV should depend on the properties of the medium.Figure 5 shows the distributions of the grain size obtained from a sieving test and an image analysis.The grain-size distribution can be obtained from Image J for the image analysis, which is provided using the function "object counter".With the exception that the image analysis area is a cube consisting of 100 voxels, the grain-size distribution obtained from the CT image fits the results of the sieving test well.The results are summarized in Table 4.The effect of the size of the reference sample is not significant.This observation supports the fact that the sample with a voxel size of 300 (or 1.5 mm) is larger than the REV when the size of one voxel is 5 µm.
4 Analysis of WRT using the proposed method

Water retention test for sand
In this study, a water retention test with a reducing elevation head method (WRT-REHM) was selected to conduct water drainage tests because it was available to measure the moisture content of identical specimens at different elevation heads during the water drainage process.The specimen used for the WRT-REHM was identical to the scanned sample.Figure 6a and b show photographs of the setup used for the water drainage test system with a suction method, and Fig. 6c shows a cross-sectional view of a mold that was tested.The mold is made of acrylic, through which an X-ray beam could be transmitted without strong beam-hardening.The dimensions of the mold were a height of 120 mm, an inner diameter of 10 mm, and a thickness of 1 mm.In order to measure the amount of drained water, a glass syringe was used with a scale of 0.01 mL.A membrane filter with a mean pore diameter of 0.2 µm was placed on a glass filter with a pore diameter between 20 and 30 µm installed on the bottom of the mold.The mold was filled with de-aired water, and an entire system with a syringe, a tube connected between the mold and syringe, a glass filter, and a membrane sheet was fully saturated in the storage mold under vacuum conditions.Toyoura sand with a porosity of 0.43 was carefully installed in the mold filled with de-aired water.The sand specimen was left for 24 h after regulation of the elevation head of the syringe to lead water drainage at each degree of saturation; then, the volume of the syringe was recorded.
The entire test system was set up in a room containing the installed microfocused X-ray CT scanner, as shown in Fig. 6d, and the temperature was controlled at 20 ± 2 • C. Owing to the change in temperature within the range of 20 ± 2 • C, the authors were concerned about the generation of condensation in the mold and extra evaporation from the specimen surface; thus, the humidity in the room was regulated.In the trial test, the authors monitored the time to reach a steady-state condition of the specimen by checking of the fluid volume in the syringe, and it was concluded that this time should be 24 h in this test.Eventually, the authors performed WRT twice in this study.

Pore structure analysis of sample for WRT
In the process described in Sect.3.2, the GIA produces a pore by overlapping many spherical elements.In this section, a pore-structure analysis based on the GIA is described.It is important to consider the three-dimensional continuity of the pore when analyzing the pore structure.In this study, a pore-structure analysis method that performs a vertical airentry simulation with the imaged pore from the X-ray CT data is proposed.This method is called the VPM.For instance, there are 13 types of spherical elements, as shown in Fig. 7.For convenience, the images in Fig. 7 are drawn in two dimensions; therefore, sphere elements should be called circle elements in the explanation of Fig. 7.In order to start the percolation flow simulation, the water-drainage process is modeled as follows in this study.First, the original image is binarized into pore space (white) and soil particles (black; see Fig. 7a).Second, the pore space is analyzed by the GIA; thus, the distribution of the labeled spherical elements can be determined (Fig. 7b).Third, the VPM finds the labeled voxel of the circular element with the largest radius from the corner of the defined side.As in Fig. 7c, only the area of the spherical element with a radius of 13 is shown as white.Fourth,VPM finds the voxels labeled 13; then, it will keep painting those voxels until it recognizes no continuous circular element (Fig. 7c-o).Likewise, the labeled number of the largest circular element is scanned to the image obtained in the second step; thus, only the voxels corresponding to the largest circular element are counted.The third step should be repeated until the smallest circular element is used.Fifth, the results in the fourth step produce the degree of saturation by the sum of the counted voxels divided by the number of voxels in the entire pore.Lastly, the capillary pressure can be evaluated by using Young-Laplace's equation with  the diameter of the circular element.In an actual analysis, the described process in the above can be processed in three dimensions.

Analyzed water retention curve
The water retention property of a soil is a typical parameter influenced by pore structure.In this section, the WRC can be reproduced by combining the GIA and VPM.On the basis of Sect.4.2, the water retention curve (WRC), h p − S r , for the drainage process can be drawn.As a first step, the VPM provides the distribution of the connected spherical elements with different sizes; then, all of the voxels with spherical elements share the same label.By sharing the same label, the behavior of the voxels seems to flow and this is percolation flow.
The diameter of the spherical element at each step contributes to the calculation of the capillary pressure head (h p ) by the Young-Laplace equation as where T is the surface tension between water and air (72.88 mN m −1 at 20 • C), θ is the contact angle (49 • ), γ w is the density of water, and d is the diameter of the tube.In this study, a WRC test for the drainage process could be performed; thus, it was simulated by the following treatments.In the first step, each labeled spherical element is categorized, and the spherical element number with the maximum diameter is recognized.Second, the number of voxels with the spherical element number corresponding to the maximum diameter from the direction of the air entry side is counted, and it should continue until the discontinuous condition is reached.Third, the first step yields the capillary pressure head using the Young-Laplace equation with the substitution of the latest perforated diameter.The second step yields the degree of saturation by dividing the number of voxels not counted by the total number of pore spaces.This can be plotted on one WRC.Fourth, once the counting process for the second step is finished, the next label of spherical element should be checked on the basis of the same process detailed in the first three steps.Lastly, the WRC can be created.In order to verify the perforated pore diameter, WRCs were obtained from an experiment at 20 • C, and an image analysis was carried out in this study.
Figure 8 shows a binary 3-D image of the pore space of Toyoura sand based on the method in Sect.3.1.In this figure, the soil particles are invisible.Figure 9a-e show binarized Xray CT images obtained from Fig. 8 in two dimensions after the GIA using 13 different spherical elements, and Fig. 9f shows the final analysis results by overlapping each result.White represents the pore space, and black represents the soil particles.Eventually, this image processing was conducted in three dimensions to obtain a 3-D map, as shown in Fig. 10.Visually, each color element is uniquely distributed, as shown in Fig. 10.Two neighboring elements from the 3-D map are found to have neighboring colors in the color bar.That is, the pore size has a continuous distribution.From the viewpoint of hydraulic behavior, the local velocities of the pore water are always different at each pore.
Figures 11a-f show X-ray CT images simulated from the VPM analysis in dimensions of 300 voxel at each capillary pressure analyzed by Eq. ( 5). Figure 11g shows the occupation ratio of the cumulative volume of the spherical elements obtained from Fig. 11a-f.In fact, the occupation ratio of the cumulative volume of the spherical element countervails the volume ratio of air in the pore structure; therefore, the degree of saturation can be evaluated by subtracting the cumulative volume of the spherical elements from the entire pore vol-ume.This is expressed as where V T is the volume of entire pore structure.The number of voxels counted produces the degree of saturation (S r ); hence, the WRC for drainage can be obtained as an imageanalyzed curve.Figure 12 presents a comparison of the degree of saturation obtained from Fig. 11 and that from the Toyoura sample; thus, the analyzed degree of saturation was verified from the results in Fig. 12. Figure 13 shows the WRC analysis results for five different voxel sizes.The authors validated the effect of the voxel size on the WRC.In Fig. 13, the test results obtained from a laboratory test are also plotted.
Focusing the effect of the voxel dimension as the REV, when S r is greater than 0.8, it is observed that a lower voxel dimension yields a 50 % underestimation of the measured data to an air entry pressure as a result of decrease in S r , remarkably.
In this case, as shown in Fig. 11, the voxel dimension of 100 is not sufficient to become the REV for the WRC evaluation, where S r is greater than 0.8; however, there is no difference between the voxel dimensions where S r is less than 0.5.Provided that the voxel dimension is greater than 300 on the basis of all of the issues discussed, the VPM could provide a reasonable pore diameter, and it is concluded that the diameter of the spherical element can be interpreted as the pore diameter.
Despite the small change in the capillary pressure between 30 and 40 cm, S r decreased from 0.9 to 0.3.This behavior indicates that the mean pore size caused by a capillary pressure head of 30-40 cm is distributed in the saturation degree between 0.2 and 0.9 as shown in Fig. 13.This behavior should be caused by sands with a value less than the uniform coefficient obtained from the sieving test.et al. (2013) and Muljadi et al. (2015) evaluated the permeability using the image data of the pores.They also discussed the REV of their sample and analyzed the pore structure using image data.The studies that consider a pore structure, in geoenvironmental engineering and petroleum engineering have to evaluate the degree of saturation and the capillary pressure; a CT image can provide the degree of saturation without the use of any model.The capillary pressure is the driving force for transferring fluids in the unsaturated condition; however, it is difficult to define the capillary pressure in porous media on the basis of both experimental work and an image analysis if the relationship between the pore size and the capillary pressure is modeled.Blunt et al. (2002) discussed the capillary pressure in porous media using Poiseuille's law to evaluate the relative permeabilities of oil and water.Blunt et al. (2013) and Iglauer et al. (2013) analyzed the basic physical properties of pores from a CT image; in particular, Blunt et al. (2013) first estimated the relative permeability and then obtained the WRC from the relative permeability distribution.The issues of how to model the migration of oil in porous media such as rocks/soils and how to inject air for remediating contaminated soil by fuels require that the water-oil flow in the soil quantitatively understands the pore structure.Mayer and Miller (1993) visualized blobs of nonaqueous phase liquids (NAPLs) in a model test and deduced the blob diameter using the Young-Laplace equation, the number of capillaries, and the number of bonds; in short, the model was used to evaluate the pore size.The MIT, SEM, and AIM (Sato et al., 1992;Kamiya et al., 1996;Uno et al., 1998) have been used to measure the pore diameter.In general, the MIT is used to evaluate the pore size in clay.Sato et al. (1992) and Kamiya et al. (1996) attempted to develop a method to evaluate the pore size in sandy soil.Uno et al. (1998) included the moisture characteristic in the results obtained from an AIM (Kamiya et al., 1996) and proposed the moisture characteristic curve method (MCCM).The measurement principle of the AIM is similar to that of the MIT, and the obtained pore size is evaluated as the diameter of a pipe; however, the water contents were not measured.Uno et al. (1998) deduced the capillary pressure head using Darcy's equation for air permeation, and then they evaluated the pore size on the basis of a pipe model (Kamiya et al., 1996).A WRC is constructed using the degree of saturation and the  capillary pressure head measured by the head method with suction or given by Young-Laplace's law with the diameter of a pipe as the representative pore diameter.The AIM gives the statistical distribution of a pore diameter model of sandy soil as a glass tube.In short, the pore diameter is defined as the diameter of the tube, as per the implicit agreement in a number of papers.

Mostaghimi
Figure 14 shows the distribution of a perforated pore diameter for Toyoura sand.A spherical element is representative of the pore size in this study.The GIA counts the number of voxels for each size of the spherical elements; therefore, we can determine the total number of voxels, which can be multiplied by the volume of one voxel to obtain the total pore volume.Then, it is possible to calculate the percentage of the perforated spherical element for each size of the spherical element on the basis of the volume of the spherical element as the percent finer by volume in Fig. 14. Figure 15 presents a comparison between the pore diameters deduced by Uno et al. (1998) and those of the authors.In particular, the measured results between 0.065 and 0.85 mm have a better fit than those between 0.03 and 0.055 mm.This indicates that the AIM overestimated between 0.03 and 0.055 mm.The AIM by Uno et al. (1998) supposes that the pore space in sand exists in a tubular form.Since air wants to intrude a large pore space, it is thought that the evaluation of the pore size by the AIM has a high precision; however, the precision decreases for the evaluation of pore space with a small size, where air does not smoothly intrude.On the other hand, the VPM can easily evaluate the size of a complicated pore space regardless of the physical interaction and does not assume that the complicated pore is a straight tube.Hence, these results imply that the pore size obtained from AIM is not a Poiseuille distribution.
The VPM also evaluates the connectivity of the pore space.The GIA provides the number of voxels of the spherical element and the spatial distribution with the VPM.The AIM can also provide the pore diameter as the inner diameter of a pipe but not the spatial distribution of the pore diameter.This issue indicates that the VPM has the advantage of being able to estimate the WRC.In fact, the distribution in Fig. 14 can provide a pore-size distribution function (PDF) with respect to the perforated pore diameter.The PDF can also provide the degree of saturation by the sum of the spherical element.Figure 16 presents the WRC analyzed by the VPM and PDF in this study.The WRC obtained from the PDF was far from the results of the VPM in terms of the measured points.The definition of the pore size could be requested based on Fig. 16.As described in Sect.4.1, the VPM considers percolation using cluster labeling based on the connectivity of the pore spaces.On the other hand, the PDF does not utilize percolation concepts.Figure 16 shows that a reasonable WRC can be obtained from the degree of saturation and the distribution of pore diameter using the properties of percolation.Therefore, it is significantly useful that the GIA and VPM can estimate the water retention property on the basis of the geometry of the pore structure without performing a WRC test.

Conclusions
In this study, a specimen of Toyoura sand was scanned using a microfocused X-ray CT scanner, and the 3-D spatial distribution of a spherical element as the pore size (d) was visualized and evaluated quantitatively by a GIA.The GIA is a useful image analysis method for evaluating the pore diameter in a 3-D CT image.Moreover, the VPM was newly developed in this study, and it was validated by comparing the analyzed WRC with the measured results.The key conclusions are summarized as follows: 1.The size of a voxel affected the results of the image analysis.When the size of one voxel was 5 × 5 × 5 µm 3 , the REV for evaluating the physical properties of granular materials, which are similar to Toyoura sand, was at least 300 voxels for the evaluation of the grain size, porosity, surface area, and perforated pore diameter.In particular, it was possible for the porosity and surface area to be evaluated with an RSD of less than 1 %.
2. The results of the GIA show that the perforated pore diameter was less than the pore diameter from the AIM and less than 0.068 mm; moreover, mostly similar pore diameters were evaluated near 0.085 mm.Hence, the AIM provided partially different pore diameters from the results of the GIA.This issue revealed that the pore diameter obtained from the AIM was not a Poiseuille distribution.
3. The AIM can estimate the pore diameter as the diameter of a pipe and the occupation ratio; however, the spatial information is not included.Therefore, it was difficult to assess the WRC based on the pore diameter and its occupation ratio.In contrast, the newly proposed VPM in this study can distinguish each spherical element by labeled number corresponding to the radius of the spherical element.As a result, the connectivity of complex pore spaces can be evaluated; therefore, the VPM provides a WRC close to the measured result.
4. It was concluded that the VPM can easily evaluate the size of a complicated pore space regardless of the physical interaction and does not assume that the complicated pore is a straight tube.Then, the capillary pressure head can be deduced on the basis of the Young-Laplace equation as long as precise image data of the pores in sand were obtained by a microfocused X-ray CT scanner.
The second, third, and fourth conclusions are based on the first conclusion.The appropriate dimension for the image analysis should be defined on the basis of the particle diameter and voxel size.When using a microfocused X-ray CT scanner, a smaller sample should be scanned if a greater resolution is required.In future work, it will be necessary to verify the appropriate dimensions (i.e., the REV) for several types of grains with round and angular shapes and wide range of grain sizes.These features will further clarify the issue of pore connectivity with respect to the aspect ratio, which affects the results of the WRC.
The Supplement related to this article is available online at doi:10.5194/se-7-929-2016-supplement.

Figure 1 .
Figure 1.Binary image of pores between grains (black indicates the pore space and white indicates the particles).
Figure 2. Explanation of features of GIA.

Figure 3 .
Figure 3. X-ray CT image of Toyoura sand with different regions of image analysis in 3-D.

Figure 4 .
Figure 4. Relative standard deviation by changing the dimension of image analysis.

Figure 5 .
Figure 5. Grain size distribution curve obtained from image analysis.

Figure 6 .
Figure 6.Water retention test apparatus with the elevation head method.

Figure 8 .
Figure 8. X-ray CT image of pore in Toyoura sand in 3-D.

Figure 10 .
Figure 10.Distribution of perforated pore size in 3-D.

Figure 11 .
Figure 11.CT images of percolated pore space as the drainage process.

Figure 12 .
Figure 12.Comparison of saturation degree measured and analyzed at each capillary pressure.

Figure 13 .
Figure 13.Water retention curves obtained from image analysis and experiment.

Figure 14 .
Figure 14.Perforated-pore size distribution for each image size.

Figure 15 .
Figure 15.Comparison of pore size obtained from AIM and VPM.

Figure 16 .
Figure 16.Water retention plots obtained from GIA and VPM.

Table 1 .
Specifications of the microfocused X-ray CT scanner.

Table 2 .
Scan conditions used in this study.