A Novel Method for Obtaining the Loess Structural Index from Computed Tomography Images: A Case Study from the Lvliang Mountains of the Loess Plateau (China) A

: The structural index is an important quantitative parameter for revealing the structural properties of loess. However, there is no a widely accepted measurement method for structural index at present. This study aims at presenting a novel method for obtaining the loess structural index (LSI), based on the application of computed tomography (CT) scanning techniques and laboratory physico-mechanical tests. The mountainous area of Lvliang in northwest China was taken as the study area, and Late Pleistocene loess samples were taken from various sites in the region. Several physical parameters were ﬁrst measured using laboratory tests, including dry density, pore ratio, and liquidity index. CT scanning was used to observe sample microstructures, and a mathematical relationship was established between CT image parameters and the physical property indices, through three dimensions (3D) reconstruction and slice porosity analysis. The results revealed that LSI can be expressed as a non-linear function related to CT image parameters, dry density, and the liquidity index of the loess. Compared with traditional calculation methods, this novel technique calculates the LSI by using an empirical formula, which is less labor-intensive. Such results indicate that the method warrants wide application in the future.


Introduction
Loess forms through the deposition of aeolian sediments [1] in arid or semi-arid environments, so it normally has unique structural properties [2] and water sensitivity [3,4]. In many cases, such properties can result in negative effects, including extensive surface erosion [5,6] and catastrophic geological hazards [7][8][9]. Covering an estimated total area of 6.4 × 10 5 km 2 , the Loess Plateau has the most complex loess topography in China and, therefore, exhibits great environmental vulnerability to human activities [10,11]. Especially in recent years, frequent and intense disturbance by external factors, including human engineering activities [12,13] and extreme climate change [14,15], have exacerbated this vulnerability. Hence, loess is regarded as one of the problematic soils that requires attention in this region [16,17]. Moreover, knowing the situation and the solutions to one of the areas of the world affected by land degradation and desertification processes is relevant to achieve the sustainable development goals of the United Nations [18,19]. This will help to achieve the land degradation neutrality challenge [20].
Defined as a combination of soil fabric and interparticle bonding [21,22], the soil structure is essentially a kind of physical property. Hence, it is usually associated with the complex mechanical or physical behavior of soil fabric [23], such as pore size [24] and distribution [25] and aggregate arrangement [26,27]. However, such work cannot be used to directly obtain or characterize the structural strength of soil, which reflects the bonding energy between soil aggregates or units [1,23]. To remedy this, laboratory physical-mechanical tests should be performed to analyze the change in the connection between particles or structural units in soil, both before and after its structure is destroyed. However, the strength change of structured soil follows a non-linear law, making the Mohr-Coulomb criterion unsuitable for this type of soil.
Despite some constitutive models for structured soils having been proposed in the literature [28][29][30], they are seldom applicable to loess due to inherent limitations. Hence, some researchers used quantitative parameters to characterize the structural properties of loess. For instance, the traditional method measured soil structural potential parameter (also called soil structural index) by comparing the different strengths of undisturbed, remolded, and saturated soil [31][32][33]. However, various states of soil must be prepared for this test, making it a rather complex and laborious process. Fortunately, the initial structure of soil has a close correlation with basic physical indicators, that is, soil particle size, density, and humidity, respectively [34]. This provides a possible method to address the limitations of complexity. Hence, in the present work, this theory was selected as the basis for analyzing the mathematical relationship between the initial soil structural index and physical indicators.
In addition to the extensive discussions around physical-mechanical laboratory tests, in recent years, high-resolution imaging has been an important tool for studying the structural properties of soil, including digital cameras [35], scanning electron microscopy (SEM) [4], mercury intrusion porosimetry [36], electronic probes [37], magnetic resonance imaging [38], and the most widely used, computed tomography (CT) scanning. Recent experimental results have demonstrated that the CT technique could offer real-time information on the variation of some physical soil parameters occurring in the test specimen, such as the process of damage evolution [39], micro-mechanical properties [40], and porosity [41]. Hence, CT scanning can be used as a bridge between micro-scale and the macro-scale geomaterials [42]. Unfortunately, the application of this technique to determine the soil structural index remains a challenge. However, since the CT image reconstruction can be used to obtain a 3D model, effectively combining laboratory and numerical tests [43], a direct connection between the micro-structural parameters and structural index of soil is possible.
This study examines an area around the Lvliang Mountains (eastern Loess Plateau) in China. Since the loess deposition in the region is a typical structured soil, the main aim of the present study is to obtain the loess structural index (LSI) using instrumental methods. The specific objectives is to establish formula which can calculate the LSI from CT scanning images to reduce the labor-intensive works (e.g., unconfined compressive strength tests) used in traditional methods.

Study Area and Loess Samples
The Lvliang Mountains were chosen as the study area, which are located in the central and northern parts of China (Figure 1a). The region spans the Shanxi and Shaanxi provinces, of which the former contains most of the mountainous area. The total length of the region is more than 400 km with a width of 40~80 km. The mountains are mainly distributed in the western part of Shanxi Province, along a northeast-southwest direction. The elevation ranges from 1000 to 2800 m above sea level, characterized by high terrain in the north and low terrain in the south. The Yellow River is the main river system in the region and lies on the western side of the study area. The climate of the area is a continental monsoon climate. The average annual rainfall is about 400~800 mm, and the rainy season occurs in June~August. Due to the hilly and mountainous location, the climate shows obvious variation with vertical differences. Loess samples were collected from five different counties in the region (Baode, Jingle, Liulin, Jiexiu, and Ji, respectively) ( Figure 1b). The samples were collected using a handheld spade and were extracted as a soil mass. The depths of the samples ranged from 3 to 16.5 m to ensure that all the samples were in the Late Pleistocene soil layer. Three typical geological layers consisting of loess soil were selected: L1 (Malan loess layer) [44], S1 (paleosol layer) [45], and L2 (Lishi loess layer) [46], respectively (Supplementary Materials S1). All these three are from the Chinese Soil Taxonomy (CST) [47]. In order to prevent moisture loss, all the samples were sealed with plastic film and tape after they were taken from the soil layer. A total of 15 undisturbed loess samples were obtained with dimensions of 30 cm × 30 cm × 30 cm.

Laboratory Test
To obtain the soil structural index, at least three types of soils are necessary, i.e., undisturbed, remolded, and saturated soils. The equation for the index is as follows [32]: where (q u ) o , (q u ) r , and (q u ) s are the unconfined compressive strengths of undisturbed soil, saturated soil, and remolded soil, respectively. Meanwhile, m u has a close connection with the comprehensive physical index (I L ρ d )/(ρ w e 0 ) [2] (with a correlation coefficient of up to 0.9 in some cases [34], where I L denotes the liquidity index, ρ d denotes dry density, ρ w is pure water density, and e 0 is soil porosity. Hence, the relationship between the LSI and comprehensive physical index can be fitted if these values are obtained. The physical properties experiments, including the moisture content testing, density testing, and limit moisture content testing, were performed on the undisturbed soil samples to obtain their physical indicators. It should be noted that all the test methods followed the relevant industry standards of China and were performed according to the given steps (Supplementary Materials S2).
The samples selected for the unconfined compressive strength test were prepared with the same size of Φ 39.1 mm × H 80 mm and the same dry density. In addition, the moisture content of the remolded samples should be the same as that of the undisturbed samples. The preparation of remolded samples can be seen in the supplementary material (Supplementary Materials S2). A unconfined pressure tester was used to measure the compressive strengths of the samples with a shear rate of 3 mm/min. When the number of manual turns reaches 18, the equipment will rise by 1 mm. The values of the axial dynamometer were recorded whenever the axial displacement reached 0.1 mm. In order to obtain a more continuous stress-strain curve, the test data were photographed and manual recorded at the same time. The equipment tested three samples with the same state at the same time, and the average values of the unconfined compressive strengths of the three samples were considered as the final result.

Computerized Tomography (CT) Scanning
Computerized tomography scanning, which was first applied in the medical field, has been widely used in recent years to solve geotechnical problems [48]. As a nondestructive characterization method, this technique is helpful for visualizing and quantitatively measuring the internal structure of samples [42]. In this study, CT scanning of undisturbed soil was performed on XTH225ST machine produced by NIKON company (https://www.nikon.com/, accessed on 20 October 2020), in the Industrial CT Scanning Laboratory of China University of Geosciences (Beijing, China). It should be noted that only four samples (samples 1-1, 6-1, 6-2, and 6-3) were selected for CT scanning owing to the limitation of the workload and the complexity of the subsequent analysis process. Because a smaller sample size can ensure a higher image resolution, the samples were prepared as cylinders with a size of Φ 20 mm × H 20 mm. The 225 kV barrel was selected as the X-ray generator, and the parameters in the scanning process were set as follows: (i) 135 kV tube voltage; (ii) 57 µA tube current; (iii) the scanning time of a single photo was 1 s; (iv) the scanning method was translation-rotation (TR) mode; (v) the image resolution was 14 µm. The CT images were obtained in TIFF format and then preprocessed with Gaussian noise-reduction filters and brightness and contrast adjustments (to improve the quality of the final model). CT Pro 3D was then used to reconstruct the images into a 3D model of each soil sample. Finally, the reconstructed 3D model was cut into many slices, each one pixel thick, using the VG Studio Max software. The output images, in Image Stack format, were then analyzed visually.
The 3D model and slices images obtained from CT scanning can clearly be used to characterize the microstructure of a sample, especially the pore size and its distribution [24,25]. Hence, the mathematical relationship between the CT image parameters and the pore index can be established through the analysis of pore volume. According to the traditional theory of soil mechanics, the density and volume of the soil can be calculated using Equations (2) and (3), respectively, as follows: where ρ, m, and V are the natural density, weight, and total volume of the soil, respectively, m s is the weight of the soil particles, W 0 is the natural moisture content, V v is the volume of the pores, and e 0 is the pore ratio. Hence, the natural density of soil can also be expressed as follows: The index V v can also be determined from microstructure images. For the 3D model of a specific soil sample, if the model was cut into q slices and the number of pores in each slice was p, the total volume of pores in the model can be calculated as follows: where A ji is the area of the i-th pore in the j-th slice (Figure 2a), and R j is the thickness of the slice, which is a fixed value in our work (i.e., one pixel). The total volume of the soil (Figure 2b) can be calculated as: where n denotes the porosity of the soil. By setting n j as the areal porosity of the j-th slice, the following expression can be obtained: In addition, the weight of soil particles m s , the total volume V, dry soil density ρ d and soil porosity n have the following relationship: By substituting Equation (7) into Equation (8), the expression can be derived as follows: If Equations (5) and (9) are substituted into Equation (4), the relationship between physical indicator (i.e., the pore ratio e 0 ) and microstructure image parameters can be obtained as follows:

Results of Laboratory Tests
As seen in Figure 3, the samples exhibited typical class I stress-strain curve, that is, with a negative slope in the post-peak phase. During the elastic deformation phase, the internal structures of the samples were continuously adjusted under an external load, along with the initial destruction. With a further increase in the axial stress, micro-cracks began to appear and extended continuously. When the peak strength was reached, the soil was destroyed, and the slope of the curve inversed immediately. The width of the cracks then increased sharply, and the curve gradually decreased to a point representing the residual strength. A typical destruction process of the undisturbed samples can be seen in Figure 4 (taking the sample 1-1 as the example), which consists of four main phases: (i) the initial state of the sample, (ii) the phase of elastic deformation, (iii) the failure of the sample, and (iv) the phase of residual deformation. In terms of the resulting unconfined compressive strength values for the different layers, L1 ranged from 55.56 to 226.39 kPa, whereas S1 and L2 had strengths of 75.46~437.82 kPa and 51.11~458.96 kPa. The largest unconfined compressive strength of all three layers was located in Baode County (i.e., the N1 sampling site). Most of the samples were completely destroyed when the strain was less than 15%. In addition, for samples in the same layer, the rank of the unconfined compressive strength was: undisturbed samples > remolded samples > saturated samples, indicating that the structural property of the loess plays an important role in soil strength. Under the effect of external loading, soil particles will slide and rotate, with microcracks developing and spreading continuously and appearing as cracks and macroscopic sample failure. As seen in Figure 5, according to the shape and direction of the cracks, the failure modes of the samples can be divided into five types: (i) linear splitting, (ii) Y-type splitting, (iii) diagonal shear, (iv) wedge shear, and (v) Y-type shear failure. Category (i), (ii), and (iii) can be seen in the undisturbed samples (Figure 5a). The remolded samples (Figure 5b) exhibited all failure modes except the (ii), whereas the saturated samples exhibited only linear splitting failure and Y-type shear failure (Figure 5c). The statistical results indicated that Y-type shear failure was the main failure mode in the undisturbed and saturated samples, while the remolded samples were dominated by wedge shear failure. The cracks on the failure surfaces of the remolded samples occurred in relatively arrangements and shapes, and the direction of the cracks were similar, with a small number of inflection points. In contrast, the cracks in the undisturbed samples had more extension directions and more inflection points. This is possibly because the preparation process of the remolded samples made their soil structures more homogeneous than that of the undisturbed samples. The undisturbed samples, influenced by the forming mechanism of loess, may have quite large heterogeneity of soil structures, causing the cracks to extend along the weak places first under loading conditions. Compared with the undisturbed samples, all the pores of the saturated samples were filled with water, which enhanced the heterogeneity of the soil structures, so more irregular cracks formed in the saturated samples. The stress-strain curves were analyzed to obtain the unconfined compressive strength of soil samples in different states, and are shown in Table 1 with the physical indicators of the samples.  Subsequently, these values were used to calculate the soil structure index according to Equation (1). The fitting curve between the LSI and the comprehensive physical index (I L ρ d )/(ρ w e 0 ) is shown in Figure 6, and their mathematical connection can be expressed as follows: where ρ w is the density of water, which is generally considered to be 1.0 g/cm 3 . Hence, Equation (11) can be simplified as follows: e 0 ≈ 2.61 I L ρ d ln 2.2979 − ln m u (12) Figure 6. The fitting curve between loess structure index and physical indicators.

Results of CT Scanning
The reconstructed 3D models of the samples were shown in Figure 7. The blue areas indicate pores while red indicates soil particles. It is clear that the pore distribution was rather uneven, which represents that the samples had high heterogeneity. VG Studio Max software was used to analyze the pore properties of the samples, including the number of pores, pore volume, and porosity of the samples. In this process, the minimum probability that the software recognized the pores was set to 5%. The whole sample was selected as the region of interest (ROI). The final results can be seen in Table 2. It should be noted that due to the size uncertainty of the ROI selected by the software, the pore characteristics in the table represent the normalized data. Overall, with an increase in the sampling depth, the 3D porosity of the samples decreased. In addition, the porosity of the L2 layer was less than that of the S1 and L1 layers. The pores were divided into four categories based on diameter (D): (i) micropores (D ≤ 100 µm), (ii) small pores (100 µm ≤ D ≤ 200 µm), (iii) coarse mesopores (200 µm ≤ D ≤ 1000 µm), and (iv) macropores (D ≥ 1000 µm). As seen in Figure 8, the percentage of micropores in the loess of Baode County (sample 1-1, 1-2 and 1-3) was slightly higher than that in other samples. The percentage of coarse mesopores in samples 3-1, 3-2, and 3-3 was much lower than in the other samples. This is considered to be associated with the content of clay particles. The pores formed between these particles were generally smaller than between other types of particles and could be filled by fine substances.   Since the pore ratio e 0 exists in both Equations (10) and (12), the relationship between the LSI and CT image parameters can be obtained by combining the two equation as follows: Based on this, the quantitative calculation formula of LSI can be determined as follows: From this equation, we can see that the loess structural index can be expressed by a non-linear function related to CT image parameters, dry density, and liquidity index of the soil. In particular, all the variables used in this formula are related only to undisturbed soil. Hence, it is not necessary to prepare the remolded and saturated samples and measure their unconfined compressive strengths. This proposed method is, therefore, more convenient and less laborious than traditional methods.

Discussion
The recognition of the pores in each slice depends on the ROI settings. As a type of supervised classification method, the selection of the ROI has a noticeable impact on the effects of classification. Generally, for samples with obvious spatial heterogeneity (like in this study), it is not possible to map the true variation of the pore indices (e.g., areal porosity, total area of the pores in each slice, etc.) without very dense sampling over the overall sample [49]. Even in a strict sense, the error or uncertainty introduced by the heterogeneity of the pores still remains a big challenge in the CT scanning tests, because true homogeneity is almost impossible in nature [50]. However, a large number of sampling operations enhance both the labor cost and the occurrence probability of errors. Fortunately, the application of visual analysis software provides a solution to address this problem. In our study, VG Studio software was used to obtain the pore parameters. The large soil samples were first cut into many slices, after which the pore and soil particles were counted. Because all the pores that were identified, according to the conditions by set the user, were included in the final result, an assessment of pore spatial heterogeneity is not necessary.
For the convenience of the scanning operations and subsequent slice analysis, most studies prepared the soil samples for CT scanning as cylinder shapes [25]. However, it should be noted that all the soil samples in this study were made manually; thus, strictly speaking, they only approximated a cylinder, due to artificial errors. The XTH225ST CT scanning machine determined the resolution of the final images as it adjusts to the size of each sample. Therefore, not all samples returned images with resolutions of 14 µm (this value can be seen in Section 2.1.). Rather, the image resolutions of the four soil samples in this study were 13.71 µm, 13.71 µm, 13.51 µm and 13.51 µm, which subsequently led to different slice thickness in VG Studio Max. Similarly, each sample was divided into a different number of slices. Another impact of the varying image resolutions on the resulting of loess structural index was reflected in the pore statistics. Even though the image resolutions were rather high, pores less than 14 µm in size still not included in the statistics. Rather, such pores were categorized as soil particles, resulting in a smaller areal porosity. The SEM results also support this assumption: the areal porosity obtained through SEM ( Figure 9) were larger than those obtained through CT scanning. An important reason for this difference is that SEM can reach nanometer-scale resolutions, which is higher than that of the CT test. Certainly, the application of the software for digital processing of CT scanning images cannot be ignored. All of the images obtained from CT tests were greyscale images, so pore and soil particle recognition was performed using VG Studio software. Although we have tried to reduce the artificial errors by selecting clearer ROIs and setting parameters (e.g., grey value) to more appropriate values, some pores were still incorrectly classified as soil particles. Hence, due to pores being smaller than the image resolution, limitations of image-processing software, and estimated formulas, the obtained results may have considerable uncertainties, especially for the samples that were not very regular in the shape, and in the ROIs containing many small pores. At the same time, we cannot effectively determine the influence of these errors on the final calculation (whether they overestimate or underestimate the results) because the LSI calculation formula associated with these parameters is a non-linear function. In this context, the standardization of the test process, including sample preparation, ROI selection, and software parameter settings, becomes very important for obtaining more accurate results. Although the method presented here can be used to obtain the soil structural index associated with the microstructure images parameters, the coefficients appearing in the estimated formula were not be proven to be universally applicable. Therefore, the purpose of this study is not to propose a widely applicable model, but to discuss the connection between the LSI and CT scanning images. Hence, a similar analysis could be adopted to determine the corresponding calculation formula applicable to other regions. Certainly, a unified formula at the regional scale is also possible if enough soil samples are prepared and a large number of experiments are carried out, which definitely calls for considerable work and cost. However, it should be noted that a widely applicable formula may only apply in a statistical sense, and more detailed experimental designs will still be necessary to obtain more accurate results. From the perspective of parameters in the formula (Equation (14)), it still requires two physical parameters (the liquidity index and dry density of the loess sample), in addition to the CT image parameters. These physical parameters are mainly related to the physical implications of the soil structure index. The traditional method [32] of calculating the LSI can explain this point: a saturated sample is required in the experiment, illustrating the influence of soil moisture content on LSI. However, as physical parameters, the liquidity index and dry density of soil also have certain morphological implications. If adequate and accurate efforts are made, it may be possible to establish a quantitative relationship between morphological parameters and soil structural index, so that physical parameters are no longer required.
Last but not least, although this study only deals with the loess structure, it would be helpful for finding the proper land management measures in fragile environments such as the loess belt in China. Many studies have identified the soil erosion and soil loss with soil structure, whereas elevated soil loss and runoff rates can reduce soil fertility [51], which subsequently cause a series of problems. Hence, to achieve the United Nations sustainable development goals, it is necessary to explore new measurement techniques of structured soils. Given that many attempts have been conducted in land management, such as the cover crop management [52] and long-term soil erosion monitoring and measurement [53], the proposed procedure in this study would be a useful supplement to improve the soil conditions and avoid soil erosion and land degradation.

Conclusions
This study proposed a novel method for obtaining the loess structural index by using the CT scanning. Based on several laboratory tests and slice porosity analysis of the reconstructed 3D model of loess samples, both the soil structural index and parameters of CT scanning images were linked with the physical indicators of the samples. Hence, a mathematical formula showing the empirical relationship between CT scanning image parameters and the soil structural index in the area could be established. The Late Pleistocene loess in the Lvliang Mountains (China) was used as a case study to clarify the application procedure. The expression of the soil structural index obtained via this method is a non-linear function associated with the CT image parameters, liquidity index and dry density of the loess samples. These results indicate that the soil structural index can be measured and calculated using two simple tests for soil density and limit moisture content, as well as CT scanning. Compared with traditional methods, the time-consuming strength tests on remolded soil and saturated soil can be avoided in the proposed method, so the workload is smaller. Moreover, the method is beneficial for the prevention of geo-hazards associated with loess, and contribute to meeting the land degradation neutrality challenge.