VZJ sign up for citetrack
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


Published online 26 January 2006
Published in Vadose Zone J 5:80-97 (2006)
DOI: 10.2136/vzj2004.0177
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
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 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 Web of Science (4)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Lehmann, P.
Right arrow Articles by Flühler, H.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Lehmann, P.
Right arrow Articles by Flühler, H.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Lehmann, P.
Right arrow Articles by Flühler, H.
Related Collections
Right arrow Tomography
Right arrow Pore-Scale Modeling
Right arrow Multiphase Fluid Flow

ORIGINAL RESEARCH

Tomographical Imaging and Mathematical Description of Porous Media Used for the Prediction of Fluid Distribution

P. Lehmann*,a, P. Wyssb, A. Flischb, E. Lehmannc, P. Vontobelc, M. Krafczykd, A. Kaestnera, F. Beckmanne, A. Gygia and H. Flühlera

a Institute of Terrestrial Ecology, Swiss Federal Institute of Technology, ETH Zurich, Switzerland
b Centre for Non-destructive Testing, Swiss Federal Laboratories for Materials Testing and Research, EMPA, Switzerland
c Spallation Neutron Source Division, Paul Scherrer Institute, PSI, Switzerland
d Institut für Computeranwendungen im Bauingenieurwesen, TU Braunschweig, Germany
5 GKSS-Research Centre, Geesthacht, Germany

* Corresponding author (peter.lehmann{at}env.ethz.ch)

Received 10 December 2004.



    ABSTRACT
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MEASUREMENTS WITH COMPUTER...
 DESCRIPTION OF SIMULATED MEDIA
 GEOMETRICAL PROPERTIES AS A...
 APPLICATION: FLUID PHASE...
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Flow and transport processes in porous media depend on the geometric properties of their pores, where the diameters typically range from micrometers to millimeters. In this study, we mapped the pore structure of glass bead and sand columns using tomography with X-rays, thermal neutrons, and synchrotron radiation. Utilizing X-rays from tubes, we mapped two 2.5- and 5.3-cm-diam. sand samples that contained particles with sizes ranging from 0.08 to 1.25 mm. The resulting voxels (i.e., the unit of a three-dimensional image, the smallest distinguishable box-shaped part of a three-dimensional image) were cubes of 60-µm length in the case of microfocus X-rays, and 70 µm in case of industrial X-rays. In the latter case, each voxel represented the material density of a rectangular parallelepiped with side lengths of 70 µm and a height of 210 µm. The material density of cubes of 70 µm was reconstructed by applying an optimized filter in Fourier space. Columns with diameters of 4.0 and 5.3 cm containing glass beads with diameters of 3.0 and 2.0 mm were scanned with thermal neutrons. The voxel size was 167 µm. Because this technique is sensitive to the presence of water, it was possible to measure the water table in a partially water-filled sample. Two sand columns were scanned with synchrotron X-rays, and the resulting voxel sizes were 11.5 and 3.5 µm. In the first case, the sample with a diameter of 15 mm contained particles of sizes ranging from 300 to 900 µm. In the second case, a sample with a diameter of 5 mm was filled with 100- to 200-µm particles. In a numerical analysis of the sphere packings, we computed various geometric properties of the porous media as a function of the resolution. The pore-size distribution and the Minkowski functionals (quantities that define the morphology of a structure) were used to describe changes in the imaged pore space as a function of voxel size. We found that the geometric properties of the mapped pore space converged to true values for a voxel size of 10 to 20% of the mean particle radius. Based on this analysis, we postulate that the resolution of a tomographic measurement must be in the range of 10% of the mean particle radius for repacked media to reconstruct the characteristic features of the pore space. This condition was fulfilled for the tomography with synchrotron light. Using the images of the sand samples measured with synchrotron light, we predicted the amount of water and air for a drainage process. For the pore space mapped with tube X-rays, it was possible to make qualitative predictions of the hysteretic water and air distribution.

Abbreviations: EPC, Euler-Poincaré characteristic • HASYLAB, Hamburger Synchrotron Laboratories, Germany • SLS, Swiss Light Source, Paul Scherrer Institute, Switzerland


    INTRODUCTION
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MEASUREMENTS WITH COMPUTER...
 DESCRIPTION OF SIMULATED MEDIA
 GEOMETRICAL PROPERTIES AS A...
 APPLICATION: FLUID PHASE...
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
FLOW AND TRANSPORT PROCESSES in soils or other porous media depend on their structure, that is, on the spatial arrangement of the soil constituents. Structure controls the storage of water, the availability of dissolved essential elements, aeration, and the transport time of hazardous substances. To understand flow of water and air, or to predict the transport of dissolved pollutants and nutrients, the interaction between soil structure and transport processes must be investigated. To do this, we face a scale dilemma. Ultimately, we are interested in processes on scales ranging from meters to kilometers, such as the transport of liquid manure constituents from the soil surface to groundwater, the dynamics of heavy metals or pesticides in the topsoil, or the amount of available water for crops within an agriculturally used area, but these processes are dominated by the geometry of microscopic structures, in particular the shape and size of the pores. Pores form a network with openings ranging from less than one to a few thousand micrometers. Air, water, and dissolved substances move through this network of fine pores. Their mobility depends on the width and connectivity of the pore space. Connectivity is probably the most pertinent property in this sense and, at the same time, the most difficult to quantify. Pore size has various effects on transport processes: the mean velocity of water and dissolved solutes in a pore increases with pore radius to the power of two and the local flow velocity in pores grows in proportion to the square of the distance from the solid wall (Hagen-Poiseuille Law). The residence time of reactive substances is determined by the ratio of the sorbing surfaces to the volume of soil solution, a ratio that hyperbolically decreases with increasing pore size. The same functional dependence exists between pore size and water retention, so the energy needed to drain water from soil increases with decreasing size of the pores (i.e., Young–Laplace equation).

To understand transport processes in porous media, we need to determine the pore structure with a voxel size of a few micrometers. To assess the size of the pores and the pore-size distribution, indirect methods are often used: the solid phase is complementary to the pore space, and some soil properties can be estimated using the measured particle-size distribution (Arya and Paris, 1981; Arya et al., 1999). However, the structure of the porous medium does not depend on the size of the particles alone. Structure is the result of soil formation processes. A living organism, for example, may form a burrow or may produce organic substances that act as cementing agents forming soil aggregates. Such structures cannot be derived or deduced from the sizes of the particles. Another indirect approach to determine the pore sizes is the measurement of the water retention curve. To drain a porous medium, a pressure difference is applied between soil water and air. The water is retained in pores that are small enough to withstand the applied suction. However, this interpretation is partially misleading because a large pore surrounded by small pores will not drain until the suction is large enough to drain the small pores. Therefore, with the applied pressure difference, only the size of the smallest pores will be detected. Vogel and Roth (1998) used a network of capillaries and investigated the effect of network topology on soil water characteristic and hydraulic conductivity to show the importance of the connectivity. Lehmann (2002) showed that the interpretation of the water retention curve as pore-size distribution to predict the hydraulic conductivity holds for soils with a wide pore- and particle-size distribution. But the prediction based on the derivation of the pore-size distribution from the water retention curve failed for soils with a narrow particle-size distribution (sand). For media with a narrow size distribution, the soil water characteristic is strongly affected by the connectivity. Also, the wetting process is affected by network topology because a tiny pore surrounded by large pores will not be refilled before the neighboring large pores are saturated. The soil water retention function may be hysteretic because small structures dominate the drainage and large pores control the wetting process (Mualem, 1974; Stauffer, 1996; Lehmann et al., 1998).

The size and shape of pores can be measured directly by thin-cut analysis with a pixel size of a few micrometers. While it is possible to deduce the connectivity of the pore space using pairs of parallel images (DeHoff et al., 1971), so-called dissectors (Sterio, 1984; Vogel and Kretzschmar, 1996), the resolution in the third dimension of consecutive thin-cuts is limited. To minimize the distance between parallel images, the sample must be impregnated with a resin. This method is time consuming, destructive and limited in the resolution of the third dimension.

To measure the three-dimensional pore network nondestructively with the same resolution in all three spatial dimensions, we must use tomographic techniques. The optimal technique depends on the size of the sample, the required resolution, and the most relevant structures of the sample. So, it makes a difference whether we are interested in the density distribution of the whole sample, the distribution of the water phase alone, or the spatial distribution of certain constituents. X-rays were used to scan the density distribution in soil columns (Hopmans et al., 1994; Tidwell et al., 2000) and to map large pores that are relevant for fast flow processes (Heijs et al., 1995; Vanderborght et al., 2002). Normally, the beam is generated with a device known as an X-ray tube, where a metallic target is bombarded with electrons from a cathode and X-rays are released. Alternatively, very intense X-rays can be emitted by very fast electrons traveling through a magnetic field. To distinguish between this synchrotron radiation and X-rays from tubes, we use the terms X-rays from tubes and X-rays from synchrotrons. The advantage of the synchrotron radiation is the wide range of the wavelength, ranging from ultraviolet light to hard X-rays, and the very intense photon flux density. To scan samples with a resolution of 1 µm or less, a high photon flux per micrometer is needed to limit the exposure time. A comparison between tomography with micro-X-ray and synchrotron radiation was given by Bernhardt et al. (2004).

To analyze fine details of the pore structure of small samples, X-rays from synchrotrons were used to measure the distribution of pore sizes in sandstone (Lindquist et al., 2000) and the porosity in chalk stone (Betson et al., 2004). To determine variations in water saturation, neutron tomography can be applied (Solymar et al., 2003). Also the nuclear magnetic resonance technique is sensitive to the amount of water in the porous media and can be used to measure the size distribution of water filled pores (Bird et al., 2005). In this study, we focus on tomography with X-rays from tubes and synchrotrons and thermal neutrons to measure the size and the arrangement of pores.

In soil physics, we want to understand and predict flow phenomena and transport processes for partially water saturated conditions. Our ultimate goal is the prediction of transport processes for a given soil type based on the information about its structure, that is the spatial arrangement of the pores and particles. In this study, we measured the structure of various porous media and calculated the water and air distribution. The manuscript has three key aspects:

  1. To make successful predictions, we need a map of the pore space that is representative of the porous material. We first describe the tomography of porous media that were scanned at different measurement facilities. In total, four sand samples and two packings of glass beads were measured.
  2. The resolution of this map must be fine enough to show all the relevant geometric properties of the pore network. To analyze the effect of resolution, we generated numerically porous media. We modeled two materials that were identical with two of the scanned materials with respect to porosity and particle size. We describe the numerical generation of the porous media, and we introduce the geometrical properties that will be analyzed. Next we discuss the properties of the images of the numerically generated porous media as a function of resolution.
  3. Finally, images of the tomographed samples are used to determine the distribution of water and air for different water unsaturated conditions. For small images we used a Lattice–Boltzmann method. For large systems we applied a pore network approach.


    MEASUREMENTS WITH COMPUTER TOMOGRAPHY
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MEASUREMENTS WITH COMPUTER...
 DESCRIPTION OF SIMULATED MEDIA
 GEOMETRICAL PROPERTIES AS A...
 APPLICATION: FLUID PHASE...
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
With different irradiation techniques we scanned columns containing glass beads or sand material. Based on the transmitted intensity of the X-rays or thermal neutrons, the three-dimensional material distribution can be calculated using backprojection algorithms. The reconstructed values do not depend on the property of the sample alone but also on the presence of random noise in the measurement (variations of the incoming beam or damaged detector cells). The artifacts resulting from these errors should be filtered before the image is analyzed with respect to the geometrical properties of the pore space. The principles of tomography and reconstruction are explained in Kak and Slaney (2001). The reconstructed pore structures of four samples are shown in Fig. 1 .


Figure 1
View larger version (78K):
[in this window]
[in a new window]
 
Fig. 1. Reconstructed geometries of scanned samples. Sand samples were scanned using X-rays from (A) tubes or (B) synchrotrons with a voxel size of 60 and 3.5 µm and a sample diameter of 25 and 5 mm, respectively. The glass bead columns with diameters of (C) 52 mm and (D) 40 mm, respectively, were scanned using thermal neutrons with a voxel size of 167 µm. In the latter case, the pore space was water filled at the bottom of the column.

 
For X-rays from tubes, thermal neutrons, and synchrotron X-rays, the attenuation of the beam by the air phase is negligible. The attenuation in dry samples is therefore proportional to the local amount and density of the solid phase. From the reconstructed three-dimensional density map with gray values proportional to the density, a segmented map of pore and solid voxels must be deduced. This classification was performed based on the frequency distribution (histogram) of the measured gray values. The shape of the histogram and the density map depend on the resolution of the tomography technique. In the best case, the histogram consists of two nonoverlapping peaks representing the pores and the solid phase, respectively. Then, a threshold with a value between the two modes can be applied to classify the voxels as pores or solids. In practice, the histogram is more ambiguous because of the existence of small pores and the influence of noise. For such a histogram, Vogel and Kretzschmar (1996) applied two thresholds to segment the image into pores, solids, and undefined voxels. The undefined voxels with a gray value between the two thresholds were determined as pores or solids in an iterative process: the undefined voxels in the neighborhood of a pore voxel were determined as pores. This process was repeated until no more gray voxels were found in the neighborhood of a pore. The remaining undefined voxels were determined as solids. Alternatively, in the case of media with well-known porosities that were scanned with a poor resolution (indicated by a histogram with a single peak), a simple threshold can be chosen to segment the image. Voxels with gray values smaller than a threshold can be designated as pores. The threshold is given by the condition that the fraction of voxels with gray values that are smaller than the threshold must correspond to the porosity.

Tomography of Sand Material with X-rays from Tubes
The porosity and particle-size distribution of the sand was determined in the laboratory. The measured porosity was 0.37 m3 m–3. The masses of the following particle-size fractions were measured with sieves: 80 to 125, 125 to 250, 250 to 400, 400 to 625, 625 to 800, 800 to 1000, and 1000 to 1250 µm. The mean radius of the measured particle-size distribution was 263 µm. The pore structure was investigated using a micro-X-ray scanner. A sand column 2.5 cm in diameter and 0.9 cm in height was irradiated with a cone beam of energy 140 keV maximum, filtered with 1-mm aluminum. The voxel size of the resulting image was 60 µm. The pore structure of the sand sample after segmentation is shown in Fig. 1A. A threshold was applied to classify the voxels as pores or particles with the constraint of the measured porosity. To scan a large (5.3-cm diam., 7-cm height) column of the same sand material, we used an industrial tomographic scanner with a 450-kV X-ray tube (focal spot size of 2.5 mm). The beam is collimated to a fan of 2-mm height. The height of the scanned sector was 5.3 cm, corresponding to 750 scans with a spacing of 70 µm. The resolution was 70 by 70 by 210 µm3, because 210 µm was the minimum opening of the collimator to receive a signal above the noise level. To enhance the resolution in the z-direction, the vertical increments between two measured planes were 70 µm. Thus, the resulting image had a grid spacing of 70 by 70 by 70 µm3, but the gray value of a voxel depended on the material density in a cuboid of size 70 by 70 by 210 µm3. Before we tried to calculate the material density in a cube of size 70 by 70 by 70 µm3, the gray-scale value images had to be filtered due to a bad signal/noise ratio. Each voxel is affected by the material in the plane below and in the plane above. We faced the problem that for the first and the last of the scanned planes, the measurements were influenced by neighboring planes whose properties are unknown. For the reconstruction of N planes, a system with N equations with N + 2 unknowns must be solved. For each position x, y of the image, a sequence of gray values as a function of the plane in height z was obtained with tomography. This sequence is denoted as g(z). From g(z), we want to derive the unknown gray values in cubic voxels, v(z). The measured gray value may thus be described as g(z) = h(–1)v(z – 1) + h(0)v(z) + h(1)v(z + 1), with a function h that defines the contribution of the planes to the measured gray value. In case of industrial X-ray tomography, the contribution of the three planes to the detected signal was equal; therefore, h(–1) = h(0) = h(1) = 1.0. In mathematical terms, this is a convolution that can be described as a multiplication after the Fourier transformation of the functions, G({lambda}) = V({lambda}) H({lambda}), with the space frequency {lambda}.

The goal was to find an optimized response function h* with its Fourier transform H*({lambda}), then to calculate the inverse Fourier transform of G({lambda})/H*({lambda}) and to use this function as an estimate for the unknown v(z). The response function h*({tau}), where {tau} is the shift with respect to the position z, is non-zero for h*(0), h*(1), and h*(N – 1). The argument N – 1 equals a shift of –1. We tested different functions h*({tau}), transformed them into the Fourier domain, divided G({lambda}) with H*({lambda}), and calculated the inverse transformation of G({lambda})/H*({lambda}). With an adequate choice of H*({lambda}), the effect of the unknown boundary values can be filtered out, and a valuable estimate for the unknown function v(z) can be calculated. Lehmann et al. (2001) estimated the function h*({tau}) by analyzing a sequence of random numbers v(z). Here, we used numerically generated sand media with the same particle-size distribution and porosity as the real one. This procedure is specified in the section "Sand Material with Wide Particle-Size Distribution." Cross sections through the images are shown in Fig. 2 .


Figure 2
View larger version (144K):
[in this window]
[in a new window]
 
Fig. 2. In the case of industrial X-ray tomography, the voxel size was 70 by 70 by 70 µm3. But the gray value of each voxel corresponded to the material density in a cuboid of size 70 by 70 by 210 µm3. We reconstructed the material density in cubic voxels of size 70 by 70 by 70 µm3 using the Fourier transformation technique. In the first row, vertical cross sections through the scanned sand sample are shown (A) before and (B) after noise filtering and inverse Fourier transformation. In the last image (C), the image is segmented into solid phase (black) and pore space (white). In the second row, the calculated density distribution in a vertical cross section of a simulated sample is shown. In the first image (D), the gray values were calculated for cubes of size 70 by 70 by 70 µm3. In the next image (E), the density in cuboids of size 70 by 70 by 210 µm3 is given. Due to the anisotropic mass contribution, the shape of the spheres is slightly distorted in the vertical direction. (F) From this asymmetric density distribution, the material density in cubes was calculated using Fourier transformation. The cross section shown in Part F is similar to the "true" structure shown in D. The same technique was applied for the measured data shown above. The numbers denote the volume that determined the gray values of the voxels.

 
Scanning of Glass Bead Columns with Thermal Neutrons
A column (5.3-cm diam., 9-cm height) containing glass beads with diameter 2 mm was irradiated with thermal neutrons at the NEUTRA radiography facility of the Paul Scherrer Institute, Switzerland. Using a slow scan Charge Coupled Device (CCD) with a 512 by 512 pixel chip, 200 projections over 180° were acquired. The resolution was 167 µm. The histogram was ambiguous because it had only one peak for the solid phase and sort of a plateau for the voxels that were (partially) in the pore space. In a first step we applied a threshold to make a preliminary segmentation. In a second step, the distance to the next pore was calculated for each voxel of the solid phase. Voxels with a distance smaller than the radius of the glass beads were assigned as pores. After this process, only the cores of the glass beads were still identified as solids. Finally, a solid sphere of 2-mm diameter was inserted in the center of each core of the solid phase. The reconstructed structure is shown in Fig. 1C.

Tomography of Sand Material using Synchrotron Radiation
The synchrotron X-rays of an electron storage ring (Swiss Light Source [SLS], Paul Scherrer Institute, Switzerland) and a positron storage ring (Hamburger Synchrotron Laboratories, HASYLAB, Germany) were used to investigate the pore space of sand samples in detail. Tomographic scans of a sand sample (1.5-cm diam.) containing 300- to 900-µm sand particles were performed at beamline W2 of the HASYLAB at the Deutsches Elektronen-Synchrotron DESY using a photon energy of about 42 keV and a voxel size of about 11.5 µm. At SLS a sand sample (0.5-cm diam.) containing 100- to 200-µm particles was scanned with a monochromatic energy of 16 keV, with a resulting voxel size of 3.5 µm. To classify the voxels as pores or solid phase, two thresholds were applied. Voxels with gray values smaller than the first were designated as pores, and voxels with gray values higher than the second threshold were designated as solid phase. In a second step, voxels with a gray value between the two thresholds were classified as a pore if one of the six neighboring voxels was already designated as pore. This second step was repeated until no further voxels with intermediate gray values were detected in the neighborhood of a pore. The remaining gray voxels were designated as solids. A section of the sample measured at SLS is presented in Fig. 1B.


    DESCRIPTION OF SIMULATED MEDIA
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MEASUREMENTS WITH COMPUTER...
 DESCRIPTION OF SIMULATED MEDIA
 GEOMETRICAL PROPERTIES AS A...
 APPLICATION: FLUID PHASE...
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
To test the dependence of the reconstructed image on voxel size, we numerically generated two porous media. The porosity and the particle sizes in the modeled and scanned media were equal. In the first case, we generated a medium that represents the sand sample scanned with X-rays from tubes. With this simulated sand we intend to analyze the pore space of a medium with a wide particle size distribution. We denote this model as the sand model. In the second case, a model of the glass beads scanned with thermal neutrons was generated. With this model, denoted as the glass model, we investigated the pore space geometry of a medium with a single particle size.

Numerical Generation of Porous Media
Medium with Particles of Different Sizes
The particle diameters ranged from 80 to 1250 µm with a mean radius of 263 µm. The mass of seven size fractions was determined in the laboratory by sieving. We summarized these mass fractions to get the measured cumulative particle size distribution function (CDF). We chose a lognormal distribution m(r) as a model for the particle size distribution as a function of the particle radius r. The cumulative distribution of this analytical function m(r) (i.e., an error function with the logarithm of the particle size in the argument) was fitted to the measured CDF. The fitted lognormal distribution m(r) was divided by the mass of a spherical particle {rho}4/3{pi}r3 to obtain the number of particles per particle size p(r), with the constant density of the particles {rho}. The function p(r) was integrated over the particle size r and normalized to obtain values between 0 and 1.0. This integral, in the following denoted as P(r), determines the probability of finding a particle of size equal to or smaller than r. To determine the size of the modeled particles, a random number between 0 and 1 was chosen. This random number n corresponds to P(r) for a certain particle radius r and a spherical particle with radius r was generated. This procedure was repeated until the volume of the generated particles was equal to the measured volume fraction of the solid phase (1 – {phi}) L3, with the measured porosity {phi} = 0.37 m3 m–3 and the size of the modeled system L. To model the media analyzed with micro-X-ray tomography, 15 139 spherical particles were randomly arranged in a cube of size L = 6.25 mm. In the case of the industrial X-ray tomography, 99 870 spheres were distributed in a cube of size L = 12.5 mm. In a second step of the numerical packing procedure, the particles were distributed randomly in space, starting with the largest particles. If one particle overlapped another, then a new random position was chosen. In literature, such an approach is denoted as a random-parking algorithm (Cooper, 1988).

Medium with Single Particle Size
In the case of a single particle size, media with small porosities with nonoverlapping particles cannot be generated with a parking algorithm. The porosity of the glass bead packing was 0.40 m3 m–3. The code we used was based on the description of the algorithm given by Jodrey and Tory (1981). The spherical particles were randomly distributed and the overlapping of the particles was allowed in a first step. In a next step, all overlapping particles were rearranged by moving the particles in a direction opposite to the center of the nearest neighbor. This procedure of redistribution was repeated until a system with nonoverlapping particles was generated. To optimize the redistribution of overlapping particles, the size of the spheres was slightly enlarged for the initial distribution of overlapping spheres. The initial size was 100.2% of the true size. When the redistribution of the particles became less effective with respect to overlapping, the size of the particles was decreased. To model the tomography of the glass bead columns, 9168 spheres of 1-mm radius were distributed in a cubic box of size 40 mm.

With this method, packings with higher densities can be generated compared with the algorithms described in the previous subsection. However, for a porosity of 0.37 m3 m–3 and the wide particle size distribution of the sand, the random parking algorithm was more efficient.

Analyzed Properties of Simulated Media
We mapped the generated structure onto a grid. The grid spacing corresponded to the voxel size of the resulting image. For each voxel of the medium we computed the volume fraction of the solid phase and determined a gray-scale value between 0 (i.e., voxel contains no solid) and 255 (i.e., voxel contains no voids). Then, a threshold was applied to obtain an image of the pores with the measured porosity. An example of the resulting density image (gray) and the corresponding image of the pore space (black and white) is given in Fig. 3 . We calculated the density map and the image of the pore space and determined the following quantities.


Figure 3
View larger version (61K):
[in this window]
[in a new window]
 
Fig. 3. In the first row, a gray-scale value image of the material density (left) and the corresponding black (particles) and white (pores) image is shown. The figures show a sector of the sand model with a side length of 50 voxels and voxel size of 60 µm. This voxel size equals the resolution of the measurement with micro-X-ray tomography. Not all voxels determined as pores are entirely within the pore space (gray value 0). The fraction of these gray voxels is given in the figure in the second row. This fraction decreases with decreasing voxel size. The large symbols denote the resolution of the tomography with X-rays from tubes (black) or thermal neutrons (gray), respectively. The positions of the arrows correspond to the resolution with 50% of voxels denoted as pores with a gray value of 0. For the tomography of the glass bead column, this requirement was fulfilled.

 
Fraction of Gray Pore Voxels
This is the fraction of voxels determined as pores that are not completely within the pore space and the mapped density value is higher than 0. In the real media, flow and transport in an element partially filled with solids is different from the processes in an element entirely in the pore space. So, these voxels may cause an error in the predictions based on the segmented image.

Minkowski Functionals
To determine the morphology of a structure, the four Minkowski functionals of the solid phase were calculated. They include the volume, surface, integral mean curvature (or mean breadth), and the Euler characteristic. The Minkowski functionals characterize the morphology of disordered media (Arns et al., 2002). The mean breadth is proportional to the integral mean curvature and depends on the curvature of the object. If the center of the curvature is within the analyzed object, the curvature is positive. In the opposite case, the value is negative. For isolated convex objects, the mean breadth corresponds to the mean diameter. The integral mean curvature equals the mean breadth multiplied with 2{pi}. For a complex network containing convex structures, the mean breadth becomes negative. The Euler characteristic, also known as the Euler-Poincaré characteristic (EPC), quantifies the connectivity of a structure. It is negative for highly connected networks with redundant connections and positive for a set of isolated elements. The EPC for the solid phase is high for spheres that are completely separated and pores that are embedded within the solid phase. In the case of an imaged structure, the structure consists of cubic voxels. Each voxel consist of one cubic body, six faces, 12 edges, and eight vertices. The Minkowski functionals of an imaged structure can be obtained by counting the total number of cubes nc, faces nf, edges ne, and vertices nv according to the following formula (Michielsen and De Raedt, 2001):

Volume = nc
Surface area = –6nc + 2nf
Mean breadth = 1/2(3nc – 2nf + ne)
Euler characteristic = –nc + nfne + nv.

If two neighboring voxels have a face, an edge, or a vertex in common, this element is counted only once. For the surface area (and the mean breadth) of spheres this formula is not exact because it is impossible to describe the smooth surface of spherical solid particles by counting the quadratic faces of the voxels. Each voxel of the solid phase has six faces, and each of the six faces is the boundary between two voxels. If one of the neighbored voxels belongs to the pore space, this face with an area of {Delta}L2 is an element of the surface of the solid phase, where {Delta}L is the length of a voxel. A more accurate estimation of the area was computed with a method described in Ohser and Mücklich (2002, p. 114–121). For each voxel of the solid phase at position x, y, z they included 26 neighbors at position x + {Delta}x, y + {Delta}y, z + {Delta}z, with the values –1, 0, or 1 for {Delta}x, {Delta}y, and {Delta}z. If the voxel at a position x + {Delta}x, y + {Delta}y, z + {Delta}z belongs to the pore space, an interface of wi is added to the total interface of the solid phase. The coefficients wi depend on the distance between the two voxels and are equal to 0.183{Delta}L2 for w1 with {Delta}x2 + {Delta}y2 + {Delta}z2 = 1, 0.105{Delta}L2 for w2 with {Delta}x2 + {Delta}y2 + {Delta}z2 = 2 and 0.081{Delta}L2 for w3 with {Delta}x2 + {Delta}y2 + {Delta}z2 = 3. For the mean breadth, the method is similar but with {Delta}L instead of {Delta}L2 and with the coefficients divided by 2. Because we generated the media with spherical particles, we will use this more accurate estimation.

Pore-Size Distribution
For each voxel determined as pore, the radius of a sphere that fits entirely into the pore space was calculated. The pore voxel is the center of the sphere. First, the radius of this sphere was attributed to the pore voxel. This result equals the distance transform (Soille, 2003), determining the minimum distance of each pore voxel from the solid wall. However, a pore voxel can be an element of a larger sphere. In such a case, the radius of the larger sphere is assigned to this voxel. The fraction of pore voxels that belongs to the different sphere sizes was computed. Another method to compute the pore size distribution was applied by Vogel (1997). He calculated the pore size distribution with a morphological filtering denoted as opening (Serra, 1982) without the preceding determination of the distance transform.


    GEOMETRICAL PROPERTIES AS A FUNCTION OF RESOLUTION
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MEASUREMENTS WITH COMPUTER...
 DESCRIPTION OF SIMULATED MEDIA
 GEOMETRICAL PROPERTIES AS A...
 APPLICATION: FLUID PHASE...
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
We calculated the material density distribution in the generated media as a function of resolution and analyzed the geometrical properties of the segmented image. To segment the image into a map of pore space and solid phase, a threshold was applied that fulfills the requirement of the porosity. In practice, the reconstructed gray values do not depend on the fraction of solid phase alone. The attenuation depends also on the variation of the mineralogical composition and on the errors of the measurement. For this reason, the segmentation of real data is more complicated and other approaches to classify the voxels may be more appropriate, but if the segmented image should fulfill the constraint of the measured porosity, a single threshold can be applied.

Sand Material with Wide Particle-Size Distribution
The density distribution in the generated media was calculated for voxel sizes ranging from 6.25 to 125 µm, corresponding to a grid with 10003 to 503 voxels. To compare the resolutions for media with different particle sizes, we define a dimensionless resolution R by dividing the voxel size with the mean particle radius.

The fraction of the solid phase was calculated for each voxel resulting in a gray value between 0 and 255. The image was segmented into a map of pore space and solid phase with the required porosity of 0.37 m3 m–3. Cross sections through the density map and the pore space image are shown in Fig. 4 for different resolutions.


Figure 4
View larger version (123K):
[in this window]
[in a new window]
 
Fig. 4. Cross sections through the density map (top row) and the segmented image (middle row) for the simulated sand media (left) and numerically generated glass bead packings (right). The numbers in the images denote the voxel sizes. In the first two rows, the images calculated with the voxel size of the tomography with X-rays from tubes or thermal neutrons, respectively, are shown. The last row gives the density map calculated with the highest resolution (10003 voxels).

 
The fraction of voxels determined as pores with a non-zero gray value increased with increasing voxel size from 11% (voxel size = 6.25 µm, R = 0.02) to 99% (voxel size = 125 µm, R = 0.48). For a voxel size that is equal to the resolution of the tomography with micro-X-ray (60 µm, R = 0.23), only 13% of the voxels determined as pores have a density value of 0. Hence, 87% of the voxels of the pore space contained a certain amount of the solid phase (Fig. 3). For resolutions better than 30 µm (R = 0.11) and 10 µm (R = 0.04), the fraction of pores with nonzero gray values was <50% and <20%, respectively.

Next, we analyzed the Minkowski functionals for different resolutions. Due to the condition that the fraction of pore space must correspond to the measured porosity, the volume of the solids is constant. The surface of the spheres in the image is compared with the analytically determined surface of all spheres. The surface increases with decreasing voxel size (Fig. 5A ) and converges to the true value. In Fig. 5B, the change of surface area with voxel size is given. For voxel sizes <30 µm (R = 0.11), the increase of surface with decreasing voxel size becomes less negative, indicating the convergence to the true value. The calculated value of the mean breadth was divided by the analytically determined value of the spheres (Fig. 6A ). At high resolutions, the mean breadth converges to the true value. The mean breadth becomes positive for resolutions better than 60 µm (R = 0.23). At lower resolutions, the solid phase is a complex network with inclusions of pore voxels. For high resolutions, the solid phase consists of separated convex spheres. For the latter case, positive breadth values are characteristic, and negative numbers result in the first case. The Euler characteristic (EPC) divided by the number of particles was computed for the solid phase (Fig. 6B). At poor resolutions, the spheres are connected, as indicated by the negative EPC. With increasing resolution, the spheres are separated and the EPC becomes positive for voxel sizes <20 µm (R = 0.08). For the best resolutions, the EPC converges to the value 1.0, with all spheres separated.


Figure 5
View larger version (15K):
[in this window]
[in a new window]
 
Fig. 5. Scaled surface area of the particles as a function of voxel size analyzed with simulated media. (A) The calculated surface area, divided by the analytically determined surface of the spheres, converges to the true value (1.0) with decreasing voxel size. The large symbols indicate the voxel size of the measurements with X-rays from tubes (black) or thermal neutrons (gray), respectively. (B) The absolute value of the change of area with decreasing voxel size is constant or decreases for voxel sizes smaller than 10 to 15% of the mean particle radius (arrows). The resolutions of the tomography did not fulfill this condition.

 

Figure 6
View larger version (16K):
[in this window]
[in a new window]
 
Fig. 6. Effect of image resolution on image properties analyzed with simulated media. (A) The mean breadth of the mapped particles, divided by the analytically determined breadth of the spheres, becomes positive for resolutions with voxel sizes smaller than 23 to 25% of the mean particle radius (arrows). The large symbols indicate the voxel size of the measurements with X-rays from tubes (black) or thermal neutrons (gray), respectively. For the tomography of the glass beads using neutron transmission technique, the voxel size was 16.7% of the mean particle radius and this requirement was fulfilled. (B) The Euler-Poincaré characteristic (EPC) of the solid phase divided by the number of particles. For high resolutions, the spheres in the image become separated and a positive EPC results. The resolutions of the tomography did not fulfill this condition.

 
Finally, the pore-size distribution was determined as a function of voxel size (Fig. 7A ). While the mean radius of the particles is 263 µm, the mean pore radius for the medium mapped with 6.25-µm voxel size is 44 µm (R = 0.17). With this high resolution, 40% of all pores are smaller than the smallest particles with a radius of 40 µm. The pore-size distribution calculated with the resolution of the X-ray tomography (60 µm) overestimates the fraction of the pores of radius 30 µm. With this resolution, 75% of the pore voxels have the radius 30 µm. In case of a voxel size of 6.25 µm, only 20% is allocated to pores ≤30 µm. Hence, the fraction of large pores is underestimated with the voxel size of the measurement, and the mapped pore-size distribution deviates from the true function. The true function (i.e., the pore-size distribution mapped with the maximum resolution) is characterized by a small volume fraction of small and large pores and large volume fraction of intermediate pore sizes. Therefore, the second derivative of the cumulative pore-size distribution must change from positive to negative values with increasing pore size. This criterion is fulfilled for resolutions ≤30 µm.


Figure 7
View larger version (19K):
[in this window]
[in a new window]
 
Fig. 7. Effect of resolution on pore-size distribution for simulated media. The cumulative pore-size distributions for (A) the sand model and (B) the generated glass bead packings were determined. The numbers in the legend are the resolution of the images in micrometers. The bold black lines with the large symbols show the size distribution computed with the resolution of the measurements with X-rays from tubes or thermal neutrons, respectively. The black curves without symbols were calculated with the highest resolution. The gray lines without symbols correspond to the largest voxel size fulfilling the criterion of a second derivative changing from positive to negative values with increasing pore size.

 
To summarize, based on the Minkowski functionals, the geometric properties turn toward the true values for voxel sizes smaller than 60 µm (R = 0.23) for the mean breadth, 30 µm (R = 0.11) for the area, and 20 µm (R = 0.08) for the EPC. To describe the shape of the pore-size distribution, a resolution of at least 30 µm is needed. Based on this analysis, we postulate that a voxel size of 10% of the mean particle radius is required to map the properties of the pore space. In addition, for a voxel size of 30 µm (R = 0.11), the majority of the voxels denoted as pores were entirely within the pore space. Therefore, compared with the mean particle radius (263 µm), the resolution R must be in the range of 0.1 (Table 1). The resolution obtained with micro-X-ray tomography (60 µm, R = 0.23) was not sufficient to reveal the details of the pore structure of this sand.


View this table:
[in this window]
[in a new window]
 
Table 1. Voxel sizes that are required to characterize the geometrical properties of the simulated sand and glass bead packings. These voxel sizes are compared with the voxel size of the tomography and with the mean particle radius determined by sieving. To compare the different media, all sizes were scaled with the respective mean particle radius.{dagger}

 
In addition, we used the generated sand medium to analyze the density distribution in cuboids of size 70 by 70 by 210 µm3. As explained above, the mass within such a volume contributes to the gray value of a voxel in the case of tomography using industrial X-rays. The calculated gray value for a position x, y as a function of height z is denoted as g(z). For the same generated medium, we calculated the density distribution in cubic voxels of size 70 by 70 by 70 µm3. These values are denoted as v(z). Then, we transformed g(z) into the Fourier domain. The inverse Fourier transformation of G({lambda})/H*({lambda}), denoted as v*(z), can be compared with v(z). We fitted the function H*({lambda}) to v(z) and applied the same procedure with the same function H*({lambda}) for the measured gray values to estimate the density distribution in cubic voxels. To optimize the function H*({lambda}), we applied thresholds to v(z) and v*(z) to obtain segmented images of the pore space with the same porosity as the measured medium. Then we counted the number of voxels with different phases in the two computed media. The smallest error was found for the values 0.70, 0.67, and 0.60 for h*(–1), h*(0), and h*(1). The error in the segmented image does not depend on the absolute values but only on the ratios between h*(–1), h*(0), and h*(1). With higher absolute values for h*, the range of the gray values becomes larger. The results for the generated and the scanned media are shown in Fig. 2.

Glass Bead Medium with Single Particle Size
To investigate numerically the effect of resolution on geometric properties of a sample with particles of equal size, 9168 spheres with a radius of 1 mm were distributed in a 40-mm cube. The fraction of glass in each voxel was calculated and expressed in terms of a gray-scale value between 0 and 255 for voxel sizes ranging from 40 to 500 µm, corresponding to images of the size 10003 and 803 voxels. Cross sections through the images are shown in Fig. 4.

The fraction of voxels determined as pores with non-zero gray values decreased with increasing resolution from 92% (voxel size 500 µm) to 11% (voxel size 40 µm). For the voxel size corresponding to the resolution obtained with thermal neutrons (167 µm, R = 0.17), 56% of the voxels determined as pores had a gray value of 0 (Fig. 3). For voxel sizes <200 µm (R = 0.20) and 70 µm (R = 0.07), the fraction of pore voxels with non-zero density values was <50% and 20%, respectively.

The calculated area divided by the analytically determined value converges to 1 with decreasing voxel size. For large voxel sizes, neighboring spheres overlapped in the image, and the surface of the intersection will be neglected. Therefore, the surface decreases with increasing voxel size (Fig. 5A). The change of surface with decreasing voxel size became constant for voxel sizes <150 µm (R = 0.15). This change of surface with voxel size differs from the function calculated for the sand. This is caused by the slightly different numerical packing procedure. In the case of the glass beads, almost all spheres are in contact with other spheres. Even for highly resolved images, the spheres are not completely separated, and the area of the intersection is neglected. In the case of the spheres of the sand model, fewer spheres are in contact with other particles. The most spheres become isolated for small voxel sizes, and the surface area is not affected by a further increase of resolution. The mean breadth of the solid phase (Fig. 6A) is positive for voxel size <250 µm (R = 0.25). Compared with the sand model, the values for high resolutions are smaller because the spheres touch each other, and these connections decrease the curvature. Then the Euler characteristic of the solid phase was determined (Fig. 6B). With resolutions better than 90 µm (R = 0.09), the spheres become separated and the EPC was positive. For high resolutions, the EPC divided by the total number of particles was higher than the expected value of 1. This deviation was caused by the inclusion of pore voxels in the vicinity of touching spheres. As mentioned above, pores enclosed in the solid phase increase the EPC. In the case of the packing algorithm used to model the glass beads, the most spheres touch other spheres because they were initially overlapping. In the discretized image, the surface of the spheres is rough, and it is possible that a pore voxel can be embedded between the voxels of two spheres that are very close. In the case of the packing algorithm used to simulate the sand media, there is only a small chance that a sphere touches another one and fewer voxels of the pore space are enclosed.

The pore-size distribution as a function of voxel size is shown in Fig. 7B. With the resolution of the tomography with thermal neutrons (167 µm), the shape of the curve corresponds to that mapped with the highest resolution. The second derivative of the cumulative pore-size distribution measured with a voxel size <200 µm is positive for small pores. The mean pore radius computed with the highest resolution is 319 µm (R = 0.32). Thus, the voxel size of the tomographic measurement was about 50% of the mean pore radius. In the case of the sand, the voxel size was 150% of the mean pore radius.

To summarize, the properties of the mapped pore space turn toward the true values for voxel sizes between 10 and 20% of the particle radius. The resolution required to map the pore space of the glass bead mixture must be smaller than 100 to 200 µm, or less than 10 or 20% of the mean particle radius (Table 1). With the resolution of the tomography with thermal neutrons (167 µm) this condition was partially fulfilled.


    APPLICATION: FLUID PHASE DISTRIBUTION
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MEASUREMENTS WITH COMPUTER...
 DESCRIPTION OF SIMULATED MEDIA
 GEOMETRICAL PROPERTIES AS A...
 APPLICATION: FLUID PHASE...
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
The purpose of scanning porous media is to analyze the effect of the pore geometry on the distribution of water and air. Instead of a direct measurement of the phase distribution, the arrangement of the water and air phase can be computed. Measurements and numerical experiments are discussed in the following subsections.

Measured Fluid Phase Distribution in Glass Bead Column Scanned with Thermal Neutrons
For thermal neutrons, the attenuation is dominated by the presence of water, so this method can be applied to measure the water distribution. Comparing the images of a dry and a wet sample, the differences due to the presence of water can be quantified. Often, heavy water (D2O) is used instead of normal water, because the attenuation of the neutrons is smaller for heavy water. The intensity of the beam attenuated by a system filled with normal water could be too small to be used for the reconstruction of the structure. We scanned a glass bead column (4-cm diam.) containing 3-mm glass beads. The bottom of the column was filled with heavy water. The reconstruction of the system and the boundary between the water saturated and dry zone is shown in Fig. 1D. For X-rays, the attenuation is density dependent and is smaller for water than for quartz sand. For this reason, the differences between the wet and a dry sample images are small, and it is difficult to distinguish between the two phases. To enhance the contrast, a potassium solution could be used instead of normal water (Hirono et al., 2003).

Fluid Phase Distribution Calculated with Lattice-Boltzmann Approach
A sample of diameter 25 mm with sand grains ranging from 80 to 1250 µm was scanned with X-rays from tubes with a resolution of 60 µm. For a reconstructed sand cube of length 3 mm (Fig. 8A ), the water and air distribution was calculated using a Lattice–Boltzmann approach. The method describes the dynamics of fluid particle distributions on a regular grid. Within each time step, particle distributions are modified by a simplified collision process and propagate to a finite set of neighboring nodes. The definition of the collision rules ensures that the dynamics of the moments of the particle distributions are equivalent to the dynamics of the variables of Navier-Stokes and the continuity equation up to second order. Thus, the Lattice-Boltzmann approach can be utilized to obtain approximate solutions of the Navier-Stokes equations. Details of this method are given, for example, in Wolf-Gladrow (2000) and Krafczyk (2001). In addition, the interface between the two-fluid phases can be computed with an extended approach (Tölke et al., 2002; Sukop and Or, 2004). For that purpose, the dynamics of two probability distributions representing the air and the water phase are computed.


Figure 8
View larger version (106K):
[in this window]
[in a new window]
 
Fig. 8. A sample of dry sand was scanned using micro-X-ray tomography with a resolution of 60 µm. (A) The particles in a cubic section of size 3 mm are shown. The water and air distributions in this cube for different pressure conditions at the lower boundary were calculated using the Lattice-Boltzmann approach. With this method, the Navier-Stokes equation for the liquid and the gas phase are solved numerically. For a drainage process, the resulting (B) air and (C) water distributions are shown. (D) For the same applied suction during a wetting process, the water distribution is different.

 
In numerical experiments, different water pressures were applied at the lower boundary of the sand cube, and the resulting water and air distributions were computed. First, the dry cube was brought in contact with water at the bottom of the reconstructed sample. Water rose due to capillary forces reaching an equilibrium distribution of the two fluid phases. Then, suction with a certain value, denoted as initial value, was applied at the bottom, which redistributed the air and water phases (Fig. 8B and 8C) and displaced some of the water out of the sample. Then, after an equilibration period with a higher suction, the suction was released to the initial value and water entered from the bottom. The resulting water (Fig. 8D) and air distribution was changed compared with the initial phase distribution because other pore structures are dominant for drainage and wetting processes. So, with this numerical experiment, the hysteresis phenomenon could be visualized.

The simulation of the full two-phase dynamics in a porous medium is a demanding task in terms of computational resources. To calculate the fluid phase distribution in larger systems, we used another approach that we denote as a "morphological pore network model." We introduce this approach in the following subsection.

Fluid Phase Distribution Calculated with Morphological Pore Network Model
In network models, the pore space consists of elements with well-defined size. Most often, the pore space is described as a composite of pore bodies connected with pore throats. Size, shape, and arrangement of these elements can be adapted to mimic the fluid distribution within a porous media. In contrast to such a conceptual representation of the porous medium, the measured image of the pore structure is used in case of morphological pore networks. The "true" pore space is subdivided into elements of well-defined size. In the approach used in this study, each pore voxel is an element of the pore space and its size is the radius of a sphere as explained. The size of each pore voxel will be needed to calculate the water and air distribution. Because the whole geometrical information is included in this network approach, we denote this type as morphological network. A similar approach to map the size of the pores and to calculate the fluid phase distribution in numerically generated sphere packings was applied by Hilpert and Miller (2001). Here, we use the morphological network to calculate the water and air distribution in sand media measured with X-rays from synchrotrons. We calculated the air and water distribution for a coarse sand with particle sizes ranging from 300 to 900 µm scanned at HASYLAB with voxel size 11.5 µm and for a fine sand with particles ranging from 100 to 200 µm scanned with a voxel size of 3.5 µm at SLS. The voxel sizes correspond to 4% of the mean particle radius (75 and 304 µm). So, the requirement of a voxel size smaller than 10% of the mean particle radius was fulfilled. In a numerical drainage experiment, we calculated the water and air distribution for a suction applied at the lower boundary. We used the three-dimensional segmented image and computed the size of each pore voxel as explained above. In Fig. 9A and 9B, the calculated pore sizes for the voxels of the pore space are shown for two scanned sand samples. For the sand shown in Fig. 9A, a small section is enlarged to explain the determination of the pore size in more detail. First, the distance from the next solid was calculated for each voxel of the pore space (Fig. 9A1). This distance corresponds to the radius of a sphere within the pore space. A voxel can be an element of different spheres and the radius of the largest sphere determined the size of the voxel. This is shown in Fig. 9A2.


Figure 9
View larger version (110K):
[in this window]
[in a new window]
 
Fig. 9. Dry sand samples were mapped with synchrotron radiation tomography (left: HASYLAB, right: SLS). Particles are shown in black. The sizes of the particles range (A and C) from 300 to 900 µm and (B and D) from 100 to 200 µm. (A and B) First, the pore sizes for all voxels of the pore space were determined. Bright values indicate large pore sizes. A small section indicated by the white border is enlarged to explain the determination of the pore size. First, the distance from the solid phase is calculated (A1). This distance corresponds to the radius of a sphere. A voxel of the pore space can be an element of different spheres and the radius of the largest sphere determines the pore size of the voxel (A2). The grid shows the size of the voxels. In the bottom row, the water (gray) and air (white) distributions computed with a morphological pore network model are shown. A head of (C) –15 cm and (D) –50 cm was applied at the lower boundary in the numerical experiment. The computed water saturation of the pore space was 56% for the coarse sand and 49% for the fine sand.

 
At the beginning of the experiment, the pore space was assumed to be completely water-filled. The lower boundary was connected to a hanging water column and air was allowed to enter at the upper boundary. A water filled pore voxel drained if the following three requirements were fulfilled: (i) a neighboring pore must be air-filled, (ii) the water-filled voxel must be hydraulically connected to water at the lower boundary, and (iii) the applied suction must be larger than the pore size dependent capillary forces. The calculated water and air distribution in a cross section for an applied hydraulic head of –50 cm for the fine and –15 cm for the coarse sand is shown in Fig. 9C and 9D. After each drainage step, the suction applied at the lower boundary was increased to withdraw more water from the sample. The volume fraction of retained water as a function of applied pressure is shown in Fig. 10 . Neglecting conditions (i) and (ii), the fluid phase distribution would be dominated by the pore size. All pores in which the size-dependent capillary forces are low enough would drain out, and the water saturations would correspond to the cumulative pore-size distribution, which are also shown for comparison in Fig. 10. Obviously, the cumulative pore-size distribution underestimates the water content. In the case of the morphological pore network model, the large pores do not drain before the surrounding small pores are air filled. In addition, with the pore network approach, the existence of residual water can be modeled.


Figure 10
View larger version (15K):
[in this window]
[in a new window]
 
Fig. 10. The effect of pore size and pore connectivity on the water retention curve. (A) The drainage was calculated with the morphological pore network model for the sand material with particles ranging from 300 to 900 µm ("network"). The dry sand sample was analyzed with X-rays from synchrotrons with a resolution of 11.5 µm. The size of the image was 900 by 900 by 300 voxels. The modeled retention curve corresponds well to the water retention curve measured in the laboratory by Ursino and Gimmi (2004). The computed and measured drainage curves are compared with the cumulated pore-size distribution ("size"). In the latter case, connectivity is neglected and all pores with low capillary forces drain out. (B) The cumulative pore-size distribution ("size") was compared with the computed drainage curve ("network") for the sand sample with fine particles of size 100 to 200 µm. The pore structure of the dry sand sample was measured with a resolution of 3.5 µm using X-rays from synchrotrons. The image size was 1000 by 1000 by 100 voxels. By neglecting the connectivity of the pore space, the water saturation was underestimated.

 
In the case of the coarse sand material, the water retention curve was determined in the laboratory (Ursino and Gimmi, 2004). The soil water characteristic computed with the morphological pore network corresponds well to the measured data. The slight underestimation of the air entry value is caused by the small height of the sample used in the numerical experiment. At the upper boundary large pores exist. The drainage of this thin layer with large pores at the top cannot be neglected for such a thin sample. For the laboratory column, the drainage of the large pores at the upper boundary is of minor importance. We conclude that the effect of pore size is not sufficient to explain the water and air distribution in porous media. However, taking into account the connectivity of the pore space, we can calculate the phase distribution. A prerequisite for a successful prediction is a structure map with high resolution on the order of 10% of the mean particle radius. By means of synchrotron radiation, such accurate structure detection was possible.


    SUMMARY AND CONCLUSIONS
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MEASUREMENTS WITH COMPUTER...
 DESCRIPTION OF SIMULATED MEDIA
 GEOMETRICAL PROPERTIES AS A...
 APPLICATION: FLUID PHASE...
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Porous media were analyzed using tomography with X-rays from tubes, thermal neutrons, and X-rays from synchrotrons. Comparing the three methods, the advantage of the tomography with X-rays from tubes is the potential to scan large samples (large with respect to the pore sizes discussed in this study). Also the existence of water within the sample does not limit the use of X-rays. However, the resolution that can be obtained is limited compared with the use of X-rays from the synchrotron. Three-dimensional structures with sizes of a few micrometers can be mapped with the synchrotron technique only. The advantage of tomography with thermal neutrons is the sensitivity with respect to the presence of water; spatial variations in the water content can be mapped.

Based on the analysis of numerically generated porous media, the voxel size must be smaller than 10 to 20% of the mean particle radius to reveal the geometrical properties of the pore space. This criterion is based on the investigation of the effect of resolution on the pore-size distribution, the surface of the particles, the breadth and the Euler characteristic of the solid phase. For a voxel size of 10 to 20% of the mean particle radius, these geometric properties turn toward true values. In addition, more than 50% of the voxels denoted as pores are entirely within the pore space for this resolution. The requirement of a voxel size smaller than 10 to 20% of the mean particle radius was not fulfilled for the tomography of a sand sample with particle diameters ranging from 80 to 1250 µm using X-rays from tubes. The resolution of the tomography divided by the mean particle size (263 µm) was 0.23. For the tomography with thermal neutrons, the resolution divided by the particle radius of 1000 µm was 0.17, and the requirements were partially fulfilled. For the tomography of two sand samples with particles of diameter 100 to 200 and 300 to 900 µm using synchrotron radiation, the resolution was about 4% of the mean particle radius (75 and 304 µm), and all requirements were fulfilled—note that this analysis is valid for the quantification of pores between particles. In the case of worm burrows, the pore structure is independent of the size of the particles, and other models of structure generation must be applied to analyze such structures and to define similar requirements.

We used the images of dry sand samples to calculate the water and air distribution for different boundary conditions. With the images obtained by X-rays from synchrotrons, we could predict the drainage of the sand sample in absolute terms with a morphological pore network model. In the case of the sand scanned with X-rays from tubes, the hysteresis phenomenon could be simulated with a Lattice Boltzmann approach.

By combining the appropriate technique of tomography with the numerical tools to model multiphase flow, new insights in the processes at the pore scale and more accurate predictions of the hydraulic properties of soil materials can be expected.


    REFERENCES
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MEASUREMENTS WITH COMPUTER...
 DESCRIPTION OF SIMULATED MEDIA
 GEOMETRICAL PROPERTIES AS A...
 APPLICATION: FLUID PHASE...
 SUMMARY AND CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
Vadose Zone JHome page
A. Carminati, D. Vetterlein, U. Weller, H.-J. Vogel, and S. E. Oswald
When Roots Lose Contact
Vadose Zone J., August 11, 2009; 8(3): 805 - 809.
[Abstract] [Full Text] [PDF]


Home page
Vadose Zone JHome page
S. E. Oswald, M. Menon, A. Carminati, P. Vontobel, E. Lehmann, and R. Schulin
Quantitative Imaging of Infiltration, Root Growth, and Root Water Uptake via Neutron Radiography
Vadose Zone J., August 13, 2008; 7(3): 1035 - 1047.
[Abstract] [Full Text] [PDF]


Home page
Soil Sci.Home page
T. Sander and H. H. Gerke
Noncontact Shrinkage Curve Determination for Soil Clods and Aggregates by Three-Dimensional Optical Scanning
Soil Sci. Soc. Am. J., August 9, 2007; 71(5): 1448 - 1454.
[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 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 Web of Science (4)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Lehmann, P.
Right arrow Articles by Flühler, H.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Lehmann, P.
Right arrow Articles by Flühler, H.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Lehmann, P.
Right arrow Articles by Flühler, H.
Related Collections
Right arrow Tomography
Right arrow Pore-Scale Modeling
Right arrow Multiphase Fluid Flow


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