Open AccessArticle

Detection of Soil Erosion Hotspots in the Croplands of a Typical Black Soil Region in Northeast China: Insights from Sentinel-2 Multispectral Remote Sensing

Lulu Qi

^1,2,

Pu Shi

^1,2,*

Klara Dvorakova

Kristof Van Oost

²,

Qi Sun

¹,

Hanqing Yu

³ and

Bas van Wesemael

College of Earth Sciences, Jilin University, Changchun 130061, China

Earth and Life Institute, UCLouvain, 1348 Louvain-la-Neuve, Belgium

Institute of Environment and Sustainable Development in Agriculture, Chinese Academy of Agricultural Sciences, Beijing 100081, China

Author to whom correspondence should be addressed.

Remote Sens. 2023, 15(5), 1402; https://doi.org/10.3390/rs15051402

Submission received: 16 December 2022 / Revised: 18 February 2023 / Accepted: 28 February 2023 / Published: 2 March 2023

(This article belongs to the Special Issue Remote Sensing for Soil Organic Carbon Mapping and Monitoring)

Download

Browse Figures

Figure 1
(a,b) Spatial distribution of the sampling points overlaid on the digital elevation model of the catchment in northeast China; (c,d) a Sentinel-2 true color image (13 May 2021) of the catchment showing the large percentage of exposed bare soils, and zoom-in areas depicting the selection of sampling points; (e,f) representative sampling locations along typical slope profiles with average length at ca. 300 m and their corresponding slope degrees. Different letters indicate significant differences among the three groups (p < 0.05). "> Figure 2
(a) Topographic position index (TPI), (b) estimated net soil erosion rates converted from 137Cs inventory, (c) soil organic carbon (SOC) content and (d) C/N ratio at Summit, Mid-slope and Foot-slope positions. Different letters indicate significant differences among the three groups (p < 0.05). "> Figure 3
Classification of three slope positions based on linear discriminant analysis (LDA) of PC scores arising from laboratory-based VNIR spectra. Subsoil spectra of Summit positions were included to investigate any spectral similarity to the topsoil of Mid-slope positions. Histograms on the right show the distribution of the first LD function’s value for the four sample groups. "> Figure 4
Laboratory-based mean spectra of topsoil in three soil erosion intensity classes. (a) Raw spectra and (b) continuum-removed reflectance. Shaded areas represent standard deviations from the mean. "> Figure 5
Boxplots of laboratory spectral indices of topsoil used for spectral separation of erosion intensity classes. (a) Mean reflectance over the visible region (400–780 nm), (b) slope between 800 and 1350 nm, and (c) continuum-removed reflectance value at 670 nm. Different letters indicate significant differences among the three groups (p < 0.05). "> Figure 6
Classification of three erosion intensity classes based on linear discriminant analysis (LDA) of PC scores derived from Sentinel-2 bare soil spectra. Histograms on the right show the distribution of the first LD function’s value for the three classes. "> Figure 7
Sentinel-2-based mean spectra of three soil erosion intensity classes. (a) Raw spectra and (b) continuum-removed reflectance. Shaded areas represent standard deviations from the mean. Sentinel-2 band width and positions are indicated in the left figure. "> Figure 8
Boxplots of Sentinel-2 spectral indices used for spectral separation of erosion intensity classes. (a) Mean reflectance over the three visible bands (B2, B3, B4), (b) slope between B8 and B11, and (c) continuum-removed reflectance value at 670 nm. Different letters indicate significant differences among the three groups (p < 0.05). "> Figure 9
Soil erosion intensity map at 10 m resolution. The classification criteria were based on the results from <a href="#remotesensing-15-01402-f008" class="html-fig">Figure 8</a>. "> Figure 10
(a) Density plots of NDVI in June 2021 for the three erosion intensity classes within the cropland extent; (b,c) zoomed-in areas with detailed soil erosion pattern and (d,e) the corresponding field-scale NDVI maps to depict the relation of crop variability with erosion-induced variations in soil productivity. ">

Versions Notes

Abstract

Global efforts to restore the world’s degraded croplands require knowledge on the degree and extent of accelerated soil organic carbon (SOC) loss induced by soil erosion. However, the methods for assessing where and to what extent erosion takes place are still inadequate for precise detection of erosion hotspots at high spatial resolution. Drawing on recent advances in multitemporal Sentinel-2 remote sensing to create a bare soil composite that reflects erosion-induced variations in soil spectral signatures, this study attempted to develop a spectra-based soil erosion mapping approach to pinpoint eroded hotspots in a typical catchment located in the black soil region of northeast China as characterized by undulating landscapes. We built a ground-truth dataset consisting of three classes of soils representing Severe, Moderate and Low erosion intensity because of their inter-class contrasts in estimated erosion rates from ¹³⁷Cs tracing. The spectral separability of different erosion classes was first tested by a combined principal component and linear discriminant analysis (PCA-LDA) against laboratory hyperspectral data and then validated against Sentinel-2-derived broadband spectra. The results show that PCA-LDA produced excellent classification accuracy (Kappa coefficient > 0.9) for both data sources and even more so for Sentinel-2 spectra, highlighting the effectiveness of the multitemporal approach to extract bare soil pixels. Further investigations into the spectral curves enabled identification of distinctive spectral features representative of shifting soil albedo and biochemical composition due to erosion-induced SOC mobilization. A classification scheme comprising the spectral features was applied to the Sentinel-2 bare soil composite for pixel-wise soil erosion mapping, in which 15.9% of the cropland area was detected as erosion hotspots, while the Moderate class occupied 65.4%. Comparing the erosion map to a NDVI map demonstrated the negative impact of soil erosion on crop growth from a spatial perspective, highlighting the potential of the proposed approach to aid targeted cropland management for food security and climate.

Keywords:

soil erosion; northeast China; Sentinel-2; multitemporal composite; mapping; SOC

1. Introduction

Land degradation affects roughly one-third of the world’s croplands [1,2], with soil erosion being arguably the most serious and widespread degradation form. Through erosion-induced lateral soil translocation and redistribution, significant truncation of soil profiles could occur in severely eroded areas, causing a spatial re-organization of carbon- and nutrient-rich topsoil materials. Such erosion hotspots often experience accelerated loss of soil organic carbon (SOC), and thus, soil fertility, thereby damaging food production in degraded croplands. Previous estimates reported an average reduction in crop yields of approximately 4% per 10 cm soil loss [3], and the extent of reduction is likely to be exacerbated in the developing world due to improper agricultural management and low levels of fertilization. Apart from the fact that accelerated SOC loss induced by erosion processes could exert detrimental effects on food production, SOC depletion in erosion hotspots could represent significant C sinks, provided that conservation measures are better targeted in those areas [4,5]. In this context, the Land Degradation Neutrality program under the framework of United Nations Sustainable Development Goals (SDGs) has explicitly adopted SOC as a key indicator to assess and monitor land degradation status [6]. Therefore, a better understanding on where and to what extent accelerated erosion and SOC loss take place could largely contribute to the global efforts on sustainable use of soil resources for food security and climate.

Ever since the release of the Global Assessment of Land Degradation (GLASOD) product in 1990, numerous attempts have been made to facilitate better quantification and spatial prediction of soil erosion, spearheaded by the widespread application of universal soil loss equation (USLE)-type algorithms to estimate gross soil erosion rates from regional to global scales [7,8]. As summarized in a comprehensive review on erosion modeling by Borrelli et al. [9], USLE-type models are by far the most widely adopted approach globally, occupying more than 40% of the erosion modeling studies conducted during the last two decades. One key advantage of USLE-type modeling can be attributed to its high degree of data accessibility, especially considering the fact that recent developments in remote-sensing and high-quality earth observation datasets could allow a more dynamic and spatially explicit erosion modeling [10]. Nevertheless, it remains an empirical method that only considers water erosion while neglecting other forms, such as tillage and wind erosion [11], and there is no simulation of soil deposition. Regarding the latter issue, a number of process-based physical models exist that simulate the interactive erosion and deposition processes during single rainfall events, mostly at the catchment scale [12,13], but the spatial predictive capability of such models is usually compromised by the imbalance between complexity in model structure and uncertainty in model parameterization [14]. Regardless of the type of modeling approaches used, the common issues hampering precise modeling and mapping of soil erosion also include: (1) the use of outdated, static and coarse resolution inputs, which are generally unable to capture the spatial and temporal variability of soil erosion at the scales where erosion processes operate; and (2) the lack of spatially distributed observations for rigorous model calibration and validation [9,15]. In addition to erosion modeling, soil erosion tracing techniques, such as the use of fallout radionuclides (e.g., ¹³⁷Cs), serve as a viable option to derive spatial estimates of net erosion, particularly at hillslope scale [16], but the costly nature of this approach indicates limited applicability across large spatial scales [17,18]. In light of the research gaps identified above, there is an urgent need for an effective soil erosion mapping approach, which enables precise detection of erosion hotspots as a result of multiple erosion processes.

Recent advances in remote sensing are well placed to address the aforementioned issues, not only because the increasing availability of high-resolution satellite images, such as Landsat-9 and Sentinel-2, allows unprecedented assessments of USLE factors at high resolution [10] and of land cover dynamics in response to soil erosion and degradation [2,19], but the rapid development in soil imaging spectroscopy also offers the potential to directly capture the erosion-induced variations in soil properties, particularly SOC [20,21]. The theoretical basis for using imaging spectroscopy for erosion detection capitalizes on the fact that soils of varying erosion intensities are characterized by distinct spectral signatures, indicative of heterogenous soil color, mineral composition, iron oxides, calcium carbonate and SOC, even over short distances [22,23,24]. These soil chemical and physical chromophores can interact with electromagnetic radiation within the visible and near-infrared (VNIR) spectral regions, producing distinct absorption features, which allows the development of chemometric approaches for efficient soil property prediction [25,26,27]. The widespread success of VNIR spectroscopic techniques in soil science has quickly expanded from laboratory to airborne and satellite platforms, and from quantification of soil properties to soil degradation assessments [21]. For example, Schmid et al. (2016) proposed a method to classify erosion stages in an agricultural region of Spain based on bare soil spectral characteristics derived from an airborne hyperspectral sensor [28]. They found that varying degree of erosion and deposition led to the emergence of different soil horizons at the surface, the morphological and chemical properties of which corresponded to distinct spectral signals, thus forming the basis for spatial classification and mapping of soil erosion intensity. In a subsequent study, Bracken et al. (2019) [29] further tested this approach against simulated satellite-based hyperspectral images (EnMAP) and achieved good results, highlighting the explanatory power of bare soil spectral imaging for erosion mapping.

Given that airborne hyperspectral imagery only provides limited spatiotemporal coverage and that spaceborne hyperspectral imaging missions are still under development [21], the current generation of satellite-based multispectral imagers, such as Sentinel-2, is well positioned to aid spatially explicit assessment of the degree and extent of soil erosion. Since the launch of Sentinel-2A in 2015, significant advances have been made in using Sentinel-2 imagery for multi-scale mapping of key soil properties, especially SOC, in croplands [20,30], thanks to its improved spatial, temporal and spectral resolution compared to its counterparts, such as Landsat-8 [31]. It is now well established that visible and near-infrared (VNIR) spectroscopic assessment of cropland soils driven by Sentinel-2 remote sensing requires (1) successful extraction of bare soil spectra representative of optimal soil conditions with minimal external disturbances (e.g., soil moisture, crop residue) [32] and (2) the employment of multitemporal compositing approaches that increase both the spatial coverage of bare soil surfaces and the robustness of the developed spectral prediction models [33,34]. However, few studies have directly related Sentinel-2-derived soil spectral information with the detection of soil erosion hotspots, except for a study by Žížala et al. (2019), who applied an unsupervised classification method against the extracted bare soil pixels from Sentinel-2 images and were able to delineate the severely eroded areas in a Chernozem region of the Czech Republic with acceptable accuracy [35]. The general applicability of this approach and the underlying mechanisms that support spectral characterization and classification among soils of varying erosion intensities remain to be further explored.

Here, we aim to build on recent progress in multitemporal Sentinel-2 remote sensing to detect erosion hotspots in a catchment located in the black soil region of northeast China, a notable breadbasket for the country, which produces more than 20% of the annual national grain output [36]. Prolonged intensive cultivation in this region has led the fertile Phaeozem and Chernozem soils to experience severe structural degradation [37], and recent reports show that this region is among the most seriously affected by soil erosion in China [7]. An effective and efficient method to detect localized erosion hotspots is therefore important for the implementation of targeted conservation measures. To this end, the goal of this study was to develop a methodological framework, which enables accurate classification and high-resolution mapping of soil erosion based solely on spectral characteristics. The specific objectives included: (1) verification of the feasibility of spectra-based soil erosion classification by first confronting laboratory VNIR soil spectra against soils in classes of different erosion intensities; and (2) development of a soil erosion classification scheme designated for Sentinel-2-derived bare soil spectra for high-resolution mapping of erosion hotspots.

2. Materials and Methods

2.1. Site Description

The catchment (44°34′–44°38′N, 125°51′–125°59′E), with an area of approximately 46.20 km², is located within the Mushi River basin of Dehui City in central Jilin Province, northeast China (Figure 1). It has a temperate continental monsoon climate, with average annual precipitation of 539 mm and average annual temperature of 4 °C. The majority of the precipitation falls from May to September during crop growth, with July typically experiencing the largest rainfall amount [38]. The “black soils” in the catchment, as they are commonly referred to by local farmers because of their dark color, are dominated by Phaeozems and Chernozems, according to the FAO World Reference Base classification [39]. The soils were developed in loess deposits of the Quaternary age with a silt loam texture. Croplands occupy nearly 90% of the catchment, with maize as the dominant crop type, mostly under a monocropping system. Conservation tillage was only introduced during the last decade to replace conventional ridge tillage, but at the time of sampling in 2021, the study area was still dominated by conventional tillage practices. This, combined with the undulating topography of long and gentle hillslopes typically found in northeast China, has led to substantial erosion-induced soil redistribution along hillslopes. Soil properties, such as SOC and soil texture, thus display a wide range of variations [40], making this area an ideal site to test the applicability of soil spectral imaging to capture erosion hotspots with accelerated SOC losses.

2.2. Soil Sampling and Analysis

To discriminate the spectral characteristics of soils under different erosion intensities, a ground-truth dataset was constructed through a comprehensive sampling campaign in October 2021. The designated sampling strategy involved three groups of soil samples taken at Summit, Mid-slope and Foot-slope positions from a set of pre-determined hillslopes (Figure 1), representing three classes of soil erosion intensity, i.e., Moderate, Severe and Low erosion intensity. Particular focus was given to the Mid-slope and Foot-slope positions, which were designated to be on the convex and concave slopes, respectively, while the Summit samples were taken from local plateaus with gentle slopes. The exact sampling locations were pinpointed in situ, accounting for any visible erosional features while also in reference to the slope and curvature maps of the catchment. The true color image shows that the sampling locations at Mid-slope were generally associated with a change of soil color due to substantial topsoil loss (Figure 1). A total of 72 sampling locations were selected, consisting of 19 at Summit, 28 at Mid-slope and 25 at Foot-slope.

To ensure a high degree of agreement between the three sample groups and their assigned class of erosion intensity, the validity of the designated sampling locations was assessed from three different aspects. First, topographic position index (TPI) was calculated from the one-arc-second SRTM digital elevation model by comparing a sampling point’s elevation to the mean elevation of its surrounding neighbors using the “spatialEco” R package. The goal was to verify whether the topographic features among the three sampling groups were distinguishable. Positive TPI values indicate convex locations, while negative values indicate concavities.

Second, at the 72 sampling locations, both topsoil (0–20 cm) and subsoil (50–70 cm) samples were collected to account for any vertical variation in SOC and nitrogen (N) content along soil profiles. It was assumed that topsoil SOC content at severely eroded areas (i.e., Mid-slope positions with significant truncation of the A-horizon) should mimic that of the subsoil at Summit with lower erosion intensity. At each sampling location, three pairs of topsoil and subsoil samples were collected within a 3 m radius and separately bulked to form two soil composites per location. Samples were then air-dried, passed through a 2 mm sieve and finely ground (<150 µm) for total C and N analysis by means of dry combustion (VarioMax CN Analyzer Elementar GmbH, Langenselbold, Hesse, Germany). SOC was considered the equivalent of the total C content, as none of the samples showed reaction upon 10% HCl treatment. The C/N ratio was also calculated to assess potential variations in SOC quality among different slope positions.

Third, the fallout radionuclide ¹³⁷Cs inventory at the sampling locations was determined and used to retrospectively estimate net erosion rates. The triplicate sampling scheme used for the C analysis above was kept the same to account for the local variability of ¹³⁷Cs inventory per location [41]. Samples were collected using a soil core sampler (38 mm diameter) at a depth of 30 cm at the Summit and Mid-slope positions, while a depth of 50 cm was sampled for the Foot-slope positions to account for excessive soil deposition. Additional soil pits were excavated at three Summit locations, and bulk soil samples were collected at a 5 cm increment up to the 50 cm depth to ensure that no additional ¹³⁷Cs activity was detected toward the subsoil layer. For the measurements of ¹³⁷Cs activity, composite soil cores (n = 3) at each sampling location were air-dried, passed through a 2 mm sieve and analyzed using a hyperpure coaxial Ge detector coupled to a multichannel analyzer. ¹³⁷Cs activity (Bq kg⁻¹) was detected at 662 keV peak with an analytical precision of ±6% requiring counting times of around 80,000 s [42] and then converted to inventory (Bq m⁻²) based on the sampling area and the weight of the soil cores. Finally, the simplified mass balance model (SMBM) for cultivated soils developed by Zhang et al. (1989) was used to convert ¹³⁷Cs inventory to net erosion rates (t ha⁻¹ yr⁻¹) [43]. Reference ¹³⁷Cs inventory from undisturbed soil required for the conversion model was taken from a 2005 study carried out in close vicinity of the catchment [44], and the value was decay-corrected to the sampling date considering the half-life of ¹³⁷Cs. Further details of the model structure and parameterization are given in the Supporting Information. Lastly, pair-wise Wilcox tests were conducted to inspect whether the differences in TPI, erosion rates and SOC among three sample groups were statistically significant at p < 0.05.

2.3. Laboratory and Sentinel-2-Based Soil Spectral Analysis

To develop a spectra-based methodological framework for the classification and mapping of soil erosion intensity, we followed a two-step procedure in which VNIR soil spectra obtained under controlled laboratory conditions were first linked with the ground-truth data to test the feasibility of using soil spectral information to discriminate the erosion classes. Then, the spectra-based classification scheme was validated against Sentinel-2-derived soil spectra in a second step, where multitemporal bare soil pixel compositing was applied to ensure spatially continuous soil erosion mapping.

2.3.1. Laboratory Soil Spectra Acquisition and Analysis

Laboratory VNIR soil spectra were obtained with an ASD Fieldspec 3 FR spectroradiometer (Analytical Spectral Devices Inc., Boulder, CO, USA). Sieved (<2 mm) and air-dried top and subsoil samples from the 72 sampling locations were placed in 9 cm diameter Petri dishes and scanned using an ASD Contract Probe with a built-in 100 W halogen lamp as the light source. All measurements were carried out in a dark room to minimize external disturbance, and the spectroradiometer was calibrated prior to the first and after every 20 measurements using a white SpectralonTM panel (Labsphere, North Sutton, NH, USA) to ensure the quality and stability of the measurements. Soil spectral reflectance was obtained from 400 to 2450 nm at 1 nm resolution for subsequent analysis. More detailed procedure is given in Shi et al. (2020) [45].

After the data acquisition, we conducted a combination of principal component analysis (PCA) and linear discriminant analysis (LDA) on VNIR spectra to investigate the spectral separability of soils under different erosion influences. PCA was first applied onto the raw spectra to reduce the dimension of the spectral information into a small number of non-collinear PCs, which explained 99% of the variance. Then, the PC scores resulting from PCA were subjected to LDA, which was used as a classifier to discriminate the three erosion classes. A confusion matrix was created by calculating the number of well and incorrectly classified classes, and the performance of the classification was assessed by Cohen’s Kappa coefficient.

Based on these analyses, further spectral characterization and processing was carried out against the raw spectral reflectance to highlight the spectral features and to ultimately develop an erosion classification scheme based on identified spectral thresholds. As adopted by Chamizo et al. (2012) [46], the first component of the classification scheme was to identify specific spectral regions, which displayed considerable differences among the sampling groups, and the slopes between those wavelengths were calculated as the difference in their reflectance divided by the spectral range. Furthermore, continuum removal (CR) baseline correction was used as an albedo normalization method to detect the prominent absorption features, which were appropriate for separating the erosion classes. CR works by fitting a convex hull to each spectrum and computing the deviations from the hull [47]. When applied against reflectance spectra, CR gives values of 1 to all parts of the spectrum that lie on the convex hull, while values less than 1 indicate the presence of absorption features.

2.3.2. Sentinel-2 Image Processing and Spectral Analysis

To test whether the above-defined classification scheme based on laboratory spectra could be translated onto the Sentinel-2 platform to facilitate spatially continuous mapping of soil erosion intensity, the key was to generate high-quality satellite-based soil spectra resembling those obtained under laboratory conditions. To achieve this, a multitemporal bare soil pixel compositing approach was adopted to optimize the stability of bare soil reflectance and maximize the spatial coverage of croplands [33]. Three Level-1C Sentinel-2 scenes (title number: T51TYK) from 29 April 2019, 13 May 2020 and 13 May 2021 were downloaded, atmospherically corrected to Level-2A using the Sen2Cor processor (standalone version 2.10) and resampled to a 10 m resolution. We selected these three scenes because of their low cloud cover percentage (<10%) and their sensing time periods in April–May, when soils in northeast China are generally prepared as seedbeds for maize sowing, thus representing an optimal bare surface condition. Further image processing procedures to extract bare soil pixels included: (1) masking of bad pixels consisting of clouds, thin cirrus and shadows based on the scene classification layer output from the Sen2Cor algorithm; and (2) using a combination of normalized vegetation difference index (NDVI: 0.10–0.25) and normalized burn ratio 2 (NBR2: 0–0.075) index thresholding to extract bare soil pixels under minimal disturbances. These NDVI and NBR2 thresholds were adopted from a recent study by Shi et al. (2022) [48], who reported the effectiveness of such index thresholding in removing pixels influenced by green vegetation, crop residue and soil moisture in the same study region. In particular, NBR2 is calculated as the normalized difference between B11 and B12 Sentinel-2 bands and has been frequently used for bare soil detection due to its sensitivity to excessive crop residue cover [32,49] and moisture content [34]. Lastly, the three processed scenes were mosaicked, and mean reflectance was calculated for pixels that had multiple bare soil occurrences.

The final multitemporal bare soil composite included 10 bands covering the visible (B2, B3, B4), red-edge (B5, B6, B7), NIR (B8, B8A) and short-wave infrared (SWIR bands, B11 and B12) regions, which served as the basis for subsequent tests on Sentinel-2-based soil erosion classification and mapping. The same PCA-LDA routine and the search for an effective classification scheme, as used against laboratory spectra, were applied to the Sentine-2-derived bare soil spectra. Classification performance and the adoption of spectral features for the classification scheme were compared between the two data sources.

2.4. Soil Erosion Mapping and Validation

The established spectral classification scheme to discriminate Low, Moderate and Severe erosion classes was applied to all bare soil pixels within the cropland extent of the catchment in order to achieve high-resolution mapping of soil erosion intensity. Particular focus was given to the percentage and spatial distribution of the severely eroded areas, namely the erosion hotspots. Finally, the soil erosion map was confronted against an NDVI map based on a Sentinel-2 image acquired on 22 June 2021, corresponding to the reflectance of the maize crop in its vegetative stage. The goal was to analyze whether the degree of erosion intensity had an impact on spatial variation in crop growth, as a measure to also validate the accuracy of the produced soil erosion map. All data analyses were carried out using R (version 4.1.3, R Core Team).

3. Results

3.1. Summary of Soil Analytical Results at Different Slope Positions

To verify the assumption that soils taken at different slope positions show differing degrees of erosion intensity, we compared TPI, estimated soil erosion rates from radionuclide ¹³⁷Cs tracing, SOC content and C/N ratio among the slope positions, as shown in Figure 2. The Summit and Mid-slope positions generally displayed positive TPI values, with significantly higher values found at Mid-slope positions, meaning that samples were collected from convexities with steeper slopes downward (Figure 1). The TPI at Foot-slope, however, showed marked differences from the other two groups. The pre-determined concave sampling locations at Foot-slope were verified by the negative TPI values, indicating a high possibility of soil deposition sourced from the upslope contributing area.

The estimated net soil erosion rates among the three slope positions were in accordance with observed significant differences in TPI values. Comparing the erosion rates converted from ¹³⁷Cs inventory data (Figures S1 and S2, Supporting Information), it can be seen that the Mid-slope positions were associated with the highest average erosion rate at 36.93 t ha⁻¹ yr⁻¹ and the widest range of variation, followed by the Summit positions with average erosion rate at 13.43 t ha⁻¹ yr⁻¹, while the Foot-slope positions had the lowest average (6.86 t ha⁻¹ yr⁻¹) and sometimes negative erosion rate, i.e., net deposition (Figure 2b). Several outliers existed, for instance the two points at Foot-slope with erosion rate unusually exceeding 20 t ha⁻¹ yr⁻¹, most likely due to human disturbances other than erosion processes, which caused significant soil loss. Overall, the fact that two-to-three-fold increases in average erosion rates were found from Foot-slope to Summit, and then to Mid-slope, corroborates the validity of using these three groups as proxies for erosion intensity classes, i.e., Low, Moderate and Severe.

This was further strengthened by the SOC measurements and C/N ratios at the three slope positions (Figure 2c,d). Between the Summit and Mid-slope positions, a clear drop in SOC content was expected and indeed observed from top to subsoils, with average SOC values of less than 10 g kg⁻¹ at Mid-slope due to erosion-induced losses of C-rich soil materials. Perhaps the more striking finding was that, although still significantly different, SOC in the topsoil of Mid-slope was close to that in the subsoil of the Summit, meaning that erosion-induced topsoil truncation at Mid-slope had at least partially exposed the B-horizon, which likely shares similar physicochemical properties to subsoils at Summit positions. This confirms the severity of accelerated SOC losses in soil erosion hotspots found at Mid-slope positions. There was no significant difference between topsoil SOC levels at Foot-slope and Summit, but the difference was that subsoil SOC content at Foot-slope positions was close to its topsoil counterparts and higher than the topsoil SOC at Mid-slope. This suggests that, although most of the Foot-slope soils still experienced net erosion (Figure 2b), the erosion-induced topsoil removal in severely eroded areas led to SOC accumulation and burial at Foot-slope positions, resulting in its relatively higher SOC content at both top and subsoils. This is supported by the fact that soils at Foot-slope were accompanied by the highest C/N ratios, on average, indicative of a higher share of particulate soil organic matter, typically released upon soil structural breakdown under erosive forces and selectively mobilized from upslope locations (Holz and Augustin, 2021), coinciding with the lowest C/N value at Mid-slope.

3.2. Laboratory-Based Spectral Discrimination of Soil Erosion Classes

3.2.1. PCA-LDA Classification

LDA was used to assess the separability of three erosion intensity classes based on a number of independent spectral features, as extracted by PCA. The scatterplot between the first and second LD showed good separability among the four groups (Figure 3). Specifically, the LD1 function could explain 82% of the spectral variance, and the distribution of its scores demonstrated clear spectral distinction between Foot-slope and Mid-slope, and to a lesser extent, between the topsoil of the Summit and Mid-slope. However, significant overlap in the LD1 existed between topsoil of the Mid-slope and subsoil of the Summit, confirming once again that severe erosion at Mid-slope led to the exposure of subsoils, which would otherwise have been protected in areas less affected by erosion. From here on forward, we will refer the three groups of soil samples taken at Summit, Mid-slope and Foot-slope positions to three classes (Moderate, Severe and Low) of erosion intensity, intended to be used as a ground-truth dataset to facilitate spectra-based soil erosion classification and mapping. In this line, PCA-LDA was further conducted against the VNIR spectra of these three groups. Assessment of the classification accuracy based on the confusion matrix comparing the predicted and observed erosion classes showed that three samples from the Low class and one sample from the Severe class were misclassified into the Moderate class, but in general, the classification accuracy was excellent, with an overall accuracy of 94% and a Kappa coefficient of 0.92 (Table 1). It should be noted that the classification accuracy refers to the performance of PCA-LDA against the entire training set without external validation, under which the rate of misclassification would likely increase.

3.2.2. Detection of Spectral Features in Support of Erosion Classification

Based on the premise that the spectral dissimilarity among erosion classes prompted an accurate classification, we further explored specific spectral features that supported such classification in an attempt to construct a classification scheme to be transferred to Sentinel-2 multispectral imagery. Inspecting the raw spectra (Figure 4), the Low and Moderate classes shared a similar pattern in the visible region (400–780 nm) but showed clear difference from the Severe class. For the NIR region, it seemed that the rate of increase, namely the slope, differed among the three classes, particularly the Low class with the steepest increase. For the SWIR region beyond the absorption peak at 1350 nm, the three classes appeared to have followed a similar spectral shape but differed in their albedo. Lastly, the continuum-removed reflectance values showed prominent absorption features at approximately 500, 670 and 2200 nm, plus two water absorption bands at 1440 and 1930 nm. In particular, the most separable absorption feature in relation to the erosion classes was found at 670 nm.

Considering the results outlined above, the mean reflectance values over the visible region, the slope between 800 and 1350 nm, and continuum-removed reflectance at 670 nm were grouped for the search of spectral thresholds to separate the three erosion classes. It can be seen in Figure 5 that mean reflectance between 400 and 780 nm could enable successful separation of the Severe erosion class but not the other two, while each of the latter two indices was found suitable to simultaneously separate all three classes. These three indices thus constitute the basis for a classification scheme to be tested against Sentinel-2 bare soil spectra below. It should be pointed out, however, that slight to significant overlap in these indices exists between Low and Moderate, hindering a complete separation for these two classes based on laboratory spectra.

3.3. Sentinel-2-Based Soil Erosion Classification and Mapping

3.3.1. Bare Soil Spectral Classification

Soil spectra with 10 Sentinel-2 bands were extracted from the 72 sampling locations from multitemporal bare soil composites, created by averaging over three scenes each acquired during the sowing season of 2019–2021, and consisting of bare soil pixels detected by the combination of NDVI and NBR2 thresholding. The same PCA-LDA classification was applied against the Sentinel-2 bare soil spectra (Table 2), and the scatterplot between LD1 and LD2 shows clear distinction among the erosion classes, as further indicated by the LD1 function values distributed in segregated ranges, explaining 94% of the variance (Figure 6).

Figure 7 shows the Sentinel-2-based soil spectra of the three erosion classes. Apart from the similarities between Sentinel-2 and laboratory-based curves in terms of the highest absolute reflectance values in the visible region (B2, B3, B4 bands) for the Severe class and the varying slopes between B8 and B11 among the three classes, notable differences in Sentinel-2-based curves include: (1) equally high separability between Low and Moderate classes in the visible region, (2) lack of spectral detail with only B11 and B12 in the SWIR region and (3) worsened separability of the absorption features (B3, B4) after CR treatment. Furthermore, the classification scheme composed of three indices was tested against the Sentinel-2 spectra. As the 1350 nm wavelength used for calculating the slope over the NIR range was not available for Sentinel-2 bands, B11 at approximately 1610 nm was used instead. It can be seen that, consistent with the interpretation of Sentinel-2 spectra above, the mean reflectance from B2 to B4 was an even better classifier of erosion intensity for Sentinel-2 than it was for laboratory-based data. After all, a complete separation between Low and Moderate classes could be achieved using a threshold of 0.09 (Figure 8). The usefulness of the slope between two NIR wavelengths to separate the Severe class was also proven for Sentinel-2 data. A threshold of 1.3 × 10⁻⁴ could allow an excellent discrimination of Severe class from the other two. Lastly, the absorption features at B3 after CR correction did not show clear separability, especially between Moderate and Severe classes.

3.3.2. Soil Erosion Mapping and Evaluation

Summarizing the results from Figure 8, the two thresholds (0.09 for the mean reflectance from B2 to B4, and 1.3 × 10⁻⁴ for the slope between B8 and B11) were applied to achieve pixel-wise erosion classification based on the bare soil composite, and a high-resolution (10 m) soil erosion intensity map with Low, Moderate and Severe erosion classes was produced (Figure 9). Specifically, the Moderate erosion class occupied the highest percentage of the cropland area with 65.4%, followed by the Low class at 18.7%, distributed mostly in concave areas, where the erosion and deposition processes interacted and led to a low net erosion and, sometimes, net deposition phenomena (Figure 2). Soils that suffered Severe erosion, also called “erosion hotspots”, occupied 15.9% of the cropland area. Comparing the NDVI map for June 2021 to the soil erosion map shows that the distribution of NDVI shifted toward higher values with decreasing erosion intensity (Figure 10), indirectly supporting the validity of the spectra-based soil erosion classification approach, as the spatial variability in crop growth responded in accordance with the predicted erosion classes. Zooms of the two areas further depicted the promising capability of such an approach to not only capture the spatial variability of soil erosion at high resolution but to also be used in assessments of soil productivity.

4. Discussion

4.1. Erosion Characteristics at the Sampled Slope Positions

Successful classification and mapping of soil erosion intensity via spectral imaging requires the establishment of a ground-truth training set for model development. For this purpose, three representative slope positions were selected, and their erosion characteristics were analyzed from the aspects of topographic features, SOC contents and ¹³⁷Cs inventories. Apart from the accelerated SOC loss in Mid-slope positions, which demonstrated the more severe truncation of topsoil materials by erosion than the other two slope positions (Figure 2), the use of radionuclide ¹³⁷Cs tracing allowed a quantitative estimation and comparison of net erosion rates, providing evidence on the validity of the selected sampling groups. A wide range of variations in the estimated erosion rates was found among the three slope positions, implying that varying the degree of erosion-induced soil redistribution would result in a highly heterogenous distribution of key soil properties at the hillslope scale, as evidenced by recent studies [38,40]. In particular, areas occupied by the Severe erosion class (36.93 ± 12.91 t ha⁻¹ yr⁻¹, equivalent to 3.1 ± 1.1 mm yr⁻¹ soil loss assuming a bulk density of 1.2 g cm⁻³), as a result of multi-faceted erosion forms consisting primarily of water, tillage and wind erosion, represent the “hotspots” that require targeted conservation measures.

Conventional soil erosion modeling approaches remain limited in their ability to detect these hotspots, especially considering the lack of high-resolution data of DEM and other USLE factors needed for such modeling framework to be most effective [10]. Additionally, different from USLE-type estimates, which give gross water erosion rates, the erosion rates converted from ¹³⁷Cs inventory represent net erosion rates as a result of multiple erosion forms, i.e., water, tillage and wind. Previous research works in various black soil regions of the world (NE China, US Corn Belt) have stressed the important role of tillage erosion in causing the topsoil mobilization and redistribution [50,51]. Hence, soil erosion mapping by imaging spectroscopy links soil spectral features to erosion-induced local variations in soil properties, thus forming the theoretical basis for a spatially explicit approach to detecting soil erosion hotspots at a high spatial resolution.

4.2. Soil Erosion Mapping Driven by Sentinel-2 Imagery

Spectra-based soil erosion classification could only be successful if there existed good spectral separability among soils under varying degrees of erosion intensity. In this study, a spectral classification scheme, consisting of three spectral indices with pre-identified thresholds, was constructed to directly link soil spectral features with different erosion classes (Figure 5 and Figure 8). More specifically, the first component of the classification scheme took advantage of the mean reflectance across the visible range, since soils depleted of SOC are typically associated with a light color and have a high albedo [52]. Therefore, soils in the Severe class that suffered extensive loss of the A-horizon could be readily detected as erosion hotspots due to their light color associated with accelerated SOC losses. This is in accordance with previous studies, which stressed the importance of visible bands in SOC prediction [53,54]. Additionally, soils in the Low erosion class showed a distinct spectral signature in the visible region of Sentinel-2 spectra but not in the laboratory spectra, suggesting that relatively higher soil moisture content generally found in concavities and/or slope bottoms led to a decrease in albedo [55]. This albedo decrease became the key factor, which contributed to the detection of hardly eroded soils. The remaining two components of the classification scheme used the spectral slope between two wavelengths in the NIR region and enhanced absorption feature after CR treatment as measures of difference in soil composition among the erosion classes. Both indices have been frequently used in spectra-based soil classification studies [46,56,57] and were also proved effective in this study. Lastly, it is worth noting that the two SWIR bands (B11 and B12) were not identified as meaningful features for separating different erosion classes, perhaps due to the reduced resolution of Sentinel-2 bands in this spectral region.

Compared to previous remote-sensing studies that used bare soil spectra for erosion classification, Schmid et al. (2016) [28] used raw airborne hyperspectral data to define the endmembers of different soil erosion stages, which were used to train a support vector machine classifier for erosion mapping. Žížala et al. (2019) [35] cross-compared raw Sentinel-2 reflectance among different erosion classes and did not find good spectral separability, except for the severely eroded class. They adopted an unsupervised classification method to map soil erosion with moderate accuracy, which was then improved by visual interpretation of erosion features from orthoimages. Here, apart from leveraging the raw reflectance, we applied spectral index thresholding to sequentially differentiate areas under varying erosion influence. This approach could further expand the potential of remote-sensing methodologies to effectively detect soil erosion and degradation.

Furthermore, the spectra-based erosion classification scheme, initially established based on distinctive spectral features from laboratory spectra, was consistent in its ability to also separate erosion classes from Sentinel-2-derived spectra. Its transferability from laboratory hyperspectral data to satellite-based multispectral images indicates that our spectral classification approach was robust and that bare soil pixels extracted from the multitemporal Sentinel-2 composite largely reflected the true soil spectral information similar to that captured under laboratory conditions. Recent studies on using Sentinel-2 imagery for soil mapping applications have facilitated methodological developments to create multitemporal composites comprising bare soil pixels of high purity [30]. In this study, NDVI and NBR2 thresholding was applied to extract bare soil pixels of minimal disturbance from green and dry vegetation for three single-date images, each of which was selected within the maize sowing period in April–May of 2019–2021, when croplands in the catchment were prepared for seedbeds, and soil surface was at its optimal condition. The multitemporal composite created from three single-date images thus ensured maximum coverage of bare croplands while also maintaining the quality of bare soil spectral information.

Pixel-wise erosion classification based on the multitemporal bare soil composite produced a 10 m resolution soil erosion map (Figure 9), which precisely located the erosion hotspots and areas that suffered less severe erosion. This type of high-resolution erosion map produced by Sentinel-2 remote sensing could be used to constrain and further develop soil erosion models, to provide spatially explicit soil degradation information for sustainable cropland management and to help assess the carbon sequestration potential in degraded soils and how soil erosion impacts crop productivity [50]. For instance, the zoomed-in figures (Figure 10) clearly depicted that localized erosion hotspots were associated with low NDVI values, demonstrating the negative impact of soil erosion on crop growth from a spatial perspective. However, it should be noted that this agreement between soil erosion pattern and crop index may not be stable from year to year, as inter-annual rainfall variability and fertilization inputs could, at least to some extent, mask the contrast in crop biomass and productivity among soils under different erosion influences [2,58].

4.3. The Way Forward

The primary goal of this study was to develop a remote-sensing method for precise detection of soil erosion hotspots at the catchment scale. The ground-truth dataset used to verify the capability of the proposed spectra-based classification scheme was only composed of three erosion classes, with particular focus given to the separability between classes. However, this ground-truth dataset consisted of three erosion classes, which were derived from samples taken from three representative slope positions. Future work should therefore develop a standardized soil erosion classification criterion, preferably based on quantitative erosion rates. Moreover, the dataset did not allow further investigation into the within-class variability in erosion intensity, particularly against the Moderate erosion class, which occupied more than 60% of the cropland area (Figure 9). The lack of detailed and accurate classification for moderately eroded areas appears to be a consistent issue in spectra-based soil erosion mapping investigations, such as the ones by Schmid et al. (2016), Žížala et al. (2017) and Žížala et al. (2019), all reporting reduced classification accuracy for moderately eroded areas in comparison to severely eroded hotspots [22,28,35]. Future studies should therefore investigate the possibility of the spectra-based classification approach to further categorize soils under moderate erosion. This requires either a more comprehensive ground-truth dataset for supervised classification or an unsupervised approach to separate the Moderate classes into sub-classes, followed by post hoc validation.

Another issue associated with the proposed approach in this study is its genericity for application in other areas or larger spatial scales. Due to the empirical nature of the spectral classification scheme and its thresholds, the application of this type of approach is likely to be area-specific, and new classification criteria will need to be calibrated based on specific training data. Notwithstanding, the principle of using soil spectral characteristics to discriminate soils of varying erosion intensities should be applicable as long as erosion-induced soil redistribution causes variations in soil spectral signature. Finally, soil erosion mapping aided by imaging spectroscopy has mostly been limited to field and catchment scales. One key reason could be that the spectral response to different soil types in a large heterogenous space would interfere with erosion-induced variations in soil spectral signals. Future upscaling studies should therefore prioritize the delineation of homogenous units and develop tailored classification models, possibly also with the incorporation of topographic [50] and vegetation metrics [58] as covariates. A spectra-based large-scale soil erosion remote-sensing approach could provide important ability to quantify the “proportion of land that is degraded over total land area” (SDG Indicator 15.3.1) in a spatially explicit manner, thereby facilitating the detection of localized soil degradation hotspots associated with accelerated SOC losses.

5. Conclusions

The potential of multitemporal Sentinel-2 remote sensing to detect erosion hotspots was tested at the catchment scale (ca. 46.20 km²) in the black soil region of northeast China. A ground-truth dataset was established to include soils collected at the Summit, Mid-slope and Foot-slope positions, corresponding to Moderate, Severe and Low erosion intensity classes due to their differences in topographic features, net erosion rates converted from ¹³⁷Cs inventory and SOC contents. Investigations into both laboratory and Sentinel-2-based soil spectral data showed that soils among the three erosion classes displayed distinctive spectral features due to erosion-induced shifts in soil albedo and biochemical composition, particularly in severely eroded areas, where substantial loss of the A-horizon was evident.

PCA-LDA demonstrated clear inter-class spectral separability to discriminate soils under varying erosion influence, thus enabling the development of a spectral classification scheme consisting of identified spectral index (mean reflectance over B2–B4, and slope between B8 and B11 for Sentinel-2 bands) thresholds for pixel-wise soil erosion mapping using the Sentinel-2 bare soil composite. The produced high-resolution soil erosion map allowed a close-up analysis of the relationship between soil erosion and crop productivity, highlighting the promising potential of our proposed approach for sustainable cropland management in the black soil region. Future studies should further test the transferability of this approach to other areas and larger spatial scales.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs15051402/s1, Figure S1: ¹³⁷Cs inventory at Summit, Mid-slope and Foot-slope positions. Red dots indicate average value per group; Figure S2: Vertical distribution of ¹³⁷Cs inventory of the three soil pits at Summit positions.

Author Contributions

Conceptualization, P.S. and B.v.W.; methodology, L.Q., P.S., K.D., K.V.O. and B.v.W.; formal analysis, L.Q. and P.S.; data curation, L.Q., Q.S. and H.Y.; writing—original draft preparation, L.Q. and P.S.; writing—review and editing, K.D. and B.v.W.; supervision, P.S. and B.v.W.; funding acquisition, P.S. and B.v.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Natural Science Foundation of Jilin Province, China, grant number 20230101079JC, and the National Natural Science Foundation of China, grant number 41807059. L.Q. is supported by the China Scholarship Council for a one-year joint Ph.D. study in Belgium (No. 202106170101). P.S. is a postdoc fellow of the Fonds de la Recherche Scientifique—FNRS (Belgium).

Data Availability Statement

Data used in this study can be made available upon reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

References

Prăvălie, R.; Patriche, C.; Borrelli, P.; Panagos, P.; Roșca, B.; Dumitraşcu, M.; Nita, I.-A.; Săvulescu, I.; Birsan, M.-V.; Bandoc, G. Arable lands under the pressure of multiple land degradation processes. A global perspective. Environ. Res. 2021, 194, 110697. [Google Scholar] [CrossRef] [PubMed]
Le, Q.B.; Nkonya, E.; Mirzabaev, A. Biomass Productivity-Based Mapping of Global Land Degradation Hotspots. In Economics of Land Degradation and Improvement—A Global Assessment for Sustainable Development; Nkonya, E., Mirzabaev, A., von Braun, J., Eds.; Springer International Publishing: Cham, Switzerland, 2016; pp. 55–84. [Google Scholar]
Bakker, M.M.; Govers, G.; Rounsevell, M.D.A. The crop productivity–erosion relationship: An analysis based on experimental work. Catena 2004, 57, 55–76. [Google Scholar] [CrossRef]
Amelung, W.; Bossio, D.; de Vries, W.; Kögel-Knabner, I.; Lehmann, J.; Amundson, R.; Bol, R.; Collins, C.; Lal, R.; Leifeld, J.; et al. Towards a global-scale soil climate mitigation strategy. Nat. Commun. 2020, 11, 5427. [Google Scholar] [CrossRef] [PubMed]
Lal, R. Soil carbon sequestration to mitigate climate change. Geoderma 2004, 123, 1–22. [Google Scholar] [CrossRef]
Cowie, A.L.; Orr, B.J.; Castillo Sanchez, V.M.; Chasek, P.; Crossman, N.D.; Erlewein, A.; Louwagie, G.; Maron, M.; Metternicht, G.I.; Minelli, S.; et al. Land in balance: The scientific conceptual framework for Land Degradation Neutrality. Environ. Sci. Policy 2018, 79, 25–35. [Google Scholar] [CrossRef]
Liu, B.; Xie, Y.; Li, Z.; Liang, Y.; Zhang, W.; Fu, S.; Yin, S.; Wei, X.; Zhang, K.; Wang, Z.; et al. The assessment of soil loss by water erosion in China. Int. Soil Water Conserv. Res. 2020, 8, 430–439. [Google Scholar] [CrossRef]
Borrelli, P.; Robinson, D.A.; Panagos, P.; Lugato, E.; Yang, J.E.; Alewell, C.; Wuepper, D.; Montanarella, L.; Ballabio, C. Land use and climate change impacts on global soil erosion by water (2015–2070). Proc. Natl. Acad. Sci. USA 2020, 117, 21994–22001. [Google Scholar] [CrossRef]
Borrelli, P.; Alewell, C.; Alvarez, P.; Anache, J.A.A.; Baartman, J.; Ballabio, C.; Bezak, N.; Biddoccu, M.; Cerdà, A.; Chalise, D.; et al. Soil erosion modelling: A global review and statistical analysis. Sci. Total Environ. 2021, 780, 146494. [Google Scholar] [CrossRef]
Alewell, C.; Borrelli, P.; Meusburger, K.; Panagos, P. Using the USLE: Chances, challenges and limitations of soil erosion modelling. Int. Soil Water Conserv. Res. 2019, 7, 203–225. [Google Scholar] [CrossRef]
Van Oost, K.; Govers, G.; de Alba, S.; Quine, T.A. Tillage erosion: A review of controlling factors and implications for soil quality. Prog. Phys. Geog. 2006, 30, 443–466. [Google Scholar] [CrossRef] [Green Version]
Van Oost, K.; Govers, G.; Desmet, P. Evaluating the effects of changes in landscape structure on soil erosion by water and tillage. Landsc. Ecol. 2000, 15, 577–589. [Google Scholar] [CrossRef]
Morgan, R.P.C.; Quinton, J.N.; Smith, R.E.; Govers, G.; Poesen, J.W.A.; Auerswald, K.; Chisci, G.; Torri, D.; Styczen, M.E. The European Soil Erosion Model (EUROSEM): A dynamic approach for predicting sediment transport from fields and small catchments. Earth Surf. Process. Landf. 1998, 23, 527–544. [Google Scholar] [CrossRef]
Jetten, V.; Govers, G.; Hessel, R. Erosion models: Quality of spatial predictions. Hydrol. Process. 2003, 17, 887–900. [Google Scholar] [CrossRef]
de Vente, J.; Poesen, J.; Verstraeten, G.; Govers, G.; Vanmaercke, M.; Van Rompaey, A.; Arabkhedri, M.; Boix-Fayos, C. Predicting soil erosion and sediment yield at regional scales: Where do we stand? Earth-Sci. Rev. 2013, 127, 16–29. [Google Scholar] [CrossRef]
Mabit, L.; Benmansour, M.; Walling, D.E. Comparative advantages and limitations of the fallout radionuclides 137Cs, 210Pbex and 7Be for assessing soil erosion and sedimentation. J. Environ. Radioact. 2008, 99, 1799–1807. [Google Scholar] [CrossRef]
Zhang, X.C.; Zhang, G.H.; Wei, X. How to make 137Cs erosion estimation more useful: An uncertainty perspective. Geoderma 2015, 239, 186–194. [Google Scholar] [CrossRef]
Walling, D.E.; Quine, T.A. Use of 137cs Measurements to Investigate Soil-Erosion on Arable Fields in the Uk—Potential Applications and Limitations. J. Soil. Sci. 1991, 42, 147–165. [Google Scholar] [CrossRef]
Vågen, T.-G.; Winowiecki, L.A. Predicting the Spatial Distribution and Severity of Soil Erosion in the Global Tropics using Satellite Remote Sensing. Remote Sens. 2019, 11, 1800. [Google Scholar] [CrossRef] [Green Version]
Tziolas, N.; Tsakiridis, N.; Chabrillat, S.; Demattê, J.A.M.; Ben-Dor, E.; Gholizadeh, A.; Zalidis, G.; van Wesemael, B. Earth Observation Data-Driven Cropland Soil Monitoring: A Review. Remote Sens. 2021, 13, 4439. [Google Scholar] [CrossRef]
Chabrillat, S.; Ben-Dor, E.; Cierniewski, J.; Gomez, C.; Schmid, T.; van Wesemael, B. Imaging Spectroscopy for Soil Mapping and Monitoring. Surv. Geophys. 2019, 40, 361–399. [Google Scholar] [CrossRef] [Green Version]
Žížala, D.; Zádorová, T.; Kapička, J. Assessment of Soil Degradation by Erosion Based on Analysis of Soil Properties Using Aerial Hyperspectral Images and Ancillary Data, Czech Republic. Remote Sens. 2017, 9, 28. [Google Scholar] [CrossRef] [Green Version]
Lin, C.; Zhou, S.-L.; Wu, S.-H.; Zhu, Q.; Dang, Q. Spectral response of different eroded soils in subtropical china: A case study in Changting County, China. J. Mt. Sci. 2014, 11, 697–707. [Google Scholar] [CrossRef]
Conforti, M.; Buttafuoco, G.; Leone, A.P.; Aucelli, P.P.C.; Robustelli, G.; Scarciglia, F. Studying the relationship between water-induced soil erosion and soil organic matter using Vis-NIR spectroscopy and geomorphological analysis: A case study in southern Italy. Catena 2013, 110, 44–58. [Google Scholar] [CrossRef]
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.; et al. A global spectral library to characterize the world’s soil. Earth-Sci. Rev. 2016, 155, 198–230. [Google Scholar] [CrossRef] [Green Version]
Nocita, M.; Stevens, A.; van Wesemael, B.; Aitkenhead, M.; Bachmann, M.; Barthes, B.; Ben Dor, E.; Brown, D.J.; Clairotte, M.; Csorba, A.; et al. Soil Spectroscopy: An Alternative to Wet Chemistry for Soil Monitoring. Adv. Agron. 2015, 132, 139–159. [Google Scholar] [CrossRef]
Soriano-Disla, J.M.; Janik, L.J.; Rossel, R.A.V.; Macdonald, L.M.; McLaughlin, M.J. The Performance of Visible, Near-, and Mid-Infrared Reflectance Spectroscopy for Prediction of Soil Physical, Chemical, and Biological Properties. Appl. Spectrosc. Rev. 2014, 49, 139–186. [Google Scholar] [CrossRef]
Schmid, T.; Rodriguez-Rastrero, M.; Escribano, P.; Palacios-Orueta, A.; Ben-Dor, E.; Plaza, A.; Milewski, R.; Huesca, M.; Bracken, A.; Cicueudez, V.; et al. Characterization of Soil Erosion Indicators Using Hyperspectral Data From a Mediterranean Rainfed Cultivated Region. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016, 9, 845–860. [Google Scholar] [CrossRef]
Bracken, A.; Coburn, C.; Staenz, K.; Rochdi, N.; Segl, K.; Chabrillat, S.; Schmid, T. Detecting soil erosion in semi-arid Mediterranean environments using simulated EnMAP data. Geoderma 2019, 340, 164–174. [Google Scholar] [CrossRef] [Green Version]
Vaudour, E.; Gholizadeh, A.; Castaldi, F.; Saberioon, M.; Borůvka, L.; Urbina-Salazar, D.; Fouad, Y.; Arrouays, D.; Richer-de-Forges, A.C.; Biney, J.; et al. Satellite Imagery to Map Topsoil Organic Carbon Content over Cultivated Areas: An Overview. Remote Sens. 2022, 14, 2917. [Google Scholar] [CrossRef]
Helder, D.; Markham, B.; Morfitt, R.; Storey, J.; Barsi, J.; Gascon, F.; Clerc, S.; LaFrance, B.; Masek, J.; Roy, D.P.; et al. Observations and Recommendations for the Calibration of Landsat 8 OLI and Sentinel 2 MSI for Improved Data Interoperability. Remote Sens. 2018, 10, 1340. [Google Scholar] [CrossRef] [Green Version]
Dvorakova, K.; Shi, P.; Limbourg, Q.; van Wesemael, B. Soil Organic Carbon Mapping from Remote Sensing: The Effect of Crop Residues. Remote Sens. 2020, 12, 1913. [Google Scholar] [CrossRef]
Vaudour, E.; Gomez, C.; Lagacherie, P.; Loiseau, T.; Baghdadi, N.; Urbina-Salazar, D.; Loubet, B.; Arrouays, D. Temporal mosaicking approaches of Sentinel-2 images for extending topsoil organic carbon content mapping in croplands. Int. J. Appl. Earth Obs. Geoinf. 2021, 96, 102277. [Google Scholar] [CrossRef]
Demattê, J.A.M.; Fongaro, C.T.; Rizzo, R.; Safanelli, J.L. Geospatial Soil Sensing System (GEOS3): A powerful data mining procedure to retrieve soil spectral reflectance from satellite images. Remote Sens. Environ. 2018, 212, 161–175. [Google Scholar] [CrossRef]
Žížala, D.; Juřicová, A.; Zádorová, T.; Zelenková, K.; Minařík, R. Mapping soil degradation using remote sensing data and ancillary data: South-East Moravia, Czech Republic. Eur. J. Remote Sens. 2019, 52, 108–122. [Google Scholar] [CrossRef] [Green Version]
Li, W.; Wang, D.; Li, Y.; Zhu, Y.; Wang, J.; Ma, J. A multi-faceted, location-specific assessment of land degradation threats to peri-urban agriculture at a traditional grain base in northeastern China. J. Environ. Manag. 2020, 271, 111000. [Google Scholar] [CrossRef] [PubMed]
Zhao, P.; Li, S.; Wang, E.; Chen, X.; Deng, J.; Zhao, Y. Tillage erosion and its effect on spatial variations of soil organic carbon in the black soil region of China. Soil Tillage Res. 2018, 178, 72–81. [Google Scholar] [CrossRef]
Hu, Y.; Li, X.; Guo, S.; Gao, X.; Ou, X.; Liu, B. On-site soil dislocation and localized CNP degradation: The real erosion risk faced by sloped cropland in northeastern China. Agric. Ecosyst. Environ. 2020, 302, 107088. [Google Scholar] [CrossRef]
IUSS Working Group WRB. World Reference Base for Soil Resources 2014, Update 2015. International Soil Classification System for Naming Soils and Creating Legends for Soil Maps. FAO: Rome, Italy, 2015. [Google Scholar]
Zhu, Y.; Wang, D.; Wang, X.; Li, W.; Shi, P. Aggregate-associated soil organic carbon dynamics as affected by erosion and deposition along contrasting hillslopes in the Chinese Corn Belt. Catena 2021, 199, 105106. [Google Scholar] [CrossRef]
Zhang, X.C.; Polyakov, V.O.; Liu, B.Y.; Nearing, M.A. Quantifying geostatistical properties of 137Cs and 210Pbex at small scales for improving sampling design and soil erosion estimation. Geoderma 2019, 334, 155–164. [Google Scholar] [CrossRef]
Li, Y.; Zhang, Q.W.; Reicosky, D.C.; Bai, L.Y.; Lindstrom, M.J.; Li, L. Using 137Cs and 210Pbex for quantifying soil organic carbon redistribution affected by intensive tillage on steep slopes. Soil Tillage Res. 2006, 86, 176–184. [Google Scholar] [CrossRef]
Zhang, X.; Li, S.; Wang, C.; Tan, W.; Zhao, Q.; Zhang, Y.; Yan, M.; Liu, Y.; Jiang, J.; Xiao, J.; et al. Study on the sediment sources in a small watershed of the Loess Plaetau using the cesium-137 method. Chin. Sci. Bull. 1989, 34, 210–213. [Google Scholar]
Yan, B.; Tang, J. Study on reference cesium-137 inventory of black soil in Northeast China. J. Soil Water Conserv. 2004, 18, 33–36. [Google Scholar] [CrossRef]
Shi, P.; Castaldi, F.; van Wesemael, B.; Van Oost, K. Vis-NIR spectroscopic assessment of soil aggregate stability and aggregate size distribution in the Belgian Loam Belt. Geoderma 2020, 357, 113958. [Google Scholar] [CrossRef]
Chamizo, S.; Stevens, A.; Cantón, Y.; Miralles, I.; Domingo, F.; Van Wesemael, B. Discriminating soil crust type, development stage and degree of disturbance in semiarid environments from their spectral characteristics. Eur. J. Soil Sci. 2012, 63, 42–53. [Google Scholar] [CrossRef]
Clark, R.N.; Roush, T.L. Reflectance spectroscopy: Quantitative analysis techniques for remote sensing applications. J. Geophys. Res. Solid Earth 1984, 89, 6329–6340. [Google Scholar] [CrossRef]
Shi, P.; Six, J.; Sila, A.; Vanlauwe, B.; Van Oost, K. Towards spatially continuous mapping of soil organic carbon in croplands using multitemporal Sentinel-2 remote sensing. ISPRS J. Photogramm. Remote Sens. 2022, 193, 187–199. [Google Scholar] [CrossRef]
Castaldi, F.; Hueni, A.; Chabrillat, S.; Ward, K.; Buttafuoco, G.; Bomans, B.; Vreys, K.; Brell, M.; van Wesemael, B. Evaluating the capability of the Sentinel 2 data for soil organic carbon prediction in croplands. ISPRS J. Photogramm. Remote Sens. 2019, 147, 267–282. [Google Scholar] [CrossRef]
Thaler, E.A.; Larsen, I.J.; Yu, Q. The extent of soil loss across the US Corn Belt. Proc. Natl. Acad. Sci. USA 2021, 118, e1922375118. [Google Scholar] [CrossRef]
Zhao, Y.; Wang, M.; Hu, S.; Zhang, X.; Ouyang, Z.; Zhang, G.; Huang, B.; Zhao, S.; Wu, J.; Xie, D.; et al. Economics- and policy-driven organic carbon input enhancement dominates soil organic carbon accumulation in Chinese croplands. Proc. Natl. Acad. Sci. USA 2018, 115, 4045–4050. [Google Scholar] [CrossRef] [Green Version]
Ben-Dor, E.; Inbar, Y.; Chen, Y. 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. 1997, 61, 1–15. [Google Scholar] [CrossRef]
Shi, P.; Castaldi, F.; van Wesemael, B.; Van Oost, K. Large-Scale, High-Resolution Mapping of Soil Aggregate Stability in Croplands Using APEX Hyperspectral Imagery. Remote Sens. 2020, 12, 666. [Google Scholar] [CrossRef] [Green Version]
Castaldi, F.; Chabrillat, S.; Chartin, C.; Genot, V.; Jones, A.R.; van Wesemael, B. Estimation of soil organic carbon in arable soil in Belgium and Luxembourg with the LUCAS topsoil database. Eur. J. Soil Sci. 2018, 69, 592–603. [Google Scholar] [CrossRef]
Nocita, M.; Stevens, A.; Noon, C.; van Wesemael, B. Prediction of soil organic carbon for different levels of soil moisture using Vis-NIR spectroscopy. Geoderma 2013, 199, 37–42. [Google Scholar] [CrossRef]
Zhang, X.; Liu, H.; Zhang, X.; Yu, S.; Dou, X.; Xie, Y.; Wang, N. Allocate soil individuals to soil classes with topsoil spectral characteristics and decision trees. Geoderma 2018, 320, 12–22. [Google Scholar] [CrossRef]
Ogen, Y.; Goldshleger, N.; Ben-Dor, E. 3D spectral analysis in the VNIR–SWIR spectral region as a tool for soil classification. Geoderma 2017, 302, 100–110. [Google Scholar] [CrossRef]
Öttl, L.K.; Wilken, F.; Auerswald, K.; Sommer, M.; Wehrhan, M.; Fiener, P. Tillage erosion as an important driver of in-field biomass patterns in an intensively used hummocky landscape. Land Degrad. Dev. 2021, 32, 3077–3091. [Google Scholar] [CrossRef]

Figure 1. (a,b) Spatial distribution of the sampling points overlaid on the digital elevation model of the catchment in northeast China; (c,d) a Sentinel-2 true color image (13 May 2021) of the catchment showing the large percentage of exposed bare soils, and zoom-in areas depicting the selection of sampling points; (e,f) representative sampling locations along typical slope profiles with average length at ca. 300 m and their corresponding slope degrees. Different letters indicate significant differences among the three groups (p < 0.05).

Figure 2. (a) Topographic position index (TPI), (b) estimated net soil erosion rates converted from ¹³⁷Cs inventory, (c) soil organic carbon (SOC) content and (d) C/N ratio at Summit, Mid-slope and Foot-slope positions. Different letters indicate significant differences among the three groups (p < 0.05).

Figure 3. Classification of three slope positions based on linear discriminant analysis (LDA) of PC scores arising from laboratory-based VNIR spectra. Subsoil spectra of Summit positions were included to investigate any spectral similarity to the topsoil of Mid-slope positions. Histograms on the right show the distribution of the first LD function’s value for the four sample groups.

Figure 4. Laboratory-based mean spectra of topsoil in three soil erosion intensity classes. (a) Raw spectra and (b) continuum-removed reflectance. Shaded areas represent standard deviations from the mean.

Figure 5. Boxplots of laboratory spectral indices of topsoil used for spectral separation of erosion intensity classes. (a) Mean reflectance over the visible region (400–780 nm), (b) slope between 800 and 1350 nm, and (c) continuum-removed reflectance value at 670 nm. Different letters indicate significant differences among the three groups (p < 0.05).

Figure 6. Classification of three erosion intensity classes based on linear discriminant analysis (LDA) of PC scores derived from Sentinel-2 bare soil spectra. Histograms on the right show the distribution of the first LD function’s value for the three classes.

Figure 7. Sentinel-2-based mean spectra of three soil erosion intensity classes. (a) Raw spectra and (b) continuum-removed reflectance. Shaded areas represent standard deviations from the mean. Sentinel-2 band width and positions are indicated in the left figure.

Figure 8. Boxplots of Sentinel-2 spectral indices used for spectral separation of erosion intensity classes. (a) Mean reflectance over the three visible bands (B2, B3, B4), (b) slope between B8 and B11, and (c) continuum-removed reflectance value at 670 nm. Different letters indicate significant differences among the three groups (p < 0.05).

Figure 9. Soil erosion intensity map at 10 m resolution. The classification criteria were based on the results from Figure 8.

Figure 10. (a) Density plots of NDVI in June 2021 for the three erosion intensity classes within the cropland extent; (b,c) zoomed-in areas with detailed soil erosion pattern and (d,e) the corresponding field-scale NDVI maps to depict the relation of crop variability with erosion-induced variations in soil productivity.

Table 1. Confusion matrix of laboratory-based spectral classification of erosion intensity using a combination of principal component and linear discriminant analysis (PCA-LDA).

Prediction	Observation
Prediction	Low	Moderate	Severe	Total	User’s Accuracy (%)
Low	22	0	0 1	22 23	100.00
Moderate	3	19	0 1	22 23	82.61
Severe	0	0	27 28	27 72	100.00
Total	25	19	27 28	27 72
Producer’s accuracy (%)	88.00	100.00	96.43
Overall accuracy (%)	94.00
Kappa coefficient	0.92

Table 2. Confusion matrix of Sentinel-2-based spectral classification of erosion intensity using a combination of principal component and linear discriminant analysis (PCA-LDA).

Prediction	Observation
Prediction	Low	Moderate	Severe	Total	User’s Accuracy (%)
Low	25	0	0 0	25 19	100.00
Moderate	0	19	0 0	25 19	100.00
Severe	0	0	28 28	28 72	100.00
Total	25	19	28 28	28 72
Producer’s accuracy (%)	100.00	100.00	100.00
Overall accuracy (%)	100.00
Kappa coefficient	1.00

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Qi, L.; Shi, P.; Dvorakova, K.; Van Oost, K.; Sun, Q.; Yu, H.; van Wesemael, B. Detection of Soil Erosion Hotspots in the Croplands of a Typical Black Soil Region in Northeast China: Insights from Sentinel-2 Multispectral Remote Sensing. Remote Sens. 2023, 15, 1402. https://doi.org/10.3390/rs15051402

AMA Style

Qi L, Shi P, Dvorakova K, Van Oost K, Sun Q, Yu H, van Wesemael B. Detection of Soil Erosion Hotspots in the Croplands of a Typical Black Soil Region in Northeast China: Insights from Sentinel-2 Multispectral Remote Sensing. Remote Sensing. 2023; 15(5):1402. https://doi.org/10.3390/rs15051402

Chicago/Turabian Style

Qi, Lulu, Pu Shi, Klara Dvorakova, Kristof Van Oost, Qi Sun, Hanqing Yu, and Bas van Wesemael. 2023. "Detection of Soil Erosion Hotspots in the Croplands of a Typical Black Soil Region in Northeast China: Insights from Sentinel-2 Multispectral Remote Sensing" Remote Sensing 15, no. 5: 1402. https://doi.org/10.3390/rs15051402

APA Style

Qi, L., Shi, P., Dvorakova, K., Van Oost, K., Sun, Q., Yu, H., & van Wesemael, B. (2023). Detection of Soil Erosion Hotspots in the Croplands of a Typical Black Soil Region in Northeast China: Insights from Sentinel-2 Multispectral Remote Sensing. Remote Sensing, 15(5), 1402. https://doi.org/10.3390/rs15051402

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu