X-ray microtomography analysis of soil structure deformation caused by centrifugation

Centrifugation provides a fast method to measure soil water retention curves over a wide moisture range. However, deformation of soil structure may occur at high angular velocities in the centrifuge. The objective of this study was to capture these changes in soil structure with X-ray microtomography and to measure local deformations via digital volume correlation. Two samples were investigated that differ in texture and rock content. A detailed analysis of the pore space reveals an interplay between shrinkage due to drying and soil compaction due to compression. Macroporosity increases at moderate angular velocity because of crack formation due to moisture release. At higher angular velocities, corresponding to capillary pressure of ψ <−100 kPa, macroporosity decreases again because of structure deformation due to compression. While volume changes due to swelling clay minerals are immanent in any drying process, the compaction of soil is a specific drawback of the centrifugation method. A new protocol for digital volume correlation was developed to analyze the spatial heterogeneity of deformation. In both samples the displacement of soil constituents is highest in the top part of the sample and exhibits high lateral variability explained by the spatial distribution of macropores in the sample. Centrifugation should therefore only be applied after the completion of all other hydraulic or thermal experiments, or any other analysis that depends on the integrity of soil structure.


Introduction
Soils, rocks and sediments are assumed to be rigid bodies in many modeling applications.Yet, the internal structure of these porous media is modified through a variety of technical and natural processes.The internal changes can either be gradual, e.g., through dissolution, biological activity or swelling/shrinking, or abrupt, e.g., landslides or tillage.Conventional laboratory methods can only provide a limited set of structural properties, such as bulk density and porosity, or provide indirect information through functional properties that are governed by the internal structure, such as gas diffusion, permeability or stress-strain relationships.Direct information on the deformation of the internal pore architecture is typically missing.X-ray microtomography has turned into a standard technique to fill this gap and measure the threedimensional internal structure of porous media (Ketcham and Carlson, 2001;Cnudde and Boone, 2013;Wildenschild and Sheppard, 2013).There is a huge variety of image processing and image analysis methods that are all tailored for the ultimate goal to quantify the complex, structural heterogeneity based on a few meaningful parameters (Kaestner et al., 2008;Vogel et al., 2010;Schlüter et al., 2014).The changes in the internal structure can be assessed statistically, e.g., by comparing the pore size distribution or pore connectivity averaged over different samples at two points in time (Jégou et al., 2002;Schlüter et al., 2011).Evidently, spatially explicit information about the internal displacement of particles or aggregates is excluded from analysis in such an approach.However, this local deformation information is of particular interest, e.g., in soil mechanics (Terzaghi et al., 1996).So S. Schlüter et al.: Soil structure deformation far there are only a few approaches to measure the deformation explicitly via imaging and image analysis.One method would obviously be to manually identify identical objects in two consecutive images and measure their separation distance.Repeating this measurement for many objects would then populate the deformation field.However, this is impractical for large, three-dimensional data sets, because of the large number of measurements required to reach an appropriate density of displacement information.Therefore, automatic methods are needed.
Automated methods to detect deformation are usually based on digital volume correlation (Bay et al., 1999;Lenoir et al., 2007;Peth et al., 2010;Hall, 2010;Son et al., 2012).The rationale of this method is to recover the displacement field by finding a geometric transformation of a deformed image that optimizes a correlation coefficient with the original, undeformed target image.The method usually comprises three steps (Bay et al., 1999): (1) the acquisition of X-ray microtomography image before and after the perturbation, (2) image registration of one image onto the other to obtain a discrete deformation vector field and (3) calculation of the strain tensor field from the displacement vector field.Another popular method for change detection is called particle image velocimetery (PIV) (White, 2003).PIV was originally developed to visualize flow paths within a fluid by tracking small, yet visible particles over time (Adrian, 1991).In this method the image is divided into a large number of sub-windows, and the displacement vector for each window is calculated via cross-correlation between two consecutive images, which results in a local velocity for the given time lag.The particles can be substituted by any moving feature such as a growing root (Bengough et al., 2010), soil displacement along an earthworm burrow (Barnett et al., 2009) or soil creeping along slopes (Baba and Peth, 2012).A serious shortcoming of PIV is that so far it can only be applied to two-dimensional sections, in which case any displacement out of plane into the third dimension is excluded.
The objective of this paper is to measure local deformations in a soil core caused by centrifugation with digital volume correlation.We put special emphasis on common pitfalls and best practices for image registration, which is the critical step for a successful application of this method.Centrifugation has been chosen for the deformation analysis in this study because it is suitable to evoke structure deformation under controlled conditions.Measuring the water loss through centrifugation of a soil is a rather old method for determining the water retention curve of a soil (Gardner, 1937;Russell and Richards, 1939;Oden, 1975;Reatto et al., 2008).It has obvious advantages over other conventional methods like multi-step outflow (Van Dam et al., 1994;Vogel et al., 2008) or evaporation (Simunek et al., 1998;Peters and Durner, 2008) in that the method is less time consuming, it captures a wide moisture range of the retention curve and provides a good reproducibility through defined experimental conditions (angular velocity, temperature, pressure).
More recently, the steady state centrifugation method was also used to measure the unsaturated hydraulic conductivity of soils (Šimůnek and Nimmo, 2005;Van den Berg et al., 2009;McCartney and Zornberg, 2010).However, a serious drawback is that soil can get compacted through the centrifugal force, i.e., the inertia that acts on the sample during rotation (Khanzode et al., 2002;Wedler and Boguslawski, 1965).The spatial variability of compaction within a sample will be investigated in this study.

Overview of core concepts
Since image registration is the backbone of our deformation detection method, we provide a brief overview over the core concepts.Optimal spatial alignment of an altered image with a target image is usually achieved by minimizing an objective function that quantifies the mismatch in terms of a predefined metric.A standard metric would be the correlation coefficient between co-localized voxels (Latham et al., 2008;Hapca et al., 2011).More advanced metrics are based on information theory like mutual information criteria (Mattes et al., 2001).Composite metrics are also possible, e.g., with an additional regularization term in case of elastic registration or an additional term for the mismatch of manually defined point pairs, so-called landmarks.The geometrical transform that produces an optimized alignment can be grouped according to the degrees of freedom by which the image can move.The simplest transform is rigid registration with six degrees of freedom, i.e., rotation around three axes and translation in three directions.A similarity transform includes isotropic scaling as one additional parameter.Affine transformations possess 12 degrees of freedom through rotation, translation, shearing in three directions and anisotropic scaling.All aforementioned transforms are global in that a single transformation matrix is uniformly applied to the entire image.Thus they do not allow for non-uniform deformation.Elastic registration with a B-spline transform, in turn, requires a regular grid of control points over the image where each performs independent, affine registrations.As a consequence, the control points move relative to each other in the course of registration.Image registration is achieved by an iterative optimization scheme with standard methods like gradient descent or more involved methods like adaptive stochastic gradient descent (Klein et al., 2010) with optimal trade-off between speed and robustness.This completes the description of the core methodology of automated image registration.However, based on our own experience a straightforward implementation can lead to a failure of the method without following some best practices.Salient objects in the altered and the target image need to have substantial overlap to identify a meaningful correspondence of features.This is facilitated by a manually edited transformation matrix, in which the resolution of the altered image is adapted to the resolution of the target image and an offset between the images is corrected by some initial translation or rotation of the altered image.This manual procedure can be avoided by exhaustive sampling, where all combinations of translations (and rotations if required) are tested at a coarse grid representation (Latham et al., 2008).A different kind of user input is provided by landmarks placed on the salient features in both images.The cumulative Euclidean distance between point pairs is then added as a second term to the objective function.In all cases, the correct orientation of the sample has to be ensured a priori, because a flipped orientation cannot be recovered by rotation.

Concatenation
It is often useful to combine several consecutive image registration steps with each other, especially when an elastic registration with thousands of degrees of freedom is involved.Then a rigid registration can place the altered image already close to the global minimum and thus anticipate a collective movement of all control points in the subsequent B-spline transform.In our example, a two-step procedure will first match the position of the cylindrical container in the original sample with that of deformed sample after centrifugation, by aligning corresponding landmarks on the core wall.Then, a B-spline transform of the centrifuged sample onto the original sample will recover the remaining internal deformations.

Pyramid schedule
As mentioned above, a coarse representation of both images at reduced resolution vastly reduces the processing time so that it may even allow for exhaustive sampling of all transform parameter combinations.At the same time the coarsening reduces the iterations necessary to achieve a certain translation in physical space.This can be achieved with a socalled pyramid schedule.That is, the registration is started at a coarse representation and when convergence is reached the registration is continued at the next pyramid level with higher resolution.Finer details, e.g., small rocks, that did not appear at a coarse level are then used to refine the registration results.

Contrast enhancement
The success of image registration depends on the existence of salient features in the two corresponding images.Here, a particular problem for structure deformation in soil may take effect that macropores are not rigid and may disappear completely due to compression.Likewise, cracks may form in dry soil that were not present at high moisture content.However, a lack of corresponding features may impair the success of image registration.To avoid this, the gray values in both images can be rescaled such that pore and soil matrix voxels have zero gray value and all rigid rocks are depicted with optimal contrast.In this way only the rock matrix is used for registration, i.e., an assemblage of bodies that change their position but not their shape.This raises the questions, what the minimum amount of rocks has to be in order to still get good image registration results.Therefore, registration results will be compared for soils with high and low rock content and also with and without contrast enhancement.

Soil sampling
Undisturbed soil was sampled from two locations using a custom-made drill for undisturbed sampling of cylindrical soil cores (UGT GmbH, Germany).One sample was taken from the upper soil layer (in 5 to 15 cm depth) of a meadow in Köllme near Halle (Saale), Germany.The second sample stems from the plow horizon of a fallow plot (in 5 to 15 cm depth) at the experimental station of the Helmholtz-Centre for Environmental Research -UFZ in Bad Lauchstädt, Germany.The samples were covered with a lid, carefully transported to the laboratory and stored in a refrigerator at 4 • C to prevent the soil from drying and to reduce biological activity.The soil samples had a mean height of 7.9 cm, a diameter of 9.4 cm.The Köllme soil had an initial bulk density of 1.24 g cm −3 and a rock content of 3 %.The bulk density in the Bad Lauchstädt soil was higher (1.32 g cm −3 ) and the rock content was much lower (0.2 %).The Köllme soil has a slightly coarser texture with 24 % sand, 58 % silt and 18 % clay as compared to 5 % sand, 71 % silt and 24 % clay for the Bad Lauchstädt soil (disregarding organic matter and rocks > 2 mm in diameter).

Centrifugation
A perforated aluminum plate covered with a filter paper (pore size 10 µm) was installed in between the lower boundary of the sample and a reservoir that collects the drained water.This setup was placed in a centrifuge (Cryofuge 6000i, Heraeus GmbH, Germany) and covered with a plastic film to prevent evaporation.The capillary pressure was initially adjusted to full saturation.Afterwards, the sample was centrifuged at increasing angular velocities.The equivalent capillary pressure of water ψ in equilibrium with the gravitational field written in differential form is as follows (Gardner, 1937;Russell and Richards, 1939): Integrating Eq. ( 1) between the inner (r i = 0.136 m) and outer (r o = 0.22 m) boundaries of the soil core one obtains where normalizing by the density of water ρ results in work per unit mass of water and hence pressure [Pa].The angular velocity ω was calculated by ω = 2π N , where N is the revolution frequency [rad s −1 ].

X-ray microtomography
The soil cores were scanned with an X-ray microtomograph (X-TEk XCT 225, Nikon Metrology) with slightly different energy settings for the Köllme soil (150 keV, 425 µA) and the Bad Lauchstädt soil (160 keV, 500 µA).In both cases we used a 1 mm copper filter to reduce beam hardening artifacts and prevent overexposure at the lateral margins of the detector panel.An entire scan comprised 2749 projections with an exposure time of 708 mS (one frame per projection) which resulted in a scan time of roughly half an hour.A detector panel (Elmer-Perkin 1620) with 1750 × 2000 pixels (200 µm resolution) captures projections with 16-bit precision.The reconstruction of three-dimensional images via filtered back projection was done with the CT Pro 3-D software package (version 3.1) at a spatial resolution of 61 µm and 8 bit grayscale resolution.

Image processing
The deformation of soil structure is analyzed in two ways.First, the changes in the macropore structure are analyzed by comparing pore size distributions as well as depth profiles of macroporosity for images before and after centrifugation.Second, the displacement of soil constituents is explored with digital image correlation, i.e., via image registration of the deformed soil to the original soil prior to centrifugation.The entire image processing workflow for pore space analysis (top row) and digital image correlation (bottom row) is summarized in Fig. 1.As a first step for the pore space analysis, noise in the raw images is removed with a non-local means filter implemented in ITK (image registration tool) (Buades et al., 2005;Tristán-Vega et al., 2012).This version of the non-local means filter requires two parameters: (1) the radius of the search window is set to four pixels and an estimate of the noise level expressed as a standard deviation of noise is set to 60; (2) vertical differences in average image intensity due to shading and cone beam artifacts are removed.To do so, the mean gray value of the soil matrix in a specific z-slice is measured and the difference to the mean gray value of the soil matrix in the entire image is subtracted from each voxel in that depth (Iassonov and Tuller, 2010).The two thresholds that separate the soil matrix from darker pore voxels and brighter rock voxels are chosen manually and do not effect the results much as long as they cover the entire grayscale range of soil without adding pores and rocks to the average.The main purpose of noise removal and intensity drift correction is to improve the robustness of automatic, histogram-based threshold detection methods.In this study the multi-level version of the popular Otsu method was used (Liao et al., 2001) to obtain the thresholds between four classes (pores, soil matrix, low density rocks and high density rocks), of which only the pore class is analyzed subsequently.Segmentation into pores and background was carried out with simple thresholding.The resulting pore space was analyzed towards two different directions.First, pore size distributions were calculated with the maximum inscribed sphere method using the BoneJ plugin (version 1.3.12) in Fiji (Doube et al., 2010).Second, depth profiles of porosity were computed in equidistant steps of 10 pixels to monitor the depth dependent changes of macroporosity caused by drying and compaction.All image processing steps except noise removal and pore size analysis were performed with the QuantIm image processing library (Vogel et al., 2010;Schlüter et al., 2014).
Digital image correlation is applied to the raw image with only little preprocessing.The images are resampled at a 3times coarser resolution of 183 µm via simple averaging in order to reduce memory consumption.Optionally, the contrast is stretched such that macropores and soil are binned to zero and the rocks are displayed with optimal contrast so that only the rock matrix is used for subsequent image registration.The lower boundary under which all gray values are set to zero is chosen manually.Image registration was carried out with elastix 1 (Klein et al., 2010), an open-source image registration tool based on ITK (Ibanez et al., 2005) which was tailored for medical imaging applications.The image of the original soil structure was registered to the deformed soil in two steps.First, we applied a rigid registration that minimizes the sum of Euclidean distances between corresponding landmarks which we placed at four notches in the container wall.The mutual information criterion (Mattes et al., 2001) between images is added to the objective function as a second term with a low weighting factor because the manual selection of landmarks is typically afflicted with an imprecision of a few pixels.This rigid registration forces the cylindrical containers of both images to be well aligned.Then we applied elastic registration with a B-spline transform that optimizes the mutual information criterion to recover internal deformations.A regularization term, called bending energy penalty, was added to the objective function which ascertains that local transitions in deformation magnitude and direction are smooth (Klein et al., 2010).The success of image registration is evaluated via visual comparison of the target image and the aligned image.The outcome of an image registration process is a parameter file that stores the global transformation matrix (for rigid registration) or the transforma- tion matrix of all grid nodes (for elastic registration).This parameter file is then applied to the input image in order to construct a new image that is perfectly aligned with the target image.This transformation matrix can also be applied to other images that were not used in the image registration process.This is essential since the image registration is optimized with auxiliary images which only contain rocks while the resulting transformation matrix is applied to the raw images that display the original soil structure.Two sample images and all parameter files to reproduce the entire workflow are available from the authors upon request.

X-ray microtomography
The three-dimensional structure of the Köllme soil before centrifugation is depicted in Fig. 2a.Rocks of different sizes are embedded in a loamy soil matrix and large macropores are present in all depths.After centrifugation down to a capillary pressure of ψ = −100 kPa the structure is markedly deteriorated (Fig. 2b).The shape and position of macropores have changed, and desiccation cracks have formed in the vicinity of macropores.A perforated, rigid plate is visible, which was mounted at the bottom of the sample to prevent soil loss during centrifugation and transport.At a capillary pressure of ψ = −500 kPa the soil is severely compacted as can been seen by the larger head space above the soil surface (Fig. 2c).Macropores are almost absent especially in the lower part of the sample and only a few desiccation cracks with vertical orientation remain in the upper part of the profile.
The low rock content of the Bad Lauchstädt soil is evident in Fig. 2d-f.After centrifugation down to a capillary pressure of ψ = −50 kPa all macropores are intact and have even slightly increased in size due to shrinkage of clay minerals in the course of drying.At a capillary pressure of ψ = −300 kPa the original macropore network is still visible but some macropores are partially compressed.
These visual observations are confirmed by quantitative analysis of the pore space (Fig. 3a).The total visible macroporosity (>61 µm) in the Köllme soil at full saturation is 6.6 %.At ψ = −100 kPa, macroporosity increased to slightly (8.1 %) due to crack formation.These new desiccation cracks due to soil drying mainly formed in that size range of 0.1-1 mm.Macropores >1 mm in diameter are less abundant in the soil at ψ = −100 kPa because larger macropores were partially compressed.The subsequent centrifugation to a capillary pressure of ψ = −500 kPa reduced visible macroporosity to 2.0 % and removed porosity in all size ranges.A similar trend is evident in the Bad Lauchstädt soil (Fig. 3b) though with lower magnitude.The visible porosity (>61 µm) at full saturation (9.4 %) increases slightly when centrifuged to ψ = −50 kPa and decreases to 7.6 % at a capillary pressure of ψ = −300 kPa.The changes in the pore size distribution are smaller in the Bad Lauchstädt soil due to less negative capillary pressures and a lower abundance of big macropores > 1 mm in diameter which are most easily compressed.The change in porosity through centrifugation is not evenly distributed across the sample (Fig. 4a).In the Köllme soil desiccation cracks at ψ = −100 kPa mainly formed in the top part of the sample, whereas the bottom part of the sample exhibits lower porosity than the reference sample due to compaction.The sample in its driest state is compacted across the entire profile.The decline in macroporosity increases with depth.A specific rock close to the soil surface served as a cut-off height for the profile.Its position changed from 68 to 58 mm and to 55 mm, respectively.In the Bad Lauchstädt soil (Fig. 4b) the same trends apply though at different magnitude.At ψ = −50 kPa macroporosity increases in the top part of the sample due to drying, whereas the bottom exhibits slightly smaller macroporosity due to compression.Further drying down to ψ = −300 kPa reduces macroporosity in all depths.The height of the sample as measured by the position of a reference rock decreases from 85 to 81 mm for ψ = −50 kPa, followed by 77 mm for ψ = −300 kPa.

Deformation
Rocks are especially suited to track internal deformations as they change in position, but not in shape.The poor spatial alignment of rocks between the saturated (green) and the centrifuged soil at a capillary pressure of ψ = −500 kPa (red) after the Euler transform is depicted in Fig. 5a (results for deformation at ψ = −100 kPa not shown).Elastic registration with a B-spline transform leads to a very good spatial alignment of rocks in all soil depths (Fig. 5b).In both images (original and deformed soil) the contrast had been optimized for rocks prior to image registration.Figure 5c shows the registration results for the resampled raw images.That is, the full grayscale range including soil matrix and macropores is used for digital volume correlation and the contrast is optimized for rocks only afterwards for visualization.This results in a less accurate spatial alignment of the rock ma-trix especially in the lower part of the sample where deformation is strongest and macropores are compressed.Even though the rock content in the Bad Lauchstädt soil is very low (0.2 %), it is sufficient to achieve a good spatial alignment when only the rock matrix is used for image registration (Fig. 5e).Results are shown for the most compacted state only (ψ = −300 kPa).The spatial alignment of the rock matrix works comparably well, if the entire grayscale range including soil and macropores is used for image registration (Fig. 5f).Hence, taking the spatial alignment of macropores also into account does not impair the spatial alignment of the rock matrix, as most of the macropores persist during the centrifugation process.A general advice whether to use the whole grayscale range or not will be discussed below.
An important result of the registration procedure is the displacement vector field.There is a clear trend towards a downward movement of soil constituents in the Köllme soil as a consequence of compaction and its magnitude increases from ψ = −100 kPa to ψ = −500 kPa (Fig. 6a-b).However, the direction and length of local displacement does not only vary with depth but also laterally.Furthermore, there is a substantial horizontal component of displacement in many locations.This lateral movement can have two different origins.First, the formation of mainly vertically aligned cracks displaces the soil normally into the crack but not along it.Secondly, regions of high macroporosity are preferential failure zones during compaction.Filling these macropores during compression with soil material from above evokes a lateral displacement component because they are not evenly distributed across the soil.The Bad Lauchstädt soil at ψ = −50 kPa exhibits a rather uniform downward displacement of soil constituents of about 3-5 mm (Fig. 6c) that is in line with the reduction of sample height shown in Fig. 4. A pressure of ψ = −300 kPa evokes lateral differences in down- ward displacement again (Fig. 6d), as already discussed for the Köllme soil.

Pore scale processes during centrifugation
The analysis of macroporosity at 61 µm resolution revealed substantial alterations of the pore space architecture during centrifugation.Without such a detailed X-ray microtomography analysis the only measurable, macroscopic changes in soil structure would have been an increase in bulk density and a decrease in sample height.The conventional, quantitative image analysis of the pore space revealed a depth-dependent increase or decrease of macroporosity that resulted from the interplay of soil shrinkage due to drying and soil compaction due to compression.The reduction of macroporosity due to compaction was most severe in the lower part of the samples.Evidently this is because the inertial force that acts on the soil in a given depth, i.e., F z = mω 2 r, increases both with increasing overburden m and increasing absolute acceleration (due to r o > r i ) (McCartney, 2007).This was corroborated by displacement vector fields obtained from digital volume correlation which showed an increase of vertical displacement with sample height.Evidently, this is due to the fact that the local displacement integrates over the distance to the lower boundary.
The deformation of the soil due to shrinkage is immanent in any drying process in presence of swelling clay minerals and capillary forces that pull unconsolidated grains closer together (Or and Ghezzehei, 2002;Stange and Horn, 2005).So it would have also occurred if drying is induced by another process, e.g., by evaporation.The compaction of the soil through centrifugal forces, however, is obviously caused by the centrifugation process and represents a severe drawback of the method.A significant breakdown of structure through centrifugation was previously reported for an equivalent capillary pressure of ψ = −100 kPa (Wedler and Boguslawski, 1965) using a loess soil with a texture (5 % sand, 77 % silt, 18 % clay) comparable to the Bad Lauchstädt soil.The deformation field analysis of the Bad Lauchstädt soil confirmed that at a capillary pressure of ψ = −50 kPa, the soil structure is still intact with fairly uniform downward displacement of all soil constituents by a few mm, whereas centrifugation at ψ = −300 kPa caused stronger lateral heterogeneity in deformation caused by the compression of heterogeneously dis-  Cartney, 2007;McCartney and Zornberg, 2010).It is clear that the susceptibility to compaction during centrifugation, just like the stress-strain relationship of any soil heavily depends on the pore size distribution and the stress history of the sample (Horn and Baumgartl, 2002;Or and Ghezzehei, 2002).The question is, how this structure deformation changes the measured water retention curve.Presumably, those macropores with the lowest mechanical stability are also the pores that drain first.That is, they have released their water before they got deformed.Comparisons between water retention curves obtained with different laboratory methods for various soils with different initial compaction states indicate that centrifugation often results in higher water content for a given suction than a pressure chamber or hanging water column test (Wedler and Boguslawski, 1965;Khanzode et al., 2002;McCartney, 2007;McCartney and Zornberg, 2010).The interpretation whether the offset is still tolerable differs among these authors.The offset emerges because soil compaction leads to a general shift of the pore size distribution, during which the absolute abundance of smaller pores grows on the expense of bigger pores (Leij et al., 2002;Assouline, 2006).However, these pores cannot be captured with X-ray microtomography.An appropriate discussion of this effect is therefore beyond the scope of this paper.
As general advice, centrifugation should not be used to measure water retention curves down to very low pressure ranges if the sample is prone to soil compaction.For practical purposes the reduction in sample height can be used as a suitable indicator to identify the critical pressure beyond which deformation has to be expected.If this method still has to be applied beyond this critical point, it should be performed at the end of all envisaged hydraulic or thermal experiments, as it causes irreversible damages to the internal soil structure.

Methodological limitations
We have developed a workflow for the automatic detection of soil structure deformation by means of free image regis- tration software and outlined best practices in order to optimize the registration results.By imposing a pyramid schedule the image registration of resampled images with roughly 400 3 voxels took approximately 1 min for the Euler transform and 150 min for the B-spline transform on a Linux work station (32 GB RAM, 12 cores with 1.2 GHz).A maximum number of iterations on each pyramid scale is the conventional stopping criterion of the program (in the range of 100-5000), which could have easily been reduced to save time, as the improvement in the last iterations is usually small.We resampled the images from 61 to 183 µm in order to reduce memory use which can reach up to 10 times the image size, depending on data type of both working copies in memory (float vs. short), interpolation type, number of pyramid scales and other internal settings of elastix.As a consequence of resampling, registration cannot be more accurate than the resolution of the final pyramid (183 µm in this case), which is still sufficient as the local deformations in our samples where in the mm range.
The benefits of focusing on the rock matrix during image registration was shown in Fig. 5. Rocks change position and orientation, but in contrast to macropores they do not change their shape.Therefore, they are more easily recovered by the B-spline transform.Taking macropores into account during image registration impaired the spatial alignment of rocks in the lower part of the Köllme soil, because many macropores disappeared completely in that depth which caused wrong feature alignment.This is corroborated by compar-ing the registration results for the full grayscale range in a representative slice (Fig. 7a-c).Fig. 7b is the result of applying the transformation matrix that was obtained from a Bspline transform of the rock matrix to the original grayscale image, whereas Fig. 7c was obtained by computing a Bspline transform with the original grayscale image directly.The pink frame highlights the perfect alignment of a set of rocks if only the rock matrix is considered for image registration, while the position of some rocks was not fully recovered without contrast adjustment.The yellow frame highlights that the disappearance of a macropore at the boundary produces a distortion of the core wall into the soil.This flaw is avoided by focusing on the rock matrix only.A critical question is the following: what is the minimum amount of rock content to guarantee a good recovery of the deformation field?In spite of a very low rock content (0.2 %) in the Bad Lauchstädt soil, the spatial alignment of the rock matrix was satisfactory both with and without contrast adjustment for rocks.This is supported by comparing the registration results in a representative slice (Fig. 7d-f).The position of rocks, e.g., in the pink frame and elsewhere, is correct.In this case taking the full grayscale range into account even improves the spatial alignment of salient features like the macropore in the pink frame.Moreover, the yellow frame highlights poor spatial alignment of the core wall after contrast adjustment because there are no rocks in the vicinity to constrain the registration result.Therefore, focusing on the rock matrix is in fact not advisable for this particular soil.
In summary, a general lack of rocks or other salient, rigid features may render elastic registration useless depending on the severity of deformation.Only additional landmarks can help to improve the registration result in that case.In turn, if the focus is not on natural soil, but on repacked substrates, then the addition of easily trackable features like rocks or metal particles is an easy way to improve the accurate detection of internal deformation.

Conclusions
Measuring the water release through soil centrifugation is a fast method to obtain soil water retention curves.Using Xray microtomography we have corroborated previous findings that the soil structure starts to deteriorate at a capillary pressure of about ψ = −100 kPa.Moreover, quantitative analysis of the pore space at 61 µm resolution revealed that the soil deformation is caused by the interplay of shrinkage and compaction.Local deformation was detected by a novel workflow for digital volume correlation based on elastic image registration.This method enables a detailed look at local soil deformation and its spatial variability.We applied this method to the measure changes in soil structure during centrifugation; however, this method has the potential to quantify the detailed mechanical deformation of soil and other materials exposed to any other type of external forcing.

Figure 1 .
Figure1.Image processing workflow for this study depicted for a small two-dimensional subset (a).Noise is removed with a non-local means filter (b).Shading and cone beam artifacts evoke differences in image intensity which are corrected by subtracting the difference between the depth-dependent gray value average of the soil matrix from the mean gray value of the entire soil matrix (c).Subsequently, image segmentation is performed with multi-level Otsu thresholding (d).For digital image correlation all images are resampled from a voxel size of 61 to 183 µm to reduce memory consumption (e).Optionally, the grayscale range is adapted such that rocks are depicted with optimal contrast and macropores vanish (f).

Figure 5 .
Figure 5.The spatial alignment of rocks between the saturated soil (green) and the soil at ψ = −500 kPa (red) before (a) and after (b) elastic registration image registration of the rock matrix in the Köllme soil.Note that the co-occurrence of rocks results in a composite, yellowish color.The spatial alignment of rocks is less accurate in the bottom of the sample, when the entire grayscale range is used for image registration (c).Same is shown for the Bad Lauchstädt soil(d-f).There, both strategies, with and without contrast enhancement for rocks, lead to equally good spatial alignment of the rock matrix.

Figure 6 .
Figure 6.Displacement vector field for the Köllme soil for the deformation (a) from ψ = 0 kPa to ψ = −100 kPa and (b) from ψ = 0 kPa to ψ = −500 kPa.Only a small percentage of all vectors is displayed to improve visibility.The vector length corresponds to the physical displacement.The two-dimensional section of the undisturbed soil is for orientation.Same is shown for the Bad Lauchstädt soil for the deformation (a) from ψ = 0 kPa to ψ = −50 kPa and (b) from ψ = 0 kPa to ψ = −300 kPa.

Figure 7 .
Figure 7. Representative slice of the Köllme soil close to saturation (a).The B-spline transformation matrix of the rock matrix in the Köllme soil at ψ = −500 kPa is applied to the corresponding grayscale image (b).B-spline transform of the original Köllme soil at ψ = −500 kPa (c).Same comparison for the Bad Lauchstädt soil close to saturation (d) and at ψ = −500 kPa (e, f).Pink and yellow frame highlight salient features.