VZJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (11)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Wildenschild, D.
Right arrow Articles by Kent, A. J. R.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Wildenschild, D.
Right arrow Articles by Kent, A. J. R.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Wildenschild, D.
Right arrow Articles by Kent, A. J. R.
Related Collections
Right arrow Laboratory Column Studies
Right arrow Pore-Scale Modeling
Right arrow Experiment Design
Published in Vadose Zone Journal 4:112-126 (2005)
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

ORIGINAL RESEARCH

Quantitative Analysis of Flow Processes in a Sand Using Synchrotron-Based X-ray Microtomography

D. Wildenschilda,b,*, J. W. Hopmansc, M. L. Riversd and A. J. R. Kenta

a Department of Geosciences, Oregon State University, Corvallis, OR 97331
b Environment and Resources, Technical University of Denmark, DK-2800 Lyngby, Denmark
c Hydrology, Department of Land, Air and Water Resources, University of California, Davis, CA 95616
d Consortium for Advanced Radiation Sources and Department of Geophysical Sciences, University of Chicago, IL 60637
* Corresponding author (wildend{at}geo.oregonstate.edu)

1 GSECARS: GeoSoilEnviroCARS; CARS: Consortium for Advanced Radiation Sources. Back


2 Simple image processing was based on cluster analysis only. Without the dry image of the sand matrix, an accurate estimation could not be determined as cluster analysis falsely detects a solid phase along air–water interfaces. This error is eliminated when the dry image is available. Back


3 The saturations in these figure are estimated from the sample porosity. Back


Received 6 April 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS
 CONCLUSIONS
 REFERENCES
 
Pore-scale multiphase flow experiments were developed to nondestructively visualize water flow in a sample of porous material using X-ray microtomography. The samples were exposed to similar boundary conditions as in a previous investigation, which examined the effect of initial flow rate on observed dynamic effects in the measured capillary pressure–saturation curves; a significantly higher residual saturation and higher capillary pressures were found when the sample was drained fast using a high air-phase pressure. Prior work applying the X-ray microtomography technique to pore-scale multiphase flow problems has been of a mostly qualitative nature and no experiments have been presented in the existing literature where a truly quantitative approach to investigating the multiphase flow process has been taken, including a thorough image-processing scheme. The tomographic images presented here show, both by qualitative comparison and quantitative analysis in the form of a nearest neighbor analysis, that the dynamic effects seen in previous experiments are likely due to the fast and preferential drainage of large pores in the sample. Once a continuous drained path has been established through the sample, further drainage of the remaining pores, which have been disconnected from the main flowing water continuum, is prevented.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS
 CONCLUSIONS
 REFERENCES
 
IN THE PAST FEW YEARS, renewed interest has emerged in the area of dynamic effects in the capillary pressure–saturation and hydraulic conductivity relationships for unsaturated or multiphase flow. Hassanizadeh et al. (2002) describe the term dynamic effect as being responsible for the dependence of the capillary pressure–saturation relationship on the rate of change of saturation; they emphasize that the dynamic effects considered in their paper differ from the flow-rate dependence of the relative permeability coefficient. Over the past three to four decades, a substantial body of evidence has materialized recognizing the nonuniqueness of the capillary pressure–saturation relationship; showing that it, in reality, depends on the flow dynamics of the system. Early experimental work was performed by, for example, Topp et al. (1967), Smiles et al. (1971), and others listed in detail by Hassanizadeh et al. (2002) and Wildenschild et al. (2001). In addition, recent experimental work by Hollenbeck and Jensen (1998), Wildenschild et al. (2001), as well as Mortensen et al. (2001) has also illustrated the nonuniqueness of unsaturated flow properties, and thus leads to questions about the unconditional use of Darcy's Law for multiphase flow problems. Mortensen and coworkers used a quantitative light transmission technique to investigate millimeter-scale phase displacement processes in two dimensions and suggested, based on their results, that the current definitions of unsaturated hydraulic properties are perhaps inadequate.

The lack of a firm theoretical foundation of the relationship between capillary pressure and saturation has been the focus of a vast body of work by Gray and Hassanizadeh over the years (e.g., Hassanizadeh and Gray, 1990, 1993; Gray and Hassanizadeh, 1998; Gray et al., 2002). In their theoretical development, which is based on macroscopic balance laws, thermodynamics, and appropriate constitutive relationships, the authors emphasize the importance of explicit inclusion of interfaces and interfacial properties, since they are known to play an important role in determining the thermodynamic state of the entire system under scrutiny. It has been hypothesized by Hassanizadeh and Gray (1993) that traditional, hysteretic, capillary pressure–saturation curves are actually a two-dimensional projection of the three-dimensional capillary pressure–saturation–interfacial area surface, and as such many of the observed nonunique or dynamic patterns can allegedly be explained by considering the full three-dimensional relationship. The concept of developing systematic conservation equations for fluid flow in porous media was first introduced with the description of phases and their properties in the 1960s by Anderson and Jackson (1967) and Whitaker (1967). There have been many other additions since then, yet, regardless of which theoretical development we are aiming to improve, and how we aim to improve it, we need experimental evidence to support these new developments.

Concurrent with the theoretical development, there has been an explosion of interest in recent years in pore-scale modeling of multiphase flow in porous media, both network models (e.g., Lowry and Miller, 1995; Celia et al., 1995; Reeves and Celia, 1996; Patzek, 2001; Dalla et al., 2002; Berkowitz and Hansen, 2001) as well as lattice Boltzmann-type models (e.g., Gunstensen and Rothman, 1993; Martys and Chen, 1996; Chen and Doolen, 1998). A detailed discussion of the physics and predictive capabilities of the various types of models was presented by Blunt (2001). Some of the above approaches have been used to predict properties such as interfacial area and common line lengths in numerical experiments. The pore-scale models are useful tools that provide important insight in lieu of experimental data, yet, there is no substitute for real experimental data. Fortunately, recent developments in microscale visualization techniques have facilitated new opportunities for imaging pore-scale flow and transport processes in three dimensions, for multiple phases, and in real porous media. Synchrotron-based X-ray microtomography has been applied to multiphase flow problems in recent years, predominantly in the petroleum industry, to describe variables such as porosity, pore topology, formation factor, electrical resistivity, and permeability in Fontainebleau sandstone (Spanne et al., 1994; Auzerais et al., 1996); for characterizing the pore network structure in a reservoir rock (Berea sandstone) for subsequent use in lattice Boltzmann models (Hazlett, 1995; Lindquist et al., 2000); and for imaging displacement of oil by water in reservoir rocks (Coles et al., 1998). However, we do not know of existing synchrotron-based microtomography work that takes on a rigorous quantitative approach to pore-scale multiphase flow problems. For instance where:

In previous work by Wildenschild et al. (2001) we listed several hypotheses that could potentially explain rate dependence or nonuniqueness of the unsaturated flow properties that was observed. However, we were unable to provide any further clue as to which process, if any of them, was most likely responsible for the observed dynamic deviations from static measurements. A natural next step, which we report on here, is to nondestructively visualize the flow process in the sample as it happens by developing minute pore-scale multiphase flow experiments that can be performed interactively while at the same time collecting high-resolution tomographic images. To accomplish this, we have used the X-ray microtomography facility at the Advanced Photon Source (APS), GSECARS1 beam-line. These first experiments only consider the drainage branch of the capillary pressure–saturation relationship. In the following, we report on the experimental set-up and procedure, derive quantitative information from the tomographic images, and we analyze the results to demonstrate a boundary condition effect on pore-scale flow. In the past we have referred to these effects as dynamic effects, however, there is currently no clear definition of dynamic effects, and we therefore attempt to report the results in terms of the controllable variable; the boundary condition. We present both quantitative and qualitative results, the latter primarily to give the reader an idea of the imaging detail possible with this technique, and thus offer insight into the possibilities of X-ray microtomography as an analytical tool for addressing multiphase flow problems.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS
 CONCLUSIONS
 REFERENCES
 
The experimental setup at the GSECARS beam-line
The bending-magnet beam-line at sector 13 (GSECARS) is equipped to perform high-resolution, three-dimensional imaging via absorption tomography. The objects can range in size up to approximately 50 mm wide, with vertical sequences 5 mm tall (Rivers et al., 1999). The stage is automated and computer controlled to rotate and translate in both the horizontal and vertical directions, and it is possible to perform fluid flow experiments in the hutch during the scans. The beam-line specifications are detailed in a previous publication (Wildenschild et al., 2002) and we will only provide the most necessary details here. The theoretical background of the tomography approach is also covered to some degree in this previous publication, and for more details we refer the reader to, for instance, McCullough (1975), and to Kinney and Nichols (1992) for the specifics of synchrotron-based tomography. All of the experiments reported here were performed at an energy level immediately above the K-shell photoelectric absorption edge of iodine, slightly above 33.17 keV. A 13% (weight) iodine-spiked solution was used as the wetting fluid in the experiments to facilitate absorption-edge imaging, whereby a very high intensity is obtained in the fluid phase, rendering the water and air phases separable in the resulting intensity maps. If the water phase is not spiked it is difficult to differentiate between the water and air intensities in the images. The desired X-ray energy level is obtained using a monochromator, which selects a single energy from the white synchrotron radiation. In the experiments reported here either an in-hutch channel-cut Si (220) crystal or a beam-line Si (311) monochromator was used for this purpose. During our use of the facility the beam was run in fill-on-fill mode, which means that the ring was refilled twice daily. Between those refills the beam decayed over a 12-h period from approximately 100 to 60 mA. Because this affects the intensity of the images slightly, the exposure time was adjusted accordingly during the decay period, to obtain approximately 2/3 absorption in any one pixel of the CCD-camera. This ensures good signal-to-noise ratios in the raw image. During a full rotation (180° with either 360 or 720 steps) a number of flat or white-fields (generally 25 or 50 angles apart) are captured, where the sample is moved out of the field of view and a baseline intensity map is acquired. This is used for normalizing the intensities before the subsequent reconstruction, correcting for beam intensity decay, nonuniformities in the incident X-ray beam, as well as nonuniform responses of scintillator and CCD-detector. Also, the camera dark current is measured and corrected for. The dark current and flat fields are used to normalize each frame so that complete absorption results in an intensity of 0.0 and no absorption produces an intensity of 1.0. The pixel values of the reconstructed images reflect absolute linear absorption coefficient in units of inverse pixel size. Dividing by the pixel size, say in millimeters per pixel, yields the linear attenuation per millimeter. The reconstruction code was written by Rivers (1998)(the online tutorial), using IDL programming (Research Systems Inc., Boulder, CO) and the Gridrec FFT (fast Fourier transform) reconstruction software developed at Brookhaven National Laboratory (Upton, NY). The imaging details are listed in Table 1.


View this table:
[in this window]
[in a new window]
  		Table 1. Tomography specifics for the two experimental runs.

 
The schematic of the basic sample-holder design is illustrated in Fig. 1a . The cell was designed for unsaturated flow experiments, using an O-ring to seal in a porous nylon membrane (1.2-µm MAGNA filter from GE Osmonics Inc., Trevose, PA) at the bottom of the cell. The air-entry pressure of this membrane was sufficient for the pressures used in the experiments (in excess of 700 cm). The sample base and top is identical for all the cells used in the experiments, whereas the centerpiece used was either a straight (Fig. 1a) or a tapered (Fig. 1b) cylinder. We refer to Wildenschild et al. (2002) for a schematic of the general beam-line setup. The centerpiece with smaller cross-section was used for high-resolution imaging. The image resolution is directly dependent on sample size. The CCD camera has a set number of pixels (1317 by 1035 for the camera used for APS-6.0 and 1300 by 1030 for APS-1.5), thus a larger sample will result in lower resolution, whereas a smaller sample allows one to "zoom in" and image at higher resolutions. This restriction on sample size is due to the requirement for quantitative tomography that air be imaged on both sides of the sample. Otherwise reconstruction of absolute linear attenuation coefficients is not straightforward. All cell components were constructed of lucite which has a low absorption relative to the sand and spiked water at 33.17 keV.



View larger version (37K):
[in this window]
[in a new window]
  		Fig. 1. (a) Sample holder with porous nylon membrane. The centerpiece is screwed into the bottom piece, sealing the O-ring against the semipermeable nylon membrane. The inside diameter is 6.0 mm. (b) Alternative centerpiece for higher resolution imaging. Height and membrane configuration is identical to Fig. 1a. The inside diameter is 1.5 mm.

 
Verification of obtained linear attenuation coefficients using absorption theory
A beam of monoenergetic photons with an incident intensity Io, penetrating a material of thickness x and density {rho} results in a reduced intensity, I, according to Lambert-Beer's Law:

[1]
which rearranges to

[2]
where µm and µ are defined as the mass attenuation coefficient (L2/M), and linear attenuation coefficient (LAC) with units of inverse length, respectively. Thus, the mass attenuation coefficient µm = µ/{rho} can be obtained from measured values of Io, I, and x, and can also be looked up in tables (e.g., Saloman et al., 1988; Hubbell, 1982) for known materials as a function of photon energy. For a mixture, the mass attenuation coefficient is obtained by simply adding the contributions from the various components:

[3]
where wi is the fraction by weight of the ith atomic constituent, and the energy specific (µm)i are obtained from tables for each constituent, or alternatively from the NIST XCOM web site http://physics.nist.gov/PhysRefData/Xcom/Text/XCOM.html (verified 9 Dec. 2004). The LAC of the mixture or compound can then be derived from the mixture–specific mass attenuation coefficient by multiplying it with the average density of the mixture. At the low X-ray energies used here (50–100 keV), X-rays interact with matter predominantly via photoelectric absorption, which is strongly dependent on atomic number. This means that the LAC depends on the mean atomic number and density of the sample material in question and consequently LAC is linearly related to sample fluid saturation. To verify the measured LAC values obtained in these experiments, we calculated the theoretical LAC for a sample at a water saturation of 81% and compared it to the measured average of all 516 slices of a 3.45 mm–tall section of the satiated 1.5-mm diameter (APS-1.5) sample. The latter was estimated from simple2 image processing, resulting in an average saturation of 81%. A few assumptions were made in calculating the theoretical value of the LAC. First, the sand was assumed to consist of 100% SiO2. Second, when calculating the fluid density it was assumed that the addition of KI did not change the fluid volume. Third, the air phase was assumed to consist of 100% N2. The constituents and their respective mass fractions were entered into the XCOM interactive spreadsheet using the 33.7 keV energy level, resulting in a total (sum of coherent and incoherent scattering) mass attenuation of 1.18 cm3 g–1. By multiplying this by the average wet density of the sample, a LAC value of µ = 2.47 cm–1 was obtained. Considering all potential error sources, this final computed value agrees well with the average measured value for the same sample with µ = 2.45 cm–1. In this calculation, the porosity is calculated from the mass of sand filled into the sample holder, the volume of the holder, and using the density for pure quartz (2.65 g cm–3). Calculations using porosity values obtained from image processing and back-calculating the corresponding sand fraction resulted in a LAC value of 2.80 cm–1, which is close considering all potential errors.

Flow Experiments
The flow experiments were performed on small cylindrical pressure cells, containing packed sandy material, with inside diameters of 6.0 and 1.5 mm, respectively. The sample holder was open to atmospheric pressure on top during the lower-pressure steps (<100 cm), and drainage (or imbibition) was induced by lowering (or raising) a graduated syringe, which was hydraulically connected to the bottom outflow port of the sample holder. For the high-pressure steps (>100 cm), air-phase pressure was applied at the top inlet port of the sample to induce sample drainage. The amount of outflow was measured with the attached syringe or pipette. Imaging was completed when outflow had ceased for the various pressure steps, generally the pressures were applied on the order of 1 to 2 h. When using the synchrotron facility, only a few days of experimental time are generally allotted for each user. Therefore, we could not let the samples drain long enough to reach equilibrium. The experiments performed are listed in Table 2. The porous material chosen for all experiments was a Lincoln sand (mean particle size of 0.17 mm), which was one of the materials found to produce dynamic capillary pressure effects in earlier work (Wildenschild et al., 2001). The material was packed to a porosity of 34%, estimated from the mass of sand and the volume of the sample holder. Saturation was achieved by raising a connected fluid reservoir above the sample and waiting for approximately 1 h. Complete saturation was never obtained using this approach because of entrapped air.


View this table:
[in this window]
[in a new window]
  		Table 2. Experimental details. The experiments are listed in the order they were performed.

 
As listed in Table 2, three different sets of data were obtained for the two runs, (i) saturated, (ii) drained slowly (multistep) and (iii) drained fast with a high air-phase pressure (one-step). The APS-6.0 sample was imaged at three steady state points for the multistep run, and the sample was imaged in five sections, which stack on top of each other to produce a vertical cross-section of 15.5 mm. Similarly, 3 sections were imaged for the APS-1.5 run to produce a 9.5 mm–tall imaged section of the sample. It should be noted that only part of the total sample height (3.88 cm) was imaged in these experiments, thus no information about boundary effects in the vertical direction was obtained.

Image Processing and Data Analyses
Linear Attenuation Coefficients
The output from the reconstruction program consists of gray-scale intensity maps based on the distribution of linear attenuation coefficients in the image. The resulting volumes consist of arrays of dimension (H x H x V), where H is the number of pixels in each of the horizontal directions of the field of view and V is the equivalent in the vertical direction. H is typically close to the maximum number of double-binned pixels on the CCD camera (658 and 650 for APS-6.0 and APS-1.5, respectively), while V varies between runs as a function of sample size and resolution and the vertical size of the incident beam. In the following, we will refer to an entire three-dimensional data set as a volume, a horizontal (or vertical) cross-section as a slice, and to several volumes imaging a larger section of a sample as a stack of volumes. A volume element, that is, the three-dimensional equivalent of a pixel is called a voxel. The linear attenuation coefficients were determined as slice averages, such that one value for each slice over the entire imaged sample height was obtained. Programs were written in IDL to derive averages and standard deviations from the presented images. The images were cropped to only consider the cylinder containing the actual sample, that is, by subtracting the lucite sample holder and surrounding air. For all of the images shown, the gray areas represent pixels dominated by sand grains, black represents the air phase, and white areas are pixels containing KI-spiked water. When images are compared, they are always scaled to the same interval before printing. No attempt has been made to detect or segment individual pores and menisci properties, however, this type of analysis is possible for high-resolution data and will be reported in a future publication.

Histogram Segmentation Supporting Spatial Analyses
Segmentation of the gray-scale images into the respective phases (air, sand, and water) is not straightforward, although it appears so to the human eye. Figure 2 shows a typical histogram for a slice in one of the volumes; the intensity values represent absolute linear absorption coefficient in units of inverse pixel size. The left peak represents the air phase (lowest intensity), while the broader peak on the right contains pixels representing both sand grains and the KI-spiked water phase (highest intensity). Obviously, there is some overlap of pixels containing a mixture of sand and water (known as the partial volume effect), making it difficult to set an absolute threshold level for the segmentation of sand and water. Clausnitzer and Hopmans (1999) reviewed the problems of partial volume effects and presented an approach for accurately determining the volume fraction of each phase within a voxel based on histogram segmentation. However, in our study where we compare images for different saturations, there is an additional complicating factor. Specifically, there are slight variations in the absolute scale of intensities among the different volumes due to factors such as the harmonic content of the monochromator and small variations in dark current, and perhaps other effects that are not understood. It is possible that we register some third harmonic beam content, which has a different absorption coefficient. This means that a specific threshold value for one volume, does not necessarily apply to another volume with a different fluid saturation. The approach by Clausnitzer and Hopmans (1999) is designed to provide volume fractions for the various phases, not to determine the physical location of the phases. Since the latter is the ultimate goal of our research, an alternate approach based on cluster analysis in conjunction with the availability of a "dry" image of the sand phase only, will eventually be used. Therefore, in the work presented here that involves separating the phases, we only perform a quantitative segmentation and a subsequent analysis of the air phase. After median filtering the images with a 3 by 3 kernel, setting a threshold to determine the volume fraction and spatial arrangement of the air phase can be done with sufficient accuracy because of the well-defined minimum between the two peaks in Fig. 2. However we had to use different thresholds for the different volumes, due to the reasons mentioned above. We used a gray-scale intensity value of 6150 for the phase separation for the slow-drainage images, and a value of 7300 for the fast drainage images. Figure 3 illustrates the effect of varying the gray-scale intensity threshold by 1.5 and 5%, resulting in average saturation3 differences of 1.6 and 5.5% for the fast-drainage images (7300 threshold), and 2.7 and 8.8% for the slow-drainage images (6150 threshold). The figures illustrate the effect of degree of saturation on the uniqueness of the threshold value; when the sample is less saturated, the threshold is easier to set because the air-phase peak in the histogram becomes larger, yet the error increases because of increased partial volume effects. Because the air and water phase are uniquely related, assuming that soil porosity remains constant, we can determine the volume fraction of one fluid phase from knowledge of the other. Selection of the threshold values also enables us to define three-dimensional surfaces, delineating the gas phase from the water phase using isosurfaces.



View larger version (14K):
[in this window]
[in a new window]
  		Fig. 2. Typical histogram for a horizontal slice of one of the partially saturated volumes.

 


View larger version (25K):
[in this window]
[in a new window]
  		Fig. 3. Sensitivity of saturation estimate to threshold of the gray-scale image. The saturation was estimated from the (mass and volume based) porosity. The analysis was performed using a gray-scale intensity of (a) 7300 as the threshold for fast drainage and (b) 6150 for slow drainage. The same axis length is used in both plots to ease comparison.

 
Spatial Analyses
The spatial distribution of attenuation values in an image can be described by the nearest neighbor distance (NND), which is a measure that, among other things, can reveal to what extent a distribution is random, clustered, or regularly spaced. As outlined in Fig. 9 and 10 of Russ (1999), a regularly spaced (or self-avoiding) distribution of values will have a distance spectrum with the largest NND, with a distinct peak in a histogram of NNDs. The clustered arrangement will have the lowest NND, but also have a distinct peak, whereas a random distribution will have a flatter (Poisson) distribution of NNDs with an intermediate average. We expect the spatial distribution of points in the presented images to resemble that of a clustered distribution, however, with a varying degree of spatial separation. To test this, scripts were written in IDL to calculate the NNDs for specific images. First, the program was tested by applying it to the three distributions presented in Russ (1999), confirming that our approach was valid. Only a squared subsection of 230 by 230 pixels of a circular horizontal slice was used in the analysis to avoid edge effects. Initially, NNDs were calculated for three slices (40, 110, and 160, see Fig. 1a) in both the slow- and fast-drained images for APS-6.0. Initially, the gray-scale images were median filtered with a 3 by 3 kernel. The air phase was then separated from the sand and water phases using histogram segmentation as described earlier. The same gray-scale intensity thresholds of 6150 and 7300 were used to create the binary images. Various erosion and dilation operations were first attempted, but eventually we settled on erosion using a 2 by 2 kernel because that best preserved all air-phase structures in the gray-scale images. For the NND analysis it is more important to preserve an air bubble than accurately calculating its actual size. All air blobs less than 4 pixels in size ({approx}68 µm) were eliminated, since such small blobs are most likely noise as opposed to true air bubbles. For a sand with a d50 of 0.17 mm this seems a reasonable resolution for eliminating air bubbles. Also, beyond 4 pixels the search routine (described below) did not work on the images because at 4 pixels or less the blobs are predominantly arranged in strings instead of as a closed blob for which the centroid can be calculated. The bubbles or blobs were subsequently individually labeled, and a two-dimensional search routine was used to obtain the x- and y-centroids, as well as the area and perimeter for each blob. Finally, the minimum Euclidean distance for each blob was calculated from which a nearest neighbor distance versus frequency histogram was obtained. The calculated blob size is naturally dependent on the selection of gray-scale threshold; however, we found that variations were small, even for the lower saturation multistep case, if the threshold value was varied by less than ±100 intensity units: Using a threshold value of 6150, the average blob size was 111.5 pixels (or {approx}0.19 mm, assuming square blobs for simplicity). Using threshold values of 6250 and 6050, average blob size values of 102.5 and 117.2 pixels, respectively, were computed, corresponding to errors of 8 and 5%. At a perturbation of 50 intensity units the error is less than 4% for thresholds above and below the mean threshold value.


   RESULTS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS
CONCLUSIONS
 REFERENCES
 
Qualitative Results
An example of the kind of information obtained with this technique is illustrated in Fig. 4 , which shows drainage across a horizontal cross-section of the 6-mm diameter APS-6.0 sample. The images illustrate that drainage is heterogeneous at first, with larger pores draining first before the smaller pores. In the last image, however, at a final capillary pressure of 490 cm, the pore matrix is drained uniformly across the sample. The sequence of images in Fig. 4 also demonstrate that individual sand grains are observed at the same spatial location in all four images, thereby confirming that no significant soil particle movement occurred during the drainage process. Since a dry soil sample was not imaged for these first experiments, we were not able to calculate absolute water saturation values. However, because the LAC values depend only on the mass fractions of the sample's atomic constituents and total mass density, that is, water saturation, we interpret resulting variations in LAC as representative of relative changes in saturations. Thus, in the following discussions, the term LAC is used interchangeably with saturation.



View larger version (164K):
[in this window]
[in a new window]
  		Fig. 4. Drainage of APS-6.0 sample during multistep experiment (a) p = 0 cm; (b) p = 24 cm; (c) p = 79 cm; (d) p = 490 cm. Horizontal slice at 13.73 mm from the bottom. In these and all other gray-scale images the white or light areas represent the water phase, black is air, and the gray areas are sand grains.

 
Visible Effects of Applied Boundary Condition
The uniform drainage patterns observed in Fig. 4d were not obtained when the sample was drained with a single step using a high air-phase pressure. A sequence of images for a vertical cross-section show the drainage process from "saturation" (Fig. 5a) through several multisteps (Fig. 5b–d), as well as the phase distribution after a single pressure step (Fig. 5e). We find that the distribution of water and air is very different in Fig. 5d and 5e although the two images were drained to the same final capillary pressure value of 490 cm. In Fig. 5e the sand matrix remains mostly saturated, while only the larger pores are emptied at this pressure. In comparison, the matrix is uniformly drained under the slow-drainage scenario of Fig. 5d. The amount and distribution of fluid drainage obtained using the 490-cm one-step approach (Fig. 5e) is comparable to the drainage obtained at a capillary pressure of only 24 cm (Fig. 5b), representing slow drainage. Similar results are illustrated for the entire imaged section of the sample in Fig. 6a and 6b , which show a complete vertical profile exposed to slow and fast drainage to a final capillary pressure of 490 cm. As in Fig. 4, individual sand grains can be observed in both cross-sections to confirm that we are indeed looking at the same cross-section for the two images. In Fig. 7 , we have isolated the air phase and rendered it as a continuous surface to illustrate the effect of the different drainage scenarios in three dimensions. The yellow (slow drainage) and blue (fast drainage) surfaces depict the air phase in the sample following the two different drainage scenarios.



View larger version (71K):
[in this window]
[in a new window]
  		Fig. 5. Comparison of multistep (slow) and one-step (fast) drainage for the same 3.6-mm vertical slice of the APS-6.0 sample. The section shown is located at the bottom of the imaged stack of volumes. Multi-step: (a) p = 0 cm; (b) p = 24 cm; (c) p = 79 cm; (d) p = 490 cm. One-step: (e) p = 490 cm.

 


View larger version (129K):
[in this window]
[in a new window]
  		Fig. 6. Difference in water and air distribution over the entire imaged section (15.5 mm) of the sample for (a) fast and (b) slow drainage (APS-6.0) to a capillary pressure of 490 cm. Fig. 6c. Cut-out of a small section of the sample illustrating in more detail the difference in fluid distribution between the fast (left) and slow (right) drainage processes.

 


View larger version (78K):
[in this window]
[in a new window]
  		Fig. 7. Three-dimensional distribution of the air phase for slow (yellow) and fast (blue) drainage for APS-6.0.

 
From Fig. 6a, 6b, and 7 it is evident that the difference in drainage pattern and resulting residual water-phase saturation among the two flow conditions is dramatic. Under high pressure, that is, high initial flow rates, a couple of processes can explain possible mechanisms controlling the observed drainage patterns:
  1. Initially drainage takes place in the most well-connected and largest pores; subsequent drainage of the smaller pores, that have a lower hydraulic conductivity is less likely because they are no longer hydraulically connected to the bulk fluid phase. This process would fit under the Category A mechanisms of Wildenschild et al. (2001):
  2. [A] Entrapment of water. This is a plausible mechanism at high flow rates. We hypothesize that water entrapment occurs through hydraulic isolation of water-filled pores by draining surrounding pores. The larger the drainage rate, the less opportunity exists for all pores to drain concurrently. It may occur throughout the soil sample and will increase water retention, but will decrease unsaturated hydraulic conductivity, as the mobile portion of the soil water is decreased. The highest flow rates occur for one-step experiments and for coarse-textured soils, and hence water entrapment is likely to be more prevalent in coarse soils with large head gradients.
  3. Some of the smaller pores are drained immediately as these are exposed to the high air-phase pressure, leaving many relatively large pores disconnected from other flowing pores, and thereby remaining undrained. This would likely have a tortuosity component since it would occur as soon as the first pore is drained continuously from the top to the bottom of the sample. This would be particularly likely to take place at the outflow boundary because of the large head gradient there. It would create a barrier to flow and isolate the pore system above the outflow boundary as described in Mechanism B of Wildenschild et al. (2001):
  4. [B] Pore water blockage. When applying a sudden large pressure step to a near-saturated or saturated soil sample, the large matric potential head gradients near the porous membrane result in faster drainage of the pores at the bottom of the sample than in the overlying soil, thereby isolating the conductive flow paths and impeding further drainage and equilibration of the soil. Consequently, the sample's unsaturated hydraulic conductivity will be reduced. Pore blockage is likely to occur for materials with a uniform pore size distribution (like the Lincoln sand). It is assumed that air is available to replace the draining water at the bottom of the sample, hence it can only occur if there is macroscopic air continuity across the whole sample (air entry value of soil must be exceeded).

A similar result related to the outflow boundary was also found by Mortensen et al. (2001). Likely, the effects of either of the two processes are similar; that is, a large number of pores end up containing residual water, resulting in less overall drainage at high flow rates than at slow flow rates.

The measured outflow for the two experiments illustrated in Fig. 5a and 5b is also significantly different: 0.22 mL for slow drainage and only 0.16 mL for fast drainage, corresponding to residual volumetric water contents of 0.16 and 0.21, respectively. This is in agreement with our earlier observations for Lincoln sand (Wildenschild et al., 2001, Fig. 3), where a higher residual water content was observed for fast drainage in one-step capillary–pressure saturation measurements, as compared to multistep or curves measured quasistatically using a syringe pump.

Studying the images in Fig. 6a, 6b, and 7 more closely, it appears that for slow drainage (multistep) the majority of the pore space is drained quite uniformly, whereas for fast drainage (one-step), only a limited number of generally larger pores are drained. This is illustrated in even more detail in Fig. 6c; the large pores drain and the smaller pores remain saturated for fast drainage. For slow drainage the smaller pores tend to drain along with the larger ones. It should be noted that in Fig. 6a, it appears that some of the largest pores drained after fast drainage cannot be identified after slow drainage (Fig. 6b); most likely these pores are formed when small solid particles are being moved by the fast-flowing water. However, this only pertains to a few of the very large pores, and the general trend in the images is definitely toward larger pores being drained and smaller ones remaining saturated under fast drainage, whereas under slow drainage the entire pore space is uniformly drained, see Fig. 6c. Even with some grain movement for fast drainage, which establishes some larger pores that were not available during slow drainage, likely, the drainage pattern observed would be analogous: For slow drainage these larger pores would drain first, but there would still be opportunity for the smaller pores to drain into the larger pores, and the end result would be the same; a much more uniform drainage of the entire pore matrix.

Based on the visible number of large drained pores, it appears that Mechanism 1 is dominant for the fast-drained images, see Fig. 6c. The largest, well-connected pores are drained rapidly as the large pressure gradient is applied, whereas the smaller pores with lower conductivity are unable to keep up with the drainage process; they are therefore left without an opportunity to drain, except through slow film and corner flow. In contrast, under slow drainage, the pores drain in a more piston-type pattern, whereby all pores stay connected. Larger pores drain first and smaller pores are drained at the subsequent higher imposed pressure, but without the disconnection of water pathways. It should be noted that tortuosity and connectivity considerations may be important when extrapolating into the third dimension. In the following section we support our qualitative interpretation through quantitative analyses.

Quantitative Results
Quantifying Effects of Applied Boundary Condition
Horizontal slice–average linear attenuation coefficients for both experimental runs are presented in Fig. 8 and 9 for the APS-6.0 (6-mm diameter sample) and APS-1.5 (1.5-mm diameter sample) data, respectively. Figure 9 also shows vertical image profiles for both fast and slow drainage that can be compared with the LAC values.



View larger version (26K):
[in this window]
[in a new window]
  		Fig. 8. Linear attenuation coefficient profiles for APS-6.0 sample. The initial "saturation" is the primary saturation of the dry sand; its measurement was limited to the top section only. The "satiated" profile is the secondary saturation following the slow primary drainage. Both were obtained by raising the water level to the top the sample (no vacuum or CO2 venting was used).

 


View larger version (90K):
[in this window]
[in a new window]
  		Fig. 9. Linear attenuation profiles for (a) fast drainage (one-step) and (b) slow drainage (multistep) for the APS-1.5 sample.

 
Figure 8 shows a large variation in LAC or water saturation for the different pressure steps, as well as significant variations over the height of the imaged region. The two "saturated" or satiated conditions have the highest LAC values, particularly for the very first saturation; however, only one volume section was imaged for this saturation. In both cases, the sample was saturated by raising a connected water reservoir above the sample and waiting for approximately 1 h. We did not use CO2 flushing or vacuum saturation to attempt full water saturation. The first imbibition was performed on the dry sample (before the multistep experiment), and therefore a higher initial saturation was obtained than for the secondary imbibition performed before the one-step experiment. The LAC values obtained in the two cases are listed in Table 2 along with other pertinent information. The average LAC values of the particular volume section, located at the top of the imaged region for both satiated cases were 2.39 and 2.27 cm–1, for the first and second saturation, respectively. This corresponds to a difference of approximately 5%. The standard deviations for the two saturations are fairly similar with values of 0.052 and 0.035 cm–1, respectively. The remaining curves in Fig. 8 illustrate the effect of the boundary condition and resulting drainage rate on water saturation. The saturation profiles for the first (p = 24 cm) and second (p = 79 cm) multistep experiments show that the sample tends to drain at the bottom and in the region around 10-mm height at these low pressures. This is seen in Fig. 8 as regions with low linear attenuation coefficients. We note that because of the higher initial saturation, the subsequent first multistep leaves some sections of the sample at a higher saturation than measured at the second full "saturation." As the pressure is increased in the last multistep (p = 490 cm), the sample drains over the entire height, however, the region of lower saturation at the bottom and around 10-mm height are still visible. The average LAC value for this last step (p = 490 cm) is 1.76 cm–1, Table 2. In comparison, the fast drainage using the single pressure step of 490 cm results in a much higher residual saturation with an average LAC value of 2.15 cm–1, which is about 22% higher. The corresponding residual saturation is equal to that attained between the first (2.27 cm–1) and second (2.05 cm–1) multistep, or for a capillary pressure between 24 and 79 cm. This result is in agreement with the qualitative observations of Fig. 5, and also supports the observations made by Wildenschild et al. (2001) for the Lincoln sand. In addition to an overall difference in residual saturation obtained for the two boundary conditions, a distinct difference in drainage profile between the two flow regimes can be observed near the 10-cm height, yet, the preferential drainage at the bottom of the sample prevails for the fast drainage profile. The fact that the sample was at an initially higher saturation before the multistep experiment than the subsequent one-step experiment only enhances the difference between the two drainage scenarios: if this was not the case an even larger difference in residual saturation (or drained volume of water) between fast and slow drainage would be anticipated.

A different result was obtained for the small diameter sample in Fig. 9. Practically no effect of the applied boundary pressure, and resulting drainage rate, was observed in this case, as the two resulting average LAC values of 1.99 and 2.03 cm–1 for the multistep and one-step, respectively, differed only by about 2%. A very slight difference in saturation is observed at the 5.5-mm height, but otherwise the residual saturation profiles are almost identical. Compared to the saturated profile, the sample drains mostly in the region around 7.5-mm height for both the fast- and slow-drainage scenario, corresponding with a higher visible porosity.

The lack of difference in fluid distribution for the two drainage scenarios is also readily seen in the vertical gray-scale profiles of Fig. 9. Only a slight variation in fluid distribution is observed near the bottom of the imaged section of the sample. This lack of a boundary condition effect is also well illustrated in the three-dimensional images of Fig. 10 , where the yellow (slow drainage) and blue (fast drainage) surfaces that describe the air-phase distribution are very similar. So, it raises the question of why there is no effect of boundary condition or resulting drainage rate for this sample. Although the final capillary pressure value for this APS-1.5 test was somewhat lower with a value of 210 cm, previous experiments for the Lincoln sand (Wildenschild et al., 2001) showed that it was sufficiently large to cause a boundary pressure effect when applied as a single pressure step. We speculate that the sample diameter could be too small to establish a broad range of pore sizes and question whether a sufficient number of larger pores are present. In other words, we question whether the sample is large enough to adequately represent the same flow process as observed in the larger sample.



View larger version (66K):
[in this window]
[in a new window]
  		Fig. 10. Three-dimensional distribution of the air phase for slow (yellow) and fast (blue) drainage for APS-1.5.

 
To address this question, we conducted a representative elementary volume (REV) analysis on the LAC for the two samples, of which the results are shown in Fig. 11 . The linear attenuation coefficient was calculated for a cube of increasing voxel size centered about the center of the volume, with the largest volume defined by a cube that fits inside the diameter and height of the cylindrical sample. The LAC values are presented as a function of the grown voxel volume in cubic millimeters to facilitate comparison of the two different-sized samples containing the same porous medium. The results show that the 1.5-mm diameter of the APS-1.5 sample is not sufficiently large to provide a sufficient volume for a representative LAC or saturation value. A stable value is not attained, whereas this requirement is nicely fulfilled for the larger APS-6.0 sample.



View larger version (13K):
[in this window]
[in a new window]
  		Fig. 11. Representative elementary volume analysis for linear attenuation coefficient of the two different samples.

 
We would like to offer a related explanation for the lack of effect of applied boundary condition for the smaller sample. Because the sample has such a small diameter, it is possible that the lack of tortuosity, or pore connectivity, which we hypothesized to influence the results for the larger sample, plays a role here. Perhaps there is simply not enough space for the three-dimensional component of the pore drainage network to develop, and therefore, no effective short-circuiting pathways are established from top to bottom of the sample. In comparison to our three-dimensional experiments on the APS-6.0 sample, and the two-dimensional experiments of Mortensen et al. (2001), we are observing a flow process that is more or less restricted to one dimension in this very small sample. It is also possible that a wall effect exists for the smaller sample. This effect would become increasingly important as the sample size decreases because irrespective of the applied pressure, the large pores at the wall control drainage.

Spatial Analyses
The nearest neighbor distance is calculated for the distribution of air in the sample following the two different drainage processes, thus the binary images shown in Fig. 12 depict only the air-phase present in the drained samples. Keep in mind that segmentation of gray-scale images is not an exact science, Russ (1999) and includes a certain level of subjectivity. The thresholds were selected to optimally represent the air-phase structures observed in the gray-scale images. This becomes a trade off between eliminating too many blobs that originate in noise and eroding away some that are actually small air bubbles. As seen in Fig. 12b in particular, there are a few air blobs that are connected in the gray-scale images, which have been disconnected in the segmented images. Also, some air blobs were very hard to detect by our routine because the image slice cuts through the very top or bottom of an air blob, and therefore it appears more gray than black in the images. To the human eye, these might still be easily distinguishable as an air bubble, however, an automated computer algorithm does not "see" it that way. We feel the values chosen are the best ones that can be obtained without actually including a reference material in the image (implemented experimentally) to scale all other gray-scale values to. More importantly, for the NND analyses, the most important criterion is to preserve the spatial location of all air bubbles, not necessarily their exact size.



View larger version (100K):
[in this window]
[in a new window]
  		Fig. 12. Gray-scale images and corresponding binary representation of the air-phase distribution for a horizontal slice of APS-6.0 following (a) fast and (b) slow drainage. The analyzed sections are 230 by 230 pixels corresponding to 3.91 by 3.91 mm2.

 
The results of the nearest neighbor analysis described earlier are shown in Fig. 13 for the three slices (40, 110, and 160), and the two flow conditions for the APS-6.0 sample. Obvious differences between fast and slow drainage can be readily observed in this figure. Under fast drainage the air bubbles tend to be few and far apart, whereas vastly more bubbles, that are more closely spaced, are obtained under slow-flow conditions. For instance, the fast-drained image shown in Fig. 12a (slice 110) has an average NND of 1672 pixels, corresponding to a distance of approximately 0.70 mm, whereas the slow-drained image in Fig. 12b has an average NND of only 372 pixels, corresponding to approximately 0.33 mm. The large average NND of the fast-drained image confirms the fact that air bubbles are far apart under this flow condition compared to slow drainage where air is introduced in pores (i.e., water is drained) in a much more uniform manner. The same pattern was obtained (not shown here for brevity) when considering the entire sample height (all slices in a volume) for the two flow conditions; large distances between air bubbles or air blobs were seen for the fast drainage process, whereas very short distances were observed for slow drainage, confirming our qualitative results presented earlier and shown for instance in Fig. 6a and 6b.



View larger version (22K):
[in this window]
[in a new window]
  		Fig. 13. Nearest neighbor distances for three slices (40, 110, and 160) for the two drainage conditions for APS-6.0.

 

   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS
 CONCLUSIONS
REFERENCES
 
A previous publication by the authors (Wildenschild et al., 2001) left us with some unanswered questions regarding the cause of observed "dynamic" flow phenomena. To address these questions further, minute pore-scale multiphase flow experiments were developed to nondestructively visualize the flow process in a sample of porous material. The flow experiments were performed concurrently with collection of high resolution tomographic images. We used Lincoln sand which had shown rate dependent flow effects in our previous research. The samples were exposed to similar flow rates and boundary conditions as in the previous experiments; one-step and multistep drainage experiments where the water was drained from the sample using either one large or several small pressure steps, respectively.

Qualitative examination of the obtained microtomographic images shows that the difference in drainage pattern, and resulting residual water phase saturation, for the two flow conditions is quite dramatic. For the high initial flow rates a couple of processes could be responsible for the observed drainage patterns:

  1. Drainage taking place in the most well-connected and largest pores first, meaning that subsequent drainage of the finer pores is unlikely because they are no longer hydraulically connected to the main flowing continuum; or
  2. Some of the smaller pores are drained the instant they are exposed to the high air-phase pressure, particularly near the outflow boundary, leaving many relatively large pores disconnected from other flowing pores, and therefore undrained.

On closer inspection of the images, we found that Mechanism 1 appeared to be dominant.

To support this observation we also performed quantitative analyses on the microtomographic data. A nearest neighbor analysis was performed on the air-phase distribution found in the images and showed characteristic differences between fast and slow drainage. Under fast drainage the air bubbles tend to be few and far apart, whereas vastly more and smaller bubbles, that are more closely spaced, are obtained under slow flow conditions. This corroborates our qualitative findings regarding the preferential drainage of larger pores as the mechanism for drainage at the higher initial flow rates. As mentioned previously, this means that mechanism A of Wildenschild et al. (2001) is indeed the most likely to be responsible for the dynamic effects that were observed at high flow rates. During the fast-drainage one-step experiments we found that some of the larger pores appeared to expand as grains moved slightly under the high pressure of the invading air, and thus altered the pore geometry. This appeared to be limited to the very largest pores, yet, it adds to our concern about using the one-step technique for measuring capillary pressure–saturation curves. The one-step technique was widely used for a number of years, but most researchers have now switched to the less problematic multistep technique where the sample is drained slowly, and capillary forces hold the grains together once the higher air-phase pressures are applied. Hopmans et al. (1992) pointed out the problems with starting one-step experiments at full saturation. Our experiments were not started at complete saturation, as the images show, yet, the initial capillary pressure might not have been large enough to prevent the slight movement of small grains.

In addition to these findings, the work presented here illustrates the imaging details possible with the X-ray microtomography technique, and thus offers requisite information about the possibilities of X-ray microtomography as an analytical tool for addressing multiphase flow problems.


   ACKNOWLEDGMENTS
 
The Danish Technical Research Council is gratefully acknowledged for financial support. Part of this research was carried out while two of the authors (Wildenschild and Kent) were employed at Lawrence Livermore National Laboratory and thus part of this work was performed under the auspices of the U.S. Dep. of Energy by Lawrence Livermore National Laboratory under contract No. W-7405-ENG-48 and supported by the Environmental Management Science Program. Use of the Advanced Photon Source was supported by the DOE Basic Energy Sciences, Office of Science, under Contract No. W-31-109-ENG-38.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
 Soil Sci.Home page
L. G. Tumlinson, H. Liu, W. K. Silk, and J. W. Hopmans
Thermal Neutron Computed Tomography of Soil Water and Plant Roots
Soil Sci. Soc. Am. J., September 1, 2008; 72(5): 1234 - 1242.
[Abstract] [Full Text] [PDF]

Home page
 Soil Sci.Home page
A. Tuli, J. W. Hopmans, D. E. Rolston, and P. Moldrup
Comparison of Air and Water Permeability between Disturbed and Undisturbed Soils
Soil Sci. Soc. Am. J., August 4, 2005; 69(5): 1361 - 1371.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (11)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Wildenschild, D.
Right arrow Articles by Kent, A. J. R.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Wildenschild, D.
Right arrow Articles by Kent, A. J. R.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Wildenschild, D.
Right arrow Articles by Kent, A. J. R.
Related Collections
Right arrow Laboratory Column Studies
Right arrow Pore-Scale Modeling
Right arrow Experiment Design

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome