Vis-NIR spectroscopic assessment of soil aggregate stability and aggregate size distribution in the Belgian Loam Belt

Vis-NIR spectroscopic assessment of soil aggregate stability and aggregate size distribution in the Belgian Loam Belt

Geoderma 357 (2020) 113958 Contents lists available at ScienceDirect Geoderma journal homepage: www.elsevier.com/locate/geoderma Vis-NIR spectrosco...

2MB Sizes 0 Downloads 25 Views

Geoderma 357 (2020) 113958

Contents lists available at ScienceDirect

Geoderma journal homepage: www.elsevier.com/locate/geoderma

Vis-NIR spectroscopic assessment of soil aggregate stability and aggregate size distribution in the Belgian Loam Belt

T



Pu Shia,b, , Fabio Castaldia, Bas van Wesemaela, Kristof Van Oosta a b

Georges Lemaître Centre for Earth and Climate Research, Earth and Life Institute, Universite Catholique de Louvain, 1348 Louvain-la-Neuve, Belgium College of Earth Sciences, Jilin University, Changchun 130061, China

A R T I C LE I N FO

A B S T R A C T

Handling Editor: Alex McBratney

Soil aggregate stability (AS) reflects a soil's resistance to external erosive forces and is an indicator that varies with changing elementary soil properties across space and time. However, the quantification of AS via conventional wet-sieving is too resource-demanding a task to be carried out at large scales. We explored the possibility of using laboratory Visible-Near infrared (Vis-NIR) spectroscopy to predict aggregate mean weight diameter (MWD), and three aggregate size fractions (i.e. clay + silt, microaggregate (63–250 μm) and macroaggregate (> 250 μm)) resulting from wet-sieving. Two spectra-based approaches, one (SPF approach) that built direct linkage between soil spectra and four AS indexes via partial least squares regression and the other (SPF + PTF approach) that established pedotransfer functions based on spectroscopically predicted elementary soil properties, were developed on a total of 83 topsoil samples collected in the Belgian Loam Belt. These two approaches were later compared to a third approach (PTF approach) that built pedotransfer functions using measured elementary soil properties. Results show that, with the PTF approach consistently giving the best performance, all three approaches produced good prediction models (R2: 0.62–0.85; RPD: 1.59–2.46) for MWD, microaggregate and macroaggregate fractions, while no correct model was developed for the clay + silt fraction using the two spectra-based approaches, due to the rather homogenous soil texture in the study area. Soil organic carbon was a major factor that controlled the variations in AS, and a critical threshold of 2% SOC content was found to be the level that separated the “stable” and “unstable” AS classes. This study demonstrated that Vis-NIR spectroscopy is a promising technique that enables large-scale assessment of AS and aggregate size fractions, especially considering that future space-borne hyperspectral imagers will provide high-resolution Vis-NIR spectral data at unprecedented scales.

Keywords: Soil aggregate stability Vis-NIR spectroscopy Aggregate size distribution Belgian Loam Belt

1. Introduction Soil aggregate stability (AS) affects a wide array of physical and biogeochemical processes in agricultural environments and reflects a soil's ability to retain its structure against external forces (Amezketa, 1999). During erosion events, AS determines the susceptibility of soil aggregates to the kinetic impact of raindrops and the shear stress of surface runoff. This in turn largely controls the rates of erosion through its influences on the availability of fragmented materials and the level of surface crusting (Le Bissonnais, 1996; Legout et al., 2005). Also, the rainfall- and/or runoff-induced aggregate breakdown products, that are

redistributed in the course of soil erosion, contain large amounts of carbon, nutrients and pollutants due to their small size classes and high specific sorption capacity (Shi and Schulin, 2018), thereby making aggregate stability an influencing factor of lateral element fluxes. As such, AS is frequently used to represent the soil erodibility parameter in previous attempts to model sediment pathways as well as source and sink areas of eroded materials (De Roo et al., 1996; Gumiere et al., 2009); aggregate size distribution has also been used as an input to simulate the multi-class sediment transport (Van Oost et al., 2004). However, the spatial predictions of such models are usually associated with large uncertainties, due, at least partly, to the inadequate

Abbreviations: AS, aggregate stability; SOC, soil organic carbon; TN, total nitrogen; MWD, mean weight diameter; SWC, soil water content; Vis-NIR, Visible-Near infrared; MIR, mid-infrared; PTF, pedotransfer function; MLR, multiple linear regression; VIF, variance inflation factor; AIC, Akaike Information Criterion; CV, coefficient of variation; PLSR, partial least squares regression; RMSE, root mean squared error; RPD, ratio of performance to deviation; RPIQ, ratio of performance to interquartile range. ⁎ Corresponding author at: Mercator B345, Place Louis Pasteur 3, Georges Lemaître Centre for Earth and Climate Research, Earth and Life Institute, Universite Catholique de Louvain, B-1348 Louvain-la-Neuve, Belgium. E-mail address: [email protected] (P. Shi). https://doi.org/10.1016/j.geoderma.2019.113958 Received 8 May 2019; Received in revised form 26 June 2019; Accepted 2 September 2019 0016-7061/ © 2019 Elsevier B.V. All rights reserved.

Geoderma 357 (2020) 113958

P. Shi, et al.

Fig. 1. Spatial distribution of the sampling points in the Belgian Loam Belt.

et al., 2016). Previous attempts have successfully predicted SOC, total nitrogen, texture, pH, soil water content, and some metal concentrations using multivariate statistics and ever-expanding soil spectral libraries (Soriano-Disla et al., 2014). Furthermore, secondary soil physical properties that are controlled by these soil variables can be estimated as well, but with various degree of success. Using Vis-NIR spectra, Askari et al. (2015) could predict the often visually-assessed soil structural quality reasonably well, due to its close relations to SOC and aggregate size distribution. Canasveras et al. (2010) combined VisNIR with mid-infrared (MIR) spectral data to assess AS in southern Spain and pointed out the potential of soil spectroscopy for distinguishing different AS classes. Within a similar Mediterranean environment, Gomez et al. (2013) used laboratory Vis-NIR spectra and achieved a relatively modest model to predict the AS index, as represented by the aggregate mean weight diameter (MWD). They concluded that a better model performance could be expected if the lack of correlation between MWD and elementary soil properties such as SOC could be improved by expanding the range of variation in those soil properties. As of yet, successful applications of soil Vis-NIR spectroscopy to predict soil AS and aggregate size distribution in other environments remain to be scarce. The Belgian Loam Belt is an important agricultural region with 65% of the land cover area being cropland (Evrard et al., 2008). Numerous

representation of the spatial variability of AS and aggregate size distribution (Jetten et al., 1999; Shi et al., 2017). Subtle variations in soil organic carbon (SOC) content, clay content, soil water conditions, and tillage operations can cause substantial changes in aggregate stability (Amezketa, 1999). Numerous studies have shown the large spatial variability of AS due to different land use/ cover and topographical characteristics (Annabi et al., 2017) and in response to land management practices (Dimoyiannis, 2009). This highlights the necessity to consider AS as a spatially dynamic variable that evolves with seasons and agricultural practices, rather than an inherently constant property for a given soil type. However, conventional methods (e.g. Le Bissonnais (1996)’s method) used to determine AS usually involve wet sieving to separate fragment size classes of distinct properties. The resources needed to carry out such laborious measurements mean that repeated sampling and analysis at relevant scales is challenging. Efficient and affordable methods are thus needed to enable the quantification of AS across large spatial and temporal scales. Visible-Near infrared (Vis-NIR) soil reflectance spectroscopy has been widely used in recent years to resolve the trade-off between large scale soil information and high costs (Nocita et al., 2015) and as an alternative to wet chemistry building on recent improvements in accuracy and standardised protocols (Ben Dor et al., 2015; Viscarra Rossel 2

Geoderma 357 (2020) 113958

P. Shi, et al.

measurement, soil samples were placed in Petri dishes and their spectral reflectance were collected over the 350–2500 nm wavelength range, with a spectral resolution of 3 nm in the 350–1000 nm range and 10 nm for the 1000–2500 nm range. The output data were sampled to a 1 nm resolution. Reflectance values from 350 to 399 nm and from 2451 and 2500 nm were removed prior to any processing because these spectral ranges had low signal to noise ratios.

studies have reported the seriousness of soil erosion in this area (Beuselinck et al., 2000; Vandaele and Poesen, 1995), and a study by Meersmans et al. (2011) suggested that, in the cropland of the Loam Belt, the SOC content suffered an average 15–30% decrease from 1960 to 2006, which could lead to a substantial decline in AS, thus increasing the risk of surface crusting and soil erosion. It is therefore of particular importance to develop a method that enables fast and efficient investigation of AS status in this area. To this end, we explored the possibility of using Vis-NIR reflectance spectra to predict aggregate stability (as represented by the aggregate MWD) and aggregate size distribution. More specifically, we used two approaches (i.e. direct vs. indirect) to predict MWD and different aggregate size fractions using soil spectra. The direct approach is to develop a spectrotransfer model to directly predict the AS indexes using laboratory Vis-NIR spectral data, while the indirect approach is to spectroscopically predict elementary soil properties (i.e. SOC, soil pH, and soil texture), which were then used to predict the AS indexes via multiple linear regression (MLR). Lastly, the performances of the two approaches were statistically compared to a third approach, which is an independent MLR model developed based on the measured elementary soil properties.

2.4. Aggregate stability analysis Aggregate stability was measured following the fast-wetting test proposed by Le Bissonnais (1996). In brief, approximately 7 g of 3–5 mm aggregates were firstly immersed in deionized water for 10 min, before being transferred onto a 63-μm sieve for the initial separation of fine and coarse fragments through wet-sieving. The wetsieving was carried out in 95% ethanol in order to preserve the aggregate structure during sieving and to prevent re-aggregation upon drying. The coarse fragments (> 63 μm) were oven-dried for 48 h at 40 °C, and then dry-sieved through a column of sieves (125 μm, 250 μm, 500 μm, 1 mm and 2 mm) and weighed to give the mass percentage for each size fraction. For brevity, the aggregate size fractions were reaggregated to form three aggregate size fractions, i.e. < 63 μm, 63–250 μm and > 250 μm, which are referred to as clay + silt, microaggregate, and macroaggregate fractions respectively. Moreover, aggregate mean weight diameter (MWD), which is the sum of the mass percentage of each size class multiplied by the mean diameter of the respective size class, was used to represent aggregate stability (AS). Hence, four indexes, namely MWD, clay + silt, microaggregate and macroaggregate, were used as AS indexes for spectroscopic predictions. In addition, three AS classes were adopted from Le Bissonnais (1996) according to the MWD values: highly unstable (< 0.4 mm), unstable (0.4–0.8 mm) and stable (> 0.8 mm), as shown also in the normalized international standard (ISO/FDIS, 10930, 2012).

2. Materials and methods 2.1. Soil sampling and preparation The study area is located in the central part of Belgium, where the main soil type is the Niveo-eolian loess-derived Luvisol (IUSS Working Group WRB, 2014). The region is dominated by cropland with a rolling topography, and the main crops are winter wheat, winter barley, sugar beet, maize and potatoes. The climate in this region is temperate oceanic with mean temperatures between 2.3 °C (January) and 17.8 °C (July) and a mean annual precipitation of 790 mm (Castaldi et al., 2019). The sampling points were randomly selected within a 9.7 kmwide and 40 km-long area (SW corner: 50.60 N, 4.65E; NE corner: 50.70 N, 5.06E) from Gembloux to Lincent, and the spatial distribution of the sampling points are shown in Fig. 1. A total of 83 soil samples were taken from the upper 10 cm of the bare soil surface in October 2018. Within 12 h after sample collection, initial soil water content (SWC) was determined for each sample gravimetrically by oven-drying an aliquot (ca. 100 g) at 105 °C. Then, all the remaining fresh samples were air-dried and divided into two subsets. The first subset was passed through a 2-mm sieve before basic soil analysis and Vis-NIR spectral measurements, while the second subset was pre-sieved using 3-mm and 5-mm sieves to separate the aggregates in the 3–5 mm fraction for aggregate stability analysis.

2.5. Methodological framework 2.5.1. Modelling approaches Three approaches, based on the work of Gomez et al. (2013), were adopted to develop multivariate statistical models that predict MWD and mass percentages of clay + silt, microaggregate and macroaggregate fractions. – The first approach, later referred to as the PTF approach, derived pedotransfer functions to estimate the four indexes from the measured elementary soil properties via MLR. – The second approach, called the SPF approach, built spectrotransfer functions to directly link the soil Vis-NIR spectra with the AS indexes. Partial least-squares regression (PLSR) was used to develop SPF models. – The third approach, called the SPF + PTF approach, combined the first two approaches by firstly predicting elementary soil properties using spectral data and then establishing pedotransfer functions based on the predicted soil properties.

2.2. Basic soil analysis The analysis of elementary soil properties included SOC, total nitrogen (TN), soil pH, and soil texture. SOC and TN were analyzed by means of dry combustion using a VarioMax CN Analyzer (Elementar GmbH, Germany). Samples containing inorganic C (6% of all the samples) were pre-treated before combustion by adding 10% HCl to remove inorganic C. Soil pH (H2O) was measured at 1:2.5 soil/water ratio by a pH meter (PHS-3E, Leici, China). Soil texture was measured by means of laser diffraction (LS13320, Beckman Coulter, USA) after removing organic matter by adding 35% H2O2 to each sample.

For each of the three predictive modelling approaches, the 83 samples were randomly divided into calibration and validation datasets at a 3:1 ratio. Then, predictive models were developed using the calibration dataset and evaluated against the independent validation dataset. To assess the model robustness and to enable statistical comparisons among three approaches, such random data subsetting plus subsequent model development was replicated 100 times for each modelling approach and the validation accuracy was quantified each time.

2.3. Soil Vis-NIR reflectance measurement Vis-NIR reflectance spectra of sieved (< 2 mm) soil samples were obtained with an ASD Fieldspec 3 FR spectroradiometer (Analytical Spectral Devices Inc., USA). The measurements were carried out using an ASD contact probe, equipped with a built-in light source (100 W halogen reflectorized lamp). All measurements were performed in a dark room to avoid disturbance by external light sources. During the

2.5.2. Multivariate regression Multiple linear regression (MLR) was used to build pedotransfer 3

Geoderma 357 (2020) 113958

P. Shi, et al.

functions, which predict AS indexes by accounting for the role of elementary soil properties in controlling AS. In the PTF approach, SWC, SOC, TN, soil pH, and clay, silt and sand content were included in the regression models. However, high levels of multicollinearity between variables could lead to over-fitting of the model. For this purpose, the Variance Inflation Factor (VIF) of each variable was calculated and variables that had high VIF values (> 10) were removed. The remaining soil properties were kept as predicators in the regression model, and Akaike Information Criterion (AIC) was used to select the best combination of predictors for each model. For the SPF + PTF approach, the same soil properties selected in the PTF approach, except for SWC, were used. Partial least-squares regression (PLSR) was applied to build spectrotransfer functions for the prediction of the elementary soil properties (i.e. SOC, TN, soil pH, and clay, silt and sand content) in the SPF + PTF approach and the MWD and aggregate size fractions in the SPF approach. Three pre-processing algorithms (i.e. log(1/Reflectance), 1st derivative with Savitzky–Golay third polynomial smoothing, and Standard Normal Variate with detrending) were used to pre-treat the soil spectra and the algorithm that gave the best performance was adopted. Tenfold cross-validation was performed to assess the predictive capability of the PLSR model for the calibration dataset, and then the model was evaluated against the validation dataset. In order to detect the main spectral regions that affect MWD estimation, we calculated the Variance Importance Projection (VIP) index, which is a weighted sum of squares of the PLS weights. Spectral bands with VIP values greater than one are considered significant for the PLSR model because the average of VIP scores is equal to one (Chong and Jun, 2005).

Table 1 Summary of the measured soil properties and aggregate stability indexes. Soil properties

Mean

Minimum

Maximum

Standard deviation

SOC (g/100 g) TN (g/100 g) SWC (g/100 g) Clay content (g/100 g) Silt content (g/100 g) Sand content (g/100 g) Soil pH (H2O) MWD (mm) < 63 μm (g/100 g) 63–250 μm (g/100 g) > 250 μm (g/100 g)

1.30 0.13 16.55 10.62 68.55 20.82 6.76 0.47 18.15 42.30 39.55

0.67 0.07 10.62 7.03 35.77 9.13 4.68 0.20 6.30 15.96 18.34

2.47 0.24 30.56 14.63 78.38 57.17 7.99 1.44 33.35 72.58 72.68

0.42 0.04 2.96 1.92 11.16 12.14 0.69 0.30 5.56 12.45 14.46

SOC, soil organic carbon; TN, total nitrogen; SWC, soil water content; MWD, mean weight diameter.

2.5.3. Model evaluation The performances of the three approaches were evaluated by examining the coefficient of determination (R2) of observed against predicted values, root mean squared error (RMSE), Ratio of Performance to Deviation (RPD) and Ratio of Performance to Interquartile Range (RPIQ). The threshold values of RPD defined by Chang et al. (2001) are widely used in soil spectroscopy literature to assess the model accuracy. The authors considered models with an excellent prediction capability to have RPD > 2, intermediate capability to have RPD values of 2–1.4 and unreliable model to have RPD < 1.4. Moreover, it was suggested by Bellon-Maurel et al. (2010) that RPIQ could provide a robust statistic to describe the model performance by correctly describing the range of variation for data with a non-normal distribution. Hence, both RPD and RPIQ were reported for model evaluation. To compare the model performances of the three approaches, the RMSE and RPD values resulting from the 100 model validations were used for significance tests. Pairwise t-test followed by a Bonferroni correction was performed to detect statistically significant differences (P < 0.01) among the three approaches. All the data analyses in this study were carried out in R version 3.4.1 (R Development Core Team, 2017).

Fig. 2. Pearson correlation coefficient between elementary soil properties and aggregate mean weight diameter (MWD). Blank cells indicate that (partial) correlation was not significant at P < 0.01. SOC, soil organic carbon; Total.N, total nitrogen; SWC, soil water content.

large variations (CV < 20%). The largest CV of all the soil properties was found for MWD, which spanned from “highly unstable” (< 0.4 mm) to “stable” (> 0.8 mm) according to the AS classification proposed by Le Bissonnais (1996). Among different AS classes, significant differences in mass percentages were detected for each of the three aggregate size fractions, except for the clay + silt fraction, in which no significant difference was found between “unstable” and “stable” classes (Fig. 3a). Also, with the clay + silt fraction remaining at around 20% across different AS classes, the “stable” soils in the study area were characterized by a high percentage of macroaggregates and low percentage of microaggregates, and vice versa for the “highly unstable” soils. A highly positive correlation between MWD and SOC (r = 0.84) was found, indicating a major control of SOC on AS level (Fig. 3b). There seemed to be a critical SOC level (1.5%) below which a majority of the samples (87%) fell in the “highly unstable” AS class. A second SOC threshold (2%), although less prominent, appeared to an important line that separates the “unstable” and “stable” soils. 67% of all the samples between 1.5% and 2% SOC were in the “unstable” class, while all the samples that had SOC

3. Results 3.1. Summary of soil properties The SOC content was rather variable (coefficient of variation (CV): 32%) and was on average below 1.5% (Table 1). Textural analysis show that the selected soils are mainly characterized by a silt loam according to the USDA classification system. A large range of variation occurred in the silt and sand content, whereas the clay content remained at a low level for all the samples. This is further evidenced by the correlation analysis (Fig. 2), in which the highly negative correlation (r = −0.99) between silt and sand implied a “seesaw relationship” between the two variables. For the other elementary soil properties, TN was positively correlated with SOC (r = 0.97), and SWC and soil pH did not display 4

Geoderma 357 (2020) 113958

P. Shi, et al.

Fig. 3. (a) Comparisons of aggregate size distribution resulting from wet-sieving test among different aggregate stability classes (i.e. highly unstable, unstable and stable) and (b) scatterplot of aggregate mean weight diameter (MWD) versus soil organic carbon content.

models for the four AS indexes, but the VIF calculations showed high multicollinearity for SOC, TN, Silt and Sand. TN and Sand were therefore excluded in the models due to their respective correlations with SOC and Silt. The remaining soil variables were used to derive the best MLR models, and regression equations that gave the best performances according to AIC are shown in Table 3. Excellent models were achieved for MWD and microaggregate (RPD > 2), followed by macroaggregate (RPD = 1.91), and a correct model was also achieved for clay + silt fraction (RPD = 1.52) (Fig. 4). Regarding the importance of each soil variable in controlling AS indexes, it can be seen that in none of the four equations was Clay a significant variable, while SOC and pH were found to be significant variables in controlling all the AS indexes. According to the Relative Importance of each soil variable in the regression equations, SOC was the most important factor for MWD, microaggregate and macroaggregate fractions, implying that the positive control of SOC over AS (in the case represented by MWD) was reflected in its positive relationship with macroaggregate fraction and negative relationship with microaggregate fraction. Other significant variables included Silt and SWC for the predictions of MWD, clay + silt and microaggregate fractions, with Silt being the most important variable for clay + silt fraction.

Table 2 Vis-NIR spectroscopic prediction of elementary soil properties using PLSR. Values are the mean values resulting from 100 predicting in the validation dataset. Soil variables

Pre-processing

R2

RMSE

RPD

RPIQ

SOC

1st derivative + Savitzky–Golay smoothing log(1/reflectance) log(1/reflectance) 1st derivative + Savitzky–Golay smoothing

0.77

0.215

2.01

2.37

0.76 0.59 0.82

0.021 0.480 4.996

2.08 1.49 2.28

2.23 1.69 2.50

TN pH Silt

concentration larger than 2% were in the “stable” class. Other notable correlations between elementary soils properties and MWD (Fig. 2) include the positive correlation between TN and MWD, which can be attributed to the propagation from the correlation between SOC and TN, and the negative correlation between pH and MWD. 3.2. Spectroscopic prediction of elementary soil properties Before introducing the results of the prediction models for AS indexes, PLSR-based spectrotransfer models were established to link soil spectra with elementary soil properties. As shown in Table 2, excellent models were achieved for SOC, TN and Silt (R2 > 0.75, RPD > 2), followed by an intermediate model for soil pH (RPD > 1.4). In addition, it should be noted that for the three textural fractions, only silt was included in the development of prediction models, because of the high correlation between silt and sand (r = −0.99, Fig. 2) and the fact that clay content in the study area remained at a low level (7%–14%, Table 1) and thus was an insignificant factor in determining the variability of AS.

3.3.2. SPF approach Using Vis-NIR spectra to directly predict AS indexes, the SPF approach provided good models for MWD, microaggregate and macroaggregate fractions, with RPD ranging from 1.59 to 1.82 (Fig. 4). However, the prediction model for the clay + silt fraction was inaccurate, with R2 and RPD below 0.5 and 1.4 respectively. In order to further examine the relationship between spectral signatures and the level of AS, a VIP figure was produced to depict the important spectral bands that determine the performance of the spectrotransfer function for MWD (Fig. 5). The most important spectral features for MWD prediction were found in the Vis region at 400–600 nm, and then in the NIR region at 1800–1900 nm and around 2300–2400 nm.

3.3. Prediction of aggregate stability indexes 3.3.1. PTF approach MLR models were built to explore the contribution of various elementary soil properties to the level of AS. Initially, all the measured elementary soil properties were included as predictors to establish MLR

3.3.3. SPF + PTF approach The SPF + PTF approach used the predicted elementary soil properties from soil spectra to build pedotransfer functions for the

Table 3 Multiple linear regression equations developed using PTF and SPF + PTF approaches. Soil variables included in the equations (P < 0.05) from left to right are sorted according to their respective Relative Importance in determining the aggregate stability indexes. Aggregate stability indexes

PTF

SPF + PTF

MWD < 63 μm 63–250 μm > 250 μm

=1.20 + 0.527 ∗ SOC − 0.174 ∗ pH − 0.015 ∗ SWC =0.26 + 0.843 ∗ Silt + 0.356 ∗ pH − 3.216 ∗ SOC − 0.486 ∗ SWC =48.40–20.598 ∗ SOC − 0.382 ∗ pH + 5.541 ∗ Silt + 0.569 ∗ SWC =50.73 + 23.730 ∗ SOC − 6.226 ∗ pH

=0.91 + 0.50 ∗ SOC − 0.16 ∗ pH =3.05 + 0.36 ∗ Silt − 0.27 ∗ pH − 5.95*SOC =64.80–20.83 ∗ SOC + 4.35 ∗ pH − 0.36 ∗ Silt =32.55 + 26.81 ∗ SOC − 4.11 ∗ pH

5

Geoderma 357 (2020) 113958

P. Shi, et al.

Fig. 4. Observed versus predicted aggregate mean weight diameter (MWD), clay + silt (< 63 μm), microaggregate (63–250 μm) and macroaggregate (> 250 μm) in the validation datasets using the three different approaches. Symbols are the average values resulting from 100 replicated simulations after randomly splitting the calibration and validation datasets each time, and error bars are standard deviations. Dashed line indicates the 1:1 line.

which could not be predicted using soil spectral information due to the fact that air-dried soils were used for soil spectra acquisition. As shown in Table 3, the regression equations resulting from the SPF + PTF approach followed a similar structure as those from the PTF approach for the four AS indexes, with the same order of relative importance for the included soil variables, excluding SWC.

estimation of AS indexes. Similar to the model performances using the SPF approach, good models were achieved for MWD, microaggregate and macroaggregate fractions, while no correct model was developed for the prediction of clay + silt fraction either. For the establishment of pedotransfer functions via MLR, same soil variables selected in the PTF were used in the SPF + PTF approach, with the exception of SWC, 6

Geoderma 357 (2020) 113958

P. Shi, et al.

stressed that the decline in SOC content was accompanied by continuous degradation in AS, until passing another critical SOC level (1.5%), below which the soils became highly susceptible to erosion. This highlights the importance of maintaining an optimal SOC level in the context of preventing soil degradation induced by surface crusting and soil erosion. The effect of SOC on soil aggregation was also reflected in its control on the mass percentages of different aggregate size fractions. SOC was found to be the most important soil variable in explaining the variance in microaggregate and macroaggregate fractions. The increase in SOC content promotes the formation of macroaggregates by cementing microaggregates and primary particles (Le Bissonnais, 1996; Tisdall and Oades, 1982), thereby resulting in the increase in macroaggregate fraction and decrease in microaggregate fraction (Fig. 3a). The variation in clay + silt fraction among different AS classes was less pronounced and was mostly explained by the silt content (Table 3), indicating that soil texture, which did not display large variation due to the rather stable clay content in the study area with homogenous soil type, could be the reason for the small variation in clay + silt fraction. From the pedotransfer functions developed using the PTF approach, although less important, soil pH and SWC were also significant soil variables in determining AS indexes, corroborating with previous studies on the effects of soil pH (Chappell et al., 1999) and moisture condition (Shi et al., 2017) on AS. Overall, because of the correlations of elementary soil properties with AS, the MLR models established using the PTF approach were able to predict AS indexes reasonably well, in particular for MWD, microaggregate and macroaggregate fractions (Fig. 4). The relatively less accurate model for clay + silt fraction was possibly due to the small variation in soil texture caused by the homogenous soil type.

Fig. 5. Variable importance projection (VIP) in partial least squares regression (PLSR) model for MWD using the SPF approach. Spectral bands with VIP values greater than one were considered significant for the PLSR model. The grey area represents the standard deviation of the VIP values (n = 100) at each wavelength.

3.3.4. Model comparisons Apart from the fact that all the three modelling approached produced acceptable results in terms of predicting AS indexes, each modelling approach gave a different predictive performance. Comparing the three different approaches based on RMSE and RPD, the predictive performance for all the four AS indexes was consistently higher in the PTF models built upon the measured elementary soil properties, as indicated by the lowest RMSE and the highest RPD values (Fig. 6). In addition, the PTF approach also produced the most robust estimates, as indicated by a smallest standard deviation of R2 from the 100 model simulations (see error bars in Fig. 5). Between the two approaches that used Vis-NIR soil spectra for the prediction of AS indexes, no significant differences in RPD was detected for MWD, clay + silt and macroaggregate fractions (Fig. 6), while the SPF + PTF approach outperformed the SPF approach for the prediction of the microaggregate fraction, with significantly higher RPD and lower RMSE. For the prediction of the macroaggregate fraction, RMSE was significant lower for the SPF + PTF approach and although there was no significant difference, the average RPD was higher for the SPF + PTF approach than for the SPF approach. Lastly, it should be noted again that the clay + silt fraction could not be well predicted in this study using spectral data.

4.2. Applicability of Vis-NIR soil spectra to aggregate stability prediction Both the direct SPF approach and indirect SPF + PTF approach produced good predictive models for MWD, microaggregate and macroaggregate fractions, but not for clay + silt fraction. For the SPF approach, since AS indexes are only secondary soil properties that do not directly react to Vis-NIR radiation, the prediction of such properties are carried out through their correlations with Vis-NIR-reactive elementary soil properties (Soriano-Disla et al., 2014). Not surprisingly, the most important spectral bands identified in the SPF model for MWD prediction were in the visible region at 400–600 nm, which is known to be associated with the darkness of the soils and thus the quantity of SOC (Ben-Dor et al., 1997; Castaldi et al., 2018), and in the NIR region at 1800–1900 nm and 2300–2400 nm, which relate to the features of organic compounds like lignin and cellulose and phenolic, amide and aliphatic functional groups (Ben-Dor et al., 1997). This confirms once again the primary control of SOC on AS prediction. For the SPF + PTF approach, SOC, TN, pH and Silt were predicted first via PLSR-based spectrotransfer functions before being used for subsequent AS predictions. Comparing the model performances to other studies, as summarized in Soriano-Disla et al. (2014), the R2 for Silt was higher than the reported median value of previous predictions using Vis-NIR spectroscopy, while the R2 for SOC, TN and pH in this study was slightly lower than their reported median values. The reason could be attributed to the relatively narrow range of variation in these properties in the study area. Nonetheless, MLR models developed using the predicted elementary soil properties gave good performances. Between the two approaches that both used pedotransfer functions to predict AS indexes, the less accurate models from the SPF + PTF approach could be in part because of the error propagation from the initial spectroscopic prediction of elementary soil properties to subsequent MLR models, and in part because of the exclusion of SWC from the SPF + PTF models. Numerous studies have demonstrated the effect of SWC on AS (Algayer et al., 2014; Dimoyiannis, 2009; Perfect et al., 1990), and failing to account for such an effect could thus possibly lead

4. Discussion 4.1. Control of SOC and other variables over the prediction of aggregate stability It is well-documented that SOC is a positive controlling factor over the level of AS (Abiven et al., 2009; Le Bissonnais et al., 2007). Loveland and Webb (2003) reported that the value of 2% SOC is widely used in soil science literature as a threshold below which soils become physically less stable and more susceptible to surface crusting and erosion. They also concluded that quantitative evidence of such a threshold was still lacking. In addition, Meersmans et al. (2011) stated in a modelling study that the temporal SOC concentration decline in the Belgian Loam Belt could lead to serious degradation in soil structure, but little quantitative information is available to enable direct assessment on the relationship between SOC and AS. Here in our study, the existence of a positive correlation between SOC and MWD and the fact that two critical SOC thresholds (i.e. 1.5% and 2%) separated the “Stable”, “Unstable” and “Highly unstable” AS classes provided quantitative evidence that supports the empirical 2% SOC threshold. It also 7

Geoderma 357 (2020) 113958

P. Shi, et al.

Fig. 6. Boxplots of RMSE and RPD values resulting from the 100 replicated simulations using each of the three approaches. Paired t-test with Bonferroni correction was performed for each AS index to detect significant differences at P < 0.01 among three approaches.

unsatisfactory results. This means that (1) predicting secondary soil properties (such as AS) via Vis-NIR soil spectroscopy rely on the relation between the target variable and soil chromophores; and (2) spectrotransfer models developed by virtue of correlations between elementary soil properties and the target variable are likely to be sitedependent and should be treated with caution when applied to areas different from where they were developed. But in areas where Vis-NIR soil spectroscopy can be used to investigate AS more efficiently, especially considering the upcoming satellite-based hyperspectral imagers that will provide continuous fine-resolution Vis-NIR spectral data across large spatial scales, monitoring AS in space and time to provide more informative data for better land management and erosion control measures becomes possible. The capability of predicting different aggregate size fractions (i.e. microaggregate and macroaggregate) with distinct SOC composition using Vis-NIR spectroscopy could also offer opportunities to assess SOC stability and C sequestration at unprecedented scales.

to biased estimations of AS. The latter has important implications for the employment of laboratory soil spectra, which in most cases are obtained with dry soils, to predict moisture-sensitive soil properties like AS. Future studies should explore the possibility of using field spectra obtained directly on fresh soils, which enable the estimation of soil water content while at the same time offer the opportunity to account for the effect of varying soil moisture on field soil spectra (Nocita et al., 2013; Wijewardane et al., 2016). Between PTF and SPF approaches, the less accurate models for the SPF approach could be because that the variations in the important controlling factors of AS indexes were not sufficiently reflected in the variability of soil spectra. Nevertheless, the SPF and SPF + PTF approaches both produced good models for not only the lumped index MWD but also for the microaggregate and macroaggregate fractions. The use of soil spectra thus proved to be effective in predicting AS indexes, as also evidence by Erktan et al. (2016) who used MIR to successfully predict MWD. On the contrary, Gomez et al. (2013) and Canasveras et al. (2010) used Vis-NIR and Vis-NIR-MIR spectroscopy to predict AS and did not achieve satisfactory models. They concluded that a lack of relation between SOC or any other elementary soil property with AS was the reason for the 8

Geoderma 357 (2020) 113958

P. Shi, et al.

5. Conclusion

and sub-erosional landforms with aggregate stability variations in a tropical Ultisol disturbed by forestry operations. Soil Tillage Res. 50 (1), 55–71. Chong, I.-G., Jun, C.-H., 2005. Performance of some variable selection methods when multicollinearity is present. Chemom. Intell. Lab. Syst. 78 (1), 103–112. De Roo, A.P.J., Wesseling, C.G., Ritsema, C.J., 1996. LISEM: a single-event physically based hydrological and soil erosion model for drainage basins.1. Theory, input and output. Hydrol. Process. 10 (8), 1107–1117. Dimoyiannis, D., 2009. Seasonal soil aggregate stability variation in relation to rainfall and temperature under Mediterranean conditions. Earth Surf. Process. Landf. 34 (6), 860–866. Erktan, A., Legout, C., De Danieli, S., Daumergue, N., Cecillon, L., 2016. Comparison of infrared spectroscopy and laser granulometry as alternative methods to estimate soil aggregate stability in Mediterranean badlands. Geoderma 271, 225–233. Evrard, O., Vandaele, K., Bielders, C., Wesemael, B.v., 2008. Seasonal evolution of runoff generation on agricultural land in the Belgian loess belt and implications for muddy flood triggering. Earth Surf. Process. Landf. 33 (8), 1285–1301. Gomez, C., Le Bissonnais, Y., Annabi, M., Bahri, H., Raclot, D., 2013. Laboratory Vis-NIR spectroscopy as an alternative method for estimating the soil aggregate stability indexes of Mediterranean soils. Geoderma 209, 86–97. Gumiere, S.J., Le Bissonnais, Y., Raclot, D., 2009. Soil resistance to interrill erosion: Model parameterization and sensitivity. Catena 77 (3), 274–284. Jetten, V., de Roo, A., Favis-Mortlock, D., 1999. Evaluation of field-scale and catchmentscale soil erosion models. Catena 37 (3–4), 521–541. Le Bissonnais, Y., 1996. Aggregate stability and assessment of soil crustability and erodibility.1. Theory and methodology. Eur. J. Soil Sci. 47 (4), 425–437. Le Bissonnais, Y., Blavet, D., De Noni, G., Laurent, J.Y., Asseline, J., Chenu, C., 2007. Erodibility of Mediterranean vineyard soils: relevant aggregate stability methods and significant soil variables. Eur. J. Soil Sci. 58 (1), 188–195. Legout, C., Leguedois, S., Le Bissonnais, Y., 2005. Aggregate breakdown dynamics under rainfall compared with aggregate stability measurements. Eur. J. Soil Sci. 56 (2), 225–237. Loveland, P., Webb, J., 2003. Is there a critical level of organic matter in the agricultural soils of temperate regions: a review. Soil Tillage Res. 70 (1), 1–18. Meersmans, J., van Wesemael, B., Goidts, E., Van Molle, M., De Baets, S., De Ridder, F., 2011. Spatial analysis of soil organic carbon evolution in Belgian croplands and grasslands, 1960–2006. Glob. Chang. Biol. 17 (1), 466–479. Nocita, M., Stevens, A., Noon, C., van Wesemael, B., 2013. Prediction of soil organic carbon for different levels of soil moisture using Vis-NIR spectroscopy. Geoderma 199, 37–42. Nocita, M., Stevens, A., van Wesemael, B., Brown, D.J., Shepherd, K.D., Towett, E., Vargas, R., Montanarella, L., 2015. Soil spectroscopy: an opportunity to be seized. Glob. Chang. Biol. 21 (1), 10–11. Perfect, E., Kay, B.D., Vanloon, W.K.P., Sheard, R.W., Pojasok, T., 1990. Factors influencing soil structural stability within a growing-season. Soil Sci. Soc. Am. J. 54 (1), 173–179. Shi, P., Schulin, R., 2018. Erosion-induced losses of carbon, nitrogen, phosphorus and heavy metals from agricultural soils of contrasting organic matter management. Sci. Total Environ. 618, 210–218. Shi, P., Van Oost, K., Schulin, R., 2017. Dynamics of soil fragment size distribution under successive rainfalls and its implication to size-selective sediment transport and deposition. Geoderma 308 (Supplement C), 104–111. Soriano-Disla, J.M., Janik, L.J., Rossel, R.A.V., Macdonald, L.M., McLaughlin, M.J., 2014. The performance of visible, near-, and mid-infrared reflectance spectroscopy for prediction of soil physical, chemical, and biological properties. Appl. Spectrosc. Rev. 49 (2), 139–186. Tisdall, J.M., Oades, J.M., 1982. Organic-matter and water-stable aggregates in soils. J. Soil Sci. 33 (2), 141–163. Van Oost, K., Beuselinck, L., Hairsine, P.B., Govers, G., 2004. Spatial evaluation of a multi-class sediment transport and deposition model. Earth Surf. Process. Landf. 29 (8), 1027–1044. Vandaele, K., Poesen, J., 1995. Spatial and temporal patterns of soil erosion rates in an agricultural catchment, central Belgium. Catena 25 (1), 213–226. Viscarra Rossel, R.A., Behrens, T., Ben-Dor, E., Brown, D.J., Dematte, J.A.M., Shepherd, K.D., Shi, Z., Stenberg, B., Stevens, A., Adamchuk, V., Aichi, H., Barthes, B.G., Bartholomeus, H.M., Bayer, A.D., Bernoux, M., Bottcher, K., Brodsky, L., Du, C.W., Chappell, A., Fouad, Y., Genot, V., Gomez, C., Grunwald, S., Gubler, A., Guerrero, C., Hedley, C.B., Knadel, M., Morras, H.J.M., Nocita, M., Ramirez-Lopez, L., Roudier, P., Campos, E.M.R., Sanborn, P., Sellitto, V.M., Sudduth, K.A., Rawlins, B.G., Walter, C., Winowiecki, L.A., Hong, S.Y., Ji, W., 2016. A global spectral library to characterize the world's soil. Earth-Sci. Rev. 155, 198–230. Wijewardane, N.K., Ge, Y., Morgan, C.L.S., 2016. Moisture insensitive prediction of soil properties from VNIR reflectance spectra based on external parameter orthogonalization. Geoderma 267, 92–101.

We explored the possibility of using Vis-NIR soil spectroscopy to predict aggregate MWD and mass percentages of clay + silt, microaggregate and macroaggregate fractions. Among the three adopted approaches, the PTF approach, which used the measured elementary soil properties to build pedotransfer functions, consistently produced the best results, while all three approaches, including two spectra-based approaches (i.e. SPF and SPF + PTF approaches) provided satisfactory prediction models for MWD, microaggregate and macroaggregate fractions, primarily due to the positive correlation between AS and the SOC content. We demonstrated that soil Vis-NIR spectroscopy is a promising technique that enables fast and efficient large-scale assessment of secondary soil properties (e.g. aggregate stability) that are correlated with chromophores such as SOC, especially taking into account the fact that upcoming satellite hyperspectral imagers will provide high-resolution Vis-NIR spectral data across large spatiotemporal scales. Acknowledgements P. Shi is a Swiss National Science Foundation Early postdoc. Mobility Fellowship holder (PEEZP2_178494). This study also received funding from National Natural Science Foundation of China via grant no. 41807059. Additional support was provided by the Science and Technology Development Plan Program of Jilin Province, China (20190103108JH). References Abiven, S., Menasseri, S., Chenu, C., 2009. The effects of organic inputs over time on soil aggregate stability — a literature analysis. Soil Biol. Biochem. 41 (1), 1–12. Algayer, B., Le Bissonnais, Y., Darboux, F., 2014. Short-term dynamics of soil aggregate stability in the field. Soil Sci. Soc. Am. J. 78 (4), 1168–1176. Amezketa, E., 1999. Soil aggregate stability: a review. J. Sustain. Agr. 14 (2–3), 83–151. Annabi, M., Raclot, D., Bahri, H., Bailly, J.S., Gomez, C., Le Bissonnais, Y., 2017. Spatial variability of soil aggregate stability at the scale of an agricultural region in Tunisia. Catena 153, 157–167. Askari, M.S., Cui, J.F., O’Rourke, S.M., Holden, N.M., 2015. Evaluation of soil structural quality using VIS-NIR spectra. Soil Tillage Res. 146, 108–117. Bellon-Maurel, V., Fernandez-Ahumada, E., Palagos, B., Roger, J.-M., McBratney, A., 2010. Critical review of chemometric indicators commonly used for assessing the quality of the prediction of soil attributes by NIR spectroscopy. TrAC Trends Anal. Chem. 29 (9), 1073–1081. Ben Dor, E., Ong, C., Lau, I.C., 2015. Reflectance measurements of soils in the laboratory: standards and protocols. Geoderma 245, 112–124. Ben-Dor, E., Inbar, Y., Chen, Y., 1997. The reflectance spectra of organic matter in the visible near-infrared and short wave infrared region (400–2500 nm) during a controlled decomposition process. Remote Sens. Environ. 61 (1), 1–15. Beuselinck, L., Steegen, A., Govers, G., Nachtergaele, J., Takken, I., Poesen, J., 2000. Characteristics of sediment deposits formed by intense rainfall events in small catchments in the Belgian Loam Belt. Geomorphology 32 (1), 69–82. Canasveras, J.C., Barron, V., del Campillo, M.C., Torrent, J., Gomez, J.A., 2010. Estimation of aggregate stability indices in Mediterranean soils by diffuse reflectance spectroscopy. Geoderma 158 (1–2), 78–84. Castaldi, F., Chabrillat, S., Chartin, C., Genot, V., Jones, A.R., van Wesemael, B., 2018. Estimation of soil organic carbon in arable soil in Belgium and Luxembourg with the LUCAS topsoil database. Eur. J. Soil Sci. 69 (4), 592–603. Castaldi, F., Hueni, A., Chabrillat, S., Ward, K., Buttafuoco, G., Bomans, B., Vreys, K., Brell, M., van Wesemael, B., 2019. Evaluating the capability of the Sentinel 2 data for soil organic carbon prediction in croplands. Isprs J Photogramm 147, 267–282. Chang, C.-W., Laird, D.A., Mausbach, M.J., Hurburgh, C.R., 2001. Near-infrared reflectance spectroscopy–principal components regression analyses of soil properties. Soil Sci. Soc. Am. J. 65 (2), 480–490. Chappell, N.A., Ternan, J.L., Bidin, K., 1999. Correlation of physicochemical properties

9