Published online 29 May 2008
Published in Soil Sci Soc Am J 72:897-907 (2008)
DOI: 10.2136/sssaj2007.0130
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SOIL PHYSICS
Three-Dimensional Quantification of Intra-Aggregate Pore-Space Features using Synchrotron-Radiation-Based Microtomography
S. Peth*,
R. Horn,
F. Beckmann,
T. Donath,
J. Fischer and
A. J. M. Smucker
Dep. of Plant Nutrition and Soil Science Christian-Albrechts-Univ. Kiel Germany
* Corresponding author (s.peth{at}soils.uni-kiel.de).
 |
ABSTRACT
|
|---|
Pore network geometries of intra-aggregate pore spaces are of great importance for water and ion flux rates controlling C sequestration and bioremediation. Advances in non-invasive three-dimensional imaging techniques such as synchrotron-radiation-based x-ray microtomography (SR-µCT), offer excellent opportunities to study the interrelationships between pore network geometry and physical processes at spatial resolutions of a few micrometers. In this paper we present quantitative three-dimensional pore-space geometry analyses of small scale (
5 mm across) soil aggregates from different soil management systems (conventionally tilled vs. grassland). Reconstructed three-dimensional microtomography images at approximate isotropic voxel resolutions between 3.2 and 5.4 µm were analyzed for pore-space morphologies using a suite of image processing algorithms associated with the software published by Lindquist et al. Among the features quantified were pore-size distributions (PSDs), throat-area distributions, effective throat/pore radii ratios as well as frequency distributions of pore channel lengths, widths, and flow path tortuosities. We observed differences in storage and transport relevant pore-space morphological features between the two aggregates. Nodal pore volumes and throat surface areas were significantly smaller for the conventionally tilled (Conv.T.) aggregate (mode
7.9 x 10–7 mm3/
63 µm2) than for the grassland aggregate (mode
5.0 x 10–6 mm3/
400 µm2), respectively. Path lengths were shorter for the Conv.T. aggregate (maximum lengths < 200 µm) compared with the grassland aggregate (maximum lengths > 600 µm). In summary, the soil aggregate from the Conv.T site showed more gas and water transport limiting micromorphological features compared with the aggregate from the grassland management system.
Abbreviations: Conv.T., conventional tillage DU, down to upward flow IK, indicator kriging LR, left to right flow PSD, pore-size distribution SOC, soil organic carbon SR-µCT, synchrotron radiation microtomography TSD, throat-size distribution
 |
INTRODUCTION
|
|---|
Assemblages of solids and associated voids within soil aggregates offer numerous microsites that become interconnected to form a plethora of storage and conductive passages important for maintaining essential soil biogeochemical and biophysical processes. To date, discrete functions of these interconnected pathways have been limited primarily to mathematical expressions of externally applied solutions, ions, gases and thermal gradients (Hillel, 1998). Statistical variations of relative fractions of these multiple phases in soils, ignoring their spatial arrangement is not sufficient to understand the mechanisms controlling soil functions. Past research has shown that soil structure is a key contributor to the intensity and rates of reaction and transport processes within soils (Boone, 1988; De Wever et al., 2004; Horn and Smucker, 2005). Soil structure effects are evident at pedon and interaggregate scales, for example, by preferential flow phenomena (Gerke, 2006). Recently, since aggregate interiors have been shown to be the repositories of soil C (Six et al., 2004) discussions on climate change and the potential of soils as a C sink have drawn the attention of soil scientists to the smaller intra-aggregate scale. Variations in C storage between soils are attributable to differences in aggregate structure formation and their stability on one hand and to different microbial activities within aggregates on the other hand (Horn et al., 1994; Chenu et al., 2001; Smucker et al., 2007). Other studies have reported a spatial heterogeneity of C storage sites within a range of aggregates that controls the distribution of microbes (Van Gestel et al., 1996).
In addition to the direct relationships between soil C sequestration by aggregates and climate change, soil aggregates and associated microbes turned out to have a great potential in mitigating soil contaminants such as organic compounds (Bouwer and Zehnder, 1993) and heavy metals (Tokunaga et al., 2003). Both bioremediation and carbon sequestration involve processes that to a great extent are controlled by the complex nature of the intra-aggregate pore network (Smucker et al., 2007). Pore-space geometry determines the transport of water, oxygen, and nutrients to the microsites populated by soil microbes and hence the efficiency of biogeochemical reactions. Despite the importance of the microhabitat structure for many microbial processes few details are known about the topology and spatial geometries of soil pore networks at the aggregate scale (Nunan et al., 2006). Six et al. (2004) pointed out that research linking physical processes in (micro)aggregates with soil biota and soil organic matter dynamics could greatly benefit from three-dimensional quantitative spatial information on the aggregate level. Kutilek and Nielsen (2007) argued that there is a need to quantify micromorphological soil characteristics in future soil science research as the configuration and arrangement of pores is a key to understand soil physical processes. First steps in relating physical and biological process in the microenvironment have been reported by Young and Crawford (2004) and Tokunaga et al. (2003).
In this paper we present an investigation of the complex nature of intra-aggregate pore spaces of two contrasting soil aggregates down to a resolution of a few microns using the SR-µCT at the Deutsches Elektronen Synchrotron (DESY) in Hamburg/Germany. Our objective was to visualize and quantify microscale morphological features of the pore-space network by three-dimensiaonal image analysis to quantitatively compare the internal pore-space architecture of soil aggregates that developed under different management systems.
 |
MATERIALS AND METHODS
|
|---|
Soil Sampling and Basic Soil Characteristics
Due to limited availability of synchrotron beam time we focused our study on two selected soil aggregates. The aggregates have been sampled at the experimental research farm Rotthalmünster (near Passau)/southern Germany from the plowed horizon of a Conv.T site and from an unplowed horizon which has been under grassland use since 1961 (hereafter referred to as grassland). Soils from both sites are classified as Alfisol (Soil Survey Staff, 2006) derived from loess with almost identical texture (10% sand, 70% silt, 20% clay). The soil aggregates had a diameter of approximately 3.3 (Conv.T.) and 4.6 mm (grassland) and were air-dried prior scanning. Additional details of soil aggregate properties are reported by Jasinska (2006).
SR-µCT Data Acquisition and Image Reconstruction
Aggregates were scanned at the SR-µCT facility operated by the GKSS Research Center at HASYLAB (Hamburger Synchrotron Strahlungslabor) of DESY (Deutsches Elektronen Synchrotron) in Hamburg/Germany. The microtomography was completed at beamline Sector BW2. Samples were mounted on a rotary stage and scanned at a photon energy of 21.5 (Conv.T.) and 24.0 keV (grassland), respectively. Photon energies and magnifications were optimized to the sample diameter to obtain maximum resolution. Given the different diameter of the soil aggregates, the achieved voxel length was slightly different with
3.2 µm for the Conv.T. aggregate and
5.4 µm for the grassland aggregate. During the microtomographical scan the sample was rotated at different projection angles between 0 and 180°(
), at 0.5° degree intervals. The absorption radiographs for the different viewing angles were projected onto a CdWO4 fluorescent screen and recorded by a 1536 by 1024 pixel CCD camera.
One advantage of using synchrotron radiation for x-ray microtomography is that due to the highly intense radiation with low divergence a monochromatic x-ray beam is obtained with the aid of a double crystal monochromator (Beckmann et al., 2005). Due to the low divergence and the large distance between the source and the experiment the sample is projected by a nearly parallel beam. The parallel beam in conjunction with the CCD camera allows for simultaneous attenuation measurements of multiple stacked planes to give three-dimensional maps of attenuation values with isotropic spatial resolutions. Image reconstruction was accomplished using a filtered back-projection algorithm in IDL (Research Systems Inc., Boulder, CO). Attenuation coefficients are expressed as greyscale values ranging from [255] for the highest attenuation coefficient to [0] for the lowest attenuation coefficient reached.
Further information on the microtomography system at the HASYLAB is provided in Beckmann et al. (1999). Details on different tomography techniques and error sources are discussed in Wildenschild et al. (2002), Stock (1999) and Ketcham and Carlson (2001). An overview on the principles of transmission tomography, three dimensional x-ray microtomography and image reconstructions is provided by Flannery et al. (1987).
Image Transformation Algorithms
The process of image transformation and the geometrical analysis of reconstructed three-dimensional voxel images are facilitated by algorithms that are based on the theoretical framework of mathematical morphology outlined by Serra (1982) and Soille (2003). A summary of the fundamental ideas and basic concepts underlying mathematical morphology and their applications to soil science was presented by Horgan (1998).
We have processed the image data sets with a suite of algorithms assembled into a software package called 3dma (Lindquist et al., 2000). It covers a comprehensive set of successive image transformations and measurements developed for the analysis of three-dimensional mineral porous media. The performance and accuracy of the algorithms has been tested on both images from natural porous media and simulated images of hexagonal closed packed spheres with relative errors <5% (Lindquist et al., 2000). A complete description of the used algorithms is beyond the scope of this paper. Therefore only principles of the involved image transformations and pore-space quantifications are outlined in this report. A schematic representation of the workflow is shown in Fig. 1
(Please note color figures are available in the online version of the article).
Extraction of Subvolume
Image analysis in 3dma is limited to regularly shaped volume objects (cylinders or cubes). Therefore a 400 x 400 x 400 voxel sized cube was extracted from the original microtomograms from near the center of the aggregate (Fig. 2
). With voxel resolutions of 3.214 x 10–6 m for the Conv.T. aggregate and 5.403 x 10–6 m for the grassland aggregate the analyzed cubes have an edge length/volume of 1.286 x 10–3 m/2.13 x 10–9 m3 for the Conv.T. aggregate and 2.161 x 10–3 m/10.1 x 10–9 m3 for the grassland aggregate, respectively. The subvolumes have been placed such that in three dimensions a representative domain of the aggregate was covered. Although it seems that the extracted cube of the Conv. T. aggregate could be larger it was restricted in size by the aggregate boundary normal to the image slice shown in Fig. 2.

View larger version (46K):
[in this window]
[in a new window]
|
Fig. 2. Tomograms and segmented images of analyzed soil aggregates. A1: Cross-section of aggregate from conventional tillage site (Conv.T.); B1: Cross-section of aggregate from conservation tillage site (grassland). Dashed white line indicates the two-dimensional representation of the extracted and analyzed cubical subvolume. A2/B2: segmented image slice of the subvolume using a simple (global) thresholding technique; A3/B3: segmented image slice of the subvolume using thresholding by indicator kriging.
|
|
Segmentation
The first step in image processing involves the conversion of each voxel gray scale value (proportionally expressing the local attenuation coefficient for x-rays) into a bimodal image that distinguishes between void and solid phase. This transformation is referred to as image segmentation and is a key step toward the quantitative analysis of image data as all following operations use these transformed datasets. The quality of the segmentation therefore determines the quality of the final results. The most commonly employed method is to choose a global (single) threshold value to segment an image (Pal and Pal, 1993). The value selected usually is based on a priori information of the sample such as bulk porosity of the imaged subvolume. When irregularly shaped objects are investigated, for example, soil aggregates, where the analyzed subvolume does not match the total volume of the object for which bulk porosity data may be available, some bias is introduced in selecting the correct threshold value. Capowiez et al. (1998) have used a simple rule to determine the threshold value based on the gray level histogram. By adding 1/3 of the distance between the pore peak and the matrix peak to the pore peak they identified the approximate minimum of the distribution function between the two peaks. This seems intuitively right when well-separated peaks with a low frequency at the minimum indicating a clear distinction between void and grain phase attenuation coefficients are encountered. However, for cases with a strong peak overlap in the intensity histogram, additional approaches are required. Sometimes there is higher variability of the attenuation coefficients that result from the greater variation of material densities associated with organic matter and different minerals. Additional reasons for misclassifications of phases include partial volume effects (voxels sharing void and solid at varying fractions due to the finite volume resolution) and noise in the tomographic reconstruction.
Alternatively, local thresholding schemes have been suggested where the voxel classification depends on the gray scale values of its surrounding voxels instead of taking a global gray level value as threshold (Pal and Pal, 1993). Oh and Lindquist (1999) developed a local thresholding method based on the Mardia-Hainsworth spatial thresholding algorithm (Mardia and Hainsworth, 1988) where the spatial covariance of a random variable (e.g., the attenuation coefficient µ) is used in conjunction with indicator kriging (IK) to determine object edges. Indicator kriging is a geostatistical interpolation procedure usually applied when the distribution function of a random variable is highly skewed or has more than one peak. This frequently is the case for tomography images due to variations in solid densities, partial volume effects and image noise. In image segmentation IK provides a means for estimating the probability that an unclassified voxel at location x0 belongs to the population void (
0) or grain (
1). The probability estimate depends on the proportion of data z(x
) valued above a threshold zc in the local neighborhood at locations x
,
= 1, ..., n, expressed by a set of indicator values (Oh and Lindquist, 1999):
 | [1] |
A linear estimate of the conditional probability that the unknown z(x0) is lower than or equal to a threshold Ti, i.e., P(Ti; x0|n)
Prob{z(x0)
Ti}, is then given by:
 | [2] |
The n weighting coefficients 
in (2) are determined by an ordinary kriging system (Journel, 1989) accounting for the proximity of x
to x0. In the kriging step the population is then assigned for all unclassified voxels x0 according to:
 | [3] |
In practice the local thresholding algorithm of Oh and Lindquist (1999) assigns voxel population in two consecutive passes. Two cut-off values, one for voids "T0" and one for solids "T1" were chosen from the gray scale histogram based on an initial estimate of the threshold value by visually comparing the segmented image with the original tomogram. In a first pass of the algorithm voxels with gray level values below T0 are assigned void, those above T1 are assigned solid grain. The remaining unclassified voxels (gray scale values between T0 and T1), representing the problematic segmentation range, were allocated in a second pass by IK. This procedure proves to significantly reduce misclassification errors (Oh and Lindquist, 1999) and resulted in our case in less noisy bimodal images. Image noise of segmented images is an important criterion for the performance of the pore network skeletonization algorithm.
As mineral composition is considered to be similar in both aggregates we have chosen a common upper gray level threshold T0 of 85 (above this value voxels are considered to belong to the grain phase) and a lower gray level threshold T1 of 75 (below this value voxels are considered to belong to the void phase). Voxels having gray levels between these two bounding values were classified by IK. For comparison simple thresholding with an intermediate global threshold value of 80 was used to segment both aggregates. However, all subsequent image transformations were done using the segmented images processed by the local thresholding scheme. This was found to significantly reduce the extent of post segmentation clean up and improved succeeding image analyses.
Post Segmentation Clean Up
Small isolated clusters of void or grain voxels may strongly distort the skeletonization process. While disconnected grain voxels are in most cases unphysical (except at image boundaries) disconnected void voxels may correspond to small isolated pores or to noise effects. Disconnected grain voxel clusters should be all removed from the image prior further analyses by converting them into void voxels. In the case of isolated void voxels, however, it is difficult to distinguish between physical and unphysical origin (resulting from segmentation artifacts). It is assumed that the probability of misidentifying voxel clusters as unphysical rapidly increases with increasing cluster size (Lindquist et al., 2000). Hence the frequency distribution of isolated void cluster sizes may be used as a decision criterion for removing clusters that are probably unphysical. Based on this distribution function it is up to the user to specify a filter size for removing disconnected isolated voxel clusters. We have followed the strategy to remove only very small isolated void clusters to not significantly change the porosity. Volume thresholds for void voxel cluster conversion have been set to 10 voxels in both cases. Although, respectively, 41 (Conv. T.) and 38% (grassland) of the isolated void voxel clusters have been removed corresponding total porosity changes were very small (Table 1
). In case of isolated grain voxel clusters, all clusters have been removed from the images.
View this table:
[in this window]
[in a new window]
|
Table 1. Volume thresholds for conversion of isolated void and grain clusters with corresponding changes in porosity.
|
|
Sealing Image Boundary
A lack of information beyond the imaged region may cause difficulties in the construction of pore throats since channels may be cut by the image boundary. This complication can be avoided by sealing the selected subvolume with a layer of grain voxels before skeletonization, as outlined by Prodanovi
(2005).
Medial Axis Construction (Skeletonization)
The geometry of the pore network in three-dimensional porous media is difficult to visually examine due to grain phases partially hiding the pore space. To circumvent this problem a lower dimensional representation of the object is desired to enhance the visualization of the pore-space geometry. Algorithms providing this are known as skeletonization or medial axis algorithms (Lee et al., 1994; Lindquist et al., 1995). The skeleton or the medial axis is a network of interconnected one-dimensional paths that has a strict geometrical relationship to the solid surface (Lindquist et al., 2000). The medial axis preserves important geometrical features of the solid void interface. The principle operation of the medial axis algorithm may be described by following process. Assume a fire is ignited simultaneously at the solid-void interface in the entire image domain and that it travels with a constant expansion rate into the pore space perpendicular to the void-grain boundary. The set of voxels where two approaching fire fronts meet (here the fire is considered to be extinguished) is defined as the medial axis. Each medial axis voxel is associated with a burn number which corresponds to the shortest distance in units of number of voxels from the medial axis toward the solid–void interface (Lindquist et al., 1995).
Medial axis construction is sensitive to surface noise which occurs (i) as irregularities of the digitized void-grain surface that otherwise would be smooth and (ii) as disconnected clusters of misclassified void or grain voxels (Lindquist and Venkatarangan, 1999). While disconnected void clusters will have its own isolated skeleton disconnected grain clusters prevent the skeleton to reduce to a union of one dimensional curves (Lindquist et al., 2000). Another problem in medial axis construction is related to the irregularity (roughness) of the grain-void interface which is particularly pronounced in soils. This generates numerous dead-end medial axis paths that may be "real" or the result of surface noise. Removal of dead-end paths is justified to some extent since they have little relevance for transport processes. However, medial axis modifications will affect throat constructions and hence the calculation of PSDs and should be done carefully. The employed image analysis software 3dma provides various options for the removal of dead-end paths and isolated medial axis pieces from the skeleton. The modification of the medial axis is controlled based on user specified trimming criteria. We have found the following criteria to be the most suitable in our case which have been applied to the medial axis for both soil aggregates: (i) Isolated voxels with a burn number less than 2 were removed; (ii) Branch-leaf paths (paths connected only at one end to a branch cluster) were removed when their length was less than the maximum cluster burn number; and (iii) Needle-eye paths (paths connected to the same branch cluster at both ends) were removed when their length was less than the maximum cluster burn number.
Using above specifications 3.3 (Conv.T.) and 5.2% (grassland) of the medial axis voxels were removed from the skeleton, respectively. Once the medial axis is constructed and pruned it is used as an embedded search path of the object to find specific sites for further analysis of geometrical features. Three-dimensional visualization of the medial axis was done with the free interactive viewing software Geomview 1.8 for Unix (http://www.geomview.org, verified 7 Apr. 2008)
Pore Throat Finding and Pore Partitioning
The network of porous media can be described by nodal pores connected by channels with each channel containing a local flow bottleneck. These bottlenecks or pore throats represent the minimum cross-sectional area of the corresponding channel. Pore throats are searched by uniformly dilating a topological cylinder from each medial axis path (single channel) in the radial direction perpendicular to its length (Lindquist et al., 2000). With continuing dilation the cylinder increases and eventually hits the grain surface of the pore pathway. The points of contact first identified are mapped and compared with progressively increasing cylinder radii which finally form into a closed loop. This loop corresponds to the perimeter of the minimum cross-sectional area and therefore represents the minimum throat perimeter. The throat surface is then calculated by constructing a triangulated surface from the medial axis path voxel at the throat center (defined by the smallest average distance to the perimeter voxels) to the respective centers of the loop of perimeter voxels (Prodanovi
, 2005). Once throats have been identified the pore space is partitioned into interconnected regions (nodal pores) separated by throat and grain voxels.
Pore Network Statistics
Pore partitioning by constructing the pore-throat-network as describe above is the basis for a quantitative analysis of the three-dimensional pore space. One of the most obvious properties of interest is the PSD. The pore volumes of all nodal pores are calculated by counting the number of corresponding voxels of the irregularly shaped pore and multiplying it with the resolution dependent voxel volume. To account for the volume occupied by throat voxels, the volume of each pore is corrected by adding to it half the volume of all of its delimiting throat barriers. Since each imaging technique has a finite resolution the derived PSD only accounts for pores larger than the voxel resolution of the image. The throat size (surface area) distribution (TSD) between connected nodal pores is derived by summing the triangulated surfaces from the throat surface construction described in the previous section. Both nodal pore volume and throat surface area are then converted into effective pore radii (radius of a sphere having the same volume) and effective throat radii (radius of a circle having the same area), respectively. An effective throat to pore radius ratio is calculated for each throat using the average of the effective pore radii of the two pores separated by the throat (Lindquist et al., 2000). This dimensionless number is an indicator for the fluid flow impedance between interconnected pore volumes through bottle-necks. The lower the ratio the more restricted is the permeability by bottle-necks.
Pore channel length is determined by summing the medial axis voxels between the centers of any two adjacent (i.e., connected) nodal pores and multiplying it with the unit length of a voxel. In this way the distance along the curved backbone of the pore channel rather than the straight line distance between nodal pore centers is measured. Note that employed clean up and pruning routines, for example, removal of isolated void/grain clusters and dead-end paths, influence the obtained statistical results. Channel width is indicated in terms of burn number, which is allocated to each medial axis voxel during medial axis construction. The higher the burn number the wider is the flow channel. Physically it is the radius in units of number of voxels of the largest sphere centered at the medial axis which just fits inside the void space (Lindquist et al., 1995). As a measure of the tortuosity of the pore networks a so-called tortuosity index (
1) is calculated for each connecting channel. It is defined as the ratio of the curvilinear channel length and the Euclidian distance between the two end voxels of the channel. A tortuosity index of one denotes a straight channel while higher numbers indicate higher channel tortuosities.
 |
RESULTS AND DISCUSSION
|
|---|
Image Segmentation and Total Porosity
Total porosities calculated from binary images segmented by global thresholding were higher for the Conv.T. aggregate (15.7%) compared with the grassland aggregate (11.1%). Note that the true threshold value is unknown, given that a priori bulk porosity information of the analyzed subvolume is not available. However, we may assume a common material cut-off value since both soil aggregates have a more or less identical solid phase composition (Jasinska, 2006). Porosities for indicator kriged images where in both cases lower with 8.4 (Conv.T.) and 9.4% (grassland), respectively. The strong decrease in porosity for Conv.T. is attributed to a higher fraction of voxels estimated by IK (Table 2
). This is related to the higher number of void-grain interface voxels for Conv.T. where population identification is impaired due to partial volume effects (Oh and Lindquist, 1999). Therefore, the number of false positives of pore identification may be higher for the Conv.T. aggregate than for the grassland aggregate.
Differences in image segmentation results using global and local thresholding procedures are shown in Fig. 2. Binary images, segmented by global thresholding, are more mottled compared with local thresholding using IK. Reduced mottling for the grassland aggregate particularly in the neighborhood of larger macro-/biopores suggests a higher soil matrix density. This is reasonable since the formation of biopores, for example, by penetrating root tips is associated with a compression of the surrounding soil material (Dexter, 1991; Dorioz et al., 1993). Consequently, given the increasing uncertainty for accurate image thresholding of pores close to the resolution limit and implying that larger pore channels are the most relevant conduits for fluid transport all further analysis were conducted using images segmented by IK. Statistical information on the pore geometry reported below refers to the "main" pore network (percolation pore network) excluding the unresolved pore space of the soil matrix. Comparison of bimodal and gray level images in Fig. 2 shows that larger connected pore volumes are reasonably identified accurately for both aggregates.
Quantification of the Three-Dimensional Pore-space Geometry
Relative frequency distributions of nodal pore volumes and throat surface areas and their corresponding effective pore and throat radii are distinctly different between the two aggregates (Fig. 3
). Note that cumulative frequencies are normalized to one in both cases. Higher relative fractions of smaller pores are observed for the Conv.T. aggregate compared with the grassland aggregate. The dominating nodal pore volume class is approximately one order of magnitude smaller for the Conv.T. aggregate than for the grassland aggregate. This result reflects well the higher number of macropores visible in Fig. 2.

View larger version (15K):
[in this window]
[in a new window]
|
Fig. 3. Relative frequency distribution of a) nodal pore volumes, b) effective pore radii, c) throat surface area and d) effective throat radii for the Conv.T. and grassland aggregate. Note that cumulative frequencies are one in both cases.
|
|
From the distribution of effective pore radii in Fig. 3b the water retention function was estimated. This was done by calculating the capillary suction
c for the mean of each effective pore radius class using:
 | [4] |
where
is the surface tension of water (72.7 x 10–3 N m–1),
is the wetting angle (assumed to be equal to 0°, i.e., perfect wetting) and rw is the radius of the capillary. Note that this procedure has limits since effective pore radii calculated from 3dma are not quite the same as effective pore radii underlying Eq. [4]. Applying Eq. [4] to soils it is assumed that the pore network consists of a bundle of interconnected capillary tubes while 3dma calculates effective pore radii based on spheres with the same volume as the irregularly shaped pore body. From a geometrical point of view the pore network may be considered as a superposition of capillaries and spheres. In other words neither capillaries nor spheres are adequate representations of the pore size but if at all a combination of both. However, from this standpoint it should at least be possible to approximate the water retention function from "spherical" effective radii with a similar quality as PSD are derived from Eq. [4] assuming straight and uniform capillaries.
From Fig. 2 we would expect a higher air entry pressure for the Conv.T. aggregate than for the grassland aggregate. This corresponds well to the calculated retention functions in Fig. 4
where the air-entry pressure is found at pF values (pF = common logarithm of matric suction in hPa) of
2.1 for the Conv.T. aggregate and
1.7 for the grassland aggregate, respectively. Note that at pF values
2.8 water contents drop to zero. This is due to the finite image resolution, which limits PSD calculations to pores with diameters
3 to 5 µm. Mind that retention curves represent merely the "main" pore network excluding unresolved matrix pores. Air capacities at matric potentials between pF 1.8 to 2.5 (field capacity) range from 0 to 6.5% for the Conv.T. aggregate and 1.0 to 9.4% for the grassland aggregate. This suggests a better aeration of the grassland aggregate under wet conditions while aeration of the Conv.T. aggregate is limited to higher matric suctions.

View larger version (6K):
[in this window]
[in a new window]
|
Fig. 4. Estimated water retention curves for the nodal pores of the pore network. Note that the retention curves represent the percolating pore network only excluding matrix pores 3 to 5 µm due to the final resolution of the tomographic images.
|
|
Fluid flow between interconnected nodal pores is predominantly controlled by the size and frequency of bottlenecks. The throat-size distribution (TSD) therefore reflects the overall resistance for fluid movement within the pore network. Prodanovi
et al. (2007) have found a positive correlation of throat permeabilities with throat surface areas computed from Lattice-Boltzmann simulations. With respect to throat sizes we observe larger throat surface areas and effective throat radii for the grassland aggregate compared with the Conv.T. aggregate (Fig. 3). This suggests stronger flow impedance within the Conv.T. aggregate, while for the grassland aggregate channels are significantly wider at the pore necks and therefore more permeable. An estimate of the number of bottleneck type of pores in the pore network is given by the frequency distribution of throat to pore radius ratios (Fig. 5
). The Conv.T. aggregate reveals a higher number of low throat to pore radius ratios (bottlenecks) while ratios for the grassland aggregate are shifted toward higher values suggesting higher hydraulic conductivities between interconnected pores. This pore connectivity phenomenon partially explains the two orders of magnitude reduction of saturated hydraulic conductivity values that paralleled porosity reductions for aggregates from long-term conventionally tilled soils compared with nearby forest soils, reported by Park and Smucker (2005). A similar shape of the distribution function shown in Fig. 5 was observed for low porosity sandstone by Lindquist (2002). Similar to the Conv.T. aggregate the maximum of the distribution function was found at throat to pore ratios between 0.3 and 0.4. This may be typical for pore systems dominated by primary pores. In contrast the distribution function for the grassland aggregate had its maximum at a throat to pore radius ratio of
0.5. The shape of the distribution function was flatter and more positively skewed suggesting a higher fraction of larger secondary pores.

View larger version (9K):
[in this window]
[in a new window]
|
Fig. 5. Frequency distribution of effective throat/pore radius ratios indicating ink bottle type of pores. The lower the ratio the higher is the resistance to fluid flow between adjacent nodal pores.
|
|
Path length frequency distributions reveal an exponential decline of number of paths with increasing channel length for both aggregates (Fig. 6
). A similar exponential decline has also been observed for sandstone samples of various porosities (Lindquist et al., 2000). Maximum channel lengths are <200 µm for the Conv.T. aggregate and >600 µm for the grassland aggregate, respectively. The decrease in frequency with increasing path length is more pronounced for the Conv.T. aggregate than for the grassland aggregate underlining the predominance of short flow channels in the Conv.T. sample. The existence of a view elongated paths > 400 µm for the grassland aggregate in Fig. 6a is also evident from the three-dimensional medial axis plot in Fig. 7d
. Such continuous channels may represent preserved biopores formed by roots. Note that the paths shown in Fig. 7d assemble continuous channels even exceeding the longest paths calculated in Fig. 6a. A visualization of the throats revealed that narrow sections exist where throats intersect the channels (data not shown). Frequency distributions for flow paths widths expressed as burn numbers are shown in Fig. 6b. Higher burn numbers corresponding to wider channels are more frequent in the grassland aggregate, while there is a strong increase in number of very narrow channels (burn number < 4) for the Conv.T. aggregate (note the logarithmic scale in Fig. 6b). Path tortuosities, larger numbers indicating larger tortuosities, are higher, on average for the Conv.T. aggregate than for the grassland aggregate (Fig. 6c). Generally the frequency decreases with increasing path tortuosity except at a tortuosity index of 3.4 where the frequency rapidly increases again in both aggregates. This increase is assumed to be related either to strongly convoluted channels enveloping individual grains or micro-aggregates and/or to branching tributary channels in irregularly shaped complex pore spaces. Latter could explain the strong contribution of branch-branch channels to the high tortuosity index, which was observed for both aggregates (data not shown).

View larger version (6K):
[in this window]
[in a new window]
|
Fig. 6. Relative frequency distributions for a) path length of flow channels, b) channel widths (burn number) and c) path tortuosity for the two analyzed soil aggregates.
|
|
Visualization and Connectivity of the Three-Dimensional Pore Network
Some of the above quantified pore-space features are also reflected by three-dimensional plots of the medial axis path networks (Fig. 7a-d). By removing isolated and dead-end paths with length
100 µm from the medial axis path network the continuity, connectivity and tortuosity of major flow channels is more easily recognizable (Fig. 7c and 7d). From a visual inspection of the three-dimensional plots we observe that the Conv.T. aggregate consists of numerous tortuous, sometimes strongly convoluted channels whereas the grassland aggregate is represented by an interconnected network of more continuous flow channels. In most cases also a greater channel diameter is indicated for the continuous flow paths shown by the color coding of medial axis voxels representing their burn number. Narrow channels are shown in red while wider channels are represented by green, blue, and purple colors (Fig. 7).
The transport of fluids over longer distances through the pore network depends not only on path continuity and width but also on the tortuosity of shortest connecting flow channels between opposite boundaries of the flow region (distance tortuosity). To get an idea of the distance tortuosity of the pore network we have calculated shortest possible flow channels connecting opposite faces of the analyzed subvolume. Further, anisotropy of the geometry of the flow network is evaluated by comparing channel tortuosities for different orientations of flow. This is done in two directions for the cubical subvolume shown in Fig. 7: for left to right flow (LR) and for down to upward flow (DU). Table 3
shows that tortuosity indices depend on flow direction with lower values for LR flow compared with DU flow suggesting that the pore network is anisotropic. Mean tortuosities are higher for the Conv.T. aggregate than for the grassland aggregate in both flow directions. This suggests in combination with the wider and more continuous flow channels (Fig. 6 and 7) a more enhanced fluid transport through and into the intra-aggregate pore space for the grassland aggregate compared with the Conv.T. aggregate. Park and Smucker (2005) identified a significant reduction in saturated hydraulic conductivity through aggregates of a similar size fraction for a conventional-tilled soil compared with a nearby non-tilled soil. Our results suggest that this could be explained by a more homogeneous structure of the intra-aggregate pore space following long term tillage which leads to a change in pore architecture mainly attributed to a reduction of continuous secondary pores.
View this table:
[in this window]
[in a new window]
|
Table 3. Distance tortuosities for the two analysed aggregates depending on direction of flow. LR denotes flow from left to right direction and DU denotes flow from down to upward direction in Fig. 7.
|
|
Some additional basic pore statistical information of the two analyzed soil aggregates is summarized in Table 4
. The total number of pores and attached throats is substantially higher for the Conv.T. aggregate corresponding to the smaller pore sizes (also note Fig. 3). About 6 to 7% of the pores are intersecting the image boundary representing a volume fraction of 25.8% for the Conv.T. aggregate and, owing to the larger pore sizes, 35.0% for the grassland aggregate. Higher fractions of boundary connected pore volumes enhance the exchange of water and gases between interaggregate and the intra-aggregate pore volumes and thus also trigger biophysical processes in the aggregate interior (Smucker et al., 2007). On the other hand the ratio of volume to surface area of the solid–void interface (representing only the percolating pore network) is a factor of 4.3 higher for the grassland aggregate compared with the Conv.T. aggregate. This value could be an important measure concerning sorption and desorption intensities along the percolating network with lower values indicating more accessible surface area per volume.
Pore Network Geometries and Soil Functions
Pore network geometries determine the biophysical and biogeochemical environment within soil aggregates controlling distribution, sorption, retention and release of soil nutrients, soil organic C (SOC), and contaminants (Smucker et al., 2007). Particularly transport properties and connectivity of flow channels are of paramount importance to understand microscale mass transfer and storage functions in natural soils. It has been shown in previous studies that the percolating pore space is anisotropic at soil core scales (Doerner, 2005; Kutilek and Nielsen, 2007). In our study we could show that also on the aggregate scale such anisotropy is evident and that it is most pronounced at a higher degree of structural development (grassland aggregate). Aggregate scale anisotropy not only affects water transport but also gas diffusion rates where path architectures have a stronger effect on the relative diffusion coefficient (ratio between diffusion coefficients in soil to that in free air) than air-filled porosity (Liu et al., 2006). Consequently, oxygen diffusion into the aggregates interior region should be enhanced for the well-structured grassland aggregate while oxygen depleted zones in the aggregate center of the Conv.T. aggregate are to be expected. This may have a strong impact on microbial environments and biogeochemical processes. Horn and Smucker (2005) pointed out that respiration rates are significantly greater for exterior soil layers of aggregates than their interior regions. They also found, however, that C respiration for soil aggregates under conventional tillage were 56% greater on the aggregate surfaces compared with their interior, while for no till sites respiration rates where only 28% greater at the aggregate surface than in the interior of the aggregate. Their results are likely to be related to a more continuous and interconnected intra-aggregate pore network for the no till aggregates which supports the supply of oxygen and distribution of substrate into the central parts of well- structured aggregates (Smucker et al., 2007). Concurrently, soil solutions are more effectively transported to the aggregates interior by the continuous, less tortuous and wider flow paths of the grassland aggregate. This could promote, for example, C fluxes into intra-aggregates regions via continuous macropores where SOC may diffuse into micropores (Park et al., 2007) and be stabilized by forming organo–mineral complexes (Priesack and Kisser-Priesack, 1993). This biogeochemical mechanism strengthens intra-aggregate structures that could enhance C sequestration (Park et al., 2007). In summary effects of soil disturbance by tillage on intra-aggregate pore network geometries, exemplarily demonstrated in this study, are in line with the conceptual understanding of the development of intra-aggregate pore domains as a function of land use (Smucker et al., 2007).
 |
CONCLUSIONS
|
|---|
In this study we have characterized and quantified transport relevant morphological features of the pore-space network of two contrasting soil aggregates on the microscale. We found distinct differences particularly in PSD and TSD, in channel length and width as well as in tortuosity of connecting flow paths. Our results suggest that the Conv.T. aggregate could suffer from impeded gas diffusion into the aggregate while gas and water transport is probably less restricted in the grassland aggregate. This may help explaining management induced differences in microbial soil functions and C sequestration potentials reported in the literature.
In further studies we will combine non-invasive x-ray microtromography and three-dimensional image analysis with physical measurements on the undisturbed imaged sample (e.g., imbibition studies, unsaturated hydraulic conductivity, oxygen partial pressure distribution and oxygen diffusion rates). This should facilitate linking pore-space geometry to soil functions. Growing knowledge on three-dimensional microscale pore-space features of structured natural soils subject to different management will not only improve our understanding of the soil microbial environment and its physical properties but is an inevitable step for developing models and verifying simulation results on water and gas transport within the intra-aggregate pore space. The quantification of these relationships is clearly needed to enhance our ability to predict changes in soil ecosystems due to management and global change (Six et al., 2004).
 |
ACKNOWLEDGMENTS
|
|---|
The authors thank HASYLAB (Hamburger Synchrotron Strahlungslabor) at DESY (Deutsche Elektronen-Synchrotron) and the Helmholtz-Association for supporting the use of the radiation source under contract no. I-05-020. We also like to acknowledge the useful comments of the anonymous reviewers.
 |
NOTES
|
|---|
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher. Permission for printing and for reprinting the material contained herein has been obtained by the publisher.
Received for publication April 6, 2007.
 |
REFERENCES
|
|---|
- Beckmann, F., U. Bonse, and T. Biermann. 1999. New developments in attenuation and phase-contrast microtomography using synchrotron radiation with low and high photon energies. p. 179–187. In Proc. Int. Soc. Opt. Eng., Conf. on developments in x-ray tomography II. Denver, CO. SPIE, Billingham, WA.
- Beckmann, F., J. Vollbrandt, T. Donath, H.W. Schmitz, and A. Schreyer. 2005. Neutron and synchrotron radiation tomography: New tools for materials science at the GKSS-Research Center. Nuclear Instr. Meth. Phy. Res. A 542:279–282.[CrossRef]
- Boone, F.R. 1988. Weather and other environmental factors influencing crop responses to tillage and traffic. Soil Till. Res. 11:283–324.[CrossRef]
- Bouwer, E.J., and J.B. Zehnder. 1993. Bioremediation of organic compounds-putting microbial metabolism to work. Trends Biotechnol. 11:360–367.[CrossRef][Web of Science][Medline]
- Capowiez, Y., A. Pierret, O. Daniel, M. P., and A. Kretzschmar. 1998. 3D skeleton reconstructions of natural earthworm burrow systems using CAT scan images of soil cores. Biol. Fertil. Soils 27:51–59.[CrossRef]
- Chenu, C., J. Hassink, and J. Bloem. 2001. Short-term changes in the spatial distribution of microorganisms in soil aggregates as affected by glucose addition. Biol. Fertil. Soils 34:349–356.[CrossRef]
- De Wever, H., D.T. Strong, and R. Merckx. 2004. A system for studying the dynamics of gaseous emissions in response to changes in soil matric potential. Soil Sci. Soc. Am. J. 68:1242–1248.[Abstract/Free Full Text]
- Dexter, A.R. 1991. Amelioration of soil by natural processes. Soil Tillage Res. 20:87–100.[CrossRef]
- Doerner, J. 2005. Anisotropie von Bodenstrukturen und Porenfunktionen in Böden und deren Auswirkungen auf Transportprozesse im gesättigten und ungesättigten Zustand. PhD Thesis. Christina-Albrechts University of Kiel, Germany.
- Dorioz, J.M., M. Robert, and C. Chenu. 1993. The role of roots, fungi and bacteria on clay particle organization. An experimental approach. Geoderma 56:179–194.[CrossRef][Web of Science]
- Flannery, B.P., H.W. Deckman, W.G. Roberge, and K.L. D'Amico. 1987. Three-dimensional x-ray micro-tomography. Science 237:1439–1444.[Abstract/Free Full Text]
- Gerke, H.H. 2006. Preferential flow description for structured soils. J. Plant Nutr. Soil Sci. 169:382–400.[CrossRef]
- Hillel, D. 1998. Environmental soil physics. Academic Press, London.
- Horgan, G.W. 1998. Mathematical morphology for analysing soil structure from images. Eur. J. Soil Sci. 49:161–173.[CrossRef]
- Horn, R., and A.J.M. Smucker. 2005. Structure formation and its consequences for gas and water transport in unsaturated arable and forest soils. Soil Tillage Res. 82:5–14.[CrossRef]
- Horn, R., W. Stepniewski, T. Wlodarczyk, G. Walenzik, and E.F.M. Eckhardt. 1994. Denitrification rate and micorbial distribution within homogeneous soil aggregates. Int. Agrophys. 8:65–74.
- Jasinska, E. 2006. Management effects on carbon distribution in soil aggregates and its consequences on water repellency and mechanical strength. PhD thesis. Christian-Albrechts University of Kiel, Germany.
- Journel, A.G. 1989. Fundamentals of geostatistics in five lessons. Short Course in Geology. Vol. 8 American Geophysical Union, Washington, DC.
- Ketcham, R.A., and W.D. Carlson. 2001. Acquisition, optimization and interpretation of X-ray computed tomographic imagery: Applications to the geosciences. Comput. Geosci. 27:381–400.[CrossRef]
- Kutilek, M., and D.R. Nielsen. 2007. Interdisciplinarity of hydropedology. Geoderma 138:252–260.[CrossRef][Web of Science]
- Lee, T.-C., R. Kashyap, and C.-N. Chu. 1994. Building skeleton models via 3-D medial surface/axis thinning algorithms. Graph. Mod. Im. Proc. 56:462–478.[CrossRef]
- Lindquist, W.B. 2002. Network flow model studies and 3D pore structure. Contemp. Math. 295:355–366.
- Lindquist, W.B., S.-M. Lee, D.A. Coker, K.W. Jones, and P. Spanne. 1995. Medial axis analysis of three dimensional tomographic images of drill core samples. J. Geophys. Res. 101B:8297–8310.[CrossRef]
- Lindquist, W.B., and A. Venkatarangan. 1999. Investigating 3D geometry of porous media from high resolution images. Phys. Chem. Earth (A) 24:593–599.[CrossRef]
- Lindquist, W.B., A. Venkatarangan, J. Dunsmuir, and T. Wong. 2000. Pore and throat size distribution measured from synchrotron X-ray tomographic images of Fontainebleau sandstones. J. Geophys. Res. 105B:21508–21528.
- Liu, G., B.G. Li, K.L. Hu, and M.T. van Genuchten. 2006. Simulating the gas diffusion coefficient in macropore network images: Influence of soil pore morphology. Soil Sci. Soc. Am. J. 70:1252–1261.[Abstract/Free Full Text]
- Mardia, K.V., and T.J. Hainsworth. 1988. A spatial thresholding method for image segmentation. IEEE Trans. Patt. Anal. Mach. Int. 10:919–927.[CrossRef]
- Nunan, N., K. Ritz, M. Rivers, D.S. Feeney, and I.M. Young. 2006. Investigating microbial micro-habitat structure using x-ray computed tomography. Geoderma 133:398–407.[CrossRef][Web of Science]
- Oh, W., and W.B. Lindquist. 1999. Image thresholding by indicator kriging. IEEE Trans. Patt. Anal. Mach. Int. 21:1–13.
- Pal, N.R., and S.K. Pal. 1993. A review of image segmentation techniques. Patt. Recogn. 29:1277–1294.
- Park, E.-J., and A.J.M. Smucker. 2005. saturated hydraulic conductivity and porosity within macroaggregates modified by tillage. Soil Sci. Soc. Am. J. 69:38–45.[Abstract/Free Full Text]
- Park, E.-J., W.J. Sul, and A.J.M. Smucker. 2007. Glucose additions to aggregates subjected to drying/wetting cycles promote carbon sequestration and aggregate stability. Soil Biol. Biochem. 39:2758–2768.[CrossRef]
- Priesack, E., and G.M. Kisser-Priesack. 1993. Modelling diffusion and microbial uptake of 13C-glucose in soil aggregates. Geoderma 56:561–573.[CrossRef][Web of Science]
- Prodanovi
, M. 2005. Fluid displacement in rock cores: A study based on three dimensional x-ray microtomography images. PhD thesis. Applied Mathematics and Statistics: Stony Brook University, New York. - Prodanovi
, M., W.B. Lindquist, and R.S. Seright. 2007. 3D image-based characterization of fluid displacement in a Berea core. Adv. Water Res. 30:214–226.[CrossRef] - Serra, J. 1982. Image analysis and mathematical morphology. Academic Press, London.
- Six, J., H. Bossuyt, S. Degryze, and K. Denef. 2004. A history of research on the link between (micro)aggregates, soil biota, and soil organic matter dynamics. Soil Tillage Res. 79:7–31.[CrossRef]
- Smucker, A.J.M., E.-J. Park, J. Dorner, and R. Horn. 2007. Soil micropore development and contributions to soluble carbon transport within macroaggregates. Vadose Zone J. 6:282–290.[Abstract/Free Full Text]
- Soil Survey Staff. 2006. Keys to soil taxonomy. 10th ed. USDA-NRCS, Washington, DC.
- Soille, P. 2003. Morphological image analysis—Principles and applications. Springer, Berlin, Germany.
- Stock, S.R. 1999. X-ray microtomography of materials. Int. Mater. Rev. 44:141–164.[CrossRef]
- Tokunaga, T.K., J.M. Wan, T.C. Hazen, E. Schwartz, M.K. Firestone, S.R. Sutton, M. Newville, K.R. Olson, A. Lanzirotti, and W. Rao. 2003. Distribution of chromium contamination and microbial activity in soil aggregates. J. Environ. Qual. 32:541–549.[Abstract/Free Full Text]
- Van Gestel, M., R. Merckx, and K. Vlassak. 1996. Spatial distribution of microbial biomass in microaggregates of a silty-loam soil and the relation with the resistance of microorganisms to soil drying. Soil Biol. Biochem. 28:503–510.[CrossRef]
- Wildenschild, D., J.W. Hopmans, C.M.P. Vaz, M.L. Rivers, D. Rikard, and B.S.B. Christensen. 2002. Using X-ray computed tomography in hydrology: Systems, resolutions, and limitations. J. Hydrol. 267:285–297.[CrossRef]
- Young, I.M., and J.W. Crawford. 2004. Interactions and Self-Organization in the Soil-Microbe Complex. Science 304:1634–1637.[Abstract/Free Full Text]