INVESTIGATING SOIL AGGREGATE PORE STRUCTURES AND THEIR
RELATIONSHIP TO BACTERIA SPATIAL DISTRIBUTION USING X-RAY
COMPUTED MICROTOMOGRAPHY
By
Wei Wang

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Crop and Soil Sciences
2012

ABSTRACT
INVESTIGATING SOIL AGGREGATE PORE STRUCTURES AND THEIR
RELATIONSHIP TO BACTERIA SPATIAL DISTRIBUTION USING X-RAY
COMPUTED MICROTOMOGRAPHY
By
Wei Wang
Soils are among the extremely complicated materials with intricate internal pore geometries.
Understanding the complex interaction between soil aggregate structure and soil processes will
greatly help us appreciate critical processes at larger scale. However, due to technical limitations
on accessing aggregate interiors, it remained relatively poorly understood until recent advances
of X-ray computed microtomography (µCT) imaging in soils. Overall, the dissertation aims to
effectively utilize soil aggregate µCT imaging, to understand soil aggregate structure difference
under contrasting land use and managements, and to investigate bacteria redistribution and
transport pattern and its relation to aggregate pore characteristics.
Synchrotron µCT provides a noninvasive approach to record three dimension (3D)
information on soil aggregates with resolution of one to several microns. Image segmentation is
the first step in µCT image analysis; however, lack of ground-truth images poses great challenge
to soil image segmentation. Proposed simulation method successfully generated soil aggregate
grey scale images from pore/solid binary images, thus provided ground-truth information.
Indicator Kriging method was found preferable to other methods in this study when the image
histogram is clearly bimodal. Region non-uniformity measure performs sufficiently well as a
segmentation criterion. The second technical difficulty involves 3D soil aggregate boundary
delineation. Image closing technique is a promising tool for boundary detection, yet it requires

selecting optimal neighborhood parameter (r) so that aggregate exteriors can be properly defined
– inclusion of all surface pores while keeping the surface sufficiently irregular.
After addressing these image processing issues, soil aggregates from long term
contrasting soil types and land use practices were compared. The studied soil managements are
conventional tillage (CT), native succession vegetation (NS), and bare soil with no vegetation
(BS). Significantly greater percent of pores > 15 µm in LTER-CT aggregates was observed as
compared to those of LTER-NS. LTER-NS aggregates had more large pores (> 90 µm) than
LTER-CT aggregates, while more medium sized pores (45-90 µm) were found in LTER-CT
aggregates. Pores were larger in the aggregate interiors than in the exterior layers; while
intermediate pores (45-90 µm) were more abundant in the aggregate exterior layers. These
results implied that a general mechanism maybe responsible for aggregates formation, but long
term land use could alter the magnitude and intensity of involved soil processes.
The last part of the dissertation addressed E.coli redistribution and transport in soil
aggregates. E.coli redistribution within soil aggregates display a significant different spatial
distribution pattern under air-dry and saturated flow condition. When E.coli was first applied to
an air-dry aggregate, its resulting spatial distribution was highly heterogeneous in aggregates of
all management practices. After saturation, equilibration, and water extraction the E.coli
redistributed within the aggregate primarily in vertical direction and became more
homogeneously spread. Only relatively small percent of E.coli has completely left the
aggregates. Redistribution was most pronounced in CT aggregates, followed by NT, and was
almost negligible in NS aggregates. E.coli spatial distribution was related to intra-aggregate pore
characteristics; however, because of high variability of the intra-aggregate pores the relationships
were weak.

I dedicated this dissertation to my wife, Ming Gu, for her love and support.

iv

ACKNOWLEDGEMENTS

I am very grateful to Dr. Sasha Kravchenko for her advice during my Ph.D. studies. Before
entering the Ph.D. program, I have very little knowledge about soil science. Her guidance and
encouragement made me a successful transition from mathematics to the cutting edge application
of X-ray imaging in soil science. She was always approachable and ready to help in every
possible way to make sure I am on the right track of the research. Doing research under her
guidance is really my pleasure and I hope I can work with her in the future.
I would like to thank Dr. Alvin Smucker, Dr. Yimin Xiao and Dr. Stephan Peth for being
on my Ph.D. committee, for their suggestions, comments, and continual support of my research.
Thanks to Dr. Smucker for teaching me theory of soil physics and introducing me to the soil Xray imaging. Thanks to Dr. Peth for detailed instruction on using the imaging softwares and
many insightful discussions on soil physics. Thanks to Dr. Xiao for many helpful suggestions
during my study and teaching several of my statistical classes which are quite helpful.
I would like to thank Tim Johnson, Sangeetha Srinivasan and Dr. Joan Rose for their help
on conducting many soil microbiology experiments which is an important part of this
dissertation.
I would also like to thank to my lab mates, Xuwen Huang, Xinmei Hao, Senthil
Subramanian, Juan David Munoz-Robayo, Kateryna Anayeva, Hyen Chung Chun, Moslem
Ladoni, Andrew Worth, Melissa Erickson and Richard Price. I am grateful to know these fellows
and enjoy the discussion in the lab. I had a good time with them during my graduate studies.

v

I also appreciate my friends in MSU – Xinxin Wang, Minjie Xu, Lei Zhang, Xiaoqin
Tang, Rui Lin, Hao Wen, Tao Feng, Ravi Krishna Shaga and Pikul Sarkar – for their friendship,
support and motivation.
My special thanks go to my wife, Dr. Ming Gu, for her unconditional love and support
during my graduate studies. We met, got to know each other, fell in love and got married at
MSU. The time we spent together at MSU pursuing out Ph.D.s was the best part in our lives.
This dissertation is also dedicated to my parents and Ming’s parents for their continuous support.

vi

TABLE OF CONTENTS
LIST OF TABLES ......................................................................................................................... ix
LIST OF FIGURES ........................................................................................................................ x
Chapter 1: INTRODUCTION......................................................................................................... 1!
BIBLIOGRAPHY ....................................................................................................................... 6
Chapter 2: COMPARISON OF IMAGE SEGMENTATION METHODS IN SIMULATED 2D
AND 3D MICROTOMOGRAPHIC IMAGES OF SOIL AGGREGATES .................................. 9!
ABSTRACT................................................................................................................................ 9!
INTRODUCTION .................................................................................................................... 10!
MATERIALS AND METHODS.............................................................................................. 13!
Soil sample collection ........................................................................................................... 13!
Synchrotron X-ray microtomography of soil aggregates...................................................... 14!
Generation of simulated grayscale images ........................................................................... 14!
Segmentation methods .......................................................................................................... 19!
Majority filtering ................................................................................................................... 25!
Segmentation performance criteria ....................................................................................... 26!
RESULTS AND DISCUSSION ............................................................................................... 28!
Parameter selection for indicator kriging .............................................................................. 30!
Comparing performance of the segmentation methods using NU and ME criteria .............. 33!
Comparing performance of the segmentation methods using pore characteristics ............... 39!
CONCLUSION ......................................................................................................................... 44!
BIBLIOGRAPHY ..................................................................................................................... 47
Chapter 3: INTRA-AGGREGATE PORE CHARACTERISTICS UNDER LONG-TERM
CONVENTIONAL TILLAGE AND NATIVE SUCCESSION VEGETATION: X-RAY
COMPUTED MICROTOMOGRAPHY ...................................................................................... 52!
ABSTRACT.............................................................................................................................. 52!
INTRODUCTION .................................................................................................................... 53!
MATERIALS AND METHODS.............................................................................................. 56!
Soil sample collection ........................................................................................................... 56!
Soil aggregate image acquisition ......................................................................................... 58!
Defining aggregate boundaries ............................................................................................. 59!
vii

Image segmentaiton .............................................................................................................. 65!
Pore characteristics ............................................................................................................... 65!
Identification of aggregate layers.......................................................................................... 67!
Statistical analyses ................................................................................................................ 68!
RESULTS AND DISCUSSION ............................................................................................... 69!
Aggregate boundary identification ....................................................................................... 69!
Pore characteristics in whole aggregates .............................................................................. 74!
Pore characteristics in aggregate layers ................................................................................ 84!
CONCLUSION ......................................................................................................................... 90!
BIBLIOGRAPHY ..................................................................................................................... 94
Chapter 4: ESCHERICHIA COLI REDISTRIBUTTION AND TRANSPORT IN INTACT SOIL
MACRO-AGGREGATES ............................................................................................................ 99!
ABSTRACT.............................................................................................................................. 99!
INTRODUCTION .................................................................................................................. 100!
MATERIALS AND METHODS............................................................................................ 102!
Soil sample collection ......................................................................................................... 102!
Soil aggregate image acquisition ....................................................................................... 103!
Aggregate image analysis ................................................................................................... 104!
Alignment of the physical and virtual cuttings on aggregates ............................................ 105!
E.coli application experiments ............................................................................................ 106!
qPCR analysis ..................................................................................................................... 109!
Statistical analysis ............................................................................................................... 111!
RESULTS AND DISCUSSION ............................................................................................. 111!
E.coli distribution and transport in the flow and slicing experiment .................................. 111!
Relationship between E.coli distribution and pore characteristics ..................................... 116!
CONCLUSION ....................................................................................................................... 125!
BIBLIOGRAPHY ................................................................................................................... 128
Chapter 5: CONCLUSION ......................................................................................................... 131!

viii

LIST OF TABLES

a

Table 2.1 : Estimated parameters of the mixed Gaussian distribution in the simulated 2D images
(images 1-8) and the 3D image. .................................................................................................... 31
Table 2.2: Assessing performance of NU in selecting parameter r for IK segmentation relative to
accuracy of image segmentation. .................................................................................................. 32
Table 2.3: Thresholds obtained from entropy based method, iterative method and Otsu’s method
for the simulated 2D (images 1-8) and 3D images. ...................................................................... 34
Table 2.4: Summary of (a) ME and (b) NU using different segmentation methods for the
simulated 2D (images 1-8) and 3D images. Numbers in bold represent best segmentation results
according to the each segmentation performance criterion. ......................................................... 38
Table 2.5: Number of images (total 9 images) where the segmentation methods produced the
‘worst’ result in terms of pore characteristics. .............................................................................. 42
Table 3.1: List of soil aggregates physical and chemical properties from top soil (0-20cm),
including texture, total C, N, total porosities, macro-porosities (pores > 14.6 µm), and porosities
of pores < 14.6 µm. ....................................................................................................................... 75
Table 3.2: Pore characteristics within exterior, transitional and interior layer of the aggregates.
(a) Macro-porosity; (b) Mean burn number on MA. .................................................................... 85
Table 4.1: Primer and probe sequence used in this study. .......................................................... 110
Table 4.2: Descriptive statistics for pore characteristics for soil aggregates with size 4-6.3mm
from conventional tillage (T1), no-till (T2) and native succession (T7). Standard deviations are
shown in parentheses. Burn number equal to 1 (BN1) is equivalent to pore size of 15-30 µm, BN
equal to 2 (BN2) is equivalent to pore size of 45-60 µm, BN equal to 3 (BN3) is equivalent to
pore size of 75-90 µm, and BN equal to or greater than 4 (BN4g) is equivalent to pore size > 105
µm. .............................................................................................................................................. 116!

ix

LIST OF FIGURES

Figure 2.1: Illustration of image simulation process. White represents the solid space and black
(gray) represents the pore space. (a) Binary image of the “true” image pixel size (14.6x14.6 µm);
(b) Binary image of the pores visible only at pixel size 7.3x7.3 µm;(c) Binary image of the pores
visible only at pixel size 3.7x3.7 µm; (d) Binary image of the pores visible only at pixel size
1.9x1.9 µm; (e) Gray scale image of the pores visible at pixel size 14.6x14.6 µm obtained by
combing a), b), c), and d); (f) Final simulated grayscale image. .................................................. 17
Figure 2.2: Examples of histograms of grayscale value distributions in two of the simulated
images. (a) Image 7 with well defined pore peak; and (b) Image 1 with poorly defined pore peak.
The blue curve represents the fitted Gaussian distribution corresponding to the pore space. The
red curve represents the fitted Gaussian distribution corresponding to the solid space. For
interpretation of the references to color in this and all other figures, the reader is referred to the
electronic version of this dissertation. .......................................................................................... 24
Figure 2.3: Simulated grayscale images with different porosities at 14.6x14.6 µm scale. ........... 29
Figure 2.4: Segmented images using different segmentation methods at 14.6x14.6 µm scale. ... 35
Figure 2.5: 3D image stack with (a) ground-truth image and segmented using (b) IK method with
r selected based on NU; (c) entropy method; (d) iterative method; (e) Otsu’s method. ............... 37
Figure 2.6: Summary of pore characteristics (a) porosity; (b) mean pore size; (c) pore surface
area obtained from images segmented using different segmentation methods............................. 40
Figure 2.7: Comparing accuracy of NU-based segmentation method selection to pore
characteristics-based selection. The relative errors are defined as the differences between the
pore characteristic value from the segmented image and that from ground-truth value divided by
the ground-truth value. .................................................................................................................. 44
Figure 3.1: Example of (a) 2D unprocessed soil aggregate grayscale image slice; (b) the same
image with the background in black and the aggregate in white where aggregate boundary has
been identified using an automated flood/seed fill algorithm; note that the intra-aggregate pore
connected to the outside is incorrectly classified as belonging to the background (red square); (c)
the same image after boundary sealing algorithm has been applied with r = 12; (d) the final
processed gray-scale image with the background in white; (e) the final processed grey-scale
image with the background in white with r = 15; note the artificial pore pixels added to the sides
of the image (red squares). ............................................................................................................ 60
Figure 3.2: An example of 2D image closing using a disk-shaped structuring element: a) the
original image; b) dilating the image a) (the portion added in a course of dilation is shown in
light blue); c) eroding the image b) (the pixels that were not present in the original image are
shown in red); d) the final image. ................................................................................................. 62
x

Figure 3.3: An example of a spherical structuring element with r = 2. a) Voxels of the structuring
element are showing in red; b) central portion of the spherical structuring element (a) is showing
in red. ............................................................................................................................................ 63
Figure 3.4: A schematic representation of medial axes and burn numbers. Black pixels represent
solid materials and white pixels represent pore space. Rainbow color scale pixels are medial axis
paths with red for burn number equal to 2, yellow for burn number equal to 3, green for burn
number equal to 4, etc. .................................................................................................................. 66
Figure 3.5: An example of the three layers of aggregates in a 2D image slice. Light, medium, and
dark gray intensities represent, respectively, exterior, transitional, and interior layers. Notice that
the layers were identified in 3D, therefore they might not look like they are of equal width when
shown in a 2D representation. ....................................................................................................... 68
Figure 3.6: (a) Total volume of the aggregates after sealing the boundaries with different r values
in four studied aggregates. Solid lines were fitted to the linear parts of the volume/r relationships;
(b) fractal dimension of the aggregate boundaries of the four studied aggregates at different r
values; (c) aggregate volume increase rate† at different r values; (d) reduction rate† of fractal
dimension at different r values. .................................................................................................... 70
Figure 3.7: Percent of pore space visible at 14.6µm resolution occupied by pores of different
burn numbers in aggregates from (a) LTER-CT and LTER-NS treatments; and (b) CERN-BS
and CERN-NS treatments. Burn number equal to 1 is equivalent to pore size of 15-30µm, BN
equal to 2 is equivalent to pore size of 45-60µm, BN equal to 3 is equivalent to pore size of 7590µm, etc. For each burn number, the letters mark the significant differences between land uses
while ‘ns’ indicates no statistically significant difference (p<0.05). Std. fraction on the y-axis is
the standardized burn number distribution. Note that for clarity purposes different scales of yaxis values were used for fractions of pores with different burn numbers. .................................. 78
Figure 3.8: Examples of 2D image slices from central portions of the aggregates from LTER-CT
(a-c), LTER-NS (d-f), CERN-NS (g-h) and CERN-BS (i-j) soils. ............................................... 81
Figure 3.9: Percent of pore space visible at >14.6 µm resolution occupied by pores of different
burn numbers in aggregates exterior (Ext), transitional (Tran), and interior (Int) layers of (a)
LTER-CT and LTER-NS treatments; and (b) CERN-BS and CERN-NS treatments. Burn number
equal to 1 is equivalent to pore size of 15-30 µm, BN equal to 2 is equivalent to pore size of 4560 µm, BN equal to 3 is equivalent to pore size of 75-90 µm, etc. For each burn number, the
lowercase letters mark significant differences among the three layers within each treatment. The
differences among the layers are not statistically significant (p<0.05) if lowercase letters are not
displayed. For each burn number, the uppercase letters mark significant differences between the
treatments within each layer, while no upper case letters are shown when the differences are not
statistically significant (p<0.05). The cases when no statistically significant difference existed
both among layers and between treatments (p<0.05) are marked with ‘ns’. Std. fraction on the yaxis is the standardized burn number distribution. Note that for clarity purposes different scales
of y-axis values were used for fractions of pores with different burn numbers. .......................... 87
xi

Figure 4.1: Schematic representation of the aggregate cutting. Numbers 1, 2, 3, 4, 5, 6 and 7
correspond to top, bottom, left, right, front, back and center slices, respectively. The sequence of
the numbers represents the cuttings in order............................................................................... 106
Figure 4.2: A picture shows the glass beads matrix chamber system – glass beads matrix
chamber for adding bacteria then flushing and vacuum extracting E.coli from soil aggregates. 107
Figure 4.3: Percent of E.coli in the effluents and glass beads in the flow experiment. .............. 112
Figure 4.4: Comparison of relative E.coli ratio in flow and slicing experiments. ...................... 114
Figure 4.5: Bi-plot of soil aggregate samples and pore characteristics from PCA for a) the whole
soil aggregates and b) middle portion of the soil aggregates including center and sides. The threedigit numbers are labels for soil aggregates with ‘1xx’ representing T1 aggregates, ‘2xx’
representing T2 aggregates and ‘7xx’ representing T7 aggregates. Soil management was
converted to indicator dummy variables CT (T1) and NS (T7). Burn number equal to 1 (BN1) is
equivalent to pore size of 15-30 µm, BN equal to 2 (BN2) is equivalent to pore size of 45-60 µm,
BN equal to 3 (BN3) is equivalent to pore size of 75-90 µm, and BN equal to or greater than 4
(BN4g) is equivalent to pore size > 105 µm. .............................................................................. 118
Figure 4.6: (a) Correlation loading plot from PLSR between pore characteristics (porosity,
tortuosity, maxflow, paths, BN1-BN4g) and E.coli spatial distribution (TOP, CENTER, SIDES,
BOTTOM). The numbers are representing aggregate ids: 1-5 represent T1 aggregates, 6-10
represent T2 aggregates and 11-15 represent T7 aggregates. (b) R2 plots for percent of variation
explained vs. number of factors in PLSR. .................................................................................. 121
Figure 4.7: The relationship between difference in E.coli concentration between top and bottom
slices and max flow/tortuosity ratio. Gray circles mark high influential points. Numbers along
with the circles are aggregates IDs. ............................................................................................ 123
Figure 4.8: 3D visualization of aggregates pore space for pieces. Pores with bright color indicate
high E.coli concentration. (a) Aggregate 230; (b) Aggregate 722; (c) Aggregate 125. ............ 124!

xii

Chapter 1: INTRODUCTION
Soils are among the extremely complicated materials with intricate internal pore geometries.
Distribution of countless pores within soil aggregates is highly spatially heterogeneous. Soil’s
hierarchical structures control various physical, chemical and biological processes at different
scales. At the soil aggregate scale, configuration of pore space plays an important role in
determining carbon (C) sequestration, water and gas flow as well as bacteria survival and
transport. Understanding the complex interaction between soil aggregate structure and these
processes will greatly help us appreciate critical processes at larger scale, e.g., C spatial
distribution, global C cycle, and subsurface water contamination caused by bacteria outbreaks.
However, due to technical limitations on accessing aggregate interiors, it remained relatively
poorly understood until recent advances of X-ray computed microtomography (µCT) imaging in
soils.
Synchrotron µCT provides a noninvasive approach to record three dimension (3D)
information on soil aggregates with resolution of one to several microns. Last decade witnessed
fast developments in soil imaging applications (Capowiez et al., 1998; Pierret et al., 1999, 2002;
Dathe et al., 2006; Peth et al., 2008; Udawatta et al., 2008; Kravchenko et al., 2009, 2011). Even
though µCT imaging is promising and powerful tool, it is only recently being applied to soil
science and there is a large number of technical difficulties that need to be resolved. Accurate
image analysis techniques are indispensable to getting meaningful and reliable results. The
overall goal of my research was to effectively utilizing soil µCT images to understand soil
aggregate structures under different land use and soil managements and their relationship to soil
1

processes, for which I focused on bacteria distribution and transport in soil aggregates. However,
prior to reaching it I had to combat a number of technical difficulties in image processing and
develop techniques for solving them.
A technique of particularly importance is image segmentation and so far a number of
methods have been published (Sezgin and Sankur, 2004; Iassonov et al., 2009). Soil aggregate
pore characteristics, such as porosity, pore-size distribution, tortuosity and connectivity can be
calculated or derived from image analysis. The accuracy of these parameters strongly depends on
the goodness of segmentation methods – the technique used to identify the pore and solid space
in soil aggregates. However, pronounced artifacts such as partial volume effect and noisy nature
of µCT pose great challenges to interpreting soil image segmentation results. Application of
segmentation on soil images was still at starting stage and different methods yielded highly
variable results (Iassonov et al., 2009; Baveye et al., 2010). One of the reasons was that soil
matrix is much more diverse in origins and densities of its solid materials than many other
natural substrates. Another difficulty is that unlike other discipline that use µCT imaging (e.g.,
sandstones (Lindquist et al., 2000); hydrate hosted sediments (Jones et al., 2007)), in soil images
there is no ground-truth information regarding about pore/solid classification. Therefore, new
approaches are needed in order to indentify best segmentation methods in soils.
The second technical difficulty I have to conquer involves 3D soil aggregate boundary
delineation. Until now, most of the aggregate image studies have been conducted by selecting a
regularly shaped region of interest within the aggregate’s central part or soil cores (Dathe et al.,
2006; Peth et al., 2008; Udawatta et al., 2008; Baveye et al., 2010). This precluded full
assessment of the aggregate’s pore structure and the differences in flow and transport properties
of the exterior and interior layers and interactions between them. Image closing technique is a
2

promising tool for boundary detection (Nunan et al., 2006), yet it requires selecting optimal
neighborhood parameter so that aggregate exteriors can be properly defined – inclusion of all
surface pores while keeping the surface sufficiently irregular.
After addressing these image processing issues, soil aggregates from long term
contrasting soil types and land use practices were compared. Long-term differences in land use
and management generate notable changes in soil physical and hydraulic properties, including
changes in soil organic matter content, porosity and hydraulic conductivity (Brye and Pirani,
2005). Surface soil disturbance through tillage alters overall soil aggregation, stability, biological
properties and distribution of C in aggregates (Lal et al., 1998, 2002; Six et al., 1999; Green et al.,
2007). Besides soil disturbance, presence of vegetation and vegetation type can also alter soil
properties through plant roots activities and their association with soil organisms (Wienhold and
Tanaka, 2001; Udawatta et al., 2008). However, precise linkage between these observations and
soil pore structures are still under development.
A well known challenge to researchers’ perception is on soil aggregates formation theory.
It should be pointed out that despite an enormous amount of work in the aggregate studies, there
is no generally accepted understanding of how the aggregates are internally assembled and what
the intra-aggregate pore structure characteristics are (Tisdall and Oades 1982; Oades et al, 1984;
Augustin et al., 1995; Horn et al., 1994; Chenu et al., 2001; Park and Smucker, 2005). Therefore,
this dissertation is aimed to fill these knowledge gaps by investigating the pore structures among
aggregates and within individual aggregate layers.
Last but not least, the dissertation studied E.coli redistribution and transport in soil
aggregates. E.coli has been regarded as an important indicator of environmental quality and
3

public health risk (CDC; Vogt and Dippold, 2005), and understanding mechanisms controlling
transport and redistribution upon precipitation in soils is crucial for agricultural practices and
food safety (Jamieson et al., 2002; Guber et al., 2005). Till now, it has been found that soil
structure and nutrient distribution at macro and micro scale are believed to be responsible for
E.coli survival and transport within soil matrix (Pachepsky et al., 2006), and precipitation
facilitates infiltration of bacteria suspensions through flow pathways into deeper ground (Jenkins
et al., 2008). Soil organic matter can hinder E.coli attaching to soil particles due to its negative
charges (Johnson and Logan (1996)), whereas high clay content in soils may promote bacteria
attachment (Huysman and Verstrate, 1993). Though findings above were mostly at field and soil
column scale, it is highly possible that E.coli travel within aggregates and further release from
the aggregates are strongly associated with pore structures they passed en route. However,
lacking of appropriate tools and modeling techniques, few past studies addressed E.coli
distribution and transport at micro scale in soils.
The dissertation is organized as follows: In Chapter 2, a simulation approach for
generating soil aggregate µCT images based on ground-truth information is proposed to compare
several segmentation methods and criteria to assess segmentation performance were discussed.
Chapter 2 studies the optimal soil aggregate boundary identification and compared soil aggregate
pore structures difference for aggregates from contrasting soil type and land use managements.
Spatial gradient of pore characteristics within individual aggregate were also discussed. Chapter
3 discusses multivariate relationship among soil pore characteristics, E.coli redistribution pattern
under saturated condition, and investigates the relationship between E.coli spatial concentration
and soil aggregate pore characteristics. Chapter 5 states the major conclusions and points out
future research to extend the notions proposed in this study.
4

BIBLIOGRAPHY

5

BIBLIOGRAPHY
Augustin, S., M. Jansen, E.Priesack, and F. Beese. 1995. Litter decomposition and matter
transport in beds of soil aggregates. In: Hartge, K.H., Stewart, B.A. (Eds.), Soil structure.
Its development and function. Advances in Soil Science. CRC Press, Inc., Boca Raton, pp.
237-256.
Baveye, P.C., M. Laba, W. Otten, D. Grinev, L. Bouckaert, P. Dello Starpaio, R.R. Goswami, Y.
Hu, J. Liu, S. Mooney, R. Pajor, S. Sleutel, A. Tarquis, W. Wang, Q. Wei, M. Sezgin.
2010. Observer-dependent variability of the thresholding step in the quantitative analysis of
soil images and X-ray microtomography data. Geoderma 157: 51-63.
Brye, K.R., and A.L. Pirani. 2005. Native soil quality and the effects of tillage in the grand
prairie region of eastern Arkansas. Am. Midl. Nat. 154: 28-41.
Capowiez, Y., A. Pierret, O. Daniel, P. Monestiez, and A.Kretzschmar. 1998. 3D skeleton
reconstructions of natural earthworm burrow systems using CAT scan images of soil cores.
Biol. Fert. Soils 27: 51-59.
CDC, Division of Bacterial and Mycotic Diseases.
http://www.cdc.gov/nczved/divisions/dfbmd/diseases/ecoli_o157h7/. (Accessed
09/26/2011).
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.
Dathe, A., A.M. Tarquis, and E.Perrier. 2006. Multifractal analysis of the pore- and solid-phases
in binary two-dimensional images of natural porous structures. Geoderma 134: 318-326.
Green, V.S., D.E. Stott , J.C. Cruz, and N. Curi. 2007. Tillage impacts on soil biological activity
and aggregation in a Brazilian Cerrado Oxisol. Soil Till. Res. 92: 114-121.
Guber, A. K., D. R. Shelton, and Y. A. Pachepsky. 2005. Transport and retention of manurebome coliforms in undisturbed soil columns. Vadose Zone J. 4:828-837.
Horn, R., W. Stepniewski, T. Wlodarczyk, G. Walenzik, and E.F.M. Eckhardt. 1994.
Denitrification rate and microbial distribution within homogeneous soil aggregates. Int.
Agrophys. 8: 65–74.
Huysman, F., W. Verstraete. 1993. Water-facilitated transport of bacteria in unsaturated soil
columns: Influence of cell surface hydrophobicity and soil properties. Soil Biology and
Biochemistry 25 (1): 83-90.

6

Iassonov, P., T. Gebrenegus, and M. Tuller. 2009. Segmentation of X-ray computed tomography
images of porous materials: a crucial step for characterization and quantitative analysis of
pore structures. Water Resour. Res. 45, W09415.
Jamieson, R. C., R. J. Gordon, K. E. Sharples, G. W Stratton, and A. Madani. 2002. Movement
and persistence of fecal bacteria in agricultural soils and subsurface drainage water: A
review. Can. Biosyst. Eng. 44:1.1-1.9.
Jenkins M. B., C. C. Truman, G. Siragusa, E.Line, J. S. Bailey, J. Frye, D. M. Endale, D. H.
Franklin, H.H. Schomberg, D. S. Fisher, and R. R. Sharpe. 2008. Rainfall and tillage
effects on transport of fecal bacteria and sex hormones 17β-estradiol and testosterone from
broiler litter applications to a Georgia Piedmont Ultisol. Science of the total environment
403:154-163.
Johnson, W.P., and B.E. Logan. 1996. Enhanced transport of bacteria in porous media by
sediment-phase and aqueous-phase natural organic matter. Water Res. 30:923–931.
Jones, K.W., H. Feng, S. Tomov, W.J. Winters, M. Prodanovic, and D. Mahajan. 2007.
Characterization of methane hydrate host sediments using synchrotron-computed
microtomography (CMT). Journal of Petroleum Science and Engineering 56: 136-145.
Kravchenko, A.N., W. Wang, A. J. M. Smucker, and M. L. Rivers. 2011. Long-term Differences
in Tillage and Land Use Affect Intra-aggregate Pore Heterogeneity. Soil Sci. Soc. Am. J.
75: 1658-1666.
Kravchenko, A. N., M. A. Martín, A. J. M. Smucker, and M L. Rivers. 2009. Limitations in
Determining Multifractal Spectra from Pore–Solid Soil Aggregate Images. Vadose Zone
Journal 8: 220-226.
Lal, R., J. Kimble, R. Follet, and C. Cole. 1998. The potential of US cropland to sequester
carbon and mitigate the greenhouse effect. Ann Arbor Press, Chelsea, MI.
Lal, R. 2002. Soil carbon dynamics in cropland and rangeland. Environ. Pollut. 116: 353-362.
Lindquist, W.B., A. Venkatarangan, J. Dunsmuir, and T.F. Wong. 2000. Pore and throat size
distributions measured from synchrotron X-ray tomographic images of Fontainebleau
sandstones. J. Geophys. Res. 105B: 21508.
Mooney, S.J., A.R. Tams, and P. Berry. 2006. A reliable method for preserving soil structure in
the field for subsequent morphological examinations. Geoderma 133: 338-344.
Mooney, S.J., C. Morris, J. Craigon, and P. Berry. 2007. Quantification of soil structural changes
induced by cereal anchorage failure: Image analysis of thin sections. J. Plant Nutr. Soil Sc.
170 (3): 363-372.

7

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.
Oades, J.M. 1984. Soil organic matter and structural stability: mechanisms and implications for
management. Plant Soil 76: 319-337.
Pachepsky, Y. A., A. M. Sadeghi, S. A. Bradford, D. R. Shelton, A. K. Guber, and T. Dao. 2006.
Transport and fate of rnanure-borne pathogens: Modeling perspective. Agric. Water
Manag. 66:81-92.
Park, E.J., and A.J.M. Smucker. 2005. Saturated hydraulic conductivity through
macroaggregates modified by tillage. Soil Sci. Soc. Am. J. 69:38-45.
Peth, S., R. Horn, F. Beckman, T. Donath, J. Fisher, and A.J.M. Smucker. 2008. ThreeDimensional quantification of intra-aggregate pore space features using synchrotronradiation-based microtomography. Soil Sci. Soc. Am. J. 72: 897-907.
Pierret, A., Y. Capowiez, C. J. Moran, and A. Kretzschmar. 1999. X-ray computed tomography
to quantify tree rooting spatial distributions. Geoderma 90: 307-326.
Pierret, A., Y. Capowiez, L. Belzunces, and C.J. Moran. 2002. 3D reconstruction and
quantification of macropores using X-ray computed tomography and image analysis.
Geoderma 106: 247-271.
Six, J., E. T. Elliott, and K. Paustian.1999. Aggregate and soil organic matter dynamics under
conventional and no-tillage systems. Soil Sci. Soc. Am. J. 63: 1350-1358.
Tisdall, J.M., and J.M. Oades. 1982. Organic matter and water-stable aggregates in soils. Journal
of Soil Science 33: 141-163.
Udawatta, R.P., S. H. Anderson, C. J. Gantzer, and H. E. Garrett. 2008. Influence of prairie
restoration on CT-measured soil pore characteristics. J. Environ. Qual. 37:219-228.
Vogt, R.L., and L. Dippold. 2005. Escherichia coli O157:H7 outbreak associated with
consumption of ground beef, June–July 2002. Public Health Rep 120 (2): 174–178.
Wienhold, B.J., and D.L. Tanaka. 2001. Soil property changes during conservation from
perennial vegetation to annual cropping. Soil Sci. Soc. Am. J. 65:1795-1803.

8

Chapter 2: COMPARISON OF IMAGE
SEGMENTATION METHODS IN SIMULATED 2D
AND 3D MICROTOMOGRAPHIC IMAGES OF
SOIL AGGREGATES

ABSTRACT
Advances in X-ray microtomography (µCT) are opening new opportunities for
examining soil pore structures. However, usefulness of µCT data for pore structure
characterization depends on how accurately the grayscale images are segmented into binary pore
and solid components. Multiple segmentation algorithms have been developed for a variety of
physical and biological substrates. However, it is not clear how applicable they are for µCT
images of soils, since soil is among the most complex porous media due to diversity in origins
and densities of composing solids. One of the difficulties in comparing the performance of
segmentation algorithms is the lack of ground-truth information in the soil samples subjected to
µCT. This means that only the criteria that do not depend on the availability of the ground-truth
data can be used in assessing performance of the segmentation methods, yet the reliability of
such criteria in soil images is unclear. In this study to resolve the problem of the lack of groundtruth information we simulated 2D and 3D soil images. The objectives of the study were (i) to
compare several commonly used segmentation methods, namely, entropy-based method, iterative
method, Otsu’s method, and indicator kriging (IK) method; and (ii) to evaluate performance of
the region non-uniformity measure (NU), the criterion that does not depend on presence of the
ground-truth image, and compare it with the ground-truth-based misclassification error criterion.
We found that though there was no single segmentation method that preserved pore
9

characteristics in all the cases, IK method yielded segmented images most similar to the groundtruth in most of the cases when histogram of image grayscale values had clearly distinguishable
peaks. For the image with poorly distinguishable histogram peaks the IK did not perform well,
while clustering methods produced acceptable segmentation results. The results indicated that
selecting a segmentation method based on NU did not always produce optimal representation of
pore characteristics. However, overall, the NU was found to be an acceptable criterion for
segmentation method selection in µCT soil images.

INTRODUCTION
Image segmentation into pores and solids is a crucial step in most soil image analyses
(Vogel and Kretzschmar, 1996; Capowiez et al., 1998; Lehmann et al., 2006; Mooney et al.,
2006, 2007; Papadopoulos et al., 2008a, 2008b; Kaestner et al., 2008; Schapp and Tuller, 2010;
Piñuela et al. 2010; Baveye et al., 2010; Vaz, et al., 2011; Vogel et al., 2010). Without accurate
segmentation, results of subsequent data analyses might be misleading.
The segmentation methods that have been applied in geosciences can be broadly
classified as global thresholding methods that generate a single threshold applied to the entire
image; and local thresholding methods that adapt the threshold(s) value depending on local
image characteristics.
Among a large number of global thresholding methods developed, the ones that appear to
be most promising for soil applications are clustering methods and entropy-based methods
(Sezgin and Sankur, 2004; Iassonov et al., 2009; Iassonov and Tuller, 2010). Sezgin and Sankur
(2004) reported the best performance of clustering methods and entropy-based methods among

10

40 segmentation methods in thermal/ultrasonic/light-microscope imaging applications. In a
recent study, Iassonov et al. (2009) evaluated 14 segmentation methods for application to
industrial and synchrotron computed tomography images of soils, sand-bentonite mixtures, and
precision glass beads, and reported Otsu’s and iterative method to be the best among the tested
global thresholding methods. However, global thresholding methods were often found to be
inferior to the local segmentation methods (Oh and Lindquist, 1999; Peth et al. 2008; Iassonov et
al., 2009).
Among many existing local segmentation methods (Pal and Pal, 1993; Oh and Lindquist,
1999; Lehmann et al., 2006; Elliot and Heck, 2007), a two-pass indicator kriging (IK) has been
shown to significantly reduce the misclassification errors (Oh and Lindquist, 1999). However,
one of the difficulties with IK application is that, unlike many automated segmentation methods,
IK requires user-specified input on its parameter values and despite popularity of IK
segmentation in soil studies (Peth et al., 2008, 2010; Udawatta et al., 2008), as of now there are
no criteria or guidelines on proper selection of the IK parameter values.
A significant difficulty in comparing accuracy of different segmentation methods for soil
images is the lack of the ground-truth binary pore/solid information. When ground-truth
information is available, the accuracy of segmentation methods can be assessed using a variety
of discrepancy based criteria, among which are: misclassification error (ME), edge mismatch,
and, Rand index and its extensions (e.g., normalized probability Rand index) (Yasnoff et al.,
1977; Sezgin and Sankur, 2004; Unnikrishnan et al., 2005). Zhang (1996) has shown that
ground-truth based segmentation criteria, which account for misclassified voxels and specific
ground-truth based features of the segmented objects, are among the most powerful tools for

11

selecting segmentation methods. However, none of these criteria can be applied to soil images
because of the impossibility of obtaining ground-truth information.
Because of the lack of ground-truth information for selecting segmentation methods in
soil images the researchers have to use criteria that do not depend on the ground-truth. A number
of such criteria have been proposed, e.g., region non-uniformity (NU) measure (Levine and
Nazif, 1985; Sezgin and Sankur, 2004), inter-region contrast (Levine and Nazif, 1985) and shape
measure of voxels’ neighborhood (Sahoo et al., 1988). Among these criteria, NU has the
potential of being particularly useful for soil X-ray imaging because of its simplicity and direct
measure of the homogeneity of the regions. Several color image segmentation studies (Liu and
Yang, 1994; Borsotti et al., 1998) had reported good matches between NU/inter-region contrast
and visual assessments. However, how well NU criterion will perform for selecting optimal
image segmentation method in soils is not known.
Information from soil imagery that is critical for modeling flow and transport processes
in soils includes various soil pore characteristics, e.g., pore size distributions, pore connectivities
and tortuosities. When working with artificial systems, e.g., glass beads or collections of
capillary tubes, some of these characteristics (porosities, solid surfaces, and pore size
distributions) can be calculated based on the properties of the system. Thus the performance of
segmentation methods can be assessed based on how closely they match the theoretically
calculated values (Iassonov et al., 2009; Porter et al., 2010). However, in an undisturbed soil
sample it is not possible to accurately assess porosities, surface areas, or distributions of the
pores visible at the X-ray image resolutions used in the study (Kasteel et al. 2000). Even though
the total porosity of the sample can be accurately measured (Black and Hartage, 1986), the total
porosity is based not only on the pores visible at the studied X-ray resolution, but on the entire
12

body of soil pores, including those invisible at the studied resolution. Soil pore size distributions
can be obtained experimentally. However, a variety of factors, including hysteresis and
experimental measurement errors, make them also unreliable as segmentation criteria.
In this study we will generate simulated 2D and 3D X-ray soil images with a wide range
of soil porosities such that they mimic features of the real soil images (Zhang, 2001; Schlüter et
al., 2010). Then, the ground-truth information available from the simulated images will be used
(i) to explore optimal parameter selection for IK segmentation; (ii) to compare the accuracy of
several commonly used segmentation methods, namely, entropy based method, iterative method,
Otsu’s method, and IK method; and (iii) to evaluate the performance of NU criterion in choosing
the optimal segmentation method for soil images when the ground-truth information is not
available.

MATERIALS AND METHODS

Soil sample collection
Soil aggregates were collected at the Long Term Ecological Research Site, Kellogg
Biological Station in southwest Michigan (85o 24` W, 42 o 24` N). The soils are well-drained
Typic Hapludalfs of the Kalamazoo (fine-loamy, mixed, mesic) and Oshtemo (coarse-loamy,
mixed, mesic) series, developed on glacial outwash. Soil blocks approximately 15 x 15 cm in
size were extracted from 0-20 cm depth using a sharp flat spade. Air-dried soil samples were
manually sieved for 30 seconds and the aggregate fractions, 4 to 6.3 mm, were retained for this
study.

13

Synchrotron X-ray microtomography of soil aggregates
Image acquisition
Aggregate image data were collected on the monochromatic bending magnet beam line,
station 13-BM-D of the GeoSoilEnvironCARS (GSECARS) at Advanced Photon Source,
Argonne National Laboratory, IL. X-ray beam line of 28 keV incident energy was used to scan
the aggregate samples. 3D image consisting of 520 slices with 696x696 pixels per slice were
generated. The voxel size was 14.6 µm. Detailed information on image acquisition is given in
Rivers et al. (1999). The scanned grayscale images were used to generate binarized ground-truth
images.

Data processing
Raw data was processed using GSECARS tomography processing software (Details are
available at http://cars9.uchicago.edu/software/idl/tomography.html). Data processing consisted
of preprocessing and reconstruction. Preprocessing of the 2D images read the raw data from
camera and wrote the data as a single 3D file. Ring artifacts were removed and a fast Fourier
transform algorithm was used to reconstruct the data (Rivers, 1998). After reconstruction, the
grayscale values (GV) ranged from 0 (black) to 255 (white) corresponding to low and high X-ray
attenuation, respectively.

Generation of simulated grayscale images
Please note that for all the following discussions we refer to ‘pixel’ for 2D images, and
‘voxel’ for 3D images.
14

2D cases
To assess the accuracy of different classification methods and compare different
segmentation criteria, we generated eight 2D binary soil aggregate images that served as the
ground-truth image information. In order to assess a wide variety of soil porosity conditions, we
generated the images with porosity ranging from low (3-5%) to high (22-28%).
Generation of the simulated grayscale images consisted of combining several image
layers. The first pore/solid image layer corresponded to the scale of the true aggregate images
(pixels of 14.6x14.6 µm) and was 100x100 pixels in size (Figure 2.1(a)). The pore/solid data of
this layer represented the soil pore structure visible at the studied scanning resolution.
In µCT images of real soil samples, the actual pore structure visible at the studied
resolution is obscured by a number of influences. These influences include partial volume effects
due to the presence of pores and solids in a single pixel (Ketcham and Carlson, 2001); variability
in densities and X-ray attenuations among the pixels of the solid material (e.g. quartz versus
organic matter); and noise generated during scanning. In generating the images we simulated
contribution of these three influences.
The partial volume effect was simulated by generating three layers representing a range
of pore sizes smaller than the scanned 14.6x14.6 µm pixel size (Figure 2.1(b-d)). The smaller
scale layers had pixel sizes of 7.3x7.3, 3.7x3.7, and 1.9x1.9 µm, respectively. To ensure realistic
representation of the soil pore structure we used manually thresholded µCT images from real soil
samples to produce the binary images for these three layers. When the size of the small scale
images exceeded that of the available aggregate images from the µCT scan, several binarized
15

µCT images were combined. For example, in order to sufficiently cover the 100x100 pixel
14.6x14.6 µm layer, the 1.9x1.9 µm layer had to consist of 800x800 pixels. To generate the
1.9x1.9 µm layer, four images from real soil samples 400x400 pixels in size were joined and
then scaled down in their physical size to match 100x100 14.6x14.6 µm layer (Figure 2.1(d)).
Simulation of the partial volume effect proceeded as follows: first, the binary pore information of
the 1.9x1.9 µm resolution image (Figure 2.1(d)) was overlaid with that of the 3.7x3.7 µm image
(Figure 2.1(c)). Then, corresponding 1.9x1.9 µm values were averaged for each 3.7x3.7 µm
image pixel. Then the 3.7x3.7 µm image was overlaid by and averaged across pixels of the
7.3x7.3 µm image (Figure 2.1(b)) and then with 14.6x14.6 µm image (Figure 2.1(a)). Their
averaging to the actual 14.6x14.6 µm pixel size produced the grayscale image that reflects the
partial volume effect (Figure 2.1(e)).

16

Figure 2.1: Illustration of image simulation process. White represents the solid space and black
(gray) represents the pore space. (a) Binary image of the “true” image pixel size (14.6x14.6 µm);
(b) Binary image of the pores visible only at pixel size 7.3x7.3 µm;(c) Binary image of the pores
visible only at pixel size 3.7x3.7 µm; (d) Binary image of the pores visible only at pixel size
1.9x1.9 µm; (e) Gray scale image of the pores visible at pixel size 14.6x14.6 µm obtained by
combing a), b), c), and d); (f) Final simulated grayscale image.

(a)

(b)

(c)

(d)

(e)

(f)

After simulating partial volume effect, the resulting image (Figure 2.1(e)) was used to
simulate variability in atomic densities and X-ray attenuations among the pixels of the solid
material. The pixels that contained no pore space (pure ‘white’ (GV = 255) pixels on Figure
2.1(e)) were regarded as pure solid space. We assumed that spatial correlation was present in the
solid material density distribution, that is, the pixels were more similar in their grayscale values
to the pixels that were closer in space than to those that were farther away. To represent the
spatial correlation, simulations based on LU decomposition technique (Christakos, 1992) were
17

performed using PROC SIM2D (SAS Inc., 2001). LU decomposition technique was employed
because it was proven to be computationally efficient in spatial simulations (Deutsch and
Journel, 1998). The mean values and covariance structures for LU simulations were determined
by analyzing actual µCT images of soil aggregates using subroutine GAM in GSLIB (Deutsch
and Journel, 1998). Pixels that contained pore space (‘gray’ pixels on Figure 2.1(e)) were kept
unchanged in the simulations.
To simulate scanning noise, after LU simulation, we added Gaussian noise. Please notice
that we manually truncated the final GV within [0, 255] to ensure its validity. The resulting
simulated grayscale image which was used for testing segmentation methods is shown on Figure
2.1(f).

A 3D case
Because of time and labor consuming process of generating simulated 3D images only
one 3D image was produced. The same procedure described in detail for 2D simulation was used
for 3D. The first pore/solid image layer corresponded to the scale of the true aggregate images
(voxels of 14.6x14.6x14.6 µm) and was 50x50x50 voxels in size. As in the case of 2D, the
pore/solid data of this layer represented the soil pore structure visible at the studied scanning
resolution. Partial volume effect was simulated by generating three image stacks with voxel sizes
of 7.3x7.3x7.3, 3.7x3.7x3.7, and 1.9x1.9x1.9 µm, respectively. These smaller scale stacks were
combined in the same fashion as was done for 2D cases. The variability in densities in solid
space was simulated using simulation annealing technique (Deutsch and Journel, 1998)

18

implemented in WinGslib V1.5 (http://www.statios.com/WinGslib/, Accessed December, 08,
2010).

Segmentation methods
Global image segmentation methods used in this study were Renyi’s entropy-based
method (Sahoo et al, 1997), iterative method (Ridler and Calvard, 1978) and Otsu’s method
(Otsu, 1979). For this study we selected several of the existing methods that have received most
positive feedback in the reviews, namely, we used entropy based and clustering based methods
and indicator kriging approach. Entropy and clustering are global threshold procedures (Sezgin
et al., 2004) that demonstrated potential for applications to soil CMT image segmentation.
Indicator kriging method is a local threshold integrated approach that considers both the
grayscale histogram and local spatial dependence structure, which has been successfully applied
in rock and fluid research (Oh and Lindquist, 1999). Here we only provide brief descriptions of
the segmentation procedures used in this study along with references to the sources where
detailed method descriptions could be found.

Entropy based method
Originated from thermodynamics, the paradigm of entropy in information theory was first
established by Shannon (1948) in order to quantify the information and uncertainty in
communication networks. Many different generalizations, such as Renyi’s entropy and Tsallis
entropy were developed thereafter by parameterizing the Shannon entropy (Renyi 1961; Tsallis,
2001). Jaynes (1957) proposed the principle of maximum entropy which stated that the
probability distribution on a proposition is the true probability distribution when it maximizes the
19

information entropy, providing that all the possible probability distribution encode this
information. This idea has been successfully applied to image segmentation field (Kapur et al.,
1985; Chang et al., 1995; Sahoo et al., 1997). The image segmentation using entropy based
method consists of finding a global threshold, such that maximizes the entropy of the two classes
of objects, i.e., pore and solid pixels.
Renyi’s entropy method (Sahoo et al., 1997) is a generalization of the maximum entropy
sum method (Kapur et al., 1985) and entropic correlation method (Chang et al., 1995) by
incorporating a parameter α in Shannon’s entropy (Shannon, 1948). The Renyi’s entropies of
the pore and solid space are calculated as

α
H Λ1 =

α
H Λ2

respectively, where

( pi %
1
ln ∑ &
& p(Λ ) #
#
1 − α i =0 '
1 $
t

( pi %
1
=
ln ∑ &
& p(Λ ) #
#
1 − α i =t +1'
2 $
255

(2.1)

α
(2.2)

pi (i = 0,1,,255) denotes the probability distribution of image GV, t is
t

the threshold,

α

p(Λ1 ) = ∑i =0 pi ,

255

p(Λ 2 ) = ∑i =t +1 pi are

the proportion of pore

and solid pixels in the whole image, and α > 0 is the parameter that generalizes Shannon’s
entropy.
Reparameterization of α ( 0 < α < 1 , α → 1 , α > 1 ) in Renyi’s entropy allows us to
obtain three distinct thresholds where two correspond to the thresholds obtained from maximum
20

entropy sum method (Kapur et al., 1985) and entropic correlation method (Chang et al., 1995).
Then a weighted average of the three thresholds is used as a threshold to segment the entire
image. The weights are determined by the relative differences among the previously found three
thresholds. This procedure was implemented using a program developed in MATLAB.

Clustering methods
This class of methods applies clustering analysis to the pore and solid pixels in order to
minimize the error of classifying a pore pixel as a solid pixel or vice versa. These methods are
helpful when it is difficult to identify the threshold due to the overlapping GV’s range of the pore
and solid pixels (Kittler and J. Illingworth, 1985 & 1986).

Iterative method
The iterative clustering method allows selecting the threshold automatically (Riddler et
al., 1978). In this approach, an initial estimate of the threshold is selected so that the image is
partitioned into “presumably” pore and solid spaces. Typically the initial estimate is the midpoint
between the minimum and maximum GV of the image. Then the new threshold is calculated as

T=

1
(µ1 + µ 2 ) ,
2

where

µ1

and

µ 2 are

the means of the two spaces. The iterations

terminate when the difference in the thresholds of the successive iteration is sufficiently small.
This procedure was implemented in MATLAB.

The Otsu’s method

21

The clustering method proposed by Otsu (1979) calculates the optimal threshold
separating the pore and solid spaces so that their GV’s intra-class variance is minimal. The GV’s
intra-class variance is defined as
2
σ w (t ) = ω 0 (t )σ 02 (t ) + ω1 (t )σ 12 (t )

Where

ω0 (t ) and ω1 (t ) are

separated by the threshold t ,

(2.3)

the proportion of pore and solid pixels in the whole image
2
σ 0 (t ) and σ 12 (t ) are

the variances for pore and solid

pixels respectively. This procedure was implemented in MATLAB by stepping through each
possible GV and finding the thresholds corresponding to the minimum

2
σ w (t ) .

Indicator kriging method
Indicator kriging is the local adaptive segmentation method where each voxel is classified
as either pore or solid based on its GV and local spatial correlation structure (Oh and Lindquist,
1999). The method involves selecting two thresholds such that the voxels with GV below the
lower threshold,

T0 , and the voxels above the upper threshold, T1 , are classified as pores and

solids, respectively. Voxels with intermediate grayscale values are classified into pores or solids
based on their proximity to previously classified voxels through indicator kriging system.
The two threshold values can be obtained by fitting the image GV histogram with mixed
Gaussian distributions. For example, if a GV histogram has a clear bimodal pattern, then we can
reasonably assume the existence of two Gaussian distributions corresponding to GV of pore and
solid space respectively. The parameters of the mixed Gaussian distribution can be estimated
using Expectation-Maximization algorithm which computes the maximum likelihood estimates
22

of the distribution parameters (Dempster et al., 1977). The histogram of GV has the form (Oh
and Lindquist, 1999) of

2

f ( x) = ∑ wi
i =1

1
−( x − µi ) 2 /( 2σ i 2 )
exp
2π σ i

(2.4)

Example of the fitted mixed Gaussian distribution for GV in a simulated image is shown in
Figure 2.2 (a). The two thresholds,

T0

T0 = min( z 0 , z1 ),

and

T1 , can be calculated as (Oh and Lindquist, 1999):

T1 = max( z 0 , z1 )

z 0 = min(µ1 + rσ 1 , µ 2 ) , z1 = min(µ 2 − rσ 2 , µ1 )
where weights ( wi ), means ( µ i ) and variances ( σ

2

(2.5)

) are the parameters estimated from

Expectation-Maximization algorithm. The user-specified parameter r controls the distance
between

T0

and

T1 . We explored the effect of parameter

r on segmentation accuracy. For that,

several r values ranging from 0 to 1.96, i.e., 0, 0.50, 1.00, 1.25, 1.50, 1.75, 1.96 were tried,
where r of 0 corresponds to setting the thresholds at the histogram peaks and r of 1.96
corresponds to thresholds separating the 5% lowest and 5% highest values in solid and pore
distributions, respectively.
The thresholds were calculated using Expectation-Maximization algorithm implemented
in R (R Development Core Team, 2006). After threshold determination, the IK segmentation was
conducted using segmentation module in 3DMA-Rock software (Lindquist, 2002).

23

Figure 2.2: Examples of histograms of grayscale value distributions in two of the simulated
images. (a) Image 7 with well defined pore peak; and (b) Image 1 with poorly defined pore peak.
The blue curve represents the fitted Gaussian distribution corresponding to the pore space. The
red curve represents the fitted Gaussian distribution corresponding to the solid space. For
interpretation of the references to color in this and all other figures, the reader is referred
to the electronic version of this dissertation.

(a)

0.008
0.006
0.004
0.002
0.000

Relative Frequency

0.010

Image 7

0

50

100

150
GV

24

200

250

Figure 2.2 (cont’d).
Image 1

0.010
0.005
0.000

Relative Frequency

0.015

(b)

0

50

100

150

200

250

GV

Majority filtering
After segmentation the images often have a number of artificial single voxel pores or
solids called ‘spots’, the removal of which benefits segmented image accuracy (Kaestner et al.,
2008; Peth et al., 2008). Majority filtering is an effective moving-window-based tool to smooth
the image by removing them. In this study we applied majority filtering following
25

recommendations of Oh and Lindquist (1999) with 3x3x3 filtering window for 3D and 3x3
window for 2D images and set the threshold for majority at 60%.

Segmentation performance criteria
In this study, we used two formal segmentation criteria, namely misclassification error
(ME) and region non-uniformity measure (NU). ME can be applied only when ground-truth
information is available and, arguably, is the most objective criterion of the segmentation method
performance (Zhang, 1996). NU can be used either with or without the ground-truth information.
In addition to using the formal criteria, we examined the effect of segmentation on
several pore characteristics, including porosity, mean pore size, and pore surface area.

Misclassification error
Misclassification error is based on enumerating pore voxels that are incorrectly classified
as solids and solid voxels incorrectly classified as pores (Yasnoff et al., 1977). It is calculated as

ME = 1 −

P+S
T

(2.6)

where P is the number of common pore voxels in the ground-truth image and the segmented
image, S is the number of common solid voxels in the ground-truth image and the segmented
image, T is the total number of voxels in the image. Misclassification error can range from 0,
corresponding to the perfect segmentation, to 1, corresponding to complete mismatch between
the ground-truth image and the segmented image.

Region non-uniformity measure
26

The region non-uniformity measure can be applied to assess segmentation quality even
when the ground-truth information is not available (Zhang, 1996, 2001). This criterion is
calculated as:
2

P σ
NU = ⋅ P2
T σ

(2.7)

where P is the number of pore voxels in the segmented image, T is the total number of voxels
in the image,

σ P2

is the variance of GV of pore space in the grayscale image,

σ 2 is the total

variance of GV in the grayscale image. The NU can range from 0, corresponding to the case
when the variance of the GV in voxels classified as pores is equal to zero, to 1, corresponding to
the case when all voxels of the image are classified as pores and thus

σ P2

is equal to

σ

2

.

Pore characteristics
Besides the criteria that describe general accuracy of pore/solid segmentation, i.e., ME
and NU, it is of interest to assess how different the segmentation methods are in representing
various characteristics of soil pore structure, such as porosity, pore size distribution and pore
surface area (Iassonov et al., 2009; Baveye et al., 2010). These characteristics are of utmost
importance for modeling soil water flow and solute transport, and such modeling is often
considered as the main goal of soil image analysis. Using simulated images with available
ground-truth information allows us to estimate how accurately segmented images obtained using
different methods reflect the pore characteristics. Note that pore characteristics discussed here
are those of the pores visible at the resolution of the scanned images used in the study (>14.6
µm).
27

Porosity is defined as the proportion of pore voxels in total number of voxels in the
image. Pore size is calculated as the maximum distance between each pore voxel and the closest
solid voxel. Mean pore size is the average of all pore size data. Pore surface area refers to the
contact surface area of pore voxels to adjacent classified solid voxels. All pore characteristics
were calculated using SCAMP V1.2 (developed in SIMBIOS Centre, University of Abertay,
Dundee, Scotland)) implemented in ImageJ (Rasband, 1997-2009).

RESULTS AND DISCUSSION
In order to compare the segmentation methods, we created 8 pore/solid binary images
with porosities ranging from low (3-5%) to high (22-28%), as well as one 3D image with
porosity 17% (Figure 2.3). The selected porosity values are representative of the porosities from
pores >14.6 µm typically observed in soil aggregates (Park and Smucker, 2005). Note that the
porosity of >14.6 µm pores is substantially lower than the total porosity, calculating using the
bulk density method, of the 4-6 mm aggregates from the studied soil, which ranges from 30% to
40% (Ananyeva et al., 2009, personal communications). The images with low, medium and high
porosities (images 1-6, Figure 2.3), contained disconnected small (~30 µm), medium (~100 µm)
and large (~200 µm) pores. The two images with the highest porosity (images 7-8, Figure 2.3)
each contained one large macropore.

28

Figure 2.3: Simulated grayscale images with different porosities at 14.6x14.6 µm scale.

Image 1

Image 3

Image 5

Image 7

Porosity = 0.048

Porosity = 0.078

Porosity = 0.165

Porosity = 0.228

Image 2

Image 4

Image 6

Image 8

Porosity = 0.036

Porosity = 0.083

Porosity = 0.158

Porosity = 0.285

Image 3D
Porosity = 0.170

29

Parameter selection for indicator kriging
The first step of the indicator kriging method (IK) application consisted of determining
the two thresholds ( T0 and T1 ). They were determined based on the parameters of the mixed
Gaussian distributions fitted to the data and based on the user-defined parameter r (Eq. (2.5)).
The mixed Gaussian distribution parameters, fitted using Expectation-Maximization algorithm
for GV in the simulated images are listed in Table 2.1. The estimated means of the Gaussian
distributions for the pore/solid space ( µ1 and µ 2 ) are quite consistent for images 2 to 8.
However, in image 1, the pore space has greater estimated mean (67.7) and larger variance
(1636.5). The main reason is that while the two peaks corresponding to pores and solids,
respectively, were clearly visible in histograms of images 2-8 and of the 3D image, they were
poorly distinguishable in image 1. The histogram shown on Figure 2.2(a) is typical of that
observed for images 2-8 and the 3D image with identifiable two peaks, while the pore peak is not
distinguishable in the histogram of image 1 (Figure 2.2(b)). The unreliable fitting of the
histogram for image 1 could significantly affect the accuracy of pore/solid classification by
providing too high value for T0 .

30

a

Table 2.1 : Estimated parameters of the mixed Gaussian distribution in the simulated 2D images
(images 1-8) and the 3D image.
Image

w2

µ1

µ2

σ 12

2
σ2

1
2
3
4
5
6
7
8
3D
a

w1
0.131
0.109
0.187
0.236
0.061
0.103
0.177
0.290
0.209

0.869
0.891
0.813
0.764
0.939
0.897
0.823
0.710
0.791

67.7
48.7
44.0
42.0
36.0
44.1
41.7
40.9
28.9

169.2
166.7
167.0
161.4
149.3
164.9
166.9
162.3
172.6

1636.5
852.4
744.5
617.1
630.8
702.7
676.8
594.4
900.2

689.3
830.8
853.1
1177.5
1050.4
913.6
893.5
1133.0
1700.4

The symbols in this table correspond to Eq. (2.4).

We examined how selection of r affects ME and specific pore characteristics in the IK
segmentation. Then we assessed how well the NU criterion performs as a tool for selecting r
value in the IK segmentation. Since NU is the measure of segmentation accuracy that does not
require ground-truth information, using it for r parameter selection in IK is of direct relevance to
soil image segmentation tasks.
In five of the images (images 1, 3, 4, 5 and 3D) the smallest ME and NU were obtained at
high r value (r >= 1.25) (Table 2.2 (a)). In three images (images 2, 6 and 8), selecting r value
based on NU would result in ME 9-14% higher than the ME that would be obtained if r value
were selected by minimizing the ME itself. In one of the images (image 7) the r selected based
on NU would produce a 25% higher ME. This indicates that even though in many cases NUbased selection can be regarded as only negligibly different from the ME-based selection, the
performances of the two criteria are not identical.

31

Table 2.2: Assessing performance of NU in selecting parameter r for IK segmentation relative to
accuracy of image segmentation.
(a) Relative error in ME;
Minimum
a
ME /r

Image

a
b
c

1
2
3
4
5
6
7
8
3D
Average

Minimum
b
NU /r

Relative errors in
c
ME

0.048/1.25
0.010/1.5
0.011/1.25
0.009/1.96
0.009/1.96
0.009/1.75
0.012/1.5
0.007/1.5
0.032/1.75

0.088/1.25
0.014/1.25
0.036/1.25
0.026/1.96
0.040/1.96
0.037/1.96
0.049/0.5
0.050/1.96
0.058/1.75

0
0.09
0
0
0
0.11
0.25
0.14
0
0.066

Minimal ME is obtained by selecting r such that it minimizes ME.
Minimal NU is obtained by selecting r such that it minimizes NU.

Relative error for ME is calculated as the difference between the minimal ME and the ME
obtained when r is selected by minimizing NU criterion. The difference is then divided by the
minimal ME.
(b) Relative errors in pore characteristics from IK segmentation with parameter r selected based
on NU and ME;
Image

Porosity

Mean pore size

Surface area

IK_ME

a

a

1
2
3
4
5
6
7
8
3D
Average

IK_NU

IK_ME

IK_NU

IK_ME

IK_NU

0.927
0.020
0.082
0.028
0.004
0.001
0.021
0.009
0.019
0.023

0.927
0.006
0.082
0.028
0.004
0.001
0
0.009
0.019
0.019

0.951
0.167
0.079
0.022
0.011
0.013
0.019
0.009
0.103
0.053

0.951
0.177
0.079
0.022
0.011
0.012
0.002
0.009
0.103
0.052

0.066
0.096
0.026
0.121
0.027
0.004
0.144
0.068
0.068
0.069

0.066
0.064
0.026
0.121
0.027
0.042
0.113
0.068
0.068
0.066

Due to relatively high segmentation error caused by lack of bimodality in the histogram, image
1 is not included in computing the averages.

32

We calculated relative differences between the pore characteristics obtained from the
images that were segmented using either NU or ME based r selection and those of the groundtruth images. Comparison between the pore characteristics from IK segmented images with r
values selected based on the NU values (IK_NU) and based on the ME values (IK_ME)
indicated that selecting r based on NU was advantageous for all the studied pore characteristics
(Table 2.2(b)). In almost all cases (except for mean pore size of image 2 and pore surface area of
image 6) the relative error from NU-based r value selection was less than or equal to that
obtained using the ME-based selection. This indicates that r selections in IK segmentation with
minimum ME would not necessarily be optimal for preserving the pore characteristics similar to
those of ground-truth images. The possible explanation is that NU measures the non-uniformity
over the pore space on the basis of GV variance of the pore space, which was evaluated at every
pore voxel. ME, on the other hand, depends on both the misclassified pore and solid voxels,
hence might not reflect the pore space properties as accurately as NU. Therefore, we recommend
NU, a pore space based measure, as a segmentation criteria for selecting r value in IK
segmentation.

Comparing performance of the segmentation methods using NU and ME criteria
The 2D segmented images obtained using IK method (with r value selected based on ME
and on NU), entropy method, iterative method, and Otsu’s method after 60% majority filtering
are presented in Figure 2.4 for images 1-8. The 3D images are presented in Figure 2.5. The
calculated thresholds for the global methods are listed in Table 2.3. For all but one image (image
2), the entropy method consistently yielded larger threshold values than both Otsu’s and iterative

33

method. This resulted in overestimated porosity in the images segmented using the entropy
method.

Table 2.3: Thresholds obtained from entropy based method, iterative method and Otsu’s
method for the simulated 2D (images 1-8) and 3D images.
Image
1
2
3
4
5
6
7
8
3D

Entropy Iterative
125
123
116
129
122
117
123
116
124
110
130
108
124
109
122
107
138
120

34

Otsu’s
122
130
116
115
109
108
108
106
119

Figure 2.4: Segmented images using different segmentation methods at 14.6x14.6 µm scale.
In row: images 1-8;
In column: (a) Ground-truth binary images; (b) IK method with r selected based on ME;
(c) IK method with r selected based on NU; (d) entropy method; (e) iterative method; (f) Otsu’s
method.
(a)

(b)

(c)

(d)

(1)

(2)

(3)

(4)

35

(e)

(f)

Figure 2.4 (cont’d).
(a)

(b)

(c)

(d)

(5)

(6)

(7)

(8)

36

(e)

(f)

Figure 2.5: 3D image stack with (a) ground-truth image and segmented using (b) IK method
with r selected based on NU; (c) entropy method; (d) iterative method; (e) Otsu’s method.

(a)

(b)

(d)

(e)

(c)

The assessment of segmentation performance using criteria ME and NU for different
methods for 2D images 1-8 and for the 3D image are listed in Table 2.4. We observed that
selecting best segmentation method based on ME and NU agreed for majority of the images.
Both ME and NU would suggest Otsu’s as the best segmentation method in image 1, and IK as
the best method in images 2-4, 7, 8 and 3D. Only in two of the studied images (images 5 and 6)
the two criteria disagreed. In these two images, the lowest ME was obtained from IK
segmentation and the lowest NU from iterative or Otsu’s methods. Visual assessment (Figure 2.4
37

(5), (6): (c) vs. (e)) suggested that the IK segmented images might be preferable to that of Otsu’s
method in these two images since they appear to better preserve similarity of pore features in
segmented images to those of ground-truth images. For example, pores marked by yellow circles
in Figure 2.4(5)(a) are mostly preserved in the IK segmented image (Figure 2.4(5)(a)), while
they are distorted in the Otsu’s segmented image (Figure 2.4(5)(a)).

Table 2.4: Summary of (a) ME and (b) NU using different segmentation methods for the
simulated 2D (images 1-8) and 3D images. Numbers in bold represent best segmentation results
according to the each segmentation performance criterion.
a) ME
Image
1
2
3
4
5
6
7
8
3D

IK_NU
0.048
0.011
0.011
0.009
0.009
0.010
0.015
0.008
0.020

Entropy Iterative
0.030
0.027
0.042
0.112
0.024
0.019
0.031
0.027
0.035
0.035
0.034
0.029
0.057
0.031
0.032
0.018
0.064
0.036

Otsu’s
0.026
0.122
0.018
0.027
0.035
0.029
0.031
0.018
0.035

b) NU
Image
1
2
3
4
5
6
7
8
3D

Groundtruth
0.014
0.013
0.018
0.023
0.032
0.032
0.039
0.045
0.031

IK_NU
0.088
0.014
0.036
0.026
0.040
0.037
0.049
0.050
0.039

38

Entropy Iterative
0.049
0.067
0.051
0.046
0.046
0.062
0.098
0.075
0.116

0.044
0.149
0.039
0.034
0.034
0.032
0.056
0.051
0.060

Otsu’s
0.043
0.157
0.037
0.033
0.033
0.032
0.053
0.050
0.056

The Otsu’s and iterative methods outperformed the entropy method in terms of both ME
and NU in all cases except image 2 (Table 2.4). The poor performance of iterative and Otsu’s
methods in image 2 can be related to the presence of a ‘dark’ region in the center of the image
which has GV below the calculated thresholds. Note that for most of the studied cases, the
iterative and Otsu’s methods generated very similar segmented results because the thresholds
calculated from the two methods were quite close (Table 2.3). Similar results were reported in
Gonzalez et al. (2004).
The IK segmentation yielded the segmented images with lowest ME as compared to the
other methods in all images except image 1 (Table 2.4). The main reason for the poor
performance of IK segmentation in image 1 is the lack of bimodal distribution in the histogram
of their GV (Figure 2.2(b)), which introduced substantial uncertainty in the ExpectationMaximization estimations of the mixed Gaussian distribution parameters. When histogram of
GV does not have a bimodal distribution, user specified values may work better for the selection
of T0 and T1 values (Oh and Lindquist, 1999). However, this introduces substantial trial and
error subjectivity in the threshold selection. Good performance of IK in all but one images of this
study are consistent with other reports that suggested great suitability of IK in studying the pore
network characteristics in various porous media (Oh and Lindquist, 1999; Peth et al., 2008).

Comparing performance of the segmentation methods using pore characteristics
Pore characteristics of the ground-truth images and of the segmented images obtained
using different segmentation methods are plotted on Figure 2.6. For porosity, IK outperformed
other methods in 7 of the 9 images (including 3D image) generating porosity values closest to
39

that of the ground-truth image (Figure 2.6(a)). For pore surface area, IK outperformed other
methods in 6 of the 9 images (Figure 2.6(b)). This result agrees well with the lowest ME values
obtained from IK segmentation (Table 2.4 (a)). However, results from ME did not match the
mean pore size very well (Figure 2.6(b)). IK produced the best mean pore size in images 5, 7,
and 8. Iterative method produced the best mean pore size in images 1 and 2. Entropy method
produced the best mean pore size in images 3, 4, and 3D. Otsu’s method produced the best mean
pore size only in image 6.

Figure 2.6: Summary of pore characteristics (a) porosity; (b) mean pore size; (c) pore surface
area obtained from images segmented using different segmentation methods.
(a)

0.3

Ground-truth
IK

0.25

Iterative

Porosity

0.2

Otsu's

0.15

Entropy

0.1
0.05
0
0

1

2

3

4

5

Image

40

6

7

8

9
3D

Figure 2.6 (cont’d).
(b)

Ground-truth
IK

Log(Surface area) (micron2)

6.55

Iterative

6.35

Otsu's

6.15

Entropy
5.95
5.75
5.55
5.35
0

1

2

3

4

5

6

Image

41

7

8

9
3D

Figure 2.6 (cont’d).
(c)

Ground-truth

165

Mean pore size (micron)

IK
140

Iterative
Otsu's

115

Entropy
90

65

40
0

1

2

3

4

5

6

7

8

9
3D

Image

We also examined the worst cases in performance of each method for the studied images
(Table 2.5). By the worst case we refer to the cases where the difference between the pore
characteristics from the segmented and the ground-truth images was the largest. Based on the
worst case assessment, IK and iterative methods were found to be preferable to other methods.

Table 2.5: Number of images (total 9 images) where the segmentation methods produced the
‘worst’ result in terms of pore characteristics.

porosity
mean pore size
pore surface area
Total

IK_NU Entropy Iterative
2
3
1
1
1
3
2
4
1
5
8
5

42

Otsu’s
3
4
2
9

Note that in most cases IK was among the methods that yielded the segmented images
with pore characteristics similar to their ground-truth counterparts, though it was not always the
best method.
To further summarize how well pore characteristics will be represented in the segmented
images when the segmentation method is selected based on NU values, we conducted relative
assessment of method performances (Figure 2.7). We calculated the relative errors as the
differences between the pore characteristic value from the segmented image and the pore
characteristic ground-truth value divided by the ground-truth value. The relative errors were
calculated for two scenarios: for the segmentation method that would be selected as the best
based on the NU criterion; and for the segmentation method that generated the pore characteristic
value the most similar to that of the ground-truth. The results showed that, NU-based selection
produced overall greater relative errors (5%-8%) than the minimal errors that were obtained
using the best method for each pore characteristic (1%-4%). Notice that high variability in the
relative errors is present in all pore characteristics for both scenarios (Figure 2.7). However, the
differences in relative errors between NU-based and the pore characteristics-based were
relatively small. This further indicates that even though inaccuracies in representation of pore
characteristics in segmented images is to be expected; in the absence of ground-truth information
the segmentation method selection using NU will produce acceptable results.

43

Relative Error, %

Figure 2.7: Comparing accuracy of NU-based segmentation method selection to pore
characteristics-based selection. The relative errors are defined as the differences between the
pore characteristic value from the segmented image and that from ground-truth value divided by
the ground-truth value.

Porosity

Mean pore size

Surface area

Pore Characteristics

CONCLUSION
Synthetic soil aggregate µCT images were simulated based on the 2D and 3D pore/solid
binary images by taking into account the partial volume effect, differences in attenuation values
of different solid materials, and noise that accompanies the scanning. Comparisons among
entropy, iterative, Otsu’s, and IK segmentation methods indicated that IK method was an optimal
segmentation tool when the histogram of studied images had a clear bimodal pattern. However,

44

for the images with poorly distinguishable histogram peaks the Otsu’s method was preferable. In
the absence of ground-truth information, choosing the value of IK parameter r using NU criteria
produced acceptable segmentation results.
The segmented images with minimal ME and NU values were not always optimal in
representing the pore characteristics, such as porosity, mean pore size, and pore surface area.
However, since the errors in pore characteristic representation when image segmentation choices
were made based on NU values were only slightly greater than the minimal errors, NU is
recommended as an acceptable criterion for selecting segmentation method.

45

BIBLIOGRAPHY

46

BIBLIOGRAPHY

Ananyeva, K., Wang, W., Smucker A.J.M., Kravchenko, A. Rivers, M., 2009. Comparison of
three approaches to measuring individual Soil Aggregate Density. Proceedings of the 19th
International Annual Meetings American Society of Agronomy, Crop Science Society of
American, and Soil Science Society of America, Pittsburgh, PA, p. 79.
Baveye, P.C., Laba, M., Otten, W., Grinev, D., Bouckaert, L., Dello Starpaio, P., Goswami,
R.R., Hu, Y., Liu, J., Mooney, S., Pajor, R., Sleutel, S., Tarquis, A., Wang, W., Wei, Q.,
Sezgin, M., 2010. Observer-dependent variability of the thresholding step in the
quantitative analysis of soil images and X-ray microtomography data. Geoderma 157, 5163.
Blake, G.R., Hartge, K.H.,1986. Bulk Density, in A. Klute, ed., Methods of Soil Analysis, Part I.
Physical and Mineralogical Methods: Agronomy Monograph no. 9 (2nd ed.), pp. 363-375.
Borsotti, M., Campadelli, P., Schettini, R., 1998. Quantitative evaluation of color image
segmentation results. Pattern Recog. Lett. 19, 741-747.
Capowiez, Y., Pierret, A., Daniel, O., Monestiez, P., Kretzschmar, A., 1998. 3D skeleton
reconstructions of natural earthworm burrow systems using CAT scan images of soil cores.
Biol. Fert. Soils 27, 51-59.
Chang, E.J., Yen, J.C., Chang, S., 1995. A new criterion for automatic multilevel thresholding.
IEEE Trans. Image Process. 4, 370-378.
Christakos, G., 1992. Random Field Models in Earth Sciences. New York: Academic Press.
Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from in- complete data
via the EM algorithm. J. Roy. Statistical Society, Series B, 39(1), 1-38.
Deutsch, C.V., Journel, A.G., 1998. Geostatistical Software Library and User’s Guide. Oxford
University Press, New York.
Elliot, T.R., Heck, R.J., 2007. A comparison of 2D vs. 3D thresholding of X-ray CT imagery.
Can. J. Soil Sci. 87, 405-412.
Gonzalez, R.C., Woods, R.E., Eddins, S.L., 2004. Digital image processing using MATLAB.
Prentice Hall, pp.405-407.
Iassonov, P., Gebrenegus, T., Tuller, M., 2009. Segmentation of X-ray computed tomography
images of porous materials: a crucial step for characterization and quantitative analysis of
pore structures. Water Resour. Res. 45, W09415.

47

Iassonov, P., Tuller, M. 2010. Application of image segmentation for correction of intensity bias
in X-Ray CT images. Vadose Zone J. 9, 1–5.
Jaynes, E.T. 1957. Information Theory and Statistical Mechanics.Physical Review. Volume 106,
#4.
Kaestner, A., Lehmann, E., Stampanoni, M., 2008. Imaging and image processing in porous
media research. Adv. Water Resour. 31 (9), 1174–1187.
Kapur, J. N., P. K. Sahoo and A. K. C. Wong. 1985. A new method for gray level picture
thresholding using the entropy of the histogram, Computing Vision Graphics Image
Process. 29, 273-285.
Kasteel, R., Vogel, H. J., Roth, K., 2000. From local hydraulic properties to effective transport in
soil. Eur. J. Soil Sci. 51, 81-91.
Ketcham, R.A., Carlson, W.D., 2001. Acquisition, optimization and interpretation of X-ray
computed tomographic imagery: Applications to the geosciences. Comput. Geosci. 27,
381-400.
Kittler, J. and J. Illingworth. 1985. On threshold selection using clustering criteria. IEEE Trans.
Syst. Man Cybern. SMC-15, 652-655.
Kittler, J. and J. Illingworth. 1986. Minimum error thresholding. Pattern Recogn. 19, 41-47.
Lehmann, P., Wyss, P., Flisch, A., Lehmann, E., Vontobel, P., Krafczyk, M., Kaestner, A.,
Beckmann, F., Gygi, A., Fluhler, H., 2006. Tomographical imaging and mathematical
description of porous media used for the prediction of fluid distribution. Vadose Zone J. 5,
80-97.
Levine, M.D., Nazif, A.M., 1985. Dynamic measurement of computer generated image
segmentations. IEEE Trans. Pattern Anal. Mach. Intell. 7, 155-164.
Lindquist, W.B., 2002. Quantitative analysis of three dimensional X-ray tomographic
images. Developments in X-ray tomography III, SPIE 4503, pp. 103-115.
http://www.ams.sunysb.edu/~lindquis/3dma/3dma_rock/3dma_rock_05_primer/primer.htm
l [Accessed December, 08, 2010].
Liu, J., Yang, Y. H., 1994. Multiresolution color image segmentation. IEEE Trans. on Pattern
Anal. Mach. Intell. 16 (7), 689-700.
Mooney, S.J., Tams, A.R., Berry, P., 2006. A reliable method for preserving soil structure in the
field for subsequent morphological examinations. Geoderma 133, 338-344.

48

Mooney, S.J., Morris, C., Craigon, J., Berry, P., 2007. Quantification of soil structural changes
induced by cereal anchorage failure: Image analysis of thin sections. J. Plant Nutr. Soil Sc.
170 (3), 363-372.
Oh, W., Lindquist, B., 1999. Image thresholding by indicator kriging. IEEE Trans. Pattern Anal.
Mach. Intell. 21, 590-602.
Otsu, N., 1979. A threshold selection method from gray-level histograms. IEEE T. Syst. Man.
Cyb. 9, 62-66.
Pal, N.R. and S.K. Pal., 1993. A review of image segmentation techniques. Pattern Recogn. 29,
1277-1294.
Park, E.J., Smucker, A.J.M., 2005. Saturated hydraulic conductivity and porosity within
macroaggregates modified by tillage. Soil Sci. Soc. Am. J. 69, 38-45.
Papadopoulos, A., Whitmore, A., White, A., Mooney, S.J., Bird, N., 2008. Combining spatial
resolutions in the multiscale analysis of soil pore-size distributions.Vadose Zone J. 8 (1),
227-232.
Papadopoulos, A., Bird, N., Mooney, S.J., Whitmore, A., 2008. Fractal analysis of pore
roughness in images of soil using the slit island method. Vadose Zone J. 7 (2), 456-460.
Peth, S., Horn, R., Beckman, F., Donath, T., Fisher, J., Smucker, A.J.M., 2008. ThreeDimensional quantification of intra-aggregate pore space features using synchrotronradiation-based microtomography. Soil Sci. Soc. Am. J. 72, 897-907.
Peth, S., Nellesen, J., Fischer, G., Horn, R., 2010. Non-invasive 3D analysis of local soil
deformation under mechanical and hydraulic stresses by µCT and digital image correlation.
Soil Till. Res. 111 (1), 3-18.
Piñuela, J., Alvarez, A., Andina, D., Heck, R.J., Tarquis, A.M., 2010. Quantifying a soil pore
distribution from 3D images: Multifractal spectrum through wavelet approach. Geoderma
155, 203-210.
Porter, M., Wildenschild, D., Grant G., Gerhard, J., 2010. Measurement and prediction of the
relationship between capillary pressure, saturation and interfacial area in a NAPL-WaterGlass bead system. Water Resour. Res. 46, W08512.
R Development Core Team, 2006. R: A Language and Environment for Statistical Computing. R
Foundation for Statistical Computing, Vienna, Austria, http://www.R-project.org
[Accessed December, 08, 2010].
Rasband, W.S., ImageJ, U. S. National Institutes of Health, Bethesda, Maryland, USA,
http://rsb.info.nih.gov/ij/ [Accessed December, 08, 2010], 1997-2009.

49

Ridler, T.W., Calvard, S., 1978. Picture thresholding using an iterative selection method. IEEE
Trans. Syst. Man. Cyb. SMC 8, 630-632.
Rivers, M. L., 1998. Tutorial Introduction to X-ray Computed Microtomography Data
Processing, University of Chicago, Chicago, IL, USA.
http://www.mcs.anl.gov/research/projects/X-ray-cmt/rivers/tutorial.html [Accessed
December, 08, 2010].
Rivers, M.L., Sutton, S.R., Eng, P., 1999. Geoscience applications of x-ray computed
microtomography. Proceedings of SPIE, the International Society for Optical Engineering,
Developments in X-ray Tomography II 3772, pp. 78-86.
Sahoo, P., Soltani, S., Wong, A., Chen, Y., 1988.A survey of thresholding techniques. Comput.
Vision Graph. 41, 233-260.
Sahoo, P., Wilkins, C., Yeager, J., 1997. Threshold selection using Renyi’s entropy. Pattern
Recogn. 1, 71-84.
SAS Institute, 2001. SAS User's Guide. Version 9.1. Statistical Analysis System Institute, Cary,
NC, USA.
Sezgin, M., Sankur, B., 2004. Survey over image thresholding techniques and quantitative
performance evaluation. J. Electron. Imaging 13(1), 146-165.
Schapp, M.G., Tuller, M., 2010. Quantitative Pore-Scale Investigations of Multiphase
Bio/Geo/Chemical Processes. Vadose Zone J., 9, 573–575.
Schlüter, S.,Weller,U.,Vogel,H.-J.,2010. Thresholding of X-raymicrotomography images of soil
using gradient masks. Comput. Geosci. 36, 1246-1251.
Shannon, C.E. 1948. A Mathematical Theory of Communication, Bell System Technical Journal,
27, 379-423 & 623-656.
Tsallis, C., 2001. In: Abe, S. and Okamoto, Y., Editors. 2001. Nonextensive Statistical
Mechanics and its ApplicationsSeries Lecture Notes in Physics, Springer, Berlin.
Udawatta, R.P., Gantzer, C. J., Anderson, S. H., Garrett, H. E., 2008. Agroforestry and grass
buffer effects on pore characteristics measured by high-resolution X-ray computed
tomography. Soil Sci. Sco. Am. J. 72 (2), 295-304.
Unnikrishnan, R., Pantofaru, C., Hebert, M., 2005. A measure for objective evaluation of image
segmentation algorithms. Proceedings of the 2005 IEEE Computer Society Conference on
Computer Vision and Pattern Recognition, San Diego, CA, USA.

50

Vaz, C.M.P, De Maria, I.C., Lasso, P.O., Tuller, M., 2011. Evaluation of an Advanced Benchtop
Micro-Computed Tomography System for Quantifying Porosities and Pore-Size
Distributions of Two Brazilian Oxisols. Soil Sci. Soc. Am. J. 75(3), 832-841.
Vogel, H.J., Kretzschmar, A., 1996. Topological characterization of pore space in soil -- sample
preparation and digital image-processing. Geoderma 73, 23-38.
Vogel, H.-J., Weller,U., Schluter, S., 2010. Quantification of soil structure based on Minkowski
functions. Comput. Geosci. 36, 1236-1245.
Yasnoff, W.A., Mui, J.K., Bacus, J.W., 1977. Error measures for scene segmentation. Pattern
Recogn. 9, 217-231.
Zhang, Y. J., 1996. A survey on evaluation methods for image segmentation. Pattern Recogn. 29,
1335-1346.
Zhang, Y. J., 2001. A review of recent evaluation methods for image segmentation. International
symposium on signal processing and its applications, Kuala Lumpur, Malaysia, 13- 16
August, 2001. pp. 148-151.

51

Chapter 3: INTRA-AGGREGATE PORE
CHARACTERISTICS UNDER LONG-TERM
CONVENTIONAL TILLAGE AND NATIVE
SUCCESSION VEGETATION: X-RAY COMPUTED
MICROTOMOGRAPHY

ABSTRACT
X-ray computed microtomography (µCT) opened new opportunities for studying internal
pore structures of individual soil aggregates. However, effective use of µCT data of irregularly
shaped soil aggregates requires properly identified soil aggregate boundaries, where image
closing technique is needed. Soil intra-aggregate pores of different land management practices
control soil physical, chemical and biological processes differentially at micro scales. Soil
aggregates from a Hapludalf soil under long-term conventional tillage (LTER-CT) and native
succession vegetation (LTER-NS) in southwest Michigan, and from an Ustochrept soil under
native succession vegetation (CERN-NT) and bare soil (CERN-BS) in northeast China were used
in this study. LTER-CT aggregates had significantly greater image based porosity than those of
LTER-NS. LTER-NS aggregates had more large pores (> 90µm) than LTER-CT aggregates,
while more medium size pores (45-90µm) were found in LTER-CT aggregates. Greater
abundance of medium sized pores in LTER-CT aggregates could be the cause of their reported
lower stability and higher macro-aggregate turnover rate. The differences in pore size
distributions between LTER-CT and LTER-NS were more pronounced in the aggregate interiors,
as compared to the exterior layers. No difference in intra-aggregate pore characteristics was
found between CERN-NT and CERN-BS aggregates. In aggregates from both treatments of both

52

studied soils pores tended to be larger in the aggregate centers while medium size pores (4590µm) were more abundant in the aggregate exterior layers.

INTRODUCTION
Soil matrix can be regarded as a hierarchical system formed by aggregates of various
sizes (Dexter, 1988) and only relatively recently we started to appreciate the role that intraaggregate heterogeneity plays in numerous physical, chemical and biological processes
(Ellerbrock and Gerke, 2004; Van Gestel et al., 1996). Pore arrangement within soil aggregates
and their connectivity to the inter-aggregate pore space are recognized to be among the main
factors controlling soil flux processes on an aggregate-scale (Revil and Cathles, 1999; Schapp
and Tuller, 2011). However, the exact mechanisms by which pore structure affects soil
functioning are still poorly understood due to poor accessibility of soil aggregate interiors to
direct measurements and lack of efficient methodology for aggregate-scale modeling.
Long-term differences in land use and management generate notable changes in soil
physical and hydraulic properties, including changes in soil organic matter content, porosity and
hydraulic conductivity (Brye and Pirani, 2005). Surface soil disturbance through tillage alters
overall soil aggregation, stability, biological properties and distribution of carbon (C) in
aggregates (Beare et al., 1994; Lal et al., 1998, 2002; Six et al., 1999; Green et al., 2007).
Besides soil disturbance, presence of vegetation and vegetation type can also alter soil properties
through plant roots activities and their association with soil organisms (Wienhold and Tanaka,
2001; Udawatta et al., 2008).

53

Soil aggregates play an important role in C sequestration (Blanco-Canqui and Lal, 2004),
and thus examining intra-aggregates pore structures in aggregates that are in land uses known for
their contrasting C sequestration potentials will help us understand the mechanisms by which
aggregates contribute to C sequestration. There were abundant literature on soil aggregate
formation, but two contrasting models exist – hierarchy concept proposed by Tisdall and Oades
(1982) and later developed by others (Oades et al, 1984) who believed soil macro-aggregates
were built of micro-aggregates with organic binding agents; and soil aggregate layer concept
(Augustin et al., 1995) where soil aggregates were formed around plant debris with soil materials
attaching to the aggregate center. Soil aggregates often exhibit spatial gradients in intraaggregate porosities, which have been assessed in a number of studies (Sextone et al., 1985;
Horn et al., 1994; Chenu et al., 2001; Park and Smucker, 2005a). These studies have
demonstrated that tillage and vegetation could affect aggregates’ exteriors and interiors
differentially. A recent study by Urbanek et al. (2011) showed that soil aggregate formation was
a result of combined hierarchical model and layer developing. It should be pointed out that
despite an enormous amount of work in the aggregate studies, there is no generally accepted
understanding of how the aggregates are built inside and of what the intra-aggregate pore
structure characteristics are.
Recent advances in X-ray computed microtomography (µCT) imaging enabled us to
quantify the intra-aggregate structures with a resolution of one to several microns (Gibson et al.,
2006; Nunan et al., 2006; Dathe et al., 2005, 2006; Peth et al., 2008). This is a significant step in
characterizing heterogeneity of soil environments at micro-scales, as it allows obtaining more
precise information regarding the physical structure of the soil aggregate matrix in three
dimensions (3D). Multi-scale soil solid-pore structures of intact soil matrix have been
54

characterized by extracting fractal dimensions or multifractal spectra (e.g., Dathe et al., 2005,
2006; Gibson et al., 2006) and by morphological characterization of soil macropore networks
(Pierret et al., 1999, 2002; De Gryze et al., 2004; Mooney et al., 2006, 2007; Mooney and Morris,
2008; Peth et al., 2008; Udawatta et al., 2008). Nunan et al. (2006) applied geostatistical tools to
quantify pore geometry for soil aggregates 1-3 mm in size at a scanning resolution of 4.4µm. De
Gryze et al. (2004) used X-ray computed tomography to study the pore morphology changes in
the soil aggregates incubated with fresh wheat residue and observed that porosity (pores > 27µm
in size) increased from 3.5% in the control aggregates to 7% in the incubated ones. Peth et al.
(2008) applied 3D µCT imaging to soil aggregates from conventional tillage and grassland at
resolution of 3.2-5.4µm, and observed morphological differences in their pore characteristics.
Papadopoulos et al. (2009) reported differences in pore shapes in 5 mm aggregates from
conventional and organic management at 8 µm resolution using µCT.
One of the difficulties in studying 3D images of the intact aggregates and pore
characteristics of the aggregate layers lies in the need for accurate delineation of the aggregate
boundaries. Until now, most of the aggregate image studies have been conducted by selecting a
regularly shaped region of interest within the aggregate’s central part. This precluded full
assessment of the aggregate’s pore structure and the differences in flow and transport properties
of the exterior and interior layers and interactions between them. Nunan et al. (2006) proposed
an innovative approach to analysis of CT aggregate images where the aggregate boundary is
identified by a sequence of image processing steps, such as binary morphological closing and
filling procedures. This is a significant advancement in soil aggregate image analysis. However,
conducting image closing and filling involves specification of certain parameters, e.g., the
neighborhood size parameter that defines the resulting ‘smoothness’ of the aggregate boundary
55

and the volume of the aggregate. The choice of the neighborhood size parameter value used
during the image processing can greatly affect the final result. How to select the value of the
neighborhood size parameter appropriate for soil aggregate image analysis remains unclear. Of
particular concern is that the pore structure characterization of the aggregate exterior layers can
be substantially affected. Note that the exterior layers are often of the greatest interest and
importance for analyses of the aggregate functioning within the soil matrix.
The first objective of this study is to investigate the selection of the optimal values of the
image processing neighborhood size parameter for defining soil aggregate boundaries using 3D
soil aggregate µCT images. Once the aggregate boundaries were identified, we (i) assessed intraaggregate pore characteristics both in the whole aggregates and in the aggregate layers at
different distances from the aggregate boundaries and (ii) analyze effects of land use on these
characteristics in the aggregates from two different soil types, namely, Hapludalfs on glacial
outwash, MI, USA and Ustochrepts on loess, Shenyang, China. The contrasting land use
practices that we will compare are conventionally tilled agricultural management, bare soil
without plant residue inputs, and native succession vegetation.

MATERIALS AND METHODS
Soil sample collection
Soil aggregate samples were collected from two long-term experiments located in the
north central part of the U.S. and in the north east part of China.

56

The first experiment is the Long Term Ecological Research site (LTER), W.K. Kellogg
Biological Station, established in 1988 in southwest Michigan (85o 24` W, 42 o 24` N). Soils are
well-drained Typic Hapludalfs of the Kalamazoo (fine-loamy, mixed, mesic) and Oshtemo
(coarse-loamy, mixed, mesic) series, developed on glacial outwash. Details on site description,
experimental design and research protocols of the experiment are provided on LTER site
(Kellogg Biological Station, 2009). The two LTER treatments sampled for this study are
conventional tillage (chisel-plowed) corn-soybean-wheat rotation with conventional chemical
inputs (LTER-CT), and native succession vegetation treatment, abandoned from agricultural
management after spring plowing in 1989 (LTER-NS). Prior to setting up the experiment in 1988
the entire area belonged to a uniformly conventionally tilled and managed agricultural field that
has been in agricultural production for at least past 100 years. Soil samples from this site were
collected in spring of 2008 from two adjacent 1-ha plots, one in LTER-CT and the other in
LTER-NS. The distance between the sampling locations of these two plots was approximately 25
m. Maintaining a short distance between the samples of the two studied treatments allowed us to
reduce effects of the local soil spatial variability on soil aggregate characteristics. Three sites
were sampled at each of the treatments.
The second experiment is located at Shenyang Experimental Station of Ecology, Chinese
Ecosystem Research Network (CERN), northeast China (41 o 31’ N, 123o 22’ W). Soils are aquic
brown soils (silty loam, Haplic-Udic Cambisols in Chinese taxonomy; Ustochrepts in USDA
taxonomy), developed on loess. The two CERN treatments used in this study were native
succession vegetation (CERN-NS) and bare soil with no vegetation (CERN-BS). The treatments
were established in 2002 in the field that has been abandoned from agricultural use in 1994. The
native vegetation in CERN- NS is dominated by Artemisia lavandulaefolia, Cephalanoplos
57

segetum, Chenopodium serotinum, Metaplexis japonica and Conyza Canadensis (Hua et al.,
2008). Soil samples were collected in summer of 2007. The distances between CERN-NS and
CERN-BS sampling locations were approximately 20 m. Five sites were sampled at each
treatment.
At each sampling site in both LTER and CERN experiments, a soil block approximately
15x15 cm2 in size was extracted from 0-20 cm depth using a sharp flat spade. Samples were airdried and manually sieved for 30 s. Aggregate 4-6.3 mm in size were retained for this study. We
decided to use the aggregates of this size since differences in porosities among the layers of the
aggregates of this size have been reported previously for a silty clay loam Alfisol (Park and
Smucker, 2005a). Thus we hypothesized that differences between the layers could also be
observed in µCT images of the LTER and CERN aggregates of this size used in this study.
Soil clay, silt, and sand contents in the aggregates were determined using Bouyoucos
hydrometer method (Gee and Bauder, 1986) and total C and nitrogen (N) were measured by dry
combustion method using Costech ECS 4010 Carbon-Nitrogen analyzer.

Soil aggregate image acquisition
Aggregate image data were collected on the bending magnet beam line, station 13-BM-D
of the GeoSoilEnvironCARS (GSECARS) at the Advanced Photon Source (APS), Argonne
National Laboratory (ANL), Illinois, USA. Data were collected with the Si (111) channel cut
monochromator tuned to 28 keV incident energy. Two dimension (2D) image slices were taken
at 1◦ rotation angle steps and combined into a 3D image consisting of 520 slices with 696x696
pixels per slice. The voxel size was 14.6 µm. Based on our preliminary results, this resolution
58

was found to be optimal for the 4-6.3 mm aggregate size fraction. A snapshot of the image on the
mirror was taken with the CCD camera representing a depth-integrated grayscale map of the
linear attenuation of the X-ray beam as it passed through the sample. After an image was
collected, the column was rotated 1º, and the process was repeated. Detailed information on the
procedure is given in Rivers et al. (1999).
The data were preprocessed by removing ring artifacts (Ketcham and Carlson, 2001) and
adjusting for backfield projections and reconstructed using Fast Fourier Transformation. After
reconstruction, the grayscale values in the image ranged from 0 (black) to 255 (white)
corresponding to low and high X-ray attenuations, respectively.
Defining aggregate boundaries
In order to assess the pore structure using images of the whole aggregates and aggregate
layers, we need to define the aggregate boundary. For that we used flood/seed fill algorithm
(Heckbert, 1990). The algorithm begins with identifying the starting point, usually the upper left
corner of the image, and then the region of the image connected to the starting point is
determined by checking whether they have the same ‘color’ (usually black), and filled with
another ‘color’ (usually white) serving as the background. One of the difficulties with the
automatic application of this process for detection of the soil aggregate boundary is that the
process may result in misclassification of the pores connected to the aggregate’s surface as the
outside space. Misclassification of these pores may significantly affect the subsequent pore
network analysis. Figure 3.1 presents an illustration of the aggregate image (Figure 3.1(a)) and
its incorrectly delineated boundary that classifies the pore connected to the aggregate’s surface as
an outside background (Figure 3.1(b)).

59

Figure 3.1: Example of (a) 2D unprocessed soil aggregate grayscale image slice; (b) the same
image with the background in black and the aggregate in white where aggregate boundary has
been identified using an automated flood/seed fill algorithm; note that the intra-aggregate pore
connected to the outside is incorrectly classified as belonging to the background (red square); (c)
the same image after boundary sealing algorithm has been applied with r = 12; (d) the final
processed gray-scale image with the background in white; (e) the final processed grey-scale
image with the background in white with r = 15; note the artificial pore pixels added to the sides
of the image (red squares).

(a)

(b)

1 mm

(c)

(d)

1 mm

(e)

1 mm

60

In order to ensure proper inclusion of such pores in the aggregate body we utilized an
image closing technique (Gonzalez et al., 2004), that allows sealing pore “necks” located at the
aggregate boundaries (Figure 3.1(c)). Image closing is a two-step procedure which involves
image dilation followed by image erosion (Gonzalez et al., 2004). An illustration of the image
closing procedure applied to a simple 2D image using a circular structuring element is shown in
Figure 3.2. Details of the image closing technique can be found at Image Processing Learning
Resources website (HIPR2, 2003). All image processing and data analyses were conducted using
3D aggregate images.
Outcomes of image closing are sensitive to the parameter called neighborhood size, r,
which defines the size of the structuring element used for image closing operation (MATLAB,
Version 7.8.0). The structuring element is a cluster consisting of a certain combination of
pixels/voxels that is used to probe the surface of the unevenly-shaped image (Gonzalez et al.,
2004). Detailed illustration of image processing using structuring elements with different sizes is
provided by Vogel et al. (2010). In this study we used 3D spherically shaped structuring
elements, from which the neighborhood size parameter r is defined by its radius. An example
shows a spherical structuring element with r equal to 2 (Figure 3.3).

61

Figure 3.2: An example of 2D image closing using a disk-shaped structuring element: a) the
original image; b) dilating the image a) (the portion added in a course of dilation is shown in
light blue); c) eroding the image b) (the pixels that were not present in the original image are
shown in red); d) the final image.
(a)

(b)
(b)

(c)

(d)

62

Figure 3.3: An example of a spherical structuring element with r = 2. a) Voxels of the
structuring element are showing in red; b) central portion of the spherical structuring element (a)
is showing in red.
(a)

Origin

Origin

(b)

63

If the selected parameter r is too small, it will lead to an insufficient closure of the pore
necks. However, large values of r will result in over-smoothing of the aggregate boundary and
inflation of the aggregate volume. In order to identify the optimal r value we examined values
ranging from 0 to 20. The r of 0 corresponds to the original image where no sealing has been
done.
The main issue is what should be used as a criterion for identifying the optimal r value. In
this study we explored the usefulness and performance of two different potential criteria for r
value selection. The criteria are the change in the overall aggregate volume as a function of r,
and the change in the fractal dimension of the aggregate boundary as a function of r values. For
analysis of the performance of these criteria in assessing the optimal value of r parameter, we
used 3D images from four aggregates, one from each treatment in each experiment.
To determine the aggregate volume we used the total number of voxels classified as those
belonging to the aggregate. When r = 0, the volume is equal to that of the original image. As r
increases, so does the size of the structuring element and, hence, the aggregate’s volume.
To determine the fractal dimension of the aggregate surface, the boundary was first
identified using Canny edge detector (Canny, 1986) implemented in MATLAB. The fractal
dimension of the boundary was then calculated using 3D box-counting algorithm (Mandelbrot,
1982). The aggregate boundary is the most irregular and has the highest fractal dimension when r
= 0. As r increases, the boundary becomes smoother and fractal dimension decreases.
Once we have identified an optimal r value, it was applied to all aggregate images used in
this study. Image closing with pre-selected r value has been conducted using Imclose in
MATLAB Image Processing Toolbox (Gonzalez et al., 2004). The “closed” pores (Figure 3.1(c))
64

were then converted into parts of aggregate bodies (Figure 3.1(d)) using Imfill in MATLAB
Image Processing Toolbox (Gonzalez et al., 2004).

Image segmentaiton
After boundary delineation, the stacks of grayscale images were segmented into pores
and solids using indicator kriging segmentation module in 3DMA-Rock software (Lindquist et
al., 2000). Multiple sets of the two thresholds in indicator kriging segmentation were chosen by
the parameters of the fitted mixed Gaussian distribution of GV within the aggregate domain
using Expectation-Maximization algorithm (Wang et al., 2010). The segmented images were
assessed using region non-uniformity measure proposed by Levine and Nazif (1985) and
evaluated by Wang et al. (2011). Detailed description of region non-uniformity calculation and
applicability in soil image analyses is provided in Wang et al. (2011). All the calculations were
conducted using 3D data sets.

Pore characteristics
The segmented soil aggregate images were used to extract pore characteristics, including
macro-porosity, medial axes, and burn numbers (BN) (Lindquist et al., 2000; Peth et al., 2008).
Macro-porosity is defined as the percent of the pore voxels within the aggregate image, and in
this study it is the porosity visible at scanning resolution of 14.6 µm, i.e. pores >14.6 µm. Medial
axis is a digitized curve which represents the skeleton of the pore space. Burn number is the
distance from the medial axis pore voxel to the closest solid voxel. A schematic representation of
medial axes and burn numbers is demonstrated in Figure 3.4. The units for BN are the numbers
65

of voxels. Note that BN values are directly related to pore sizes, e.g., in this study, BN equal to 1
is equivalent to pore size of 15-30 µm, BN equal to 2 is equivalent to pore size of 45-60 µm, BN
equal to 3 is equivalent to pore size of 75-90 µm, etc. In this study the BN distribution was
standardized by dividing the number of BN voxels of each size by the total number of pore
voxels and we will refer to it as standardized fraction for BN (std. fraction). Calculations of the
pore characteristics were conducted in 3DMA-Rock software.
Figure 3.4: A schematic representation of medial axes and burn numbers. Black pixels represent
solid materials and white pixels represent pore space. Rainbow color scale pixels are medial axis
paths with red for burn number equal to 2, yellow for burn number equal to 3, green for burn
number equal to 4, etc.

The aggregate total porosity was determined from aggregate weights and volumes (Rossi
et al., 2008). The aggregate volume was obtained using 3D aggregate images. Particle density for
total porosity determination was measured using pycnometer method following procedure
described in Bielders et al. (1990) using ~10 g samples consisting of several 4-6.3 mm
aggregates. The particle density was not different among the treatments and was equal to 2.65
g/cm3.

66

Identification of aggregate layers
Since one of the hypotheses of this study is that aggregate pore networks may differ in
their characteristics depending on the distance from the aggregate boundaries, we conducted pore
characterizations in three aggregate layers. This layer classification was consistent with and thus
allowed direct comparisons with the study by Park and Smucker (2005b). The layers were
identified based on the shortest distance between each voxel inside the aggregate and the
aggregate’s surface (Figure 3.5). The determination of the layers was conducted in 3DMA-Rock
software using the same technique as the one applied for the BN determination. For that, first the
maximum BN of each aggregate was determined. This number represented the distance between
aggregates’ surfaces and their centers. Three layers, namely, exterior, transitional and interior,
were determined as each layer constituted about 1/3 of the aggregate’s volume. The width for
each layer was calculated based on three concentric spheres:

&
& 2
2#
1#
1
W2 = $ 3 − 3 ! R , W3 = $1 − 3 ! R
W1 =
R,
$
$ 3
3!
3!
3
%
"
%
"
3

Where

(3.1)

W1 , W2 , and W3 are the radius of the interior, transitional and exterior layer of the

aggregates, R is the distance between aggregate’s surface and its center.

67

Figure 3.5: An example of the three layers of aggregates in a 2D image slice. Light, medium,
and dark gray intensities represent, respectively, exterior, transitional, and interior layers. Notice
that the layers were identified in 3D, therefore they might not look like they are of equal width
when shown in a 2D representation.
Exterior

Transitional

Interior

Statistical analyses
Statistical analyses for comparing of pore characteristics of the whole aggregates from
different experiments and treatments were conducted as a completely randomized design with
treatments in LTER or CERN experiments, respectively. For analyses of pore characteristics in
exterior, transitional and interior layers of the aggregates, layer was used as an additional factor
in a split-plot arrangement. The data analyses were conducted using PROC MIXED in SAS 9.2
(SAS Institute, 2008).

68

RESULTS AND DISCUSSION
Aggregate boundary identification
The relationship between the aggregate volume and the values of the parameter r is
presented on Figure 3.6(a). The aggregate volumes and r values were positively related.
However, while the magnitude of the relationship appeared to be inconsistent at 0-10 range of r
values, the aggregate volumes were slowly approaching plateaus at r values greater than 10. This
indicated that structuring element with small r values (r <= 10) was not sufficient to ‘close’ all
the necessary pores connected to aggregate boundaries, e.g., cracks and channels induced by
pneumatic pressures and/or plant roots and fungi hyphae. On the other hand, when the size of the
structuring element became too large (r >= 15), we observed artificial layers of pores added to
the aggregate boundary (Figure 3.1e). This result indicates that substantial caution must be
exercised when selecting optimal r value since large r values significantly increases the
porosities which are falsely taken from the background. Figure 3.6(c) showed that the aggregates’
volume increase rates were generally declining until r value of 12, after which the increase rates
went up a bit. We need to point out that there exists a large aggregate-to-aggregate variability in
the relationship between aggregate volume and r value, so the best strategy is to determine the
optimal r value based on each individual aggregate in the study. However, considering the large
number of samples in this study (49 aggregates) and significant computer processing time for
large r values (r > = 12), we decide to choose a universal r value for all the aggregates in this
study.
The complexity of the aggregate boundary as described by fractal dimension decreased
with increasing r values (Figure 3.6(b)). As expected, the aggregate boundaries became much
smoother when large structuring elements were probing the image. Notice that at r >= 15 the
69

aggregate boundary had become sufficiently smooth not to be strongly affected by further
increase in r values. The reduction in the fractal dimension for different r values mostly agreed
with the corresponding increase in aggregate volumes (Figure 3.6(b)(d)).
Figure 3.6: (a) Total volume of the aggregates after sealing the boundaries with different r
values in four studied aggregates. Solid lines were fitted to the linear parts of the volume/r
relationships; (b) fractal dimension of the aggregate boundaries of the four studied aggregates at
different r values; (c) aggregate volume increase rate† at different r values; (d) reduction rate† of
fractal dimension at different r values.

(a)

2.18
Agg1
Agg2
Agg3
Agg4

Fractal dimension

2.16
2.14
2.12
2.10
2.08
2.06
2.04
0

1

2

3

4

5

6

7

8

9

10

r

eduction in fractal dimension, %

(b)
Agg1
Agg2
Agg3
Agg4

1.2
1.0
0.8
0.6

70
0.4
0.2

2.06
2.04
0

1

2

3

4

5

Figure 3.6 (cont’d).

6

7

8

9

10

r

Reduction in fractal dimension, %

(b)
Agg1
Agg2
Agg3
Agg4

1.2
1.0
0.8
0.6
0.4
0.2
0.0
1

2

3

4

5

6

r

71

7

8

9

10

Figure 3.6 (cont’d).

(c)
Agg1
Agg2
Agg3
Agg4

Volume increase rate, %

8

6

4

2

0
1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19

r

Reduction in fractal dimension, %

(d)

6
Agg1
Agg2
Agg3
Agg4

5
4
3
2
1
0
1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19

r

† Increase rate and reduction rate are calculated as the difference between the value at
current r and the value at (r-1), and then divided by the value at (r-1).
72

`

Although there was a large variability among the volume and fractal dimension changes
as a function of r values in different aggregates, r value of 12 appeared to be suitable for
boundary sealing in all four aggregates. The r value of 12 ensured the inclusion of boundary
pores having opening less than about 350 µm (14.6x12x2 µm, where 14.6µm is the image
resolution and is multiplied by optimal image closing diameter value (2x12)) but at the same
time maintained sufficiently irregular aggregate boundaries. Once optimal r value (r = 12) was
determined based on the four aggregates, it was applied to all aggregate images of the study.
However, we need to point out that still further caution should be exercised during boundary
sealing. We observed that in some of the studied aggregates, sealing with r = 12 still resulted in
several voxels of an artificial layer of pores added to the aggregate boundaries. Removing 5
layers of pore voxels from the entire aggregate boundary remediated the problem. Please note
that there is a possible implication that doing this could result in fewer pores in the aggregates’
exterior layer, so we carefully checked several aggregates’ images to ensure that removal of 5
layers of pores from the aggregate boundary is only sufficient for removing the artificial layer.
Overall, we recommend examination of aggregate volumes as a main tool and fractal
dimension as an auxiliary tool for identifying appropriate r values. We expect that the optimal r
value will depend on the size of the studied aggregates, the scanned resolution and the average
size of pore connected to aggregate surface.

73

Pore characteristics in whole aggregates
Total porosity values for LTER-CT and LTER-NS aggregates were equal to 34.4% and
36.4%, respectively (Table 3.1). Although total porosity for LTER-NS aggregates was
numerically higher than that of LTER-CT aggregates, the difference was not statistically
significant. The tendency for lower total porosity in macro-aggregates from conventional tillage
as compared to those from no-till and forest soils was reported by Park and Smucker (2005a).
They attributed lower total porosities in the conventional tillage aggregates to more frequent
aggregate turnover rates, and accumulation of clay particles inside the aggregates. Munkholm
and Kay (2002) found that aggregates of 4-8 mm in size had lower porosity in soils under
continuous corn management than in soils from a dairy farm with 5-year crop grass rotation
management. Lower total porosities of the aggregates from conventionally tilled as compared to
native vegetation or pasture soils appear to be associated with lower organic matter content and
lower aggregate stability (Six et al., 1999; Lal, 2002). In this study LTER-CT has significantly
lower C concentration than LTER-NS (Table 1). Total aggregate porosity values for CERN-NS
and CERN-BS aggregates were 45.5% and 46.7%, respectively (Table 3.1), and also not
significantly different from each other.

74

Table 3.1: List of soil aggregates physical and chemical properties from top soil (0-20cm), including texture, total C, N, total
porosities, macro-porosities (pores > 14.6 µm), and porosities of pores < 14.6 µm.

Location Treatment Sand, % Silt, %
LTER
CERN

CT
NS
NS
BS

21.3
26.9
32.7
24.4

44.1
33.8
35.6
46.3

Clay, %
34.6
39.4
31.7
29.4

Total
C, %
0.78
1.02
1.31
1.19

Total
N, %
0.07
0.09
0.12
0.10

Total
porosity, %
34.4 (2.3)‡
36.4 (3.7)
45.5 (4.3)
46.7 (2.9)

Macroporosity, %
10.4 (3.2) a†
5.9 (1.7) b
14.2 (5.4) a
15.3 (5.8) a

Porosity < 14.6
µm, %
24.0 b
31.5 a
31.3 a
31.4 a

† Means within the same columns of the same location followed by the same lowercase letter are not significantly different (p<0.05).
‡ Standard deviations are shown in parentheses.

75

Even though there were no differences in total porosity between LTER-CT and LTERNS aggregates, pore size distributions of the two treatments were noticeably different. First, the
macro-porosity, that is µCT measured porosity based on visible pores >14.6 µm, was found to be
significantly higher in the aggregates from LTER-CT than in the aggregates from LTER-NS
soils (Table 3.1). No information on intra-aggregate porosity directly comparable to the results of
this study is available from literature. One example is the study by Peth et al. (2008) who
investigated the microtomographic images of two aggregates with resolution of 5.4 µm and 3.2
µm, representing a conventionally tilled soil and a long-term (>40 years) grassland Alfisol on a
loess parent material. They observed that the µCT measured porosity (pores~>5µm) was lower in
the aggregate from the conventionally tilled soil. Papadopoulos et al. (2009) compared soil
aggregate porosities using µCT with resolution of 8 µm for a silty clay loam soil from organic
and conventional management. They reported higher porosity values of aggregates in organic
management in one of the two studied fields. However, both studies only selected limited regions
of interest within the aggregates and only two aggregates (one from each land use) were used.
Note that in this study, the macro-porosity is based on the proportion of pores >14.6µm thus
direct comparisons with the results of Peth et al. (2008) and Papadopoulos et al. (2009) are
difficult.
Together with numerically greater total porosity values of LTER-NS aggregates, their
lower macro-porosity indicates that the small pores (<14.6µm) were more abundant in LTER-NS
than in LTER-CT aggregates (Table 3.1). The proportion of pores 15-30 µm in size (pores with
BN=1) were found to be similar in LTER-NS and LTER-CT soil (Figure 3.7(a)). However, pores
of intermediate size, 45-60 µm (pores with BN=2) and 75-90 µm (pores with BN=3), were more
abundant in LTER-CT than in LTER-NS soil (Figure 3.7(a)). The largest pores (>90 µm,
76

BN>=4) were much more abundant in LTER-NS than in LTER-CT aggregates and proportion of
pores 195-210 µm (BN=7) and >225 µm (BN>=8) were 3-4 times greater in LTER-NS than in
LTER-CT (Figure 3.7(a)). These results indicated that aggregates in long-established (>20 years)
native succession vegetation had pore size distributions different from those of the aggregates in
conventionally tilled agricultural management. Pore size distributions of the aggregates from
long term native vegetation were characterized by greater proportions of very small and very
large pores, but by lower proportions of intermediate pores.
Visual inspection of the aggregate images indicated that most of the large pores present in
LTER-NS aggregates are round pores of likely biological origin, while elongated pores and
cracks of non-biological origin were prevalent in LTER-CT aggregates. Six representative 2D
aggregate images from LTER-CT and LTER-NS aggregates are shown on Figure 3.8(a-f).

77

Figure 3.7: Percent of pore space visible at 14.6µm resolution occupied by pores of different burn numbers in aggregates from (a)
LTER-CT and LTER-NS treatments; and (b) CERN-BS and CERN-NS treatments. Burn number equal to 1 is equivalent to pore size
of 15-30µm, BN equal to 2 is equivalent to pore size of 45-60µm, BN equal to 3 is equivalent to pore size of 75-90µm, etc. For each
burn number, the letters mark the significant differences between land uses while ‘ns’ indicates no statistically significant difference
(p<0.05). Std. fraction on the y-axis is the standardized burn number distribution. Note that for clarity purposes different scales of yaxis values were used for fractions of pores with different burn numbers.

(a) 3e-2

3.5e-3
LTER-CT
LTER-NS

ns
3e-2

a

3.0e-3

Std. Fraction

2e-2
2e-2
a
1e-2

b

2.0e-3
1.5e-3

a
b

1.0e-3

a

5e-3

5.0e-4

0
1

Burn number

3

2.5e-4

78

2.0e-4
1.5e-4

a
b

4

Burn number
a

3.0e-4

b

0.0

2

3.5e-4

td. Fraction

Std. Fraction

b
2.5e-3

a

5

a

1.0e-2
1.0e-3

b

Figure 3.7 (cont’d).

b

1e-3

2

Burn number

b

0
1
3

(a) 3.5e-4

4

3e-4

2.5e-4
2.0e-4
1.5e-4

a
b

a

1.0e-4

2

Burn number
Burn number
(b) 3e-4

a

3.0e-4

ns

a

Std. Fraction

1

0.0
0.0

ns

ns

5.0e-3
5.0e-4

5e-3
0

1.5e-3

S

1e-2

Std. Fraction

St

a

3

5

Burn number

ns
ns

2e-4
2e-4

ns

1e-4

b

b

5e-5

5.0e-5
0.0

0
6

7

>=8

6

Burn number

7

Burn number

79

4

>=8

5

Figure 3.7 (cont’d).

(b)

3.5e-2
3.0e-2

4e-3
CERN-NS
CERN-BS

ns

3e-3

Std. Fraction

2.5e-2
2.0e-2
1.5e-2

ns

1.0e-2

2e-3
ns
1e-3

ns

5.0e-3
0.0

0
1

2

3

Burn number

3e-4

ns
ns

2e-4
2e-4

ns

1e-4
5e-5
0

4

Burn number

3e-4

Std. Fraction

Std. Fraction

ns

80

5

Figure 3.8: Examples of 2D image slices from central portions of the aggregates from LTER-CT
(a-c), LTER-NS (d-f), CERN-NS (g-h) and CERN-BS (i-j) soils.

(a)

(d)

1 mm

1 mm

(e)

(b)

1 mm

1 mm

(f)

(c)

1 mm

1 mm
81

Figure 3.8 (cont’d).

(i)

(g)

1 mm

1 mm

(h)

(j)

1 mm

1 mm

Marked differences in the aggregate size distributions in the LTER-CT and LTER-NS
treatments were reported by Grandy and Robertson (2006), who observed that aggregates from
LTER-CT had lower mean weight diameter, and LTER-CT soil had lower proportion of macroaggregates (2-8 mm) and higher proportion of micro-aggregates (<250 µm) as compared to
LTER-NS. De Gryze et al. (2004) reported lower stability of LTER-CT than LTER-NS
aggregates with aggregate mean weight diameter being significantly lower in soil of LTER-CT
as compared to LTER-NS. These lower stabilities of the LTER-CT aggregates could be caused
82

by the prevalence of the intermediate sized (45-90 µm) pores that we observed in this study.
Such pores might possibly serve as planes of aggregate breakage. Arguably, if a compressive
force is applied to LTER-CT and LTER-NS aggregates shown on Figure 3.8(a-f), the LTER-CT
aggregates will easily break into smaller fractions along the numerous micro-cracks, while
greater force would be necessary to break the LTER-NS aggregates. Additionally, CERES
modeling predicts surface soil aggregates from LTER-CT are exposed to 53 to 57 severe, eg, >10%
volumetric soil water contents and drying and wetting (DW) cycles annually in contract to 35 to
37 DW cycles for LTER-NS aggregates (Smucker and Park, 2006). These significantly greater
numbers and greater ranges of soil water content explain the larger number of planes of
weakness in the LTER-CT aggregates.
Greater abundance of large pores of biological origin in LTER-NS aggregates (Figs. 3.7(a)
and 3.8(d-f)) can be explained in part by greater numbers of plant roots and greater activity of
soil fauna in LTER-NS soil as compared to LTER-CT (Smith et al, 2008). Once the roots are
decomposed they leave behind pores with smoother curved walls enriched by presence of stable
1-10 µm organo-mineral microaggregate formations and root residues (Oades, 1993; Rasse and
Smucker, 1998; Six et al., 2004; Watteau et al., 2006). In this study, roots and pieces of
particulate organic matter were observed within pores of many LTER-NS aggregates. It also
should be pointed out that significant differences in soil textural percentages (Table 3.1) could
lead to significant difference in pore size distributions between LTER-CT and LTER-NS
aggregates as reported by Park and Smucker (2005ab). It could be speculated that higher clay
content in LTER-NS aggregates could be responsible for greater number of larger pores and less
number of smaller pores, since clay particles could possibly fill in the small space and create
smoother yet stable walls favoring for larger pore formation in LTER-NS.
83

No significant differences in pore size distributions were observed between CERN-NS
and CERN-BS aggregates (Table 3.1 and Figure 3.7(b)). Likewise, visual examination of the
aggregate images did not reveal any notable differences between the aggregates of the two
treatments (Figure 3.8(g-j)). This lack of differences probably reflected a relatively young age of
CERN experiment, insufficient for the differences to develop.

Pore characteristics in aggregate layers
At exterior, transitional and interior layers, macro-porosity of LTER-CT aggregates was
significantly higher than that of LTER-NS aggregates, consistently with the discussed above
macro-porosity results for the whole aggregates (Table 3.2(a)). However, in either LTER-CT or
LTER-NS aggregates there were no significant differences in macro-porosity among exterior,
transitional and interior layers. The possible explanation for lack of macro-porosity gradient in
LTER-CT is that aggregates were pieces of former compacted fragments, resulting from repeated
mechanical stress exerted by tillage, transport and harvest processes. In contrast, there was a
tendency for macro-porosity to increase from exterior to interior in both CERN-NS and CERNBS aggregates (Table 3.2(a)).
While no previous studies have examined the differences in macro-porosities of the
aggregate layers, the total porosity has been measured. E.g., Park and Smucker (2005b) reported
greater total porosity in exterior layer for 4-6.3 mm aggregates from silty clay loam soils of
Wooster (fine-loamy, mesic Typic Fragiudalf) and 6.3-9.5 mm aggregates from silt loam soils of
Hoytville (fine, illitic, mesic Mollic Epiaqualf) series under conventional tillage, and no tillage
management. Park and Smucker (2005b) attributed the more profound spatial gradient in total
84

porosity found in CT aggregates to compression of repeated tillage, which led to condensed
center of the aggregates, compared to less obvious porosity gradient in NT aggregates.
Table 3.2: Pore characteristics within exterior, transitional and interior layer of the aggregates.
(a) Macro-porosity; (b) Mean burn number on MA.
Location
LTER
CERN

Location
LTER
CERN

(a) Macro-porosity
Treatment Exterior
Transitional
CT
9.95 (0.73) a†A‡
10.76 (0.73)aA
NS
5.65 (0.76)bA
5.78 (0.76)bA
NS
12.45 (1.64)aA
13.65 (1.64)aA
BS
12.98 (1.71)aA
15.52 (1.71)aB
(b) Mean burn number on MA
Treatment Exterior
Transitional
CT
1.64 (0.05)aA
1.58 (0.05)aA
NS
1.61 (0.05)aA
1.63 (0.05)aA
NS
1.60 (0.05)aA
1.54 (0.05)aA
BS
1.63 (0.05)aA
1.60 (0.05)aA

Interior
10.45 (0.73)aA
6.27 (0.76)bA
15.38 (1.64)aB
16.56 (1.71)aB

Interior
1.59 (0.05)aA
1.73 (0.05)aB
1.62 (0.05)aA
1.65 (0.05)aA

† Means within the same columns of the same location followed by the same lower case letter
are not significantly different (p<0.05).
‡ Means within the same row followed by the same upper case letter are not significantly
different (p<0.05). Standard errors for means are shown in parentheses.

A strong spatial trend in mean burn number values was present in LTER-NS aggregates,
where the highest burn numbers were observed in the interior, followed by the transitional and
then by the exterior layer (Table 3.2(b)). This result provides further quantitative support to the
visual observations of large pores being often observed in the central parts of the LTER-NS
aggregates (Figure 3.8), reflecting the most complete absence of soil disturbance and high
biological activity for the longest period of time in this treatment.
Pore size distributions in the aggregate layers of the two LTER treatments had different
spatial gradients (Figure 3.9(a)). In LTER-NS the exterior layers had more pores of intermediate
sizes of 45-60 µm (BN = 2) and 75-90 µm (BN = 3) while noticeably fewer large pores
85

(pores >120 µm (BN>=5)) than in the interior layers. In contrast, in LTER-CT the exterior layers
were somewhat more abundant in 75-90 µm (BN = 3) and 105-120 µm (BN = 4) pores. However
no differences in terms of either pores with BN<=2 or BN>=5 were observed among the LTERCT layers. These differences point to greater spatial uniformity in pore distributions in
aggregates from this heavily disturbed conventionally tilled management as compared to the
undisturbed native vegetation.
The trends for higher proportions of 45-90 µm pores (BN=2 and 3) and lower proportions
of pores 135-150 µm (BN = 5) in LTER-CT than in LTER-NS that were previously mentioned in
the whole aggregate data also consistently appeared in the data from all three individual layers.
However, the greater abundance of the large pores, that is pores >165 µm (BN>=6), in LTER-NS
than in LTER-CT was only present in the transitional and interior layers. The differences
between the two treatments in terms of large pore presences were not significant in the aggregate
exteriors.
The tendency for more 75-90 µm (BN = 3) pores in the aggregate exteriors than in
transitional and interior layers noted for LTER-CT and LTER-NS aggregates was also present in
both CERN-NS and CERN-BS aggregates. Also similar to the native vegetation LTER-NS
treatment, large pores of 165-180 µm (BN = 6) and 195-210 µm (BN = 7) in size were more
abundant in the interiors of the native vegetation CERN-NS aggregates. These similarities point
to possible ubiquity in aggregate formation processes acting in the studied soils despite
substantial differences in their parent materials and the duration of the native vegetation presence.
No differences in pore size distribution between CERN-BS and CERN-NS aggregates were
observed in the individual layers (Figure 3.9(b)).

86

Figure 3.9: Percent of pore space visible at >14.6 µm resolution occupied by pores of different
burn numbers in aggregates exterior (Ext), transitional (Tran), and interior (Int) layers of (a)
LTER-CT and LTER-NS treatments; and (b) CERN-BS and CERN-NS treatments. Burn number
equal to 1 is equivalent to pore size of 15-30 µm, BN equal to 2 is equivalent to pore size of 4560 µm, BN equal to 3 is equivalent to pore size of 75-90 µm, etc. For each burn number, the
lowercase letters mark significant differences among the three layers within each treatment. The
differences among the layers are not statistically significant (p<0.05) if lowercase letters are not
displayed. For each burn number, the uppercase letters mark significant differences between the
treatments within each layer, while no upper case letters are shown when the differences are not
statistically significant (p<0.05). The cases when no statistically significant difference existed
both among layers and between treatments (p<0.05) are marked with ‘ns’. Std. fraction on the yaxis is the standardized burn number distribution. Note that for clarity purposes different scales
of y-axis values were used for fractions of pores with different burn numbers.

(a)

0.030
CT-Ext
CT-Tran
CT-Int
NS-Ext
NS-Tran
NS-Int

ns

Std. Fraction

0.025
0.020
0.015
A A
0.010

A
aB

0.005
0.000
1

2

Burn number

87

bB bB

Figure3.9 (cont’d).

(a)

0.004

0.003

aA

CT-Ext
CT-Tran
CT-Int
NS-Ext
NS-Tran
NS-Int

bA aB
bA
bB
bB

0.002
a

Std. Fraction

0.001

A

b
bB

bAaA
bA
B BB

0.000
3

4

5

0.0004
aA aA
0.0003

aA

b

0.0002
B

aA
aA

B

0.0001
B

B

ab

b

aB
ab b
b

0.0000
6

7

Burn number

88

>=8

Figure3.9 (cont’d).

(b)

0.035
ns

NS-Ext
NS-Tran
NS-Int
BS-Ext
BS-Tran
BS-Int

Std. Fraction

0.030
0.025
0.020
0.015

ns

0.010
0.005
0.000
1

2

Burn number

89

Figure3.9 (cont’d).

(b)

0.004

a

a
bb

NS-Ext
NS-Tran
NS-Int
BS-Ext
BS-Tran
BS-Int

ab
b

0.003

0.002
a

b

ab

Std. Fraction

0.001

ns

0.000
3

4

5

0.00030
a

a

0.00025
0.00020

ab
b
ns

0.00015

a

a
0.00010

ab

a
b

b

0.00005
0.00000
6

7

>=8

Burn number

CONCLUSION
Automatic flood fill algorithm for identifying aggregate boundaries on 3D µCT images
often fails to include pores connected to the aggregate outsides. Image closing is an effective
technique to delineate soil aggregate boundaries from µCT images. However, successful use of

90

image closing technique requires determining an optimal value of the neighborhood size
parameter, r. We found that examining the plots of the total aggregate volume, aggregate
boundary fractal dimension, and the rates of change in total volume and fractal dimension versus
r could serve as useful criteria for selecting the optimal r value. However, close attention should
be paid to the artificial layers of pores that might be added to aggregates’ boundary. The optimal
r value might depend on a number of factors including the size of the studied aggregates, the
scanned resolution and the average sizes of pores connected to aggregate surface. In this study
the r value of 12 appeared to be adequate for boundary delineation.
Twenty years of contrasting management in the LTER experiment led to pronounced
differences in the intra-aggregate pore size distributions. We observed significantly greater
percent of pores > 15 µm in LTER-CT aggregates as compared to those of LTER-NS. LTER-NS
had more large pores (pore size of > 90µm) than LTER-CT aggregates, while more medium size
pores (45-90 µm) were found in LTER-CT aggregates. Visual examination of aggregate images
indicated that large pores of biological origin were often present in the interiors of the LTER-NS
aggregates while medium-sized cracks were commonly found in the LTER-CT aggregates.
Further studies describing and quantifying the role of intra-aggregate pore size distributions in a
number of soil processes, including, but not limited to C sequestration are necessary.
The differences in pore size distributions between LTER-CT and LTER-NS were more
pronounced in the aggregate interiors, as compared to the exterior layers. For aggregates from
LTER and CERN, pores were found to be larger in the aggregate interiors than in the exterior
layers; while intermediate pores (45-90 µm) were more abundant in the aggregate exterior layers.
This trend was present in all land use and management treatments of the study and was most

91

pronounced in LTER-NS. Future research is needed to confirm spatial pattern in pore size
distributions within soil aggregates of other sizes.

92

BIBLIOGRAPHY

93

BIBLIOGRAPHY
Augustin, S., Jansen, M., Priesack, E., Beese, F., 1995. Litter decomposition and matter transport
in beds of soil aggregates. In: Hartge, K.H., Stewart, B.A. (Eds.), Soil structure. Its
development and function. Advances in Soil Science. CRC Press, Inc., Boca Raton, pp.
237-256.
Beare, M. H., P. F. Hendrix, and D. C. Coleman. 1994. Water-stable aggregates and organic
matter fractions in conventional- and no-tillage soils. Soil Sci. Soc. Am. J. 58: 777-786.
Bielders, C.L., L.W. Backer, and B. Delvaux. 1990. Particle density of volcanic soils as
measured with a gas pycnometer. Soil Sci. Soc. Am. J. 54: 822-826.
Blanco-Canqui, H., and R. Lal. 2004. Mechanisms of carbon sequestration in soil aggregates.
Critical Reviews in Plant Sciences 23(6): 481-504.
Brye, K.R., and A.L. Pirani. 2005. Native soil quality and the effects of tillage in the grand
prairie region of eastern Arkansas. Am. Midl. Nat. 154: 28-41.
Canny, J. 1986. A computational approach to edge detection. IEEE Trans. Pattern Analysis and
Machine Intelligence 8: 679-714.
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.
Dathe, A., and M.Thullner. 2005. The relationship between fractal properties of solid matrix and
pore space in porous media. Geoderma 129: 279-290.
Dathe, A., A.M. Tarquis, and E.Perrier. 2006. Multifractal analysis of the pore- and solid-phases
in binary two-dimensional images of natural porous structures. Geoderma 134: 318-326.
De Gryze, S., J. Six, K. Paustian, S.J. Morris, E.A. Paul, and R. Merckx. 2004. Soil organic
carbon pool changes following land-use conversions. Glob. Change Biol. 10:1120-1132.
Dempster, A.P., N.M. Laird, and D.B. Rubin. 1977. Maximum likelihood from in- complete data
via the EM algorithm. J. Roy. Statistical Society, Series B, 39(1): 1-38.
Dexter, A. R. 1988. Advances in characterization of soil structure. Soil Till. Res. 11: 199-238.
Ellerbrock, R.H., and H.H. Gerke. 2004. Characterizing organic matter of soil aggregate coatings
and biopores by fourier transform infrared spectroscopy. European Journal of Soil Science
55: 219-228.

94

Gee, G.W., and J.W. Bauder. 1986. Particle-size analysis. In: Klute, A. (Ed.), Methods of Soil
Analysis. Part I. Agron. Monogr. 9, second ed. ASA and SSSA, Madison, WI, USA, pp.
383-411.
Gerke, H. H., and Köhne, J. M. 2002. Estimating Hydraulic Properties of Soil Aggregate Skins
from Sorptivity and Water Retention. Soil Sci. Soc. Am. J. 66:26 -36.
Gibson, J.R., Lin, H. and Bruns, M.A., 2006. A comparison of fractal analytical methods on 2and 3-dimensional computed tomographic scans of soil aggregates. Geoderma 134: 335348.
Gonzalez, R.C., R.E. Woods and S.L. Eddins. 2004. Digital image processing using MATLAB.
Prentice Hall, 405-407.
Grandy, A. S. and G. P. Robertson. 2006. Aggregation and Organic Matter Protection Following
Tillage of a Previously Uncultivated Soil. Soil Sci. Soc. Am. J. 70: 1398-1406.
Green, V.S., D.E. Stott , J.C. Cruz, and N. Curi. 2007. Tillage impacts on soil biological activity
and aggregation in a Brazilian Cerrado Oxisol. Soil Till. Res. 92: 114-121.
Heckbert, P. A Seed Fill Algorithm. Graphics Gems, Andrew Glassner, ed. Academic Press,
1990: 275-277, 721-722.
Hypermedia Image Processing Reference 2 (HIPR2). 2003.
http://homepages.inf.ed.ac.uk/rbf/HIPR2/close.htm (verified 09 October, 2011).
Horn, R., W. Stepniewski, T. Wlodarczyk, G. Walenzik, and E.F.M. Eckhardt. 1994.
Denitrification rate and microbial distribution within homogeneous soil aggregates. Int.
Agrophys. 8: 65–74.
Hua, J., Y. Jiang, and W. Liang. 2008. Effects of vegetation coverage on the spatial distribution
of soil nematode trophic groups. Front. Biol. China 3(1): 63-67.
Kellogg Biological Station. 2009. Agronomic protocol: Long-term ecological research (LTER)
in row-crop agriculture. Available at
http://houghton.kbs.msu.edu/data/agronomic_protocol/2009_lter_agronomic_protocol.pdf
(verified 09 October, 2011).
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.
Lal, R., J. Kimble, R. Follet, and C. Cole. 1998. The potential of US cropland to sequester
carbon and mitigate the greenhouse effect. Ann Arbor Press, Chelsea, MI.
Lal, R. 2002. Soil carbon dynamics in cropland and rangeland. Environ. Pollut. 116: 353-362.
95

Levine, M.D., and A.M. Nazif. 1985. Dynamic measurement of computer generated image
segmentations. IEEE Trans. Pattern Anal. Mach. Intell. 7: 155-164.
Lindquist, W.B., Venkatarangan, A., Dunsmuir, J. and Wong, T.F., 2000. Pore and throat size
distributions measured from synchrotron X-ray tomographic images of Fontainebleau
sandstones. J. Geophys. Res. 105B: 21508.
Mandelbrot, B.B. 1982. The fractal geometry of nature. San Francisco, CA: W.H. Freeman.
MATLAB, Version 7.8.0. The MathsWorks, Inc. 2009.
Mooney, S.J., A.R. Tams, and P.Berry. 2006. A reliable method for preserving soil structure in
the field for subsequent morphological examinations. Geoderma 133: 338-344.
Mooney, S.J., C. Morris, J. Craigon, and P. Berry. 2007. Quantification of soil structural changes
induced by cereal anchorage failure: Image analysis of thin sections. J. Plant Nutr. Soil Sc.
170 (3): 363-372.
Mooney, S. J. and C. Morris. 2008. A morphological approach to understanding preferential flow
using image analysis with dye tracers and X-ray Computed Tomography. Catena 73: 204211.
Munkholm, L. J. and B. D. Kay. 2002. Effect of water regime on aggregate tensile strength,
rupture energy, and friability. Soil Sci. Soc. Am. J. 66: 702-709.
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.
Oades, J.M. 1984. Soil organic matter and structural stability: mechanisms and implications for
management. Plant Soil 76: 319-337.
Oades, J.M. 1993. The role of biology in the formation, stabilization and degradation of soil
structure. Geoderma 56: 377-400.
Oh, W., and B. Lindquist. 1999. Image thresholding by indicator kriging. IEEE Trans. Pattern
Anal. Mach. Intell. 21: 590-602.
Papadopoulos, A., N.R.A. Bird, A.P. Whitmore and S.J. Mooney. 2009. Investigating the effects
of organic and conventional management on soil aggregate stability using X-ray computed
tomography. Eur. J. Soil Sci. 60: 360-368.
Park, E.J., and A.J.M. Smucker. 2005 a. Saturated hydraulic conductivity through
macroaggregates modified by tillage. Soil Sci. Soc. Am. J. 69:38-45.

96

Park, E.J., and A.J.M. Smucker. 2005 b. Erosive strengths of concentric regions within soil
macroaggregates. Soil Sci. Soc. Am. J. 69: 1912-1921.
Peth, S., R. Horn, F. Beckman, T. Donath, J. Fisher, and A.J.M. Smucker. 2008. ThreeDimensional quantification of intra-aggregate pore space features using synchrotronradiation-based microtomography. Soil Sci. Soc. Am. J. 72: 897-907.
Pierret, A., Y. Capowiez, C. J. Moran, and A. Kretzschmar. 1999. X-ray computed tomography
to quantify tree rooting spatial distributions. Geoderma 90: 307-326.
Pierret, A., Capowiez, Y., Belzunces, L., and Mora,n C.J. 2002. 3D reconstruction and
quantification of macropores using X-ray computed tomography and image analysis.
Geoderma 106, 247-271.
Rasse, D.P. and A.J.M. Smucker. 1998. Root recolonization of previous root channels in corn
and alfalfa rotations. Plant and Soil 204: 203-212.
Revil,A., and L.M. Cathles. 1999. Permeability of shaly sands. Water Resour. Res. 35: 651-662.
Rivers, M.L., S.R. Sutton, and P. Eng. 1999. Geoscience applications of x-ray computed
microtomography. Proceedings of SPIE, Developments in x-ray tomography II, 3772: 7886.
Rossi, A.M., D.R. Hirmas, R.C. Graham, and P.D. Sternberg. 2008. Bulk density determination
by automated three-dimensional laser scanning. Soil Sci. Soc. Am. J. 72:1591-1593.
SAS Institute, 2008. SAS user's guide. Version 9.2. SAS Inst., Cary, NC.
Schapp, M.G., and M. Tuller, 2010. Quantitative Pore-Scale Investigations of Multiphase
Bio/Geo/Chemical Processes. Vadose Zone J., 9: 573–575.
Sextone, A.J., N.P. Revsbech, T.B. Parkin, and J.M. Tiedje. 1985. Direct measurement of
oxygen profiles and denitrification rates in soil aggregates. Soil Sci. Soc. Am. J. 49:645651.
Six, J., E. T. Elliott, and K. Paustian.1999. Aggregate and soil organic matter dynamics under
conventional and no-tillage systems. Soil Sci. Soc. Am. J. 63: 1350-1358.
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 Till. Res. 79:7-31.
Smith, R.G., C. P. McSwiney, A. S. Grandy, P. Suwanwaree, R. M. Snider, and G. P. Robertson
2008. Diversity and abundance of earthworms across an agricultural land-use intensity
gradient. Soil Till. Res. 100: 83-88.

97

Tisdall, J.M., J.M. Oades. 1982. Organic matter and water-stable aggregates in soils. Journal of
Soil Science 33, 141-163.
Udawatta, R.P., S. H. Anderson, C. J. Gantzer, and H. E. Garrett. 2008. Influence of prairie
restoration on CT-measured soil pore characteristics. J. Environ. Qual. 37:219-228.
Urbanek, E., Alvin J.M. Smucker, and R Horn. 2011. Total and fresh organic carbon distribution
in aggregate size classes and single aggregate regions using natural 13C/12C tracer.
Geoderma 164: 164-171.
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.
Vogel, H.J., U.Weller,S.Schluter. 2010. Quantification of soil structure based on Minkowski
functions. Computers & Geosciences 36:1236-1245.
Wang, W., A. N. Kravchenko, A.J.M. Smucker, and M. Rivers. 2011. Comparison of image
segmentation methods in simulated 2D and 3D microtomographic images of soil
aggregates. Geoderma 162:231-241.
Watteau, F., G.Villemin, G. Burtin, and L. Jocteur-Monrozier. 2006. Root impact on the stability
and types of microaggregates in silty soil under maize. Europ. J. Soil Sci. 57:247-257.
Wienhold, B.J., and D.L. Tanaka. 2001. Soil property changes during conservation from
perennial vegetation to annual cropping. Soil Sci. Soc. Am. J. 65:1795-1803.
Young, I.M., Crawford, J.W., Nunan, N., Otten, W., and Spiers, A. 2009. Microbial distribution
in soils: physics and scaling. In Donald L Sparks (ed). Advances in Agronomy, Vol 100.
Burlington: Academic Press, 2008, pp.81-121.

98

Chapter 4: ESCHERICHIA COLI
REDISTRIBUTTION AND TRANSPORT IN
INTACT SOIL MACRO-AGGREGATES

ABSTRACT
Past research had demonstrated the ability of E. coli to survive and regrow in the soil
environment, and further to infiltrate with water flow under saturated conditions. Knowledge of
the exact pore locations within soil aggregates opens new opportunities for understanding soil
microorganism movement and functioning. The first objective of this study was to assess E.coli
spatial distribution potential release within soil aggregates from three studied soil treatments,
namely, conventionally tilled (T1) and no-till (T2) corn/soybean/wheat rotation and native
succession vegetation (T7) on Typic Hapludalfs of the Kalamazoo (fine-loamy) series at NSF
Long-Term Ecological Research site, southwest Michigan. The second objective was to study the
redistribution pattern of E.coli within soil aggregates and relate E. coli movement within
aggregates to the aggregates’ internal pore structures. The results confirmed that E.coli
movement in soil aggregates were mainly driven by water flow via capillary forces. Most of the
seeding E.coli cells were remained within the aggregates (97%-99%) under the saturated
condition even when water was evacuated, however, they were found to be differentially
redistributed within the soil aggregates from three studied treatments. Redistribution was most
pronounced in T1 aggregates, followed by T2, and was almost negligible in T7 aggregates.
Principal component analysis showed that T1 aggregates are quite different from T2/T7
aggregates in pore characteristics, which was possibly responsible for E.coil redistribution
pattern. E.coli transport from top to other parts of soil aggregates were related to pore
99

characteristics of aggregates middle part, however, due to high variability in pore characteristics
among soil aggregates, the relationship was weak. E.coli redistribution was somehow related to
pore tortuoisities and connectivities of the aggregates’ center portion.

INTRODUCTION
Escherichia coli (E.coli), especially E.coli O157:H7, has been regarded as an important
indicator of environmental quality and public health risk (CDC; Vogt and Dippold, 2005;
Hellweger and Masopust, 2008). Effective strategies need to be developed to control the passage
of E.coli to waterways; yet thorough understanding of its survival and transport mechanisms is
indispensible. In agricultural fields, transport of manure-borne E.coli into soils, surface and
ground water can cause severe in-field contamination thus threatening food safety (Jamieson et
al., 2002; Guber et al., 2005; Jay et al., 2007). The soil structure and nutrient distribution at
macro and micro scale is believed to be responsible for E.coli survival and transport within soil
matrix (Pachepsky et al., 2006), and precipitation facilitates infiltration of bacteria suspensions
through flow pathways into deeper ground (Jenkins et al., 2008).
Soil matrix is among the most complex media in that it controls many physical, chemical
and biological processes at various scales. Past research had shown soil processes should be
studied at the scale where they function with applicable modeling techniques (Vogel and Roth,
2002; Vogel et al., 2003; Lin, 2003). Among various scales, it was believed that pore
arrangement within soil aggregates and their connectivity to the inter-aggregate pore space are
recognized to be among the main factors controlling soil processes (Revil and Cathles, 1999).
When entering soil matrix, E.coli will transport with water flow and/or attach to soil particles.
100

There were ample research published on E.coli movement and attachment in soils (Hagedorn et
al., 1978; Lindqvist and Bengtsson, 1991; Johnson and Logan, 1996; Ling et al., 2002; Guber et
al., 2005a&b). Soil organic matter can hinder E.coli attaching to soil particles due to its negative
charges (Johnson and Logan (1996)), whereas high clay content in soils may promote bacteria
attachment (Huysman and Verstrate, 1993).
Though findings above were mostly at field and soil column scale, it is highly possible
that E.coli travel within aggregates and further release from the aggregates are strongly
associated with pore structures they passed en route. However, lacking of appropriate tools and
modeling techniques, few past studies addressed E.coli distribution and transport at micro scale
in soils.
Recent advance of X-ray computed microtomography (µCT) imaging enabled us to
quantify the intra-aggregate structures in 3 dimensions (3D) with a resolution of several microns
(Gibson et al., 2006; Nunan et al., 2006; Dathe et al., 2005, 2006; Peth et al., 2008). Nunan et al.
(2006) investigated semivariograms and pore shape parameters of microbial habitat for soil
aggregates 1-3 mm in size at a scanning resolution of 4.4µm. Our previous work (Kravchenko et
al. 2011) showed that soil aggregates under different land use managements displayed
contrasting pore distribution patterns, thus impacting on E.coli fate and transport differentially.
The objectives of this study were 1) to compare E.coli spatial distribution within intact
soil aggregates from the same soil type but under long-term (>18 years) differential land use and
management settings, namely, conventionally plowed row crop agriculture, conservationally
managed (no-till) row crop agriculture, and native vegetation on land with previous agricultural
use under air dry and saturated flow conditions; 2) to explore the potential of E.coli release from
101

the aggregates and relate E.coli spatial redistribution under saturated flow to internal soil
aggregate pore structures quantified from µCT analysis.

MATERIALS AND METHODS
Soil sample collection
Soil samples were collected from the Long Term Ecological Research site (LTER),
Kellogg Biological Station, established in 1988 in southwest Michigan (85o 24` W, 42 o 24` N).
Soils are well-drained Typic Hapludalfs of the Kalamazoo (fine-loamy, mixed, mesic) and
Oshtemo (coarse-loamy, mixed, mesic) series, developed on glacial outwash. For details on site
description, experimental design and research protocols of the LTER experiment see
http://lter.kbs.msu.edu (Kellogg Biological Station, 2009). The three LTER treatments that were
sampled are corn-soybean-wheat rotation either conventionally plowed (T1) or no-till (T2) and
native succession vegetation (T7) treatment, abandoned from agricultural management after
spring plowing in 1989. Prior to setting up the experiment in 1988 the entire area belonged to a
uniformly conventionally tilled and managed agricultural field that has been in agricultural
production for at least past 100 years. Soil samples from this site were collected in spring of
2009 from three adjacent 1-ha plots. The distance between the sampling locations of these plots
was approximately 25 m. Maintaining a short distance between the samples of the two studied
treatments allowed us to reduce effects of the local soil spatial variability on soil aggregate
characteristics. Three sites were sampled at each of the treatments.
At each sampling site, a soil block approximately 15x15 cm2 in size was extracted from
0-20 cm depth using a sharp flat spade. Samples were air-dried and manually sieved for 30 s and
102

the aggregate fraction, 4-6.3 mm, was retained for this study. 15 scanned aggregates (5 per
treatment) were used in the flow experiment. 9 scanned aggregates (3 per treatment) were used in
the slicing experiment.

Soil aggregate image acquisition
Aggregate image data were collected on the bending magnet beam line, station 13-BM-D
of the GeoSoilEnvironCARS (GSECARS) at the Advanced Photon Source (APS), Argonne
National Laboratory (ANL), IL. Data were collected with the Si (111) channel cut
monochromator tuned to 28 keV incident energy. 2D image slices were taken at 1◦ rotation angle
steps and combined into a 3D image consisting of 520 slices with 696x696 pixels per slice. The
voxel size was 15.51 µm. Based on our preliminary results, this resolution was found to be
optimal for the 4-6.3 mm aggregate size fraction. A snapshot of the image on the mirror was
taken with the CCD camera representing a depth-integrated grayscale map of the linear
attenuation of the X-ray beam as it passed through the sample. After an image was collected, the
column was rotated 1º, and the process was repeated. Detailed information on the procedure is
given in Rivers et al. (1999).
The data were preprocessed by removing ring artifacts (Ketcham and Carlson, 2001) and
adjusting for backfield projections and reconstructed using Fast Fourier Transformation. After
reconstruction, the grayscale values in the image ranged from 0 (black) to 255 (white)
corresponding to low and high X-ray attenuations, respectively.

103

Aggregate image analysis
Soil aggregate boundaries were identified using a combination of image closing and
filling techniques and specific parameterization was applied to each individual aggregate to
ensure the best volume inclusion and boundary roughness (Wang et al., in preparation). After
boundary delineation, the stacks of grayscale images were segmented into pores and solids using
indicator kriging segmentation module in 3DMA-Rock software (Lindquist et al., 2000).
Multiple sets of the two thresholds in indicator kriging segmentation were chosen based on the
parameters of the fitted mixed Gaussian distribution of grayscale values within the aggregate
domain using Expectation-Maximization algorithm (Wang et al., 2011). The segmented images
were assessed using region non-uniformity measure proposed by Levine and Nazif (1985) and
evaluated by Wang et al. (2011). Detailed description of region non-uniformity calculation and
applicability in soil image analyses is provided in Wang et al. (2011). All the calculations were
conducted using 3D data sets.
The segmented soil aggregate images were used to extract pore characteristics, including
macro-porosity, medial axes, and burn numbers (BN) (Lindquist et al., 2000; Peth et al., 2008).
Macro-porosity is defined as the percent of the pore voxels within the aggregate image, and in
this study it is the porosity visible at scanning resolution of 15 µm, i.e. pores >15 µm. Medial
axis is a digitized curve which represents the skeleton of the pore space. Burn number is the
distance from the medial axis pore voxel to the closest solid voxel. The units for BN are the
numbers of voxels. Please note that BN values are directly related to pore sizes, e.g., in this study,
BN equal to 1 is equivalent to pore size of 15-30 µm, BN equal to 2 is equivalent to pore size of
45-60 µm, BN equal to 3 is equivalent to pore size of 75-90 µm, etc. The BN distribution is
standardized as dividing the number of BN voxels for each size by the total number of pore
104

voxels. And we will refer it to standardized fraction for BN (std. fraction). Calculations of the
pore characteristics were conducted in 3DMA-Rock software (Lindquist et al., 2000).
Based on medial axes information, the pore network connectivity was calculated as the
maximum flow in a network N = (V, E) using Edmonds-Karp algorithm (Edmonds and Karp,
1972). Medial axis voxels were used as the vertex sets V and burn numbers on medial axes were
used as capacities on edge sets E. The super-source node and super-sink node were artificially
superimposed to the pore network to avoid difficulty in multi-source and multi-sink maximum
flow calculation, where infinity capacity was assigned on each edge (Ford and Fulkerson, 1962).
The super-source node was set so that it connected to all the medial axis voxels on the first slice
of the 3D volume. The super-sink node was set so that it connected to all the medial axis voxels
on the last slice of the 3D volume. The MATLAB code for generating capacity matrix is in
Appendix 4.1. The C++ code for maximum flow calculation implementing Edmonds-Karp
algorithm is in Appendix 4.2.
Due to the very large size of the aggregate images, all the image data and codes are run at
High Performance Computing Center (HPCC) at Michigan State University.

Alignment of the physical and virtual cuttings on aggregates
Before cutting, the orientation of an aggregate was matched to that of the scanned data.
During actually cutting notes were taken on the order of cuts and sizes of resulting pieces. Then
based on the notes, the cutting planes on aggregate images were identified so as to ensure
maximal correspondence with the physical cuts within each individual aggregate. Final
assessment of quality of the cuts was made by visually examining the aggregate images along

105

with cutting planes and by calculating density of the aggregate slices. Density of each slice was
calculated based on the slice’s weight and volume as determined based on the cut location. The
density estimates could serve only as an approximate assessment tool since aggregate slices were
moist when weighted. A schematic representation of aggregates cutting was shown in Figure 4.1.
Figure 4.1: Schematic representation of the aggregate cutting. Numbers 1, 2, 3, 4, 5, 6 and 7
correspond to top, bottom, left, right, front, back and center slices, respectively. The sequence of
the numbers represents the cuttings in order.

E.coli application experiments
Two E.coil application experiments, namely, E.coli flow experiment and E.coli slicing
experiment, were conducted.
In the E.coli flow experiment, 15 scanned aggregates, 5 from each treatment, were used.
This experiment was conducted in the glass beads matrix chamber system designed for the study
(Chun et al., 2011 in preparation). Pictures of the system were shown in Figure 4.2.

106

Figure 4.2: A picture shows the glass beads matrix chamber system – glass beads matrix
chamber for adding bacteria then flushing and vacuum extracting E.coli from soil aggregates.

3

Approximately 6 mm of glass beads (1mm in diameter, with density of 2.5 g/cc) was
measured using centrifuge tubes, and the weight (Wbeads) was recorded. A single layer of glass
beads were added into the chamber, and then aggregates were carefully transferred into it. Prior
to placing an aggregate in the chamber its image was examined and then, consistent with the
image, the orientation of the aggregate (e.g., top, right, left, etc.) was marked with colored
permanent markers. This ensured that after the flow experiment relative positions of the
aggregate slices corresponding to the image could be identified for further analyses of the
relationships between intra-aggregate pores and E.coli distributions. Once in the chamber, E. coli
7

was added to the top of the aggregate at a concentration of approximately 10 CFU/ml. The
volume of E. coli culture dilution to be added to each aggregate was calculated as 30% of the

107

total porosity of the aggregate determined from the aggregate’s image. Then remaining glass
beads were added to cover the aggregate.
o

The whole system was allowed to equilibrate in the fridge (4 C) for 0.5 hour. Then
Phosphate Buffered Water (PBW) (pH 7-7.2) was added into the system. The volume of PBW
was determined as the in-between glass beads space volume:

VPBW = (VGB _ in _ tube −
where

VGB _ in _ tube

WGB

ρ GB ) + 70% × V pore

(4.1)

3

is the volume of glass beads placed in the centrifuge tube (6 mm ),

WGB

V pore

is the

is the weight of the glass beads,

ρGB

is the density of glass beads (2.5 g/cc) and

pore volume based on total porosity calculation.
o

After PBW application, the system was further equilibrated in the fridge (4 C) for 2
hours. Then the effluent was extracted from the system using a vacuum and collected into
centrifuge tubes. During the extraction process, the same suction pressure (120 mbar) was
applied to all the aggregates for approximately 8 mins. This pressure ensured the maximum
extraction of the liquid from the system. About 69% of the added PBW was extracted. The flow
system was kept in the ice box during the effluent extraction process. Effluent was then stored in
o

the freezer at -20 C for qPCR analysis.

The flow system was dried using compressed air for approximately 50 mins, while
chambers were surrounded by ice to minimize E. coli growth. Aggregates were cut into 7 pieces
108

using medical scalpels, flame-ethanol sterilizing in between each cut. Positions of the cutting
planes were recorded as described above. Aggregate slices were transferred into soil DNA kits
for qPCR analysis. Glass beads were transferred into centrifuge tubes with 5 ml PBW in. All the
o

samples were stored in the freezer (-20 C) for qPCR analysis.

In the E.coli slicing experiment, a total of 9 not-scanned aggregates, 3 from each
treatment, were used. We added 20 µL of E. coli culture dilution at concentration of
7

approximately 10 CFU/ml on top of each aggregate. Based on the assessment of scanned
aggregates, 20 µL is the average volume of the 30% aggregate pore space. Aggregates were
o

allowed to equilibrate for 2 hours after E.coli application and kept in the fridge (4 C) inside
closed containers within beakers containing water to avoid evaporation. Aggregates were cut into
7 pieces and transferred into soil DNA kits for qPCR analysis.

qPCR analysis

DNA was extracted from the seeded soil aggregates using UltraClean™ Soil DNA kit
(MoBio Laboratories, Solana Beach, CA) according to the manufacturer’s protocol. The
o

extracted DNA was stored at -20 C until further analysis. These samples were analyzed by
qPCR method at Research Technology Support Facility (RTSF), Michigan State University
(MSU) using the assay targeting uidA gene of E. coli (Srinivasan et al. 2011). The primers and
probe sequences used for the analysis is listed in Table 4.1. The final 15 µL reaction mixture for
amplification of E. coli uidA gene consisted of 10µL of Roche Light Cycler Taqman Mastermix,
109

0.5 µM of forward and reverse primers each, 0.2 µM probe, MgCl2 (3.2 mM) and nuclease free
water. Cycle threshold (Ct) was measured during each amplification and target gene
concentration was analyzed automatically by absolute quantification method by the
LightCycler® Software 4.0.

Table 4.1: Primer and probe sequence used in this study.
Bacteria/gene
E. coli/uidA

Primer/probe
Primer/probe sequence
cconcentrations
0.5 µM
5’-CAATGGTGATGTCAGCGTT-3’
0.5 µM
5’-ACACTCTGTCCGGCTTTTG-3’
0.2 µM
6FAM-TTGCAACTGGACAAGGCACCAGC-BBQ

To assess E.coli distribution within aggregate slices the number of E.coli cells observed
in each piece, i, were expressed per unit of mass and then divided by the total amount of E.coli
cells in the aggregate expressed per unit of mass:

n

obsi
Re lative _ E.coli _ ratio =
Weighti

∑ obs
i =1

n

∑Weight
i =1

The resulting relative E.coli ratio was used in all further analyses.

110

i

i

(4.2)

Statistical analysis
Descriptive statistics of pore characteristics were summarized using PROC MEANS
(SAS Inc., V9.2). The pore characteristic data were analyzed using principal component analysis
with soil management as qualitative variables using PROC PRINQUAL (SAS Inc., V9.2).
Categorical principal component analysis maximizes relationship among variables having
numerical, ordinal or nominal measurements using alternating lease squares (Young et al., 1978).
The univariate relationship between E.coli spatial distribution and pore characteristics in
aggregates pieces was analyzed using linear regression model in PROC REG (SAS Inc., V9.2).
The multivariate relationship between E.coli spatial distribution and pore characteristics in
aggregates pieces was analyzed using partial least regression (PLSR) method in PROC PLS
(SAS Inc, V9.2). Transformation of the variables was applied when normality assumption was
violated.

RESULTS AND DISCUSSION
E.coli distribution and transport in the flow and slicing experiment
The breakthrough of E.coli found in effluents and glass beads were ranging from 1.5% to
3.2% for aggregates with large variability observed from T1, T2 and T7 treatments (Figure 4.3).
Most of the E.coli cells were retained in soil aggregate matrix even though negatively charged
E.coli cell surface will probably repel clay particles in soils. Therefore, there are possibly other
mechanisms that control E.coli transport in soil aggregates. Moreover, the retention of E.coli
within soil aggregates is not a function of soil management in this study.

111

Figure 4.3: Percent of E.coli in the effluents and glass beads in the flow experiment.
E.coli in effluents and glass beads
10

Percent of E.coli, %

8
6
4
2
0
T1

T2

T7

It shows that E.coli cells distributed spatially heterogeneous within soil aggregates under
different water conditions and soil managements (Figure 4.4). Here we use differences in relative
E.coil ratio between soil aggregates pieces to compare the E.coil distribution in the flow and
slicing experiments. E.coli distribution in slicing experiment showed high degree of spatial
heterogeneity (Figure 4.4(a)). This is consistent with another set of similar experiments, where
we have shown that at air-dry condition E.coli transport in soil aggregates was mainly driven by
water flow via capillary suctions and by specific configuration of the intra-aggregate pores
(Kravchenko et al., submitted). This mechanism had been confirmed in previous studies on the
same and larger scales, such as Artz et al. (2005) and Brennan et al. (2010), who reported that
bacteria travels fast through macro-pore in soil columns with preferential flow; and Guber et al.
(2009) who found that E.coli entering air-day soil aggregates along with water flow was mainly
driven by capillary forces.

112

After saturation, equilibration, and water extraction the E.coli redistributed within the
aggregate primarily in vertical direction and became more homogeneously spread (Figure 4.4(a)).
Redistribution was most pronounced in T1 aggregates, followed by T2, and was almost
negligible in T7 aggregates. Top slices of soil aggregates retained more E.coli cells than other
parts of the aggregates for all soil aggregate managements in both the flow and slicing
experiments (Figure 4.4(a)). There were significant greater differences between top and other
portion of aggregates in relative E.coli ratio found in the slicing experiment, while no difference
among aggregates’ center, sides and bottom was observed (Figure 4.4(a)). However, the notable
differences in E.coli between top and other portion of the aggregates were found differentially in
three soil managements. There were significant greater differences between top and
bottom/center of aggregates in relative E.coli ratio found in the slicing experiment in T1 and T2,
while no such difference was found in T7 (Figure 4.4(b,d)). There were significant greater
differences between top and sides of aggregates in relative E.coli ratio found in the slicing
experiment in T1, while no such difference was found in T2 and T7 (Figure 4.4(c)). This
suggested that saturated water flow during rainfall and runoff events, could lead to significant
E.coli redistribution within soils. Moreover, the magnitude of redistribution highly depends on
soil managements; specifically conventional tilled soils will facilitate E.coli transport while soils
under native succession are more likely to retard its transport. The reason could be presence of
organic matter in T7 aggregates which will reduce their hydrophilicity. These differences
observed in the two experiments among soil managements could be attributed to pore structure
differences, such as total porosity, macro porosity and pore size distributions.

113

Figure 4.4: Comparison of relative E.coli ratio in flow and slicing experiments.

Differences in relative E.coli ratio between aggregate pieces

(a)

a

Flow
Slicing

3.0
a
a

2.5
2.0
b

1.5

b
b

1.0
NS

NS

0.5
0.0
top-bottom

top-sides

top-center sides-center center-bottom

Difference in relative E.coli ratio between aggregates' top and bottom Difference in relative E.
6

5
Flow
Slicing

4

a
a
4
Relative E.coli ratio

a
a

5

a
a

3
NS
NS

b
b

2
b
b

1

3
2

b
b

1

0

0
T1

T2

T7

T1

Difference in relative E.coli ratio between aggregates' top and

114

5
4

tio

(b)

Relative E.coli ratio

Difference in relative E.coli ratio

3.5

Flow
Slicing

a
a
a
a

Figure 4.4 (cont’d).

een aggregates' top and bottom Difference in relative E.coli ratio between aggregates' top and sides

(c)

NS
NS
3

b
b

2

Relative E.coli ratio

Relative E.coli ratio

Relative E.coli ratio

5
fference in relative E.coli ratio between aggregates' top and bottom Difference in relative E.coli ratio between aggregates' top
Flow
Flow
6
5
a
a
Slicing
Slicing
4
Flow
Flow
a
a
a
a
Slicing
Slicin
5
4
NS
NS
3
4
a
a
NS
b
b

NS
NS

1

2

b
b

1

NS

3

NS
NS
2

b
b
NS
NS

1

0

T7

T1

T2

0

T7
0

T1
T2
rence in relative E.coli ratio between aggregates' T7 and center
top

5

(d)

Difference in relative E.coli ratio between aggregates' top and center
Flow
Slicing
5

a
a

4

a
a

Relative E.coli ratio

b
b
b
b

1

a
a
3

NS
NS

2

b
b

NS
NS
b
b

1

0
T1

Flow
Slicing

a
a

4

3

2

T1

T2

T7

0
T1

T2

115

T7

T2

T7

Relationship between E.coli distribution and pore characteristics
Table 4.2 summarized pore characteristics for soil aggregates from T1, T2 and T7. Total
porosity values for T1, T2 and T7 aggregates were similar, equal to 39.3%, 35.2% and 37.9%
(Table 4.2). The macro-porosity for T1, T2 and T7 aggregates were 10.2%, 6.3% and 5.4%
(Table 4.2). Statistical analysis showed that T1 aggregates has significant greater macro-porosity
than those of T2 and T7 aggregates, where no difference found between T2 and T7 aggregates in
macro-porosity values. Pore size distributions for T1, T2 and T7 showed a declining trend with
increasing pore sizes, yet no differences were found among treatments for each pore size.

Table 4.2: Descriptive statistics for pore characteristics for soil aggregates with size 4-6.3mm
from conventional tillage (T1), no-till (T2) and native succession (T7). Standard deviations are
shown in parentheses. Burn number equal to 1 (BN1) is equivalent to pore size of 15-30 µm, BN
equal to 2 (BN2) is equivalent to pore size of 45-60 µm, BN equal to 3 (BN3) is equivalent to
pore size of 75-90 µm, and BN equal to or greater than 4 (BN4g) is equivalent to pore size > 105
µm.

Treatment

Total
porosity, %

Macroporosity, %

CT

39.3 (3.4)

10.2 (2.5)

NT

35.2 (4.7)

6.3 (2.8)

NS

37.9 (4.1)

5.4 (1.6)

BN1

BN2

BN3

BN4g

0.025
0.011
0.003
0.001
(0.007) (0.002) (0.0003) (0.0006)
0.028
0.009
0.002
0.002
(0.009) (0.001) (0.0006) (0.0007)
0.032
0.008
0.002
0.002
(0.012) (0.002) (0.0002) (0.0006)

Univariate analysis using each pore characteristics individually was suffering small
sample size issues. Thus we explored multivariate relationship between soil aggregates and their
pore characteristics. The two-dimensional solution with eigenvalues greater than 1 accounted for
73% of the variance from categorical PCA using soil pore characteristics variables were plotted
in Figure 4.5(a). The scores for each soil aggregates (soil aggregate id) were also plotted in the
116

same figure. It is clear that T1 aggregates were well separated from T2 and T7 aggregates with
the first PC, which accounts for 41.27% of the total variability (Figure 4.5(a)). However, there is
no distinction between T2 and T7 aggregates.
On the first dimension, total porosity, macro-porosity, BN2, BN3 and CT had positive
component loadings (Figure 4.5(a)). Therefore, these features are more favorable to T1
aggregates. Large projection of macro-porosity on the first principal component indicated that it
contributed a significant portion in explaining the variability in T1 aggregates. BN3 and BN2 are
the second and third largest positively contributing variable for T1 aggregates, indicating
prevalence of medium size pores (45-90 µm) in conventionally tilled soil aggregates. Similar
findings were reported in Wang et al. (2011) who studied another set of 4-6.3 mm T1 and T7
aggregates sampled from the same KBS-LTER experiment, but at different time. These results
perhaps point to a universal nature of observed tendency for CT aggregates which contain a large
amount of cracks and elongated pores, leading to low aggregate stability reported in other studies
(DeGryze et al., 2004; Grandy and Robertson, 2006). However, BN1 and BN4g were found not
correlated to neither T1 nor T7 aggregates, which demonstrate small pore size (15-30 µm) and
large pore size (> 105 µm) are probably uniformly distributed in all soil aggregates in this study.
This is inconsistent with Wang et al. (2011) for pore size distributions found in T7 soil
aggregates where more large pores were present. However, with only 5 aggregates from T7 from
this study, it is not clear how this results could be comparable to 12 T7 aggregates used in that
study. The differences in E.coli redistribution observed between T1/T2 and T7 (on Figure 4.4)
possibly are related to the differences in pore structures. It is also possible that higher organic
matter content and higher clay percentage in T7 aggregates (data not shown) could be the reason
for relatively high E.coli retention.
117

Figure 4.5: Bi-plot of soil aggregate samples and pore characteristics from PCA for a) the whole
soil aggregates and b) middle portion of the soil aggregates including center and sides. The threedigit numbers are labels for soil aggregates with ‘1xx’ representing T1 aggregates, ‘2xx’
representing T2 aggregates and ‘7xx’ representing T7 aggregates. Soil management was
converted to indicator dummy variables CT (T1) and NS (T7). Burn number equal to 1 (BN1) is
equivalent to pore size of 15-30 µm, BN equal to 2 (BN2) is equivalent to pore size of 45-60 µm,
BN equal to 3 (BN3) is equivalent to pore size of 75-90 µm, and BN equal to or greater than 4
(BN4g) is equivalent to pore size > 105 µm.

(a)

118

Figure 4.5(cont’d).

(b)

To relate E.coil transport from soil aggregates’ top to bottom part, we examine the pore
characteristics that connect these two parts, i.e. pore characteristics in the center and side portion.
Decomposing all variables in two-dimensional solutions using PCA with eigenvalues greater
than 1 accounted for 71% of the variance from soil pore characteristics in the center and sides
(Figure 4.5(b)). Similar as those of the whole soil aggregates’ pore characteristics, porosity, BN2
and BN3 are the most positive contributing variables to T1 aggregates’ center and sides.
Moreover, flow and connectivity related features, i.e., the number of paths, max flow and
tortuosity had positive loadings on T1 aggregates (Figure 4.5(b)). In general, this indicates pore
119

networks in center of T1 aggregates are more convoluted yet connected. This is in line with
aggregate pore characteristics (5 µm) from conventional tillage from another soil type reported
by Peth et al. (2008).
The relationship between E.coli spatial distribution in soil aggregates from the flow
experiment and soil pore characteristics were investigated using partial least square regression
(Figure 4.6(a)). Because of high variability of the intra-aggregate pores, only 49% of total
variation in pore characteristics and 31% of total variation in E.coil spatial distribution were
account using the first two PLSR factors. The R2 against number of factors in PLSR was plotted
in Figure 4.6(b). While three independent factors could explain 89% of total variation in pore
characteristics, only about 40% of total variation in E.coil data could be explained (Figure
4.6(b)). This indicates very high variability in E.coli data which could not be revealed by pore
configuration investigated in this study.
Among E.coli distribution variables, top was mostly correlated to pore characteristics,
followed by center, then sides and bottom (Figure 4.6(a)). At top of the aggregates, where E.coli
was applied, E.coli concentration was negatively related to max flow in the center portion of the
soil aggregates. E.coli retained in the center of the soil aggregates was somehow positively
related to porosity, medium to large pore sizes (BN3 and BN4g), and tortuosity of the pore space.
At aggregates sides, small pores (BN1) were more responsible for E.coli distribution.

120

Figure 4.6: (a) Correlation loading plot from PLSR between pore characteristics (porosity,
tortuosity, maxflow, paths, BN1-BN4g) and E.coli spatial distribution (TOP, CENTER, SIDES,
BOTTOM). The numbers are representing aggregate ids: 1-5 represent T1 aggregates, 6-10
represent T2 aggregates and 11-15 represent T7 aggregates. (b) R2 plots for percent of variation
explained vs. number of factors in PLSR.

(a)

121

Figure 4.6(cont’d).
(b)

It seems that max flow and tortuosity of the pore networks in aggregate center are mostly
related to E.coli distribution. Figure 4.7 shows that the difference in E.coli concentration
between top and bottom slices are somehow negatively related to max flow/tortuosity ratio in the
center portion of the pore space. However, the R2 is quite low mainly because of several high
influential points (marked as 125, 722 and 230). Visual assessment of aggregate 230 indicates
that there exist large pores connected from top and center portion to the bottom, which may
facilitate E.coli transport (Aggregate 230, Figure 4.8(a)). However, pores in the center is not well
connected which leads to low connectivity in aggregate 230 (Figure 4.7). Aggregate 722 shows
opposite trend in difference between top and bottom in E.coil concentration, where E.coli are
122

still mostly trapped at top with only a minor translocation to other parts of the aggregate (Figure
4.8(b)). There were disconnected medium to large pores found in the center of aggregate 722.
Abundance of randomly scattered pores of various sizes penetrating the entire aggregate led to
substantial E.coli redistribution to the aggregate’s sides as well as center and bottom (Aggregate
125, Figure4.8(c)).

Difference in E.coli distribution between top and bottom

Figure 4.7: The relationship between difference in E.coli concentration between top and bottom
slices and max flow/tortuosity ratio. Gray circles mark high influential points. Numbers along
with the circles are aggregates IDs.

5
722

4
3
2
1

R2 = 0.13

0
-1
125

-2
230
-3
0

2

4

Max flow : tortuosity ratio

123

6

8

Figure 4.8: 3D visualization of aggregates pore space for pieces. Pores with bright color indicate
high E.coli concentration. (a) Aggregate 230; (b) Aggregate 722; (c) Aggregate 125.

(a)

(b)

(c)

124

CONCLUSION
Despite high variability in enumeration of E.coli cells using qPCR analysis, E.coli
redistribution and transport within soil aggregates display a significant different spatial
distribution pattern under air-dry and saturated flow condition. When E.coli was first applied to
an air-dry aggregate, its resulting spatial distribution was highly heterogeneous in aggregates of
all management practices. After saturation, equilibration, and water extraction the E.coli
redistributed within the aggregate primarily in vertical direction and became more
homogeneously spread. However, only relatively small percent (T1: 2.5%, T2: 3.2%, T7: 1.5%)
of E.coli has completely left the aggregates. Redistribution was most pronounced in T1
aggregates, followed by T2, and was almost negligible in T7 aggregates.
T1, T2 and T7 aggregates are well separated by aggregate pore characteristics. There was
greater macro-porosity and medium sized pores (75-120 microns, BN2 and BN3) found in T1
aggregates compared to those of T2 and T7 aggregates. The differences in E.coli redistribution
observed between T1/T2 and T7 possibly are related to the differences in pore structures.
E.coli redistribution was related to intra-aggregate pore characteristics, however, because
of high variability of the intra-aggregate pores the relationships were weak. Only 49% of total
variation in pore characteristics and 31% of total variation in E.coil spatial distribution were
account using the first two partial least square regression factors. In aggregates’ top, E.coli
concentration was negatively related to max flow in the center portion of the soil aggregates.
E.coli retained in the aggregate’s centers was positively related to porosity, medium to large pore
sizes (BN3 and BN4g), and tortuosity of the pore space. Max flow to tortuosity ratio seems to
relate to difference in E.coli concentration between top and bottom slices, however, this
relationship is week (R2=0.13). It also should be pointed out that other than pore characteristics,
125

higher organic matter content in NS aggregates could be the reason for relatively high E.coli
retention. Future analyses are needed to address E.coli release and redistribution from soil
aggregates with repeated saturation and extraction processes, which simulates multiple rain
events.

126

BIBLIOGRAPHY

127

BIBLIOGRAPHY
Artz, R. R. E., J. Townend, K. Brown, W. Towers, and K. Killham. 2005. Soil macropores and
compaction control the leaching potential of Escherichia coli O157:H7. Environ.
Microbiol. 7:241–248.
Brennan, F.P., V. O’Flaherty, G. Kramers, J. Grant, and K.G. Richards. 2010. Long-Term
persistence and leaching of Escherichia coli in temperate maritime soils. Applied and
Environ. Microbiol. 76(5): 1449-1455.
CDC, Division of Bacterial and Mycotic Diseases.
http://www.cdc.gov/nczved/divisions/dfbmd/diseases/ecoli_o157h7/. Accessed 09/26/2011.
Edmonds, J. and R. M. Karp. 1972. Theoretical improvements in algorithmic efficiency for
network flow problems. Journal of the ACM 19 (2): 248–264.
Ford, L.R. and D.R. Fulkerson. 1962. Flows in networks. Princeton University Press, Princeton,
N.J.
Guber, A. K., D. R. Shelton, and Y. A. Pachepsky. 2005a. Transport and retention of manurebome coliforms in undisturbed soil columns. Vadose Zone J. 4: 828-837.
Guber, A. K., D. R. Sheiton, and Y. A. Pachepsky. 2005b. Effect of manure on Escherichia coli
attachment to soil. J. Environ. Qual. 34: 2086-2090.
Guber, A.K., Y. A. Pachepsky, D. R. Shelton, and O. Yu. 2009. Association of fecal coliforms
with soil aggregates: Effect of water content and bovine manure application. Soil Science
174 (10): 543-548.
Hagedorn, C., D. T. Hansen, and G. H. Somonson. 1978. Survival and movement of fecal
indicator bacteria in soil under conditions of saturated flow. J. Environ. Qual. 7:55-59.
Haugland, R.A., S.C. Siefring, L.J. Wymer, K.P. Brenner, and A.P. Dufour. 2005. Comparison
of Enterococcus measurements in freshwater at two recreational beaches by quantitative
polymerase chain reaction and membrane filter culture analysis. Water Research 39(4):
559-568.
Hellweger, F. L., and P. Masopust. 2008. Investigating the fate and transport of E. coli in the
Charles River, Boston using high-resolution observation and modeling. J. Am. Wat. Res.
Assoc., 44(2): 509-522.
Huysman, F., and W. Verstraete. Water-facilitated transport of bacteria in unsaturated soil
columns: Influence of cell surface hydrophobicity and soil properties. Soil Biology and
Biochemistry. 25 (1): 83-90.

128

Jamieson, R. C., R. J. Gordon, K. E. Sharples, G. W Stratton, and A. Madani. 2002. Movement
and persistence of fecal bacteria in agricultural soils and subsurface drainage water: A
review. Can. Biosyst. Eng. 44:1.1-1.9.
Jenkins, M.B., C. C. Truman, G. Siragusa, E. Line, J. S. Bailey, J. Frye, D. M. Endale, D.H.
Franklin, H. H. Schomberg, D. S. Fisher, and R. R. Sharpe. 2008. Rainfall and tillage
effects on transport of fecal bacteria and sex hormones 17β-estradiol and testosterone from
broiler litter applications to a Georgia Piedmont Ultisol. Science of the total environment
403: 154-163.
Johnson, W.P., and B.E. Logan. 1996. Enhanced transport of bacteria in porous media by
sediment-phase and aqueous-phase natural organic matter. Water Res. 30: 923-931.
Kellogg Biological Station. 2009. Agronomic protocol: Long-term ecological research (LTER)
in row-crop agriculture. Available at
http://houghton.kbs.msu.edu/data/agronomic_protocol/2009_lter_agronomic_protocol.pdf
(verified 09 October, 2011).
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.
Khan, I.U.H., V. Gannon, R. Kent, W. Koning, D.R. Lapen, J. Miller, N. Neumann, R. Phillips,
W. Robertson, E. Topp, E. van Bochove, and T.A. Edge. 2007. Development of a rapid
quantitative PCR assay for direct detection and quantification of culturable and nonculturable Escherichia coli from agriculture watersheds. Journal of Microbiological
Methods 69(3): 480-488.
Kinzelman, J.L., R.N. Bushon, R.T. Noble, and S.E. Dorevitch. 2010. WERF - PHASE I
Progress Report. Comparative Evaluation of Molecular and Culture Methods for Fecal
Indicator Bacteria for use in Inland Recreational Waters. In: WERF.
Kravchenko, A. N., H.-C. Chun, M. Mazer, W. Wang, J.B. Rose, A. Smucker, and M. Rivers.
Escherichia coli transport within soil macro-aggregates from long-term contrasting
management practices: relationships with pore-size distributions. Submitted to Journal of
Environmental Quality.
Lin, H. 2003. Hydropedology: Bridging Disciplines, Scales, and Data. Vadose Zone Journal 2: 111.
Lindqvist, R., and G. Bengtsson. 1991. Dispersal dynamics in groundwater bacteria. Microb.
Ecol. 21: 49–71.
Lindquist, W.B., Venkatarangan, A., Dunsmuir, J. and Wong, T.F., 2000. Pore and throat size
distributions measured from synchrotron X-ray tomographic images of Fontainebleau
sandstones. J. Geophys. Res. 105B: 21508.
129

Ling, T. Y., E. C. Achberger., C. M. Drapcho, and R. L. Bengtson. 2002. Quantifying adsorption
of an indicator bacteria in a soil-water system. Trans. ASAE 45: 669-674.
Michele, T. J., M. Cooley, D. Carychao, G. W. Wiscomb, R. A. Sweitzer, L. Crawford-Miksza,
J. A. Farrar, D. K. Lau, J. O’Connell, A. Millington, R. V. Asmundson, E. R. Atwill, and
R. E. Mandrell. 2007. Escherichia coli O157:H7 in Feral Swine near Spinach Fields and
Cattle, Central California Coast. Emerging Infectious Diseases, 13(12): 1908-1911.
Muirhead, R. W., R. P. Collins, and P. J. Bremer. 2006. Interaction of Escherichia coli and Soil
Particles in Runoff. Appl Environ Microbiol 72(5): 3406–3411.
Pachepsky, Y. A., A. M. Sadeghi, S. A. Bradford, D. R. Shelton, A. K. Guber, and T. Dao. 2006.
Transport and fate of rnanure-borne pathogens: Modeling perspective. Agric. Water
Manag. 66:81-92.
Rivers, M.L., S.R. Sutton, and P. Eng. 1999. Geoscience applications of x-ray computed
microtomography. Proceedings of SPIE, Developments in x-ray tomography II, 3772: 7886.
SAS Institute, 2008. SAS user's guide. Version 9.2. SAS Inst., Cary, NC.
Srinivasan, S., A. Aslan, I. Xagoraraki, E.Alocilja, and J. B. Rose. 2011. Escherichia coli,
enterococci, and Bacteroides thetaiotaomicron qPCR signals through wastewater and
septage treatment. Water Research 45: 2561-2572.
Vogel, H. J., I. Cousin, and K. Roth. 2002. Quantification of pore structure and gas diffusion as a
function of scale. Eur. J. Soil Sci. 53: 465-473.
Vogel, H. J., and K. Roth. 2003. Moving through scales of flow and transport in soil. J. Hydrol.
272: 95-106.
Vogt R.L., and L. Dippold. 2005. Escherichia coli O157:H7 outbreak associated with
consumption of ground beef, June–July 2002. Public Health Rep 120 (2): 174-178.
Wang, W., A.N. Kravchenko, A.J.M. Smucker, and M. Rivers. 2011. Comparison of image
segmentation methods in simulated 2D and 3D microtomographic images of soil
aggregates. Geoderma, 162: 231-241.
Wang, W., A.N. Kravchenko, A.J.M. Smucker, and M. Rivers. 2012. Intra-aggregate pore
characteristics: X-ray computed microtomography analysis. Soil Science Society of
American Journal, 76: 1159-1171.
Young, F. W., Y. Takane, and J. de Leeuw. 1978. The Principal Components of Mixed
Measurement Level Multivariate Data: An Alternating Least Squares Method with Optimal
Scaling Features. Psychometrika, 43: 279–281.
130

Chapter 5: CONCLUSION
Soil aggregates play an important role in several soil processes, e.g., C sequestration,
water flow and solute transport, and microbial activities. Advancement in X-ray
microtomography (µCT) imaging provides us opportunity to understand the mechanisms behind
these processes. Overall, the dissertation aims to effectively utilize soil aggregate µCT imaging,
to understand soil aggregate structure difference under contrasting land use and managements,
and to investigate bacteria redistribution and transport pattern and its relation to aggregate pore
characteristics.
Lack of ground-truth images poses great challenge to soil image segmentation. Proposed
simulation method successfully generated soil aggregate grey scale images from pore/solid
binary images, thus provided ground-truth information. Performance of several commonly used
image segmentation methods, including entropy, iterative, Otsu’s, and indicator kriging (IK)
were assessed using ground-truth based (ME) and non-ground-truth based criteria (NU). IK was
found preferable to other methods when the image histogram is clearly bimodal. NU performs
sufficiently well as a segmentation criterion, which could be used in assessing soil µCT image
segmentation results in future. Identification of irregular soil aggregate boundaries from µCT
images can be done using image closing technique with a properly determined neighborhood size
parameter, r. Examining the plots of the total aggregate volume, aggregate boundary fractal
dimension against r could serve as useful criteria for selecting the optimal r value. The optimal r
value might depend on a number of factors including the size of the studied aggregates, the
scanned resolution and the average sizes of pores connected to aggregate surface.

131

Using processed soil aggregate µCT images, soil aggregates pore structures from
contrasting land use and management as well as spatial gradient of pores in individual aggregate
was addressed. 20 years of contrasting management in the LTER experiment led to pronounced
differences in the intra-aggregate pore size distributions. Significantly greater percent of pores >
15 µm in LTER-CT aggregates was observed as compared to those of LTER-NS. LTER-NS
aggregates had more large pores (> 90µm) than LTER-CT aggregates, while more medium sized
pores (45-90 µm) were found in LTER-CT aggregates. With respect to spatial distributions,
pores were larger in the aggregate interiors than in the exterior layers; while intermediate pores
(45-90 µm) were more abundant in the aggregate exterior layers. This trend was present in all
land use and management treatments of the study and was most pronounced in LTER-NS. These
results implied that a general mechanism maybe responsible for aggregates formation, but long
term land use could alter the magnitude and intensity of involved soil processes. One implication
of aggregates’ pore size distribution on aggregate stability may serve as a warning indicator of
soil degradation on large scale. The difference in pore size distribution at aggregate scale
between contrasting land uses could possibly upscale to soil cores/fields to understand macropore flow, C sequestration mechanisms and water retention properties.
E.coli redistribution and transport within soil aggregates display a significant different
spatial distribution pattern under air-dry and saturated flow condition. When E.coli was first
applied to an air-dry aggregate, its resulting spatial distribution was highly heterogeneous in
aggregates of all management practices (CT, NT and NS). After saturation, equilibration, and
water extraction the E.coli redistributed within the aggregate primarily in vertical direction and
became more homogeneously spread. Only relatively small percent of E.coli has completely left
the aggregates. Redistribution was most pronounced in CT aggregates, followed by NT, and was
132

almost negligible in NS aggregates. E.coli spatial distribution was related to intra-aggregate pore
characteristics; however, because of high variability of the intra-aggregate pores the relationships
were weak. In aggregates’ top, E.coli concentration was negatively related to pore connectivity
in the center portion of the soil aggregates. E.coli retained in the aggregate’s centers was
positively related to porosity, medium to large pore sizes and tortuosity of the pore space. It
could be speculated that other than pore characteristics, complex processes are also involved for
E.coli transport in soil aggregates. Implication of these results may help us understand potential
in-field contamination caused by E.coli movement to groundwater under precipitation events.
Further research is necessary in several aspects: developing new approaches for
multiphase (e.g., pore, water, and solid; pore, organic matter, and solid) segmentation in soil
images; modeling pore scale flow using fluid dynamic theories and related theory to lab
measurements in order to understand water and solute transport mechanisms; describing and
quantifying the role of intra-aggregate pore size distributions in a number of soil processes,
including, but not limited to C sequestration; confirming spatial pattern in pore size distributions
within soil aggregates of other sizes; evaluating E.coli release and redistribution from soil
aggregates with repeated saturation and extraction processes, which simulates multiple rain
events.

133