Use of Triple-Exposure X-ray Computed Tomography for Nondestructive Identification of Heavy Elements in Soil Samples

ABSTRACT For the purpose of the nondestructive measurements of hazardous heavy-element-containing substances in soil samples, the emerging technique, triple-exposure X-ray computed tomography (CT), is proposed. The use of triple exposure (i.e. three different acceleration voltages for the X-ray tube) improves the accuracy of the identification of the atomic number of a heavy element. With lead as the target element, virtual X-ray microtomography simulations utilizing a polychromatic X-ray source accelerated at 70, 100, and 140 kV were performed. CT images for the three voltages were reconstructed on a computer for clayey soil samples with heavy-element-containing natural minerals, namely galena (PbS), witherite (BaCO3), stibnite (Sb2S3), sphalerite (ZnS), and magnetite (Fe3O4). Quantitative analysis of the reconstructed CT images was performed to obtain the following results. (i) Galena can be distinguished from other heavy minerals (i.e. witherite, stibnite, sphalerite, and magnetite). (ii) The lead content in the clayey soil can be estimated with reasonable accuracy. These successful results demonstrate that the triple-exposure CT technique is promising for the nondestructive characterization of soil samples contaminated with hazardous heavy elements such as lead.


Introduction
Soil and sediment contaminated with heavy elements are often hazardous (e.g. Golui et al. 2020;Li et al. 2019;Sivalingam, Al Salah, and Poté 2021;Sulaiman, Ghazali, and Ismail 2020;Zhang et al. 2020). Thus, a nondestructive measurement technique is desirable for analyzing such contaminated samples in a laboratory. X-ray computed tomography (CT) is a promising technique for nondestructively identifying heavy elements confined in samples. However, the voxel values of the reconstructed CT images strongly depend on both the molar concentration and the atomic number, Z, of the heavy element in the samples (e.g. Nakashima 2000). Thus, the use of the single-exposure (i.e. a single acceleration voltage for the X-ray tube) X-ray CT method leads to error and ambiguity (e.g. Nakashima and Nakano 2012;Van Geet, Swennen, and Wevers 2000). This problem can be reduced by employing the double-exposure (i.e. two acceleration voltages) X-ray CT method (e.g. Bourque, Carrier, and Bouchard 2014;Champley et al. 2019;Krauss, Schmidt, and Flohr 2011;Nakashima and Nakano 2020;Osipov et al. 2018;Paziresh et al. 2016;Rogasik et al. 1999;Zheng et al. 2020).
The double-exposure or dual-energy method has been a conventional standard CT imagery for element identification. However, the method still has drawbacks when applying it to the identification of heavy elements whose atomic numbers are much larger than approximately 40. The dual-energy ratio (Bourque, Carrier, and Bouchard 2014) and the dual-energy index (Krauss, Schmidt, and Flohr 2011) is no longer a one-to-one onto function of Z for such heavy elements. This makes the inversion of the experimentally obtained ratio and index into a unique value of Z impossible, yielding error and ambiguity. Moreover, the double-exposure method employing the analysis of the cupping artifact caused by beam hardening (Nakashima and Nakano 2020) requires a sample that is large enough to exhibit a strong cupping artifact (e.g. 8 cm in diameter; Nakashima and Nakano 2020). Thus, it is difficult to apply this method to the microtomography of small samples, which exhibit a weak cupping artifact.
To overcome this difficulty, the use of the emerging technology, triple-exposure X-ray CT method (e.g. Niu et al. 2020;Sato et al. 2018;Zhao et al. 2018), is proposed in the present study. This method allows more accurate identification of heavy elements through careful analyses of three CT images acquired at three different acceleration voltages. The general objective of the present study is to demonstrate the effectiveness of triple-exposure CT, when applied to the analysis of soil samples contaminated with heavy elements. The specific objective and procedure are as follows. With lead as the target heavy-element, and the atomic number (i.e. Z = 82) was determined and the lead concentration was quantified for virtual soil samples contaminated with heavy-element-containing natural minerals. Theoretical consideration and numerical simulations of the triple-exposure virtual microtomography of the soil samples were performed for the determination and quantification.

Methods
This section comprises two parts. The triple-exposure X-ray CT of single-phase samples is performed in the first part. The single-phase sample refers to a macroscopically homogeneous sample in terms of the spatial distribution of heavy elements. An analytical expression for the reconstructed CT images is available for the single-phase sample. Thus, the analytical consideration for the single-phase sample is useful for the accurate and quantitative discussion on the principle of the triple-exposure CT. In contrast, real soil samples are often of multi-phase, and comprise several regions with different heavy elements. Thus, in the second part, virtual CT simulations were performed for such multiphase soil samples to quantify the atomic number and concentration of the target element (i.e. lead) embedded in a region.

X-ray CT of single-phase samples
X-ray CT of single-phase samples with or without heavy elements was performed using the analytical expression for the reconstructed CT image of a cylindrical sample. Consider a virtual CT system for a homogeneous cylindrical sample with a radius of R. As the sample rotates, parallel X-ray beams penetrate it. The degree of attenuation of these X-ray beams is detected by a one-dimensional (1-D) array with a separation of δ between detectors ( Figure  1). The X-ray source is assumed to be polychromatic. Its energy-dependent spectra, i 0 (E), are shown in Figure 2, where E is the X-ray photon energy. Examples of the mass attenuation coefficient (MAC) spectrum, which is also energy-dependent, for some elements (as simple substances) are shown in Figure 2.
An axisymmetric profile of the reconstructed voxel value, f(r), where r is the distance from the sample center, can be obtained by inversion using the X-ray attenuation data acquired by the system in Figure 1. The Abel, onion-peeling, and filtered back-projection methods are conventionally employed for the one-dimensional inversion (e.g. Dasch 1992). An analytical expression for f(r) was obtained by the filtered back-projection method assuming the Ramachandran reconstruction filter with infinitesimal δ (Nakano and Nakashima 2018). Based on the exact solution by Nakano and Nakashima (2018), examples of the f(r) distribution for a homogeneous sample are shown in Figure 3, showing the cupping artifact (i.e. an inhomogeneous distribution of the voxel value with respect to r). This is one of the beam hardening artifacts, and is a consequence of the energy dependence of the MAC and X-ray source spectra (e.g. Nakano and Nakashima 2018;Vedula and Munshi 2008).
For the analytical expression of the cupping profile, the quantity C av calculated using Eq. (1) is important.
where E max is the maximum value of the X-ray photon energy and µ(E) is the linear attenuation coefficient (LAC) spectrum of the sample (Nakano and Nakashima 2018). E max is theoretically determined from the acceleration voltage of the X-ray tube. For example, E max is 100 keV for an acceleration voltage of 100 kV. The quantity C av is the linear attenuation coefficient of the sample averaged over the X-ray source spectrum. Figure 3 demonstrates the following aspects of the cupping artifact. (i) f(r) approaches C av as r → R. (ii) The C av value is independent of the sample size, R. (iii) The f(r) distribution becomes almost r-independent as R → 0. Therefore, for a small cylindrical  Figure 1. Configuration of a virtual X-ray CT system for a single cylindrical sample (gray) to obtain the analytical solution for the cupping profile, f(r). The cylindrical sample is homogeneous and the radius is R. The separation between detectors is δ. The r-coordinate is the distance from the sample center. No bowtie filter for beam hardening reduction is attached.
sample, it is reasonable to assume that the reconstructed voxel value is almost constant, and thus C av could be a good estimate of the voxel value. The present study focuses on the quantification of heavy elements in such small samples.
With this in mind, the following systematic calculation was performed. Consider a small homogeneous cylindrical sample, such as that in Figure 1. The sample consists of a simple substance, ranging from 11 Na to 84 Po, and has a bulk density of 1 g/cm 3 . The quantity C av was calculated using Eq. (1) for each simple substance to evaluate the dependence of C av on the atomic number. The photon energy spectra, i 0 (E), for a polychromatic X-ray source accelerated at the three voltages employed in the measurements are shown in Figure 2. They are generated by the program X-Tucker ver. 4, by Kato (2019), which is based on the model by Tucker, Barnes, and Chakraborty (1991). Three acceleration voltages, namely 70, 100, and 140 kV, were assumed for the X-ray tube made of tungsten. The target angle, Al filter thickness, and Cu filter thickness were taken to be 12°, 1 mm, and 0.1 mm, respectively. The LAC spectrum, µ(E), was calculated using the XCOM database by Berger et al. (2010). Effects of the coherent and incoherent scattering and of the photoelectric absorption are considered for the calculation of the LAC spectrum. We consider soil and sediment contamination by naturally occurring heavy minerals at sites near ore deposits (e.g. Hoshino et al. 2020;McGrath, Zhang, and Carton 2004;Noerpel et al. 2020;Obasi 2020;Ohta et al. 2005;Tabelin et al. 2018;Ujiie-Mikoshiba et al. 2006;Vreča, Pirc, and Šajn 2001;Zhang et al. 2012) in the present study. Thus, in addition to the C av values for the simple substance, those for galena (PbS), witherite (BaCO 3 ), stibnite (Sb 2 S 3 ), sphalerite (ZnS), magnetite (Fe 3 O 4 ), native sulfur (S 8 ), and quartz (SiO 2 ) were also calculated. Galena, witherite, and stibnite are hazardous minerals. Galena was chosen because lead is the target heavy element in the present study. Witherite, stibnite, native sulfur, and quartz were chosen because their constituents, namely Ba, Sb, S, and Si, yield C av ratios that are nearly equal to that of Pb, and thus may cause ambiguity. Sphalerite was chosen because it often occurs with galena in ore deposits (e.g. Lambert and Sato 1974). Magnetite was chosen because it is one of the most common heavy-element-containing minerals in the crust and yields bright voxels in reconstructed CT images (e.g. Falvard and Paris 2017; Nakashima and Komatsubara 2018), possibly causing ambiguity.

X-ray CT of multi-phase samples
Computer simulations of microtomography were performed for multi-phase soil samples to check the applicability of the triple-exposure CT method to the soil samples more realistic  Figure 1 for various values of sample radius, R. The source spectrum, i 0 (E), for 100 kV in Figure 2 was employed. The sample is montmorillonite clay with galena, whose bulk density and chemical composition are given in Table 1 (phase A1). The quantity C av calculated using Eq.
(1) is indicated as a broken line, showing that C av is independent of R and is a reasonable estimate of the f(r) value very near the sample rim.
than single-phase samples. A virtual microtomographic system for a small soil sample contaminated with heavy elements was constructed ( Figure 4). The same polychromatic X-ray source spectra as those in Figure 2 were employed. Five values (i.e. 600 nm, 1 µm, 2 µm, 5 µm, or 10 µm) were assumed for the square voxel dimension (i.e. detector separation, δ). Thus, the radius, R, of the small phases, A1 to A4, in Figure 4 was 36, 60, 120, 300, or 600 µm. Although each phase was macroscopically homogeneous, the microscopic structure was heterogeneous, as shown in Figure 5. The matrix (phase B) is a nominally dry soil consisting of stacked montmorillonite clay grains. According to a study on the dry density of sediments (Avnimelech et al. 2001;Hossain, Chen, and Zhang 2015), the bulk density of phase B was taken to be 0.4 g/cm 3 , implying that its porosity was 84.6 vol%. Galena, witherite, stibnite, and native sulfur were inserted into the pore space as phases A1 to A4, respectively, in case I; galena, sphalerite, magnetite, and quartz were inserted as phases A1 to A4, respectively, in case II. The volume fraction of the inserted minerals was 6 vol%. The grain size of all minerals was assumed to be much smaller than the voxel dimension, so that the shape of individual grains did not appear in the reconstructed CT images. The details of the physical and chemical properties of the phases for cases I and II are listed in Tables 1 and 2, respectively. The LAC spectra for each phase were calculated based on Tables 1 and 2 and the XCOM database (Berger et al. 2010). The calculated LAC spectra are shown in the Electronic Supplementary Material (Fig. ESM1). The X-ray attenuation caused by air in Figure 5 was neglected.
Original programs (Ikeda, Nakashima, and Nakano 2019; Nakashima and Nakano 2014) were used for reconstructing 2-D CT images. The details of these programs are described elsewhere (Nakano et al. 2000;Nakashima and Nakano 2012). Briefly, projections were calculated assuming parallel X-ray beams (not fan beams). The number of projections obtained with sample rotation was 805. Metal bowtie filters (e.g. Mail et al. 2009) for the 1-D array of detectors polychromatic X-rays Figure 4. Virtual X-ray microtomographic system, same as that in . Figure 1 but for four small cylindrical phases (A1 to A4) embedded in a large cylindrical one in gray (B), used to perform numerical simulations of CT image reconstruction. The diameters of the small and large phases are 120 and 450 voxels, respectively. The physical and chemical properties of the phases are given in Tables 1 and 2. beam hardening reduction are not assumed. The programs employed a convolution backprojection algorithm (e.g. Medoff et al. 1983;Uesugi et al. 1999) with the Ramachandran reconstruction filter for a finite Nyquist wavenumber of 1/(2δ). The detector separation, δ, was the same as the voxel dimension of the reconstructed CT images. The field of view of the reconstructed image was 512 2 voxels. The volume fraction of the inserted minerals was limited to 6 vol% for the microtomographic simulations mentioned above (see Tables 1 and 2). A lower fraction of heavy minerals would reduce the effective atomic number of phases A1 to A4, resulting in worse accuracy for the identification of the heavy element embedded in the clay. To investigate the effects of the fraction of heavy minerals on voxel values, C av values were calculated using Eq. (1) for various fractions of galena (phase A1), witherite (phase A2), and stibnite (phase A3). The volume fraction of montmorillonite in the phases was kept constant at 15.4 vol% in the calculation.

X-ray CT of single-phase samples
The results of the theoretical calculation for the single-phase samples are summarized in Figure 6, which shows the characteristic dependence of C av on the acceleration voltage and atomic number. The quantity C av decreases with increasing acceleration voltage because the  Figure 5. Schematic microstructure of the macroscopically homogeneous phases in Figure 4. The matrix (phase B) is a stack of montmorillonite clay grains. Its porosity is 84.6 vol%. An amount (6 vol%) of specific contaminant minerals (hexagonal grains in gray) was inserted into the pore space for phases A1 to A4. For example, the contaminant mineral is galena for phase A1 (see also Tables 1 and 2).
LAC spectrum generally (except for the jump at the K absorption edge) decreases with increasing E. As for the dependence on the atomic number, the three curves of C av in Figure  6 monotonically increase for Z approximately less than 50. This monotonical increase allows to calculate the effective atomic number of compounds using the conventional power-low dependence of C av on Z (e.g. Van Geet, Swennen, and Wevers 2000). As for Z approximately larger than 50, however, the quantity C av decreases to minima and increases again. This complex behavior causes the breakdown of the power-low dependence on Z, and the function C av (Z) is no longer a one-to-one onto function. Thus, the conversion of a C av value to a Z value is not unique, rendering the calculation of the effective atomic number of heavy-element-containing compounds such as galena and witherite impossible.
If the bulk density of the simple substance is 2 g/cm 3 , the quantity C av becomes twice that shown in Figure 6. In contrast, the ratios C av (70 kV)/C av (140 kV) and C av (100 kV)/C av (140 kV) are independent of the bulk density, depending only on the atomic number, Z. According to Figure 6, such ratios are sensitive to the atomic number, making them similar to the dual-energy ratio (Bourque, Carrier, and Bouchard 2014), the dual-energy index (Krauss, Schmidt, and Flohr 2011), and the low/high cross-section ratio (Champley et al. . Quantity C av for elements ( 11 Na to 84 Po) as a simple substance with a bulk density of 1 g/cm 3 calculated using Eq. (1), showing the dependence on the atomic number, Z. Three source spectra, i 0 (E), in Figure 2 were employed. "C av (70 kV)" refers to a C av value calculated using i 0 (E) for 70 kV. As indicated by the horizontal solid lines, the values for 82 Pb, 51 Sb, and 16 S are very similar in terms of the ratio C av (70 kV)/C av (140 kV); those for 82 Pb, 56 Ba, and 14 Si are very similar in terms of the ratio C av (100 kV)/C av (140 kV). 2019). The sensitivity to the atomic number is clearly depicted in the cross-plot of the ratios C av (70 kV)/C av (140 kV) and C av (100 kV)/C av (140 kV) in Figure 7. It shows a complicated trajectory due to the complex interaction of i 0 (E) and µ(E) via Eq. (1). Thus, when C av (70 kV)/C av (140 kV) and C av (100 kV)/C av (140 kV) can be obtained in triple-exposure CT experiments, the atomic number of the heavy element in the sample can be identified using Figure 7. The next step after the identification of the atomic number of a heavy element is the estimation of the concentration or bulk density of the element. The C av curves in Figure 6 are for a simple substance with a bulk density of 1 g/cm 3 . Thus, Figure 6 allows us to convert the measured C av value into the bulk density of the heavy element in the sample. This conversion also provides a supplementary chart that is useful for determining the atomic number of the heavy element. For example, the C av values for the simple substance of lead are 14.1540, 9.9069, and 8.1868 1/cm for 70, 100, and 140 kV, respectively, in Figure 6. From the figure, we can calculate the bulk density of each simple substance that yields the same C av values as those of lead (density, 1 g/cm 3 ). The result is plotted in Fig. ESM2, which is similar to Figure 6 in the study of Nakashima and Nakano (2020). In terms of the selfconsistency of the method, all three estimated bulk density values should be within the measurement error if the candidate element is the ground truth. Thus, in addition to its use in bulk density estimation, Fig. ESM2 can be used for checking the self-consistency of the triple-exposure CT method.  Tables 1 and 2 are superimposed. The two solid lines that include Pb in Figure 6 are re-plotted as broken lines. Cav ( Figure 9. Examples of the numerically reconstructed CT image for the sample arrangement in Figure 4 and for the chemical composition given in Table 1 (case I). The source spectra, i 0 (E), employed for (a), (b), and (c) are 70, 100, and 140 kV, respectively. The square voxel length (i.e. δ) is 600 nm, and the image dimensions are 512 2 voxels = 307.2 2 µm 2 . The same CT image as that in (c) with a different gray scale is shown in Fig. ESM6(a) to clearly show the dark phases A4 and B. The cupping profile sampled along the horizontal red line in (b) is shown in Figure 10.
Some minerals (not simple substances) are also plotted in Figure 7. Galena (PbS), witherite (BaCO 3 ), and stibnite (Sb 2 S 3 ) are located at almost the same positions as those of 82 Pb, 56 Ba, and 51 Sb, respectively. This is a consequence of the contribution of the heaviest element in the compound (e.g. Pb in PbS) dominating those of the lighter elements (e.g. S in PbS) (Nakashima and Nakano 2020). In contrast, there is a significant discrepancy between the locations of sphalerite (ZnS), magnetite (Fe 3 O 4 ), and quartz (SiO 2 ) and those of the corresponding simple substance (Zn, Fe, and Si). This is because the atomic numbers of their constituent elements are almost the same. A detailed explanation is given for Si and SiO 2 in Fig. ESM3. Figure 7 shows that the behaviors of Pb and PbS are almost the same in triple-exposure CT, which would simplify the data interpretation because the contributions of light elements such as sulfur can be neglected.
The effects of the fraction of galena (phase A1), witherite (phase A2), and stibnite (phase A3) on the C av (70 kV)/C av (140 kV) vs. C av (100 kV)/C av (140 kV) plot are shown in Fig.  ESM4. A more dense plot is shown in Figure 8 for the galena-montmorillonite mixture. Both figures demonstrate that a significant shift of data points occurs with increasing montmorillonite fraction, making triple-exposure CT inaccurate.  Table 1 is shown as a broken line.

X-ray CT of multi-phase samples
Examples of the numerically reconstructed microtomographic images for a voxel dimension of 600 nm are shown in Figure 9, ESM5, and ESM6. The LAC value decreases with E except for the contribution of the K absorption edge (Fig. ESM1). Thus, the brightness of the reconstructed images in Figure 9 and ESM5 decreases with increasing acceleration voltage. Examples of the profile, f(r), for three reconstructed images are shown in Figure 10. Each profile shows a cupping artifact (see also Fig. ESM7). Unlike the exact solution in Figure 3, however, the profiles for the reconstructed images are not smooth but fluctuate, probably due to (i) the discretization of the cylindrical sample shape using a finite voxel size (Goertzen, Beekman, and Cherry 2002;Nakashima and Nakano 2020) and (ii) the interaction among phases A1 to A4 and B. The voxel values are also modulated near the sample rim (i.e. r ≈ R) due to the use of a reconstruction filter with a finite Nyquist wavenumber (Nakano and Nakashima 2018). This fluctuation/modulation could reduce the accuracy of simple substance (theoretical) mineral without clay (theoretical) mineral with clay (theoretical) mineral with clay (simulation, 600 nm) mineral with clay (simulation, 1 µm) mineral with clay (simulation, 2 µm) mineral with clay (simulation, 5 µm) mineral with clay (simulation, 10 µm) Cav(70 kV)/Cav(140 kV) and

µ70/µ140
Cav(100 kV)/Cav(140 kV) and µ100/µ140 84Po 11Na Figure 11. Same chart as that in Figure 7. Results of the numerical CT simulations for phases A1 to A4 with montmorillonite clay (Figure 4)  the calculation of C av , that is, the f(r) value at r ≈ R. Thus, as a reasonable substitute for C av , a new quantity was employed for the CT image analysis, namely a voxel value averaged over a circular region of interest (diameter: 116 voxels) within circular phases A1 to A4 (diameter: 120 voxels), denoted as µ 70 , µ 100 , and µ 140 for 70, 100, and 140 kV, respectively. The quantities µ 70 , µ 100 , and or µ 140 were calculated for each simulated CT image. The resultant ratios µ 70 /µ 140 and µ 100 /µ 140 were plotted in Figure 11. Magnifications of Figure 11 are shown in Figs. ESM8 and ESM9, depicting the dependence of the ratios on voxel dimension. Figure 11, ESM8, and ESM9 are the most important results for tripleexposure CT because they allow the identification of the atomic number of a heavy element, as described in the next section.

Discussion
The applicability of triple-exposure CT to the quantification of heavy elements in soil is discussed using the virtual microtomographic results described in the previous section. The problem to be solved is as follows. Consider the multi-phase clayey soil sample in Figure 4 with four regions, A1 to A4. The four regions are contaminated with unknown minerals. The candidates of the minerals are galena, witherite, stibnite, sphalerite, magnetite, native sulfur, and quartz. With microtomographic images acquired at three acceleration voltages (70, 100, and 140 kV), we should determine which region is contaminated with lead and estimate the lead content in the region.

Double-exposure X-ray CT
Before the triple-exposure CT data analysis, the double-exposure CT data analysis was performed for a pedagogical purpose. Either C av (70 kV)/C av (140 kV) or C av (100 kV)/C av (140 kV) curve in Figure 6 would be provided by double-exposure or dual-energy CT experiments. Each curve is independent of the bulk density, depending only on the atomic number, Z. Thus, the atomic number of the heavy element in the sample may be identified using one of the curves in Figure 6. Once the atomic number is identified, it is straightforward to estimate the density or molar concentration of the heavy element using Figure 6 and ESM2. However, the double-exposure CT method has the following drawback. For example, as indicated by the horizontal solid line in Figure 6, the C av (70 kV)/C av (140 kV) value for

Triple-exposure X-ray CT
This technical difficulty described above can be solved by employing triple-exposure CT using 70, 100, and 140 kV as follows. The cross-plot of the ratios C av (70 kV)/C av (140 kV) and C av (100 kV)/C av (140 kV), is shown in Figure 7 for the microtomography of singlephase samples (Figure 1). The theoretically calculated cross-plot shows a complicated trajectory due to the complex interaction of i 0 (E) and µ(E) via Eq. (1). Fortunately, however, the candidates ( 51 Sb, 16 S, 56 Ba, and 14 Si) that cause ambiguity are located significantly far from the ground truth, 82 Pb, in Figure 7, ensuring their elimination. This is the basis for the heavy element identification using triple-exposure X-ray CT technique.
Results of the CT simulations of multi-phase samples (Figure 4) are shown in Figure 11. Figure 11 and its magnifications in Figs. ESM8 and ESM9 are useful for determining the lead-containing region from among phases A1 to A4. Although slight effects of the voxel dimension can be seen, the data points for µ 70 /µ 140 and µ 100 /µ 140 with respect to phase A1 are located at almost the same positions as those for C av (70 kV)/C av (140 kV) and C av (100 kV)/C av (140 kV) with respect to lead as a simple substance and pure galena mineral. These data points for galena are located far from those for other candidates (e.g. witherite and stibnite). Thus, it can be concluded that triple-exposure CT successfully identified phase A1 as the galena-containing region. It should be noted that the discrimination of Pb from other candidates, such as Ba, Sb, S, and Si, was demonstrated, a task at which conventional double-exposure CT would fail.  Figure 12. Estimate of the lead content in phase A1 for the CT image in Figure 9 (voxel dimension = 600 nm). This chart is essentially the same plot as that in Fig. ESM2, but the bulk density of each simple substance was re-calculated to reproduce the µ 70 , µ 100 , and µ 140 values obtained from the image analysis of Figure 9. The ground truth (0.395 g/cm 3 ) of the lead content as a simple substance in phase A1 is shown by dotted lines.
Once the atomic number of a heavy element is identified as that of lead, the next step is the determination of the lead content in phase A1. This could be done using the procedure described with respect to Fig. ESM2. According to the image analysis of the CT images in Figure 9 and ESM5, µ 70 , µ 100 , and µ 140 were 5.7769, 4.0535, and 3.3578 1/cm, respectively, for the lead-containing phase A1. Approximating that C av (70 kV) ≈ µ 70 , C av (100 kV) ≈ µ 100 , and C av (140 kV) ≈ µ 140 , the element content in the soil that reproduces these obtained values of µ 70 , µ 100 , and µ 140 was calculated using the C av curves in Figure 6. The results are shown in Figure 12 for a voxel dimension of 600 nm.
According to Figure 12, the lead content as a simple substance in the soil that reproduces the µ 70 , µ 100 , and µ 140 values obtained from the analysis of the CT images was 0.408 g/cm 3 for the 70 kV data, 0.409 g/cm 3 for the 100 kV data, and 0.410 g/cm 3 for the 140 kV data. These three bulk density values are very similar, confirming the self-consistency of the method. Although the chart is omitted, the lead bulk density for the microtomography conducted with a voxel dimension of 10 µm was estimated to be 0.354 g/cm 3 for the 70 kV data, 0.349 g/cm 3 for the 100 kV data, and 0.358 g/cm 3 for the 140 kV data, again confirming self-consistency. According to Table 1, the galena content in phase A1 is 53.3 wt% and the bulk density of the galena-clay mixture is 0.856 g/cm 3 , implying that the ground truth of the lead content as a simple substance is 0.395 g/cm 3 . This ground truth value agrees well with the values estimated in the CT image analysis based on the microtomographic simulations noted above, demonstrating the reasonable performance of the triple-exposure CT technique.
The undesirable effects of an increase in the clay fraction are shown in Figure 8 and ESM4. As described, the relatively large galena content of 6 vol% in phase A1 ensured the successful identification of the atomic number and content of lead. However, the data points for the galena-montmorillonite mixture shift to the position of pure montmorillonite in Fig. ESM4, suggesting that the detection of lead in the lead-poor clay mixture is difficult. Figure 8 also shows that the contribution of lead to the C av value is almost negligible for the soil with a bulk lead content of 150 mg/kg (environmental quality standard for soil pollution in Japan; Ministry of the Environment, Government of Japan 2003). Thus, it would be difficult to apply triple-exposure CT to soils with low lead concentration. However, this difficulty can be avoided by performing high-resolution microtomography of galenacontaining samples (e.g. Jussiani and Appoloni 2015). The voxel dimension was assumed to be much larger than the mineral grain size in the present study ( Figure 5). It can be reduced to be as small as the galena grain size in high-resolution X-ray microtomography. The local lead content increases for small voxels containing large galena grains, which ensures the detection of lead by triple-exposure CT.
It should be noted that the theoretical C av values calculated using Eq. (1) are only valid for the sample with a single phase (Figure 1). If the sample is a composite consisting of several phases (e.g. Figure 4), the spectrum of the incident X-ray beam on phase A1 is not exactly equal to that in Figure 2 due to the beam hardening effect by surrounding phases A2 to A4 and B. This yielded a slight but systematic underestimate of the voxel values obtained in the tomographic simulations compared with those obtained from theoretical calculations ( Figure 10). Fortunately, the successful determination of lead in clay samples demonstrates that this undesirable composite effect was almost negligible for the case considered in the present study. However, Figures 7 and 11 may need to be remade based on the sample structure when analyzing composite samples with a more complicated structure than that in Figure 4. Figures 7 and 11 may also need to be remade for large samples. The present study focuses on the quantification of heavy elements in samples as small as R = 36, 60, 120, 300, and 600 µm for the phases, A1 to A4, in Figure 4. It should be noted that the approximation that C av (70 kV) ≈ µ 70 , C av (100 kV) ≈ µ 100 , and C av (140 kV) ≈ µ 140 is only valid for such small samples (see Figure 10 for 100 kV). A cupping artifact more severe than that in Figure 10 and Fig. ESM7 would occur for samples much larger than those employed in the present study. Such a severe cupping artifact would reduce the accuracy of the above approximation, rendering the discussion based on Figures 7 and 11 less reliable. The possible solution for the large-sample problem is to remake Figures 7 and 11 based on the large sample geometry. The method analyzing voxel values at the sample rim and center (Nakashima and Nakano 2020) would be helpful because it is applicable to large-sample CT images with serious cupping artifacts.

Conclusions
The use of triple-exposure X-ray CT was presented for the nondestructive measurements of hazardous heavy-element-containing substances in soil samples. The employment of three voltages for the X-ray tube acceleration improves the accuracy of the identification of the atomic number of a heavy element. Computer simulations of X-ray microtomography with a polychromatic X-ray source accelerated at 70, 100 and 140 kV were performed. Microtomographic images for the three voltages were reconstructed for clayey soil samples in which mineral grains (galena, witherite, stibnite, sphalerite, magnetite, native sulfur, and quartz) were embedded. An analysis of the CT images utilizing the ratios of the reconstructed voxel values for the three acceleration voltages (i.e. µ 70 /µ 140 and µ 100 /µ 140 ) enabled us to distinguish lead-containing galena from other minerals, which is difficult to achieve using conventional double-exposure X-ray CT. The estimate of the lead content in the clayey soil was also achieved with a reasonable accuracy. The results show that triple-exposure CT is promising for the nondestructive characterization of soil and sediment samples contaminated with hazardous heavy elements such as lead.