A QUANTITATIVE EVALUATION OF THE WATER DISTRIBUTION IN A SOIL SAMPLE USING NEUTRON IMAGING

This paper presents an extension of the empirical method by Kang et al. recently proposed for correcting two-dimensional neutron radiography for water quantification in soil. The method was tested on data from neutron imaging of the water infiltration in a soil sample. The raw data were affected by neutron scattering and by beam hardening artefacts. Two strategies for identifying the correction parameters are proposed in this paper. The method has been further developed for the case of three-dimensional neutron tomography. In a related experiment, neutron imaging is used to record ponded-infiltration experiments in two artificial soil samples. Radiograms, i.e., two-dimensional projections of the sample, were acquired during infiltration. A calculation was made of the amount of water and its distribution within the radiograms, in the form of two-dimensional water thickness maps. Tomograms were reconstructed from the corrected and uncorrected water thickness maps to obtain the 3D spatial distribution of the water content within the sample. Without the correction, the beam hardening and the scattering effects overestimated the water content values close to the perimeter of the sample, and at the same time underestimated the values close to the centre of the sample. The total water content of the entire sample was the same in both cases. The empirical correction method presented in this study is a relatively accurate, rapid and simple way to obtain the quantitatively determined water content from two-dimensional and three-dimensional neutron images. However, an independent method for measuring the total water volume in the sample is needed in order to identify the correction parameters.


Introduction
Neutron imaging (NI) is a modern non-destructive, non-invasive material structure visualization technique.The neutron itself was discovered by James Chadwick in 1932 [2], and neutrons have been utilized in many scientific fields.In the geosciences, neutrons have frequently been used for imaging hydrogen-rich fluids in geomaterials and in engineered porous media (as documented e.g. by [3]), in material research (as reviewed e.g. in [4]), in archaeology for investigations of the composition and microstructure of ancient materials (e.g. in [5]), and in industry (e.g. in [6]).Due to the high contrast between water and other soil constituents, neutron imaging is able to map the water content distribution in soils [7].Neutrons are also highly efficient in mapping root structures [6,8].
Neutron imaging is based on the ability of flying neutrons to interact with materials.The degree of neutron beam interaction in a sample depends on its thickness and composition.Light elements (e.g.hydrogen, the main component of water) scatter neu-trons more than some heavier elements (e.g.silicon, the main component of soil).Although scattering is the main type of interaction between the sample and the neutrons, it is convenient to use the concept of the linear attenuation coefficient to evaluate the behaviour of neutrons in imaging experiments.Thus, the neutron imaging method uses the Lambert-Beer law to find the relation between the neutron beam intensity before and after passing through the sample that reaches the detector [9]: where I 0 is the intensity of a neutron beam before passing through the sample, and I is the intensity of the neutron beam after it has passed through the sample, Σ is the attenuation coefficient, and d is the thickness of the material through which the beam passes.
The result of neutron imaging is either a twodimensional image (i.e., a radiogram) or a threedimensional image (i.e., a tomogram), which is ob-tained by reconstruction from a series of radiograms acquired at different angular positions of the sample.
Recently, neutron imaging has been employed to investigate the behaviour of soil water.For example, Lehmann et al. [10] demonstrated the ability of neutrons to investigate the dependence of flow and transport processes in porous media from the geometry of the pores.Snehota et al. [11] reported on how neutrons helped in an investigation of the influence of trapped air redistribution on the quasi-saturated hydraulic conductivity of a composed heterogeneous sample.Heavy water (D 2 O) is sometimes substituted for normal (light) water (H 2 O) in experiments, because it attenuates neutrons approximately seven times less than normal water, making it possible to investigate thicker samples.The final goal when imaging variably saturated porous media is to make precise measurements of the spatio-temporal distribution of the water content.This can be done when a reference image of dry samples is available, but the results may be biased by various imaging artefacts.The effect of beam hardening occurs if the neutron beam is polychromatic, since low-energy neutrons are more attenuated than high-energy neutrons, which results in variations in the measured attenuation.Some neutrons that are scattered at small angles may still strike the detector.As a result, they add an extra image intensity in the vicinity of the pixel that they would have struck if they had continued in the incident trajectory.This extra intensity is falsely interpreted as more neutron transmission, and therefore as less thickness of the material [12].Studies [13,14] have noted the difference between the NI-derived water content and the results of gravimetric measurements of water content.The discrepancies are caused by scattering associated with the large water thicknesses of the saturated column.The bias was suppressed by calculating the relative water saturations from the difference between two consecutive tomograms.
The method of Point Scattered Functions (PScF) calculated by Monte-Carlo simulations was recently proposed for correcting the scattering artefacts in water-containing samples [15], and has been implemented in the Quantitative Neutron Imaging (QNI) tool [16].PScF considers the neutron energy spectrum, the sample material, the thickness of the material layer, the type of detector and the sample-detector distance to simulate the statistical distribution of the neutrons on the detector.The use of PScF implemented in QNI is not straightforward for a user from outside the field of particle physics, since not all combinations of detectors and material data are currently available in QNI [17].The Monte Carlo N-Particle Transport Code System (MCNPX) [18] can prepare the input data for PScF, but it is not generally accessible.A simpler empirical approach [1] for correcting radiograms was recently proposed for quantifying the water content in porous media.
The aim of this paper is to quantify the water con-tent based on the empirical method of radiogram correction.The experiment investigates the water infiltration in two composed soil samples.Specifically, various optimization methods parameters are proposed for identifying the correction.The main goal was to extend the use of the empirical correction method to three-dimensional neutron imaging.

Experiment
Two infiltration experiments were conducted on two soil samples packed in a cylindrical container.The experimental method and the imaging set-up have been described in detail in our previous publications [19,20].
A soil sample composed of fine and coarse fractions of sand (Sample 1) was used in the first experiment, while in the second experiment the sample was made of coarse sand and ceramic (Sample 2).Both samples were packed into quartz glass tubes.For the first sample, the inner diameter of the sample was 34.0 mm and the height was 87.5 mm.For the second sample, the inner diameter was 29.3 mm and the height was 30.7 mm.A ponded infiltration experiment was conducted on each sample by delivering heavy water (a mixture formed by 10 % H 2 O and 90 % D 2 O) with a peristaltic pump to the top of the sample and maintaining a constant water level, so that the water infiltrated under slight positive water pressure.The water was allowed to flow freely by gravity through the perforated disc at the bottom of the sample [21].
The samples were imaged in two sessions with different configurations.Sample 1 was scanned with a field of view (FOV) of 901 × 1451 pixels, producing a pixel size of 0.102 × 0.102 mm.The exposure time for one radiogram was 8 s.The interval between shutter openings was 17 s.Sample 2 was imaged with a FOV of 1300 × 1500 pixels and a corresponding pixel size of 0.045 × 0.045 mm.The acquisition time for one radiogram was 10 s.The interval between shutter openings was 13 s.The actual spatial resolution of the images, as measured by gadolinium-based Siemens star, was better than 0.243 mm for Sample 1 and better than 0.135 mm for Sample 2.
To record the rapid changes in water distribution at the beginning of each infiltration experiment, the series of radiograms was acquired on a motionless sample.Then, after steady state flow was reached, neutron tomography of the sample was performed.The samples were rotated in a 180°angular range to acquire a series of 201 radiograms.Imaging was performed on the NEUTRA beam line [22] of the Paul Scherrer Institute in Switzerland.

Image processing
The raw images were normalized for background noise and for spatial non-homogeneities of the detector, for fluctuations of the neutron flux in time, and for spatial non-homogeneities of the neutron beam.We used the following steps, which are described in detail in [23,24].The radiograms were normalized to open beam (OB) and dark current (DC) images.An OB image was obtained by performing image acquisition without the sample in the beam, while the dark current image (DC) was acquired without the active neutron beam.The effect of the temporal non-homogeneity of the neutron beam on the images was corrected by multiplying each image by rescale factor f , calculated as a reciprocal of the mean intensity value in the region-of-interest (ROI) selected outside the projected sample.Equation (2) details the correction of the raw images: where I n (a.u.) is the normalized image, I Raw , I DC and I OB (a.u.) are the intensities of the raw image, the dark current image and the open beam image, respectively, and f is the rescale factor performing the dose correction.
If the sample contains more than one material, the exponent in (1) can be calculated as the sum of the multiples of the individual attenuation coefficients and the corresponding material thicknesses.In the experiment presented here, the unknown thickness of the water d (cm) and the attenuation coefficient of the water Σ (cm −1 ) were determined for each image pixel representing the sample, according to [25]: where and Σ Wet and Σ Dry (cm −1 ) are the mean attenuation coefficients of the sample in the dry and wet state, and d Wet is the thickness of the wet sample and d Dry (cm) is the thickness the dry sample.
When beam hardening and the neutron scattering effects are neglected, the water thickness in (3) is the unknown variable, and the attenuation coefficient is assumed to be constant.The validity of this assumption was verified in the region of the images that represented the layer of ponding water above the top of the sample, where the thickness of water is known and Σ can be calculated easily.If the beam hardening artefacts cannot be neglected, the empirical formula developed by Kang [1] can be utilized to adjust the attenuation coefficient according to the water thickness, as follows: where Σ (cm −1 ) is the total attenuation coefficient of the water, Σ w (cm −1 ) is the attenuation coefficient of the water, d w (cm) is the thickness of the water, and β (cm −2 ) is the beam hardening correction parameter.The water thickness corrected for beam hardening artefacts is then obtained by substituting (3) and ( 4) into (1): The total water volume in the sample was calculated from each radiogram as the sum of the water volumes below each pixel: where V (cm 3 ) is the total water volume in the image, d w (cm) is the water thickness, a (cm) is the pixel size, D (pixels) is the diameter of the sample, and H (pixels) is the height of the radiogram.

Reconstruction of tomograms
To obtain three-dimensional images, the samples were rotated around their vertical axis between angles 0°a nd 180°with an angular step of 0.9°.One radiogram per step was acquired.The resulting series of 201 radiograms was used to reconstruct one tomogram and subsequently to derive the water content in each voxel.
In this study the tomograms were reconstructed from projections representing the uncorrected water thickness (3) and the corrected water thickness (5).The MuhRec program [26], developed at the Paul Scherrer Institute, was utilized to reconstruct the images using the filtered back-projection algorithm.The result of the reconstruction is then directly the three-dimensional matrix of the water contents.

Results and discussion
The Σ profiles measured for both samples in the cylindrical region of heavy water above the sample (ponding water) are shown in Figure 1a.The attenuation coefficients calculated with use of (1) for known water thickness and corrected I/I 0 are clearly underestimated close to the center of the sample.Figure 1b then shows the relationship between the known water thickness and the attenuation coefficient.
Figure 1b clearly indicates that the attenuation coefficient is linearly dependent on water thickness.The correction method (4) by Kang was therefore employed to calculate the corrected water thicknesses in the sample.
Figure 2 compares the gravimetrically determined volume of water (continuous line) and the volume of water derived from the corrected radiograms based on ( 5) and the parameters obtained from linear fitting of Σ w and β in the ponding water above the sample (blue symbols).The resulting parameters were β = −0.11cm −2 and Σ w = 0.89 cm −1 for Sample 1, and β = −0.10cm −2 and Σ w = 0.86 cm −1 for Sample 2. The agreement was not so good, due to the presence of sand and its influence at scattering.The difference between the gravimetric water volumes and the radiography-determined water volumes is +4.3 % in Sample 1 and −9.8 % in Sample 2 at the end of radiography imaging.In the second step, the coefficients of (5) were determined from the water volume according to (6) by fitting to the gravimetrically determined volume of water  (Figure 2, red symbols).A non-linear method was used, where the value of coefficient β was fixed at a value of 0.105 cm −2 (the average from the previous method), while parameter Σ was optimized by the least squares method in order to minimize the difference between the actual water volume in the sample and the neutron derived water volume.The resulting Σ w was 0.92 cm −1 for Sample 1 and 0.76 cm −1 for Sample 2.
The numerical limitation of (5) that requires Σ w ≥ 4β ln I Dry I Wet (7) was not fulfilled in the region of ponding water above the sample.The amount of water accumulated above the sample was therefore calculated manually from the known water layer thickness (measured from the images) and the cross-section area of the water above the sample.The influence of nonlinear correction on the noise in the radiograms is shown in Table 1.The value of the corrected water thickness in the radiograms and its standard deviation (SD) was calculated for three regions.The first region is outside the sample, the second region contains ceramic only, and the last region contains sand.The values are very similar for a corrected tomogram and for an uncorrected tomogram.
The effect of the correction of the attenuation coefficient on the water content distribution in the to-

Background Ceramic Sand
Mean SD Mean SD Mean SD (10 −5 cm) (10 −2 cm) (10 0 cm) (10 −1 cm) (10 −1 cm) (10  1.The mean water thickness and its standard deviation in three regions: outside the sample, in the region containing ceramic, and in the region containing sand.mogram is demonstrated by the example of Sample 2. The tomogram of the sample in a saturated state was reconstructed from radiograms corrected by a nonlinear method and from uncorrected radiograms.Figure 3 compares the water content distribution reconstructed for Sample 2 from corrected and uncorrected radiograms.Slice 1 is directed through the ceramic disk with a small cylindrical element of coarse quartz sand in the middle, while Slice 2 is directed through the layer of coarse quartz sand only.Figure 3 shows that the non-constant attenuation coefficient causes a significant bias of the reconstructed image.The uncorrected tomogram overestimated the water content values close to the wall of the sample and underestimated the values close the centre of the sample.A better interpretation of the same effect is given in Figure 4, which shows the mean water content in voxels that lie at each radial distance from the axis of rotation in two selected slices that correspond to the slices selected in Figure 3.The water content profile of the ceramic in Slice 1 (Figure 4a) was signif-icantly biased without correction, ranging from 0.32 to 0.43, while after correction it was almost constant at a value of 0.37.The radial water content profile for the sand layer (Figure 4b) shows that the correction did not affect the water content close to the centre of rotation.While the water content at the perimeter and in the profile was significantly equalized it was not constant, as it would be assumed to be for a fully water-saturated sample.The mean water content of the whole sample, 0.34, was the same for a corrected tomogram and for an uncorrected tomogram.In comparison with the water content calculated directly from the radiograms using (6) and from the known total volume of the sample, the reconstruction tends to underestimate the water content by 7 %.In our future work, we will investigate this bias, which originates from the reconstruction process.This paper has presented our work on a correction method that removes the bias of the water content distribution within the sample, rather than the total bias of the water content.

Conclusions
Two ponded-infiltration experiments were monitored by the neutron imaging method.Radiograms were normalized and corrected for the beam hardening effect by a simple empirical method.Tomograms were created from corrected water thickness maps, and the water distribution in three dimensions was calculated.Without the correction, the beam hardening effect overestimated the water content values close to the walls of the sample and underestimated the values close the centre of the sample, but the mean water content was the same.This artefact has to a large extent been corrected.The empirical method used in this study has a mathematical limitation that prevented the correction of the parts of the images containing the thickest layer of water.Linear fitting of the parameters successfully corrected the images from beam hardening and neutron scattering effects, but the total NI determined water volume is different from the gravimetrically determined water volume.The independently determined water volume in the sample is therefore needed in order to identify correction parameters β and Σ correctly.The method presented here is relatively simple and easy to use for corrections of three-dimensional quantitative neutron images of the infiltration process in cylindrical soil samples.

Figure 1 .
Figure 1.(a) The attenuation coefficient profiles across the water layer calculated from radiograms for the two samples under study where the zero is in the middle of the samples, and (b) the measured relationship between the attenuation coefficient of the water and the thickness of the water in the area of the image representing the water layer above the sample.

Figure 2 .
Figure 2. A comparison of the temporal development of the water volume present in the sample calculated from corrected neutron radiograms and the actual volume of water in the sample (including ponding water) measured gravimetrically (a) for Sample 1 and (b) for sample 2. Results are shown for correction parameters obtained by linear fitting of the attenuation coefficients in the ponding (blue symbols) and by nonlinear fitting of the water volume in the sample (red symbols).The stepwise increase in volume in Figure 2b is caused by the dosing of the water.

Figure 3 .
Figure 3.The three-dimensional distribution of the water content obtained from tomograms.Horizontal slices of the uncorrected tomogram show the effect of a seemingly higher water content close to the perimeter, while the corrected image is free of this artefact.Slice 1 contains two materials -coarse quartz sand in the middle and ceramic surrounding the sand.Slice 2 contains just sand.

Figure 4 .
Figure 4.The mean water content of homogeneous material (a) for Slice 1 and (b) for Slice 2, for each radial distance from the centre of the sample.