Combination of Imaging Infrared Spectroscopy and X-ray Computed Microtomography for the Investigation of Bio- and Physicochemical Processes in Structured Soils

Soil is a heterogeneous mixture of various organic and inorganic parent materials. Major soil functions are driven by their quality, quantity and spatial arrangement, resulting in soil structure. Physical protection of organic matter (OM) in this soil structure is considered as a vital mechanism for stabilizing organic carbon turnover, an important soil function in times of climate change. Herein, we present a technique for the correlative analysis of 2D imaging visible light near-infrared spectroscopy and 3D X-ray computed microtomography (μCT) to investigate the interplay of biogeochemical properties and soil structure in undisturbed soil samples. Samples from the same substrate but different soil management and depth (no-tilled topsoil, tilled topsoil and subsoil) were compared in order to evaluate this method in a diversely structured soil. Imaging spectroscopy is generally used to qualitatively and quantitatively identify OM with high spatial resolution, whereas 3D X-ray μCT provides high-resolution information on pore characteristics. The unique combination of these techniques revealed that, in undisturbed samples, OM can be found mainly at greater distances from macropores and close to biopores. However, alterations were observed because of disturbances by tillage. The correlative application of imaging infrared spectroscopic and X-ray μCT analysis provided new insights into the biochemical processes affected by soil structural changes.


INTRODUCTION
Soils are heterogeneous and complex mixtures of organic matter (OM), mineral particles and pore space. Their inherent functionality emerges from the spatial arrangement of these constituents on various spatial scales (Ritz et al., 2004;Portell et al., 2018;Wanzek et al., 2018). Important soil physical properties and processes such as conductivity, water retention, gas exchange and root penetration are defined by soil structure, i.e., the spatial arrangement of solids and voids (Bronick and Lal, 2005;Rabot et al., 2018). The interplay of soil chemical and biological properties through physical protection has been considered a key factor in stabilizing soil carbon (Dungait et al., 2012;Mueller et al., 2012;Stockmann et al., 2013;Wiesmeier et al., 2019). This, in turn, is influenced by soil structure, which defines the microenvironmental conditions for microbial carbon turnover. OM is not evenly distributed in soils. For instance, the walls of biopores have been often reported as OM rich compared with bulk soil (Banfield et al., 2017;Hoang et al., 2017;Hobley et al., 2018). It is assumed that the heterogeneous location of OM in soil may be a significant factor controlling the extent of OM mineralization (Dungait et al., 2012;Steffens et al., 2017). Therefore, to understand the relationship between soil structure, carbon dynamics and biogeochemical processes, it is crucial to characterize undisturbed soil samples. However, so far this relationship is often investigated by techniques based on soil disturbance, e.g., by extracting aggregate fractions (Young et al., 2008) from the bulk soil, which indeed provides insights into microscale processes, but the complete image of the soil architecture is lost .
Visible light near-infrared spectroscopy (Vis-NIR) is an established method for the simultaneous analysis of various soil properties (Stenberg et al., 2010;Soriano-Disla et al., 2014). However, soil samples have to be ground and homogenized for optimum performance. Imaging Vis-NIR (imVNIR) is an emerging technique that allows the spatial analyses of intact soil samples (Buddenbaum and Steffens, 2011;Steffens and Buddenbaum, 2013). Steffens and Buddenbaum (2013), Steffens et al. (2014) and Hobley et al. (2018) demonstrated the potential of imVNIR by improving our understanding in soil classification, chemical composition and carbon storage on the micrometer scale for whole pedons.
The combinatory application of 2D and 3D imaging of chemical composition and soil structure, respectively, may provide new insights for the physical properties and the spatial extent of soil structures along with their chemical information. X-ray computed microtomography (X-ray µCT) is nowadays considered as a commonly used technique for describing soil structure with physical parameters such as pore size distribution and connectivity (Rabot et al., 2018). As 2D imVNIR and 3D X-ray µCT give different but corroborating information on the same material, the combination of the two techniques can therefore be used to uncover and quantify biogeochemical processes in intact soil samples. The combination of different biogeochemical imaging methods, designated as correlative microscopy or correlative imaging, is increasingly used in life sciences (Caplan et al., 2011;Handschuh et al., 2013;Guyader et al., 2018). A few recent applications have demonstrated that there is a great potential for the application of correlative imaging in soil sciences (Hapca et al., 2011;Juyal et al., 2019;Kravchenko et al., 2019;Schlüter et al., 2019). Schlüter et al. (2019) highlighted that a combination of different 2D techniques for chemical and microbial imaging, such as scanning electron microscopeenergy dispersive using X-ray (SEM-EDX) analysis and nanoscale secondary ion mass spectrometry (nanoSIMS) or fluorescence microscopy and X-ray µCT, can be used to reveal specific soil microenvironments and thus biogeochemical and physical processes in structured soil. Compared with these 2D techniques, the spatial resolution of Vis-NIR imaging spectroscopy is lower (several µm instead of nm) but allows the investigation of larger samples up to a complete pedon. Therefore, adding imVNIR to the toolbox of correlative analyses of 2D imaging techniques and 3D X-ray µCT for fast imaging of several samples may reveal patterns of soil biochemical processes at the pedon level of intact soils.
In this study, we developed a procedure for the correlative image analysis of 2D imVNIR and 3D X-ray CT on intact soil cores. Our approach consisted of X-ray µCT scanning, followed by impregnation with resin, slicing, imVNIR and classification of the same undisturbed soil core sample. After image registration, the classified imVNIR images and the X-ray µCT images were analyzed with correlative image analysis, e.g., by analyzing the OM distribution as a function of distance to pores and biopores. We used cores sampled from a spacefor-time reclaimed chronosequence. These displayed both structural and chemical differences within the same substrate Pihlap et al., 2019) and thus provided perfect conditions for validating the correlation analysis method on undisturbed field samples.

Sampling and Sample Preparation
Intact soil cores were obtained from reclaimed soils in the Garzweiler open pit mine area, 40 km west of Cologne (Germany). A detailed site characterization together with a description of the reclamation procedure and crop rotation of the space-for-time chronosequence is described in the study of Pihlap et al. (2019). Briefly, an approx. 20-meter-thick layer of loess, including approx. 2.2 m of developed Luvisol, was excavated during coal mining, mixed to afford a Luvisol/loess ratio of 1:10 and deposited. This mixture was used for reclamation with alfalfa (Medicago sativa) in a pioneering phase for 3 years. Between the 4th and 6th year, the reclamation was based on a crop rotation including Triticum aestivum L. (wheat) and Hordeum vulgare L. (barley). In the years after, a variation of T. aestivum L. (winter wheat and summer wheat), H. vulgare L. (winter barley), Brassica napus L. (rapeseed) and Zea mays L. (maize) were grown. About 30 t/ha of organic fertilizer (compost or manure) were added in the 4th and 7th year of the reclamation procedure.
Three cores with different soil structures were selected on the basis of the study of Lucas et al. (2019) for correlative analysis with imVNIR, representing different soil development levels (3, 6, and 12 years after reclamation) and management stages (no tillage/tillage) ( Table 1). The selected undisturbed soil cores (height 20 cm and diameter 10 cm) were obtained from fields with different ages, in depths of 0-20 and 40-60 cm, using a custom-made drill for undisturbed sampling of cylindrical soil cores (UGT GmbH, Germany). To overcome the trade-off between resolution and sample size, Macroporosity and biopore length density are shown for the exact same sample , whereas bulk density, organic carbon (OC) concentration and pH values represent a mean value from the plot where the sample was taken (Pihlap et al., 2019). three subsamples of 3 cm diameter and 3 cm height were collected per core. The no-tilled topsoil is an immature soil with low OM content and dense soil structure that was not affected by soil tillage. After 6 years of reclamation, the tilled topsoil was characterized by high macroporosity due to plowing and high OM content owing to the organic fertilizer amendment. The subsoil represents a sample with many macropores due to high biological activity, namely, high biopore density, which developed for 12 years under a tilled topsoil.

X-ray µCT Scanning at 19 µm Resolution
The cylindrical subsamples with a diameter of 3 cm were scanned using an X-ray microtomograph (X-TEk XCT 225, Nikon 162 Metrology). Reconstruction was achieved with a spatial resolution of 19 µm, and the images were processed and segmented into solids and pores, as described in detail by Lucas et al. (2019). Specifically, all biopores were segmented in the binary images using the Tubeness plugin in Fiji to separate all tubular objects. All biopores, including filled biopores with roots or earthworm cast, were separated from the remaining, irregularly shaped, pore network in the binary images. Biopores in the resulting images were therefore empty, and the rough pore wall was excluded. These biopores were dilated five times in all 3D directions (kernel diameter of 3 voxels), and then the image of the segmented pores was subtracted in Fiji. This allowed us to include the walls of the biopores and solids, like earthworm cast within the biopores, to the resulting image.
imVNIR Scanning at 53 µm Resolution After X-ray µCT analysis, the same intact soil cylinders were dehydrated gradually with acetone [series from 30-100% (v/v)], and subsequently, the soil cores were impregnated with polyester resin (PALATAL P 6-01, BÜFA, Germany). After 5 weeks of polymerization, the cores were placed in an oven at 40 • C for 48 h to cure the surface of the impregnated samples. All cores were cut in two slices of 1 cm thickness, polished and scanned with a hyperspectral camera (VNIR-1800, Norsk Elektro Optikk Ås, Norway). The camera lens distance from the sample surface was set to approx. 30 cm, and the sample was illuminated with two light sources in front of and behind the camera (angle about 45 • ). An image of 53 × 53 µm 2 /pixel (1800 pixel/line) and spectral range with 196 bands between 400 and 990 nm was thus obtained. A certified reflectance standard (reflectance 50%) was adopted and used to normalize the reflectance on the target image in order to adjust the differences in illumination: where L obj is the radiance of the recorded sample, L ref the radiance of the certified reflectance standard and ρ ref the reflectance of the certified reflectance standard (Peddle et al., 2001;Steffens and Buddenbaum, 2013) Hyperspectral images were processed in ENVI Classic (Version 5.2, Exelis Visual Information Solutions, Boulder, Colorado, United States) following the workflow in Figure 1. In particular, we used a principal component analysis (PCA) to concentrate the spectral information in images and remove the correlation between neighboring bands (Rodarmel and Shan, 2002;Kavzoglu et al., 2018). PCA images show spectral spatial variability of different components on the soil surface and helps us to define number of endmembers (classes) (Steffens et al., 2014). Each endmember represents pure spectra from a selected component (class) and it is a basis for determining the composition of each pixel with linear spectral unmixing (LSU). The first 10 principle components were used to separate polyester resin and pores filled with polyester from the soil matrix and to identify pure spectral endmembers of observed components for the subsequent image processing steps. The main criteria for the endmember collection were to select the most abundant features (from PCA analyze) and to cover pixels from different spatial segments in the image. According to the PCA analyze in total we selected six different components as endmembers, such as loess parent material, OM particles, Feoxide (FeOx) dominated pixels and iron-manganese concretions. For the tilled topsoil and subsoil, two additional endmembers were selected because of their specific features. The tilled topsoil was affected by plowing, which created voids on top of the loess. As a polyester resin is a transparent material, in tilled topsoil we captured also a reflectance "noise" from the loess that was underneath the shallow tillage void. This gave different reflectance than the voids created by roots and earthworms extending vertically through the entire sample. Thus, additional endmembers for tillage voids were necessary. An additional endmember was also selected from the subsoil to describe fresh, OM-rich pixels, which displayed the most defined OM spectrum compared with that of other OMdominated areas. The statistical significance of the selected endmembers was tested with the Jeffries-Matusita distance test (Dabboor et al., 2014). For detecting and mapping the FIGURE 1 | Workflow of imVNIR image processing. After X-ray µCT analysis, the intact soil cylinders impregnated with polyester resin were cut in two slices of 1 cm thickness. These slices were scanned with a hyperspectral camera. In the first image processing step, scanned images were used to compute a principal component analysis (PCA), which concentrates the spectral information on images and removes correlations between neighboring bands. Principle components were used to separate polyester resin, pores filled with polyester from the soil matrix, and to identify pure spectral endmembers of observed components. In the example of the workflow, we selected five different classes as endmembers (ROI-s) and they were used to calculate relative abundance in each pixel with linear spectral unmixing (LSU). With LSU, abundance maps of each endmember is created, where the most intensive color refers to the highest contribution of selected endmember to the pixel. Obtained abundance maps were subsequently used for a supervised classification using maximum likelihood, and the area coverage of each class can be calculated.
abundance of each spectral endmember, LSU was conducted using the manually selected endmembers (Supplementary Figure 1). LSU generates abundance maps, where the relative contribution of the selected endmembers to each pixel is depicted (Ravel et al., 2018). The obtained abundance maps (Supplementary Figures 2-4) were subsequently used for a supervised classification using maximum likelihood classifier (Borra et al., 2019). For the supervised classification, new additional ROIs were selected along with the Jeffries-Matusita distance evaluation. The subsequent area coverage of the classes was calculated in the classified images. Maximum likelihood classification is created by results of the abundancy maps (LSU), of which we know every class contribution to every pixel in 2D space. This means that for example a pixel classified as Fe is dominated by iron, but may also contain fractions of other endmembers. The accuracy and the kappa coefficient of the classification were tested with a confusion matrix processed in ENVI Classic software. For the confusion matrix we selected random pixels that were not part of the ROIs used for the classification itself and compared them with the classification results. The kappa coefficient is a basis to estimate the agreement of classification and it is computed in the ENVI Classic software as follows (L3Harris Geospatial Solutions, 2020): where i is the number of classes, N is the total number of classified values compared to truth values, m i,i is the number of values that belongs to the truth class i and have also been classified as class i, C i is the total number of predicted values belonging to class i, G i is the total number of truth values belonging to class i.

Image Registration
The images obtained from X-ray µCT and imVNIR had different resolution and orientation. To overcome these constraints, an image registration was needed prior to correlative analysis. The aim of the registration is to find a geometric transformation to map the pixels of the first (moving) image into the second (target) image. Therefore, 2D/3D image registration was performed using the elastix software (Klein et al., 2010;Shamonin et al., 2013) along with a protocol similar to that described by Schlüter et al. (2019). In the current study, the moving image was always the filtered X-ray CT image and the imVNIR image was the target image. Both images were used as gray scale images, however, the information transported by a certain gray scale value differs with the source. We also used a combination of two metrics, the mean of the Euclidean distance between the corresponding landmark points and the mutual information criterion (Mattes et al., 2001). The mutual information is a statistical measure of the degree of dependences of two random variables A and B. It is closely related to the entropy, i.e., the amount of uncertainty, of the two variables and hence also to the joint entropy. More precisely, this means it quantifies the "amount of information" given by the random variable A about the variable B and vice versa (Maes et al., 1997). Here, it defines the objective function that has to be optimized and such allows the entropy to be matched with the corresponding image pairs For example, homogenous features, such as macropores in µCT images, should be matched with homogenous features visible in imVNIR (Figure 2). A second metric, which uses landmarks (at least three manually selected points), was also added. This ensured that the optimization of the moving image did not result in local minima. A pyramid schedule was used to achieve a fast registration, starting with a coarse resolution and moving toward the next finer scale. The best results were obtained by using a similarity transform, defined firstly by the objective function of the two metrics. Afterwards, the resulting parameters were used as the initial transform for affine transformation, for which only one metric was used to define the objective function. A first rough registration was thus performed and was mainly guided by the landmark points, and a second step, using only the mutual . For this, the two metrics Euclidean distance and mutual information were used for optimization. The Euclidean distance metric aims to minimize the distance between the manually set points (red crosses). In this example, it is minimized from a mean of 237 to 4 pixel. The second metric represents the extent of mutual information between the registered µCT image and the imVNIR image, i.e., it measures the "amount of information" given by the µCT image about the imVNIR and vice versa. The scatter plots of gray values from imVNIR and CT (B) at the beginning of the registration reveals two clusters, which are formed by the µCT image (pores and matrix) but are not shared by the corresponding imVNIR pixel values. This results in a low mutual information of 0.19. After the final iteration, pore and matrix clusters were formed by both variables, which is also reflected by the isolines showing high kernel densities. Thus at the end of the registration mutual information is high (0.86). information metric, resulted in a more accurate result that did not rely on manually set points. Additionally, the similarity transformation involved rotation, translation and scaling and therefore had 7 degrees of freedom in 3D. Instead, the affine transformation added aspect ratio and shear to the allowed transformation, resulting in 12 degrees of freedom in 3D (Klein et al., 2010). Corresponding elastix script and example images can be found in the Supplementary Material.

Correlative Image Analysis
X-ray µCT analysis allowed the characterization and the localization of pores >38 µm (2 pixel), whereas imVNIR enabled the determination of the soil's chemical composition. From imVNIR we derived OM and FeOx abundance maps obtained with LSU, whereas from CT images we used the binary (pore) images and the images of Euclidean distances to pores and biopores. The Euclidean distance transformation (EDT) in Fiji was computed to obtain the straight-line distance of each pixel to the next pore/biopore in 3D. In order to compare the data from the two sources (correlative image analyses) the CT images were transformed with transformix, a subroutine of elastix, using the transformation parameters of the original CT image to assure that the images match the corresponding imVNIR image (see section "Image Registration"). Thereafter, multichannel images of the Euclidean distances, pores and results of LSU were created in Fiji, from which we could calculate the relative contribution of OM-rich pixels (obtained with LSU) with distance to pores or biopores (obtained with EDT). In addition, the solid and pore volume was also measured with distance to biopores. This made possible to take into account density shifts around biopores by normalizing the relative contribution of OM-rich pixels by the solid volume fraction. Only distances with a certain number of pixels (100), equal to an area of about 0.3 cm 2 , were taken into account, to avoid high standard deviation for smaller quantities. It should be noted that correlative imaging does not necessarily include correlation in a statistical sense. Figure 3 displays the maximum likelihood classification images for the three soil cores and the respective coverage of the five endmembers. The maximum likelihood classification was conducted through relative abundance of each component obtained from LSU. The overall classification accuracy was 92.0% (k = 0.89) in the no-tilled topsoil, 91.6% (k = 0.90) in the tilled topsoil and 88.3% (k = 0.85) in the subsoil. Our technique enabled the clear separation of voids and soil matrix in all samples. The soil matrix was characterized irregular distribution of loess and Luvisol that were mixed during the soil reclamation procedure. Images recorded at later reclamation stages (tilled topsoil and subsoil images) also showed soil management effects on the soil chemical structural composition. Specifically, the no-tilled topsoil was characterized by a dense soil matrix with a higher area coverage of loess soil and spherical aggregates, as e.g., shown in the areas E5,F, L1-2, and M3 in Figure 3. In two tilled topsoil slices, we observed an increase in OM content and aggregate shape, due to the organic fertilizer amendment (Figure 3, B,C11 and B12, K11-12, L,M10, and L,M11). The area coverage of OM in the slices of the subsoil was comparable with that of the tilled topsoil, although the total organic carbon of the topsoil was rather low (Table 1). Although Fe/Mn concretions could only be found in the no-tilled topsoil, all samples contained FeOx-rich areas. The most distinguishable was slice one from the no-tilled topsoil sample, which had exceptionally high percentages of FeOx-dominated areas.

Image Registration
Only 4-5 common points were sufficient to align the µCT images with the imVNIR images by a similarity transformation. This initial registration was suitable to ensure a good result of the affine transformation, which relied only on mutual information. Figure 4 depicts all 2D slice cuts corresponding to original µCT images, which were mapped according to the imVNIR images. The similarity could be distinguished especially in large pores, such as the central biopore in the subsoil sample or in large cracks. An example for a large crack can be found in A-E14,15 and G-J14,15, whereas a large biopore is visible in N-P15,16 and S-U15,16 (Figure 4). These characteristics were less pronounced in other denser or more mixed samples. However, other features such as aggregates also indicated correct registrations. For example, in the tilled topsoil sample (e.g., Figure 4, areas G10,11 and U9), the identified aggregates were rich in OM and therefore had lower reflectance in the imVNIR images. The corresponding regions in the µCT image were also darker because of the low electron density of the OM compared with that of the other solid substances. However, in some areas, slight differences could be observed because of the pre-treatment of the intact soil cores during resin impregnation. It should be noted that some pores in the imVNIR images are not black. In shallow pores, solid material from deeper layers was visible through the transparent resin (Figure 4, areas I9, H10, I10, H12, and T9), hence the reflectance was affected by the signal from the material below.

Correlative Analysis of µCT and imVNIR Detectable Features
The EDT was computed individually for segmented pores and biopores, and the corresponding images show the Euclidean distance from any non-pore/non-biopore voxel to the next pore/biopore. Figure 5 illustrates the results of the correlative analysis by observing the OM intensity distribution over the distance to the macropores and biopores. Correlative analysis with macropore distances in the no-tilled sample demonstrated a small increase in OM content with increasing distance. In contrast, the tilled topsoil sample revealed a decreasing trend of OM content with increasing distance to macropores. Furthermore, in the tilled topsoil sample, there were no soil matrix pixels more than 0.5 mm away from the next macropore. The subsoil sample (slice 1) had a relatively small abundance of OM-rich pixels near the pores, which then increased to about FIGURE 3 | Classification images and relative coverage of each class. Classification was computed on images of linear spectral unmixing (LSU). LSU images contain information about each endmember contribution to each pixel in the 2D space. When we classify images according to their abundancy, classes are calculated based on the highest domination of selected component. It means that for example a pixel classified as Fe is dominated by iron, but may also contain fractions of other endmembers.
100% at a distance of 0.75 mm. This relatively high distance was observed only in the subsoil sample.
The density of biopores in samples was much lower compared with that of the total amount of macropores, leading to much longer distances (>4 mm). Compared with the distance to pores, the OM trends around the biopores were less pronounced (Figure 5). In both slices of the tilled topsoil sample and in one slice of the subsoil, a small increase in OM was observed with increasing distance. Moreover, the biopores (0 mm distance) of the tilled soil samples exhibited low OM content. Figure 6A illustrates the macroporosity distribution relative to the distance to biopores. In all samples, macroporosity decreased with increasing distance (up to roughly 0.2 mm), but this gradient was predominant in the subsoil (Figure 6A). The relative abundance of OM in the vicinity of the subsoil biopores increased after normalizing the local OM contribution by the corresponding fraction of dense soil material, i.e., taking the shift in macroporosity into account ( Figure 6B). Figure 5) did not significantly change in any sample with the distance to pores and biopores. Interestingly, the FeOx abundance next to biopores increased in the no-tilled topsoil sample. The deviation of the FeOx distribution directly in the pores from zero can be attributed to the high amount of pores, which are smaller than the spatial resolution of the hyperspectral camera (53 µm). In this case, smaller pores are not detectable in the imVNIR images, but they are visible in µCT (19 µm resolution), resulting in a mixed spectrum of pores and FeOx-rich solid volume. This can explain the strong signal of FeOx in the LSU images at locations that are classified as pores in the µCT image instead.

DISCUSSION
ImVNIR is a non-destructive analysis to characterize physicochemical properties in soil. However, when slicing field moist samples, several artifacts such as smearing soil can FIGURE 4 | Image pairs of registered 2D slices of µCT (gray) and imVNIR images (visible spectra). Gray values of µCT represent dense areas (white) and pores (black). The pores in the imVNIR images can be black or light colored, since the resin used is transparent.
be created, thus destroying the macropore system and creating surface topography that leads to shadows. In our study, we successfully applied, for the first time, resin impregnation with polyester followed by cutting in combination with subsequent imVNIR. The developed impregnation procedure allowed us to sustain the integrity of the soil microstructure and avoid the destruction of the macropore system (Rasa et al., 2012;Jangorzo et al., 2014;Gutiérrez Castorena et al., 2016;Mueller et al., 2017;Schlüter et al., 2019). In rare instances, we observed small changes occurring especially in looser areas, due to sample preparation artifacts. In the acetone dehydration pre-treatment, a large pore collapsed in the tilled topsoil (Figure 4, areas C11 and H,I11) during the impregnation, whereas some cracks that were visible in the no-tilled topsoil (Figure 4, areas O6 and T6) diminished. Although there are different dehydration methods available (e.g., drying in the oven), they can create other specific artifacts, such as crack formation, that could complicate the registration and interpretation (Jongerius and Heintzberger, 1975). Despite the occurrence of small changes, the currently developed sample pre-treatment method allowed the horizontal cut of the soil in larger slices without smearing (Figure 4), and thus, a good registration of µCT and imVNIR images was achieved (Figure 4).
The impregnated soil slices provided a perfect basis for the segmentation of different compositional and structural features, such as OM-and FeOx-dominated areas (Figure 3), with imVNIR. The slopes in the spectrum (Supplementary Figure 1) of selected endmembers were in line with data reported by Steffens and Buddenbaum (2013), as their selected regions, such as Mn-and FeOx-dominated areas in a loess-derived stagnic Luvisol, exhibited comparable spectra. As soil cores originated from the same substrate, imVNIR image processing revealed that the chemical and structural features were distributed differently within the samples due to the soil development FIGURE 5 | Area based relative contribution of OM in linear spectral unmixing images with the 3D Euclidean distance to pores (upper graphs) and biopores (bottom graphs). The pictures on the right display the same cutout of macropores and biopores (gray) and the corresponding distances (small distance = purple, high distance = orange).
after reclamation. In the no-tilled topsoil and subsoil, we observed features that were influenced by the reclamation process. Moreover, µCT revealed structural pores inside the OM-rich aggregates of these samples. These aggregates, known as "rolled aggregates", are typically formed on the conveyor belt during transportation of the reclamation soil in the area and can be rich in OM due to the admixture of old topsoil (Luvisol) (Pihlap et al., 2019). In contrast, the sample from the tilled topsoil was characterized by a tillage-induced structure and contained distinctively shaped, OM-rich aggregates (Figure 3, areas B,C 11,12 and K11,12). Segmentation and analysis of the soil composition in undisturbed soil cores defined their 2D spatial distribution and their heterogeneity in 53 µm resolution from samples of 3 cm diameter.
Correlative imaging allowed the enhancement of the information obtained by imVNIR with data provided by µCT. Hence, we were able to correlate the compositional information of the solid phase with the soil structure, such as the OM distribution as a function of the 3D Euclidean distance to pores and biopores ( Figure 5). Nowadays, besides imVNIR, several techniques are available for chemical mapping of 2D slices, such as SEM-EDX, near edge X-ray absorption fine structure spectroscopy, nanoSIMS, or matrix-assisted laser desorption/ionization . These techniques could be also applied for correlative analyses to obtain additional chemical information on the 3D pore structure derived from the µCT, as reported in the study of Schlüter et al. (2019) for SEM-EDX and nanoSIMS analysis. Most of these techniques can even have a resolution on the nm scale and thus exceed the resolution of both imVNIR and µCT. However, high-resolution results in small sample sizes, thus limiting these techniques to microaggregate analyses . In contrast, the combination of imVNIR and µCT can be applied to scan soil cores up to 10 cm width and infinite length. Thus the linkage between soil structure and soil composition can be investigated in larger soil volumes up to a complete pedon. In addition, relative fast imaging (µCT ∼30 min, imVNIR ∼5min), the fast and parallelized elasitx registration method (∼5 min) and the possibility to embed all samples at the same time make it possible to use the method on multiple soil cores in a full field study, even if samples are sliced several times (every 0.5 cm).
The methods employed in the present correlative imaging had different spatial resolutions. The resolution of the 3 cm scanned samples in X-ray µCT images was higher (19 µm) than that of the Vis-NIR images (53 µm), thus making only the pores with >38 µm (2 pixel) detectable, whereas local differences within the micropores (diameter <38 µm) could not be addressed. Alternatively, sub-resolution features from greyscale µCT images could have been taken into account by using the gray values around macropores to normalize the OM contribution , instead of using only visible pore and solid fraction (Figure 6). However, the normalization we used indicated a closer linkage between the soil composition and the soil structural features, visible in µCT, i.e., macropores >38 µm. Potentially, an even broader spectrum of pore sizes could be investigated using subsamples for µCT images and combining the information of pore size distribution from different sample sizes (Vogel et al., 2010). In the current study, we demonstrated that registering higher-resolution µCT images and combining them with imVNIR images expands the information level on the sub-resolution features in larger samples.
Correlative imaging with imVNIR and µCT revealed patterns in FeOx and OM distribution through the distance to pores and biopores (Figure 5). The different trends in OM, as a function of distance to macropores, demonstrated that OM accumulated with increasing distance to macropores, most probably due to physical protection, as the absence of macropores may lead to a deficient aeration but may also hinder microorganisms of reaching OM (Dungait et al., 2012). The comparison between the tilled topsoil and subsoil samples revealed how tillage affects the OM distribution, as shown by previously observed homogenization in OM distribution (Yang and Wander, 1999;Kay and VandenBygaart, 2002). In contrast to our expectations, we did not identify higher OM quantities in biopore walls; in some samples even lower intensities were detected (Figure 5), perhaps due to local differences in macroporosity that were higher close to biopores (Figure 6). When the OM contribution was normalized according to the density shift from µCT, the OM in the no-tilled topsoil and subsoil was higher next to the walls of biopores compared with that in bulk soil (Figure 6).
Consequently, based on the OM quantity per gram soil, the OM content would be higher close to pore walls. This shows the possibility and strength of the new method to explain spatial OM distribution obtained by imVNIR image analysis and additionally explain possible contradictory results in OM distribution by incorporating the shifts in density revealed by µCT. Considering that the impact of soil management on carbon stabilization has been particularly investigated (Follett, 2001;Chivenge et al., 2007;Cárcer et al., 2019), correlative analysis with imVNIR and X-ray µCT may broaden our understanding of OM distribution, associated with physical protection due to its location in the porous soil.

CONCLUSION
In this study, we present an approach for the identification and correlation of soil chemical composition and pore structure in intact soil cores. To that end, we developed a novel combination of 2D imVNIR and 3D µCT through registration and correlative analysis. As a proof of principle, we analyzed three intact soil samples, each with two slices that differed in soil structure and OM content (no-tilled topsoil, tilled topsoil and subsoil). Correlative analysis along with 3D soil structural information about pores can give new insights into the linkage between the biopore network and the soil chemical composition including OM distribution. The analysis indicated that, when the soil is not disturbed by tillage, OM accumulates and is more protected as the distance to macropores increases. µCT data gave valuable input not only by locating pores and biopores but also through local shifts in bulk density. As both techniques give complementary and partly corroborating information on the same material, the development of the applied method improves our understanding on the spatial distribution of soil components and their correlation to soil structure.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS STATEMENT
All authors listed have made a substantial, direct and intellectual contribution to the work and approved it for publication. ML conducted image registration and correlative imaging; EP conducted imVNIR image analysis.

ACKNOWLEDGMENTS
We gratefully acknowledge Manuel Endenich (Department of Reclamation, RWE Power AG) for the access to the study fields and his support during sampling campaigns. We are thankful to Carsten W. Mueller and Gertraud Harrington for their valuable advice and assistance with sample preparation and to Steffen Schlüter for his advice on elastix software. We are grateful to reviewers CS and HZ, as their comments helped to improve the manuscript.