2 MATERIALS AND METHODS
2.1 Experimental design and soil sample collections
Two of the experimental sites used for this study were located in
Wisconsin (Oregon and Hancock) and three in Michigan (Lux Arbor, Lake
City, and Escanaba), USA. The soils of Oregon, Lux Arbor, and Escanaba
sites are Alfisols, and of Hancock and Lake City sites are Entisol and
Spodosol, respectively. At each site a randomized complete block design
experiment with 3 (Hancock) or 4 (the rest of the sites) replications
has been established in 2013. Details on the research sites and soil
descriptions have been reported by Kasmerchak and Schaetzl (2018) and
Lee et al. (2023).
The two studied vegetation systems were: (1) non-fertilized monoculture
switchgrass (Panicum virgatum L.; Cave-in-Rock variety); and (2)
restored prairie, which consisted of 18 plant species of grasses
(including switchgrass), forbs, and legumes. Soil sampling was conducted
in 2020 (Oregon site) and 2019 (the other 4 sites). Two types of soil
samples were collected from each replicated plot. First, three intact
soil cores (5 cm in height and 5 cm Ø) were collected from 5 to 10 cm
depth for X-ray computed micro-tomography (µCT) scanning. Then, the
loose soil surrounding the cores was also collected for measurements of
soil texture and mineralogy. All samples were stored at 4 °C until
scanning and measurements.
2.2 Measurements of soil porosity, texture, and mineralogy
Soil porosity was calculated from bulk and particle densities of the
collected samples (i.e., total porosity), and texture was determined
using the hydrometer method (Gee and Or, 2002) in all replicated plots
of all sites. Soil mineralogy composition was measured using X-ray
powder diffraction (XRD) carried out at the Illinois State Geological
Survey (Champaign, Illinois, USA). Because of the high costs of XRD
analyses only three replicated samples from each system in each site
were subjected to these measurements. Prior to XRD analyses, the samples
were cleaned, dried in a vacuum oven, and ground to < 44 µm.
Then, one subset of the prepared sample was powdered by the McCrone mill
(MBP) (McCrone Accessories & Components, Westmont, IL, USA) for
quantification of non-clay minerals (quartz, clay, K-feldspar,
P-feldspar, calcite, dolomite, siderite, and pyrite/marcasite), and the
other was fractionized into < 2 µm powder for clay minerals
(smectite, mica, kaolinite, and chlorite). The prepared two types of
powders were then spread on a glass slide and analyzed using a
Siemens/Bruker D5000 X-ray Powder Diffraction instrument (Billerica, MA,
USA). The JADE™ software was used to identify percentages of
constituents in each powdered sample from XRD patterns.
2.3 Soil core scanning and image analysis
Pore structure assessments were performed via X-ray µCT. Soil cores were
drained at -28 kPa using a 5-bar pressure plate extractor (Soilmoisture
Equipment Corp., Goleta, CA) prior to μCT scanning to remove water from
pores Ø >10 μm and increase the contrast between the solids
and air on X-ray μCT images. Then, the cores were scanned using X-ray
μCT machine (North Star Imaging, X3000, Rogers, MV, USA) at the
Department of Horticulture facility, Michigan State University. The
energy settings were 75 kV and 450 μA. The scanning resolution of 18 μm
was achieved using the Subpix-mode of the scanner, combining four
individual scans shifted half pixel in vertical and horizontal
directions. Scanned images of 3014 projections were reconstructed by efX
software (North Star, Rogers, MN, USA).
A schematic summary of the steps involved in the image processing for
this study is outlined in Figure 1. The image pre-processing was
conducted using ImageJ-Fiji software (Schindelin et al., 2012) to remove
artifacts and noises. First, to exclude sampling artifacts near the soil
core walls, images were centered and cropped into prisms (1500 × 1500 ×
2240 pixels corresponded to 2.7 cm in length, 2.7 cm in width, and 4.1
cm height). Then, ‘Remove Background’ tool in Xlib/Beat plugin was used
to remove shadowing effects from the images, followed by the removal of
ring artifacts on the image’s polar domain using a stripe filter of the
Xlib/Beat plugin. After that, a 3D non-local mean filter (σ = 0.1)
implemented in scikit-image (Walt et al., 2014) was used to reduce the
noise (Buades, Coll, & Morel, 2011; Darbon, Cunha, Chan, Osher, &
Jensen, 2008). The pre-processing steps dropped the resolution of images
from 18 to 36 μm.
Root residues, which we will refer to as particulate organic matter
(POM), were segmented from the filtered images with Ilastik software, a
machine learning-based tool (Berg et al., 2019). A random forest
classifier was used on a multi-dimensional feature space of the filtered
gray scale images. The classifier was trained using two cores from each
combination of different vegetation and sites (20 of total 114 cores)
and then applied on entire cores. The training dataset produced
out-of-bag error rate estimates less than 1.8% in overall, and all
segmented POM images were visually inspected to ensure the accuracy and
integrity of the segmentation process. The outcome of POM segmentation
was denoised by removing objects smaller than 4 voxels in diameter from
the images.
Segmentation of the filtered grayscale images into pore and solid binary
images was performed to identify the pores visible at the image
resolution, referred further on as image-based pores. For each sample
the segmentation threshold was estimated as an ensemble of six
segmentation methods (i.e., Otsu, Triangle, Huang, IsoData, Li, and
Moments). The global thresholds for the stack of images in each
individual core, estimated using the six segmentation methods, were
averaged and applied to that stack to separate the solid and air-filled
voxels in the images using SimpleITK in Python (Beare, Lowekamp, &
Yaniv, 2018; Lucas et al., 2022). Obtained images were used to compute
pore size distributions using ‘Local Thickness’ tool, an approach based
on the maximum inscribed sphere method (Hildebrand & Rüegsegger, 1997;
Vogel, Weller, & Schlüter, 2010).
Biopores were identified as described by (Lucas et al., 2022; Lucas,
Schlüter, Vogel, & Vetterlein, 2019b). Specifically, to employ
tubular-shaped features of biopores in differentiating them from other
irregularly shaped pores, we used the Tubeness plugin in ImageJ-Fiji for
shape detection. As rising σ-values significantly increased the
computational time, binary images were scaled down to 50% and 20% for
Tubeness filtering with σ-values ranging from 1-4 and 2-30, for each
scale respectively, with a step size of 1. Gaussian blurring was applied
to the entire binary image with varying σ-values in order to efficiently
identify biopores of various diameters. The resulting tubular channels
were slightly smaller than the root channel itself due to the exclusion
of rough surface on biopore walls. Thus, to better capture the actual
width of biopores, 3D dilation steps were employed as a postprocessing
measure. After combining all elongated objects, misclassified objects
were removed (Phalempin, Lippold, Vetterlein, & Schlüter, 2021a).
Proportions of biopores in the entire pore system were calculated. After
that, proportions of biopores occupied by POM were computed by first
calculating the volume of POM located in biopores and then dividing this
volume by the entire volume of biopores.
Nine masks corresponding to interval regions nine distances away from
the segmented POM (0-0.25, 0.25-0.5, 0.5-1.0, 1.0-1.5, …,
3.0-3.5, and 3.5-4.0 mm) were created using 3D distance transform in
ImageJ-Fiji. Then, masks of interval regions were applied to the
pore-solid segmented image and to the pore size distribution image of
the entire sample to calculate the porosities and the size distributions
individually for each interval region. Contributions of pores of
different size classes to image-based porosity of the distance interval
regions were expressed as pore fractions (%). We considered three pore
size classes, namely 36-150 μm, 150-300 μm, and > 300 μm Ø.
The 36 μm Ø corresponded to the smallest pore size that could be
reliably detected on the studied images. Pores < 150 μm Ø are
known to have especially high microbial activity and strongly contribute
to the C processing (Kravchenko & Guber, 2017; Kravchenko et al., 2019;
Strong et al., 2004), and pores < 300 μm Ø function as the
secondary pathways for water and nutrient supplies to resident
microorganisms (Franklin et al., 2021).
A Connectivity tool of BoneJ plugin in the ImageJ-Fiji was used for the
pore connectivity calculations: first, Euler numbers (χ) were computed,
and the numbers were divided by the total volume of corresponding
regions (V ) (Odgaard & Gundersen, 1993; Vogel & Roth, 2001):
\(\chi_{V}=\ \frac{N-C+H}{V}\) (1)
Where N is the number of isolated objects, C is the number
of redundant connections or loops, and H is the number of
completely enclosed cavities, which are typically negligible in soil
pore system (Lucas, Vetterlein, Vogel, & Schlüter, 2021; Vogel, 2002).
The minimum size of the object was 2 × 2 × 2 voxel. Higher, e.g.,
positive, χV values calculated via Eq. (1)
correspond to lower connectivity, while lower, e.g., negative, to the
higher connectivity. To simplify the presentation of the connectivity
data we report the results as negative values of
χV , that is, the high values of
-χV correspond to high connectivity while the low
values to the low one.
Since pore connectivity can be affected by the volume of the soil in
which it is calculated, we did not assess it at the same distance
intervals as those that were used for the image-based porosity and pore
size distributions described above. Instead, we only calculated it in
immediate vicinity of the residue, i.e., the region of 0-0.25 mm away
from the POM, and for the entire soil volume. The resultant two
estimates of the connectivity were used for comparisons among the five
experimental sites and plant systems.
2.4 Statistical analysis
The data were analyzed using SAS 9.4 (SAS Institute Inc., NC, USA)
procedures of PROC MIXED and PROC GLIMMIX. Since we did not expect that
6-7 years of disparate vegetation influenced soil texture and
mineralogy, the statistical models for texture and mineralogy
characteristics included only the fixed effects of the experimental
sites. For the other soil properties, the statistical models included
the fixed effects of sites, plant systems, and their interaction.
Statistical models for analyses of image-based porosity data at
different distances from POM additionally included the same fixed
effects as the soil properties and individually tested by distance
intervals, as the interaction among the sites, plant systems, and
distances was significant. All models included the random effects of
experimental blocks nested within the sites and, when necessary, the
random effects of cores nested within the blocks, plant systems, and
sites. The latter were used as an error term for testing the plant
system effects. The assumptions of normality and variance homogeneity
were assessed using normal probability plots, plots of residuals vs.
predicted values, and Levene’s tests for equal variances.
Additionally, we grouped the sites into two soil texture classes for
comparing pore size distributions and connectivity between
finer-textured soils and coarser-textured soils. The first group
included Oregon, Lux Arbor, and Escanaba sites, the three soils with
< 66 % sand content, and the second group consisted of
Hancock and Lake City sites with > 82 % sand content
(Table 1). Models for analysis of pore size distribution data within
each distance interval and of connectivity data within 0-0.25 mm
distance and entire image stack included the fixed effects of plant
systems, soil groups, and their interaction, and random effects of
experimental sites nested within soil groups, blocks nested within the
sites, and cores nested within the blocks, plant systems, sites, and
groups.
Linear relationships among soil texture variables, mineralogical and
clay mineralogy variables, proportions of biopores and POM in biopores,
and the distance-based porosities and connectivity were assessed using
Pearson’s correlation coefficients.