Soil Organic Carbon Mapping Using LUCAS Topsoil Database and Sentinel-2 Data: An Approach to Reduce Soil Moisture and Crop Residue E ﬀ ects

: Soil organic carbon (SOC) loss is one of the main causes of soil degradation in croplands. Thus, spatial and temporal monitoring of SOC is extremely important, both from the environmental and economic perspective. In this regard, the high temporal, spatial, and spectral resolution of the Sentinel-2 data can be exploited for monitoring SOC contents in the topsoil of croplands. In this study, we aim to test the e ﬀ ect of the threshold for a spectral index linked to soil moisture and crop residues on the performance of SOC prediction models using the Multi-Spectral Instrument (MSI) Sentinel-2 and the European Land Use / cover Area frame Statistical survey (LUCAS) topsoil database. The LUCAS spectral data resampled according to MSI / Sentinel-2 bands, which were used to build SOC prediction models combining pairs of the bands. The SOC models were applied to a Sentinel-2 image acquired in North-Eastern Germany after removing the pixels characterized by clouds and green vegetation. Then, we tested di ﬀ erent thresholds of the Normalized Burn Ratio 2 (NBR2) index in order to mask moist soil pixels and those with dry vegetation and crop residues. The model accuracy was tested on an independent validation database and the best ratio of performance to deviation (RPD) was obtained using the average between bands B6 and B5 (Red-Edge Carbon Index: RE-CI) (RPD: 4.4) and between B4 and B5 (Red-Red-Edge Carbon Index: RRE-CI) (RPD: 2.9) for a very low NBR2 threshold (0.05). Employing a higher NBR2 tolerance (higher NBR2 values), the mapped area increases to the detriment of the validation accuracy. The proposed approach allowed us to accurately map SOC over a large area exploiting the LUCAS spectral library and, thus, avoid a new ad hoc ﬁeld campaign. Moreover, the threshold for selecting the bare soil pixels can be tuned, according to the goal of the survey. The quality of the SOC map for each tolerance level can be judged based on the ﬁgures of merit of the model.


Introduction
Since the onset of agriculture, soils have lost their organic carbon to such an extent that the soil functions of many croplands are threatened. There is strong demand for mapping and monitoring critical soil properties and, in particular, soil organic carbon (SOC). After all, SOC stock is one of the three sub-indicators for the Sustainable Development Goal (SDG) 15.3.1 defining the proportion of land and airborne data for SOC mapping without spectral transfer function between laboratory and remote spectral data. While this approach proved to be successful in a relatively small pilot area restricted to bare cropland soils, the main challenge for the SOC estimation from an entire Sentinel-2 tile (i.e., 100 by 100 km) is the minimization of disturbing factors, such as vegetation, residues, roughness, and soil moisture. Clearly, the search for ideal soil conditions (i.e., croplands in seedbed condition) reduces the time window to get a useful image in temperate regions. Rogge et al. [12] used temporal series of Landsat and Sentinel-2 scenes. Their approach masks photosynthetically active vegetation using normalized difference vegetation index (NDVI) as a vegetation index. This map produces, among other things, an average reflectance of bare soil pixels in croplands. Given the position of most of the MSI/Sentinel-2 bands, such an approach is effective in removing the disturbing effect of growing crops, but is less suitable to mask wet soils or crop residues. Although, the two wavelengths used in hyperspectral imagery to calculate the Normalized Soil moisture index [13] are very close to the two SWIR bands of the Sentinel-2, the capability of the multi-spectral data for dry vegetation detection and soil moisture quantification still has to be proven. In this regard, Demattê et al. [14] applied different thresholds based on Landsat spectral data in order to select bare soil during the dry season in Brazil and remove soil spectra affected by straw or vegetation residues.
In this case, we aim testing the effect of the threshold for a spectral index linked to soil moisture and crop residues on the performance of SOC prediction models using the Sentinel-2 and the European LUCAS topsoil database. On the one hand, such an approach will allow deriving a qualified compromise between prediction quality and spatial coverage of the map. On the other hand, exploiting the LUCAS database reduces the need for calibration data and a new ad hoc field campaign. We tested this approach within a large area (100 by 100 km) in North-eastern Germany characterized by different soil types.

Study Area and Sentinel-2 Data
The study area measures 10,000 km 2 , which corresponds to the coverage of the Sentinel-2B image acquired over North-eastern Germany on 30 [15], which is managed by the German Research Center for Geosciences Potsdam (GFZ), both institutions have cooperated in the region of Demmin. The primary objective of TERENO-NE is long-term monitoring (>15 years) and analysis of environmental change. The DEMMIN test site is located in the intensively-used agricultural state Mecklenburg-Western Pomerania, consisting of 80.3% cropland, 19.5% grassland/pastureland, and 0.2% other agrarian land use [16]. The site is within the young morainic soil-landscape of northern Germany [17] characterized by geomorphological features, including (sandy-) loamy morainic parent material of the late Pleistocene, kettle holes, and slightly undulating relief. Stagnosols, Luvisols, and their transitional types evolved from unconsolidated glacial till and colluvial soils from eroded loamy material or deposits covering former bogs [17]. The hummocky region shows a complex variability at a short distance of eroded summits to colluvial foot slopes [18]. The truncated soil profiles are mainly classified as Luvisols and Calcaric Regosols, while the colluvial soils are Gleyic-Colluvic Regosols and Terric Histosols. Tillage erosion is the dominant process of soil redistribution due to industrial farming on the very large fields. The large Borentin field (107 ha) illustrates the complex topography and its relation with the spatial patterns in soil properties (Figure 1c,d). The soils of the undulating ground moraine are truncated, while, in the center of the field, a flat depression occurs with a few kettle holes and low moor areas that are partly covered by colluvium [16].

SOC Calibration Models, Validation, and Mapping
The LUCAS topsoil dataset is a large European soil spectral library (about 20,000 soil samples from 25 countries), which was collected in the framework of the European land use/cover area frame statistical survey in 2009 [20]. The topsoil sampling locations were selected using a latin hypercube-base stratified sampling design from the LUCAS master sample grid of 2 km by 2 km. The LUCAS dataset consists of laboratory spectra and 12 chemical and physical variables for each soil sample. The spectra were acquired with an XDS Rapid Content Analyzer (FOSS NIR Systems Inc., Laurel, MD, USA) spectrometer between 400 and 2500 nm. For the calibration dataset, we selected from the entire LUCAS dataset all samples (35) collected in croplands within the study area ( Figure  1, Table 1). The laboratory spectra of the selected LUCAS soil samples were re-sampled according to the spectral bands of the MSI/Sentinel-2B sensor ( Figure 2).  [6]) and the elevation map obtained from a digital terrain model (d). All panels are in the Demmin area (Terrestrial Environmental Observatory Northeast test site) in North-eastern Germany. Coordinates are in WGS 84/UTM zone 33N.

SOC Calibration Models, Validation, and Mapping
The LUCAS topsoil dataset is a large European soil spectral library (about 20,000 soil samples from 25 countries), which was collected in the framework of the European land use/cover area frame statistical survey in 2009 [20]. The topsoil sampling locations were selected using a latin hypercube-base stratified sampling design from the LUCAS master sample grid of 2 km by 2 km. The LUCAS dataset consists of laboratory spectra and 12 chemical and physical variables for each soil sample. The spectra were acquired with an XDS Rapid Content Analyzer (FOSS NIR Systems Inc., Laurel, MD, USA) spectrometer between 400 and 2500 nm. For the calibration dataset, we selected from the entire LUCAS dataset all samples (35) collected in croplands within the study area ( Figure 1, Table 1). The laboratory spectra of the selected LUCAS soil samples were re-sampled according to the spectral bands of the MSI/Sentinel-2B sensor ( Figure 2).   The soil albedo decreases with increasing SOC content ( Figure 2). This effect is particularly developed in the MSI/Sentinel-2 NIR and SWIR bands. Consequently, in order to have a proxy of the soil brightness, we re-sampled the spectra to MSI/Sentinel-2 bands and calculated the correlation of all combinations of two bands within the SOC content (Equation (1)). These band (B) combinations were used as SOC indices (CI) to fit an exponential regression model to estimate SOC.
These models were later applied to real MSI/Sentinel-2 spectra extracted from pixels in the satellite image. The prediction accuracy of the SOC maps was evaluated on a completely independent validation dataset. The validation dataset was retrieved from the GASI+Demmin (German Agricultural Soil Inventory) dataset, that consists of 253 soil samples ( Figure 1, Table 1). The ratio of performance to deviation (RPD), i.e., the ratio between the standard deviation and the root mean square error (RMSE), was used as a measure for the estimation accuracy. Additionally, 78 samples of the GASI+Demmin dataset were collected in the framework of the German Agricultural Soil Inventory [21]. These samples were collected from a soil profile in cropland sites at 0-10 cm depth using three 200 cm 3 soil cores. The samples were dried at 40 °C and sieved at 2 mm. The SOC content was determined for a sample aliquot that was finely ground and analyzed via a dry combustion (TRUMAC, LECO, Saint Joseph, USA). If inorganic carbon was present in the samples, a thermos gradient combustion was used and the organic carbon was determined until a combustion temperature of 550 °C.
The remaining 175 GASI+Demmin samples were taken in the north-eastern part of the study area. The soil samples were collected from the topsoil (0 cm to 10 cm) using a gouge auger in 2017. Each sample consisted of five sub-samples taken within an area of a 5-m radius. The samples were air-dried and sieved at 2 mm and the total carbon was measured by dry combustion using a Carbon Nitrogen (CN) analyzer (VarioMax, ElementarGmbh, Hanau, Germany). If present, inorganic carbon The soil albedo decreases with increasing SOC content ( Figure 2). This effect is particularly developed in the MSI/Sentinel-2 NIR and SWIR bands. Consequently, in order to have a proxy of the soil brightness, we re-sampled the spectra to MSI/Sentinel-2 bands and calculated the correlation of all combinations of two bands within the SOC content (Equation (1)). These band (B) combinations were used as SOC indices (CI) to fit an exponential regression model to estimate SOC.
These models were later applied to real MSI/Sentinel-2 spectra extracted from pixels in the satellite image. The prediction accuracy of the SOC maps was evaluated on a completely independent validation dataset. The validation dataset was retrieved from the GASI+Demmin (German Agricultural Soil Inventory) dataset, that consists of 253 soil samples ( Figure 1, Table 1). The ratio of performance to deviation (RPD), i.e., the ratio between the standard deviation and the root mean square error (RMSE), was used as a measure for the estimation accuracy. Additionally, 78 samples of the GASI+Demmin dataset were collected in the framework of the German Agricultural Soil Inventory [21]. These samples were collected from a soil profile in cropland sites at 0-10 cm depth using three 200 cm 3 soil cores. The samples were dried at 40 • C and sieved at 2 mm. The SOC content was determined for a sample aliquot that was finely ground and analyzed via a dry combustion (TRUMAC, LECO, Saint Joseph, USA). If inorganic carbon was present in the samples, a thermos gradient combustion was used and the organic carbon was determined until a combustion temperature of 550 • C.
The remaining 175 GASI+Demmin samples were taken in the north-eastern part of the study area. The soil samples were collected from the topsoil (0 cm to 10 cm) using a gouge auger in 2017. Each sample consisted of five sub-samples taken within an area of a 5-m radius. The samples were air-dried and sieved at 2 mm and the total carbon was measured by dry combustion using a Carbon Nitrogen (CN) analyzer (VarioMax, ElementarGmbh, Hanau, Germany). If present, inorganic carbon Remote Sens. 2019, 11, 2121 6 of 15 was obtained by measuring the pressure of CO 2 emitted as a result of hydrochloric acid addition. The SOC content was determined by subtracting the inorganic carbon to total carbon.
The SOC distribution of the GASI+Demmin and calibration dataset is right skewed and most samples have an SOC content below 30 g kg −1 , while only a few samples show a very high SOC content ( Figure 3).
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 15 was obtained by measuring the pressure of CO2 emitted as a result of hydrochloric acid addition. The SOC content was determined by subtracting the inorganic carbon to total carbon. The SOC distribution of the GASI+Demmin and calibration dataset is right skewed and most samples have an SOC content below 30 g kg -1 , while only a few samples show a very high SOC content ( Figure 3). The SOC prediction models based on Sentinel-2 indices were applied to the exposed bare soil pixels of the Sentinel-2 image, i.e., the pixels that meet all the following criteria ( Figure 4).

•
Zero clouds probability. For this purpose, the clouds probability layer provided by European Space Agency (ESA) was used to exclude pixels with a cloud probability higher than zero.
• NDVI (normalized difference between B8 and B4) lower than 0.35 to exclude green vegetation, growing crops, and mixed pixels.
• Differences between B3 and B2 (Green Vegetation Index 1 = GVI1) and between B4 and B3 (Green Vegetation Index 2 = GVI2) both higher than 0. The application of these filters improves the soil mask [14].
• Lastly, different Normalized Burn Ratio 2 (NBR2) index thresholds (from 0.025 to 0.35 in steps of 0.025) were tested in order to exclude spectra affected by high soil moisture content or crop residues. This index was computed as the normalized difference between B11 and B12 (Equation (2)). These two bands are strongly correlated with soil moisture [22] and their differences allow discriminating between dry soil spectra, moist soil spectra, and spectra influenced by straw, other residues, or dry vegetation [14]. The difference between the reflectance around 1600 nm and that around 2100 nm is close to 0 for the soil, while the presence of straw increases this difference ( Figure  2).
The choice of the NBR2 index threshold affects the number of pixels for which SOC can be mapped and, consequently, the number of GASI+Demmin samples that can be used for the validation of the SOC models. The SOC prediction models based on Sentinel-2 indices were applied to the exposed bare soil pixels of the Sentinel-2 image, i.e., the pixels that meet all the following criteria ( Figure 4).

•
Zero clouds probability. For this purpose, the clouds probability layer provided by European Space Agency (ESA) was used to exclude pixels with a cloud probability higher than zero. • NDVI (normalized difference between B8 and B4) lower than 0.35 to exclude green vegetation, growing crops, and mixed pixels.

•
Differences between B3 and B2 (Green Vegetation Index 1 = GVI1) and between B4 and B3 (Green Vegetation Index 2 = GVI2) both higher than 0. The application of these filters improves the soil mask [14].

•
Lastly, different Normalized Burn Ratio 2 (NBR2) index thresholds (from 0.025 to 0.35 in steps of 0.025) were tested in order to exclude spectra affected by high soil moisture content or crop residues. This index was computed as the normalized difference between B11 and B12 (Equation (2)). These two bands are strongly correlated with soil moisture [22] and their differences allow discriminating between dry soil spectra, moist soil spectra, and spectra influenced by straw, other residues, or dry vegetation [14]. The difference between the reflectance around 1600 nm and that around 2100 nm is close to 0 for the soil, while the presence of straw increases this difference ( Figure 2).
The choice of the NBR2 index threshold affects the number of pixels for which SOC can be mapped and, consequently, the number of GASI+Demmin samples that can be used for the validation of the SOC models.

SOC Indices
We explored the linear correlation (Pearson correlation coefficient) between the reflectance indices using the laboratory spectra of the LUCAS dataset and the SOC content. The results showed a strong negative correlation between all indices and SOC. The best values were observed for the indices involving B11 and B12 (Figure 5a), while the combination of bands using B2 and B3 provided the weakest correlations (−0.71). Although the correlation results of Figure 5a refer to a linear correlation, they are consistent with the cross-validation RPD values of the SOC exponential regression models (Figure 5b and 6). The highest RPD values were obtained using B12 in combination with other bands (RPD > 5, Figure 5b).

SOC Indices
We explored the linear correlation (Pearson correlation coefficient) between the reflectance indices using the laboratory spectra of the LUCAS dataset and the SOC content. The results showed a strong negative correlation between all indices and SOC. The best values were observed for the indices involving B11 and B12 (Figure 5a), while the combination of bands using B2 and B3 provided the weakest correlations (−0.71). Although the correlation results of Figure 5a refer to a linear correlation, they are consistent with the cross-validation RPD values of the SOC exponential regression models (Figures 5b and 6). The highest RPD values were obtained using B12 in combination with other bands (RPD > 5, Figure 5b).

SOC Indices
We explored the linear correlation (Pearson correlation coefficient) between the reflectance indices using the laboratory spectra of the LUCAS dataset and the SOC content. The results showed a strong negative correlation between all indices and SOC. The best values were observed for the indices involving B11 and B12 (Figure 5a), while the combination of bands using B2 and B3 provided the weakest correlations (−0.71). Although the correlation results of Figure 5a refer to a linear correlation, they are consistent with the cross-validation RPD values of the SOC exponential regression models (Figure 5b and 6). The highest RPD values were obtained using B12 in combination with other bands (RPD > 5, Figure 5b).

SOC Maps and Validation
Furthermore, 12.5% of the pixels within the Sentinel-2 image were characterized by cloud cover, and only 26.9% of the non-cloud pixels were classified as 'without green vegetation' using the NDVI, GVI1, and GVI2 thresholds (Figure 7). Thus, 34 out of the 253 GASI+Demmin samples were localized within crop pixels or cloudy areas (Tables 1 and 2). An NBR2 threshold from 0.35 to 0.25 did not affect the number of selected pixels (26.8%) and the validation samples remaining in a suitable condition to apply the SOC models (219). The selected pixels decreased markedly when an NBR2 threshold of 0.1 was reached (58 remaining validation samples). For the most severe NBR2 threshold (0.05), only 10% of the area was selected (Figure 7) and only 17 validation samples remained (Table 2). Clearly, the NBR2 threshold affects the number of the pure pixels selected and, consequently, the surface that can be mapped. Without applying NBR2, we can map 26.9% of the investigated area, while using NBR2, we reduce the mapped area to 10% of the entire Sentinel-2 image. The mapped area in one restricted framework of the Sentinel-2 image decreases when the number of pixels to which an SOC estimation model is applied is restricted using an NBR2 threshold of 0.05 compared to a threshold of 0.075 ( Figure 8). On the one hand, the decrease of the NBR2 threshold reduces the area for which the SOC is mapped (Figure 7) and, on the other hand, a decrease in NBR2 improves the SOC validation accuracy ( Table 2). In Table 2 and Figure 9 (black solid line), the best RPD value among all the band combinations for different NBR2 thresholds is displayed. Moreover, we showed the results for two indices: the average value between B5 and B6 (Red-Edge Carbon Index: RE-CI) and the average between B4 and B5 (Red-Red-Edge Carbon Index: RRE-CI). The SOC validation accuracy, in terms of RPD, decreases with the NBR2 threshold ( Figure 9). Thus, the highest RPD (4.4) was detected for an NBR2 value of 0.05 using the RE-CI ( Table 2). The RRE-CI calibration model showed RPD values larger than 1.5 for NBR2 values below 0.2. Overall, RE-CI performed less well and already reached an RPD of 1.5 for the NBR2 value larger than 0.075. Increasing the NBR2 threshold from 0.05 to 0.1 yielded an underestimation of the additional samples involved in the validation (Figure 10).

SOC Maps and Validation
Furthermore, 12.5% of the pixels within the Sentinel-2 image were characterized by cloud cover, and only 26.9% of the non-cloud pixels were classified as 'without green vegetation' using the NDVI, GVI1, and GVI2 thresholds (Figure 7). Thus, 34 out of the 253 GASI+Demmin samples were localized within crop pixels or cloudy areas (Tables 1 and 2). An NBR2 threshold from 0.35 to 0.25 did not affect the number of selected pixels (26.8%) and the validation samples remaining in a suitable condition to apply the SOC models (219). The selected pixels decreased markedly when an NBR2 threshold of 0.1 was reached (58 remaining validation samples). For the most severe NBR2 threshold (0.05), only 10% of the area was selected (Figure 7) and only 17 validation samples remained (Table 2). Clearly, the NBR2 threshold affects the number of the pure pixels selected and, consequently, the surface that can be mapped. Without applying NBR2, we can map 26.9% of the investigated area, while using NBR2, we reduce the mapped area to 10% of the entire Sentinel-2 image. The mapped area in one restricted framework of the Sentinel-2 image decreases when the number of pixels to which an SOC estimation model is applied is restricted using an NBR2 threshold of 0.05 compared to a threshold of 0.075 ( Figure 8). On the one hand, the decrease of the NBR2 threshold reduces the area for which the SOC is mapped (Figure 7) and, on the other hand, a decrease in NBR2 improves the SOC validation accuracy ( Table 2). In Table 2 and Figure 9 (black solid line), the best RPD value among all the band combinations for different NBR2 thresholds is displayed. Moreover, we showed the results for two indices: the average value between B5 and B6 (Red-Edge Carbon Index: RE-CI) and the average between B4 and B5 (Red-Red-Edge Carbon Index: RRE-CI). The SOC validation accuracy, in terms of RPD, decreases with the NBR2 threshold ( Figure 9). Thus, the highest RPD (4.4) was detected for an NBR2 value of 0.05 using the RE-CI ( Table 2). The RRE-CI calibration model showed RPD values larger than 1.5 for NBR2 values below 0.2. Overall, RE-CI performed less well and already reached an RPD of 1.5 for the NBR2 value larger than 0.075. Increasing the NBR2 threshold from 0.05 to 0.1 yielded an underestimation of the additional samples involved in the validation ( Figure 10).

Discussion
Several differences exist between a laboratory spectrum of a dried and sieved soil sample and the spectrum of the satellite sensors representing the pixel where the sample was taken. One of the most evident differences is related to the distance between the sensor and the target. The atmospheric disturbance and lower signal-to-noise ratio of the remote sensor affect the spectral output, and, consequently, the estimation accuracy of the soil property. The B12 (2185 nm) and B11 (1610 nm) of the MSI/Sentinel-2 sensor, in combination with other bands, showed the highest correlation with SOC using the LUCAS spectra acquired in the laboratory ( Figure 5). We observed a weakening of the correlations removing the two calibration samples with the highest SOC values (Pearson's coefficient between 0.5 and 0.62). This decrease in correlation especially concerns the combination of bands including B11 and B12. These two bands proved to be very important for calibration sets with a large SOC variability [6]. In other words, the differences in terms of soil brightness become more evident.

Discussion
Several differences exist between a laboratory spectrum of a dried and sieved soil sample and the spectrum of the satellite sensors representing the pixel where the sample was taken. One of the most evident differences is related to the distance between the sensor and the target. The atmospheric disturbance and lower signal-to-noise ratio of the remote sensor affect the spectral output, and, consequently, the estimation accuracy of the soil property. The B12 (2185 nm) and B11 (1610 nm) of the MSI/Sentinel-2 sensor, in combination with other bands, showed the highest correlation with SOC using the LUCAS spectra acquired in the laboratory ( Figure 5). We observed a weakening of the correlations removing the two calibration samples with the highest SOC values (Pearson's coefficient between 0.5 and 0.62). This decrease in correlation especially concerns the combination of bands including B11 and B12. These two bands proved to be very important for calibration sets with a large SOC variability [6]. In other words, the differences in terms of soil brightness become more evident.

Discussion
Several differences exist between a laboratory spectrum of a dried and sieved soil sample and the spectrum of the satellite sensors representing the pixel where the sample was taken. One of the most evident differences is related to the distance between the sensor and the target. The atmospheric disturbance and lower signal-to-noise ratio of the remote sensor affect the spectral output, and, consequently, the estimation accuracy of the soil property. The B12 (2185 nm) and B11 (1610 nm) of the MSI/Sentinel-2 sensor, in combination with other bands, showed the highest correlation with SOC using the LUCAS spectra acquired in the laboratory ( Figure 5). We observed a weakening of the correlations removing the two calibration samples with the highest SOC values (Pearson's coefficient between 0.5 and 0.62). This decrease in correlation especially concerns the combination of bands including B11 and B12. These two bands proved to be very important for calibration sets with a large SOC variability [6]. In other words, the differences in terms of soil brightness become more evident.
However, when we applied these indices involving the SWIR bands (B11 and B12) to real Sentinel-2 data, the accuracy was always unsatisfactory and it decreased with an increasing NBR2 value. The loss of prediction capability of the SWIR bands from the laboratory to the Sentinel-2 data is mainly due to three factors: i) the large bandwidth of the two bands (141 nm for B11 and 238 nm for B12) and their lower spatial resolution (20 m) as compared to the other Sentinel-2 bands, ii) the lower signal-to-noise ratio (SNR) as compared to that of the visible and infrared bands [6], iii) the mixed influence of several disturbing factors on the spectral signal of the SWIR bands due to variable soil texture and mineralogy (clay or carbonate content), soil moisture [23], or dry vegetation cover [24].
Generally, the best validation results were obtained using a combination of red edge bands (RE-CI) or red/red edge (RRE-CI, Table 2, Figure 11). The exponential calibration models of Figure 6a,b showed how a small variation of the index values between 0.10 and 0.12 corresponds to abrupt changes in SOC content, while the two indices are less sensitive for SOC values lower than 100 g kg −1 . Many authors reported a wide absorption feature related to SOC around 664 nm [4,25], which includes the spectral range of B4 and B5. The MSI/Sentinel-2 band centered at 740 nm (B6) is spectrally close to the absorption feature associated with the overtone of the N-H bond [25,26]. This is a chemical bond that characterizes some compounds of soil organic matter, e.g., amino acids [27]. Castaldi et al. [6] analyzed the relative variable importance for estimating SOC in the Demmin area using the MSI/Sentinel-2 bands as variables in a random forest regression model. Their results showed the most important bands to be the B4 and B5 as well as the two SWIR bands (B11 and B12).
The structural complexity of the organic matter entails a large heterogeneity of the components, which, in turn, involves a large variability of spectral responses and, consequently, the need to calibrate 'local' SOC estimation models based on geographical or spectral proximity. After all, these empirical models can hardly be applied in areas characterized by different soil types. The SOC indices mainly exploit the link between SOC content and soil brightness. This link can be masked or weakened as a result of variation in mineralogy that affects the soil color as well as the type of clay minerals and iron oxide contents. Consequently, the transferability of these spectral indices is likely to be more straightforward for young soil types, where variability in clay minerals and iron oxide content is not yet pronounced.
The availability of a sufficient number of LUCAS soil samples within the study area and their wide spatial distribution ( Figure 1) together with the large variability in terms of SOC content (Table 1) allowed calibrating reliable SOC models [9]. The only drawback in the LUCAS database is that it is based on a stratified random sampling strategy in geographical space (rather than feature space) [6]. Since the area covered by peaty topsoil in cropland is limited to the center of the kettle holes [28], the number of high SOC content samples in the LUCAS database is very limited (Table 1). It is likely that the LUCAS samples in other European regions are not numerous enough [29] or not representative of the soil variability of the investigated area. Until now, remote sensing or airborne spectroscopy was used for small areas that are mostly rather homogeneous in parent material and soil forming factors [30]. Castaldi et al. [9] demonstrated that sampling strategies based on the feature space, where the spectral bands were used as ancillary data, were most efficient when the Demmin area was stratified, according to 'soil scapes,' i.e., distinguishing sandy and clay topsoil (Figure 1). Other factors that are likely to improve the SOC prediction models are the absence of a soil crust (since it does not represent the spectral signal of the topsoil [30] and the large variability in SOC content). Gomez et al. [31] reported that, in contrast to clay and calcium carbonate content, SOC could not be predicted in Mediterranean environments due to the very low and rather constant SOC content. It should be noted that only 10% to 15% of the landscape is not affected by erosion [28]. The fact that the landscape is derived from relatively young glacial sediments explains that differences in pedogenic oxides, leaching of calcium carbonate, and differentiation in clay minerals have not yet had enough time to develop and are, therefore, unlikely to greatly influence the soil spectra. Consequently, the SOC indices can effectively exploit the link between SOC content and soil brightness in young soil types where this link is not masked or weakened from variation in mineralogy due to the pedogenesis age. Straw, other crop residues, or dry vegetation are disturbing factors for soil imaging spectroscopy. Currently, several methods exist to map dry vegetation fractional cover based on hyperspectral imagery. Such methods can be used to separate pure soils from soils partly covered by vegetation [24]. We selected green vegetation pixels using NDVI and two other indices based on the visible region and applied the NBR2 index to detect dry vegetation or crop residues. Demattê et al. [14] already showed that the NBR2 values are correlated with the amount of dry vegetation at the soil surface. In their study, the index provided negative, or close to zero, values for pixels with a 100% soil signal, while the NBR2 values increased with straw cover. Therefore, it is reasonable to adopt a very low NBR2 threshold in order to select pure soil pixels for remote sensing applications ( Figure  9). We tested different NBR2 thresholds, and sought to understand how they affect the SOC estimation accuracy. The best validation accuracy was obtained using an NBR2 threshold of 0.05. While using a less restrictive threshold, the estimation accuracy sharply decreased ( Figure 9). The SOC prediction of the RE-CI and RRE-CI indices decreases with an increasing NBR2 threshold ( Figure 10). This is mainly due to an underestimation of the low SOC content. Although the statistical accuracy in terms of R 2 , RPD, and RMSE is still satisfactory, the SOC underestimation already appears for an NBR2 of 0.075 and becomes striking for a NBR2 of 0.1 (Figures 10b and 10c). Demattê et al. [14] remarked that the NBR2 index is not particularly effective for pixels with a residue cover between 80% and 20% and that a stricter NBR2 threshold (<0.05) does not improve detection of pure soil pixels. On the contrary, a less severe threshold (>0.05) will result in the selection of pixels that are almost completely covered by dry vegetation or residues.
The spectrum of a bare soil pixel can also be affected by soil moisture, which entails a nonlinear reduction of the reflectance over the entire spectrum [32]. This challenge in predicting soil moisture effects made many researchers refrain from the use of laboratory spectral libraries for the calibration of remote sensing data or incited them to find alternative solutions to the problem by spectral transfer methods [33] or by a bottom-up approach [11]. We faced the issue from another perspective, detecting and selecting only the bare soil pixels mimicking the conditions of a dry soil sample. Several methods exist that could be applied to mask moist soil pixels based on hyperspectral imagery [34,35]. The adaptation of these methods to multispectral satellites has yet to be investigated. Some authors investigated the suitability of normalized soil moisture indices to be applied for remote sensing [13,34]. Both the normalized indices calibrated by Haubrock et al. [13] and Castaldi et al. [34] exploited the combination of two SWIR wavelengths: 1800 nm and 2119 nm in Haubrock et al. [13] (NSMI = normalized soil moisture index) and between 1700 nm and 2100 nm in Castaldi et al. [34] (SMIR_A = soil moisture index for remote sensing). The SWIR bands of the MSI/Sentinel-2 sensor Straw, other crop residues, or dry vegetation are disturbing factors for soil imaging spectroscopy. Currently, several methods exist to map dry vegetation fractional cover based on hyperspectral imagery. Such methods can be used to separate pure soils from soils partly covered by vegetation [24]. We selected green vegetation pixels using NDVI and two other indices based on the visible region and applied the NBR2 index to detect dry vegetation or crop residues. Demattê et al. [14] already showed that the NBR2 values are correlated with the amount of dry vegetation at the soil surface. In their study, the index provided negative, or close to zero, values for pixels with a 100% soil signal, while the NBR2 values increased with straw cover. Therefore, it is reasonable to adopt a very low NBR2 threshold in order to select pure soil pixels for remote sensing applications (Figure 9). We tested different NBR2 thresholds, and sought to understand how they affect the SOC estimation accuracy. The best validation accuracy was obtained using an NBR2 threshold of 0.05. While using a less restrictive threshold, the estimation accuracy sharply decreased (Figure 9). The SOC prediction of the RE-CI and RRE-CI indices decreases with an increasing NBR2 threshold ( Figure 10). This is mainly due to an underestimation of the low SOC content. Although the statistical accuracy in terms of R 2 , RPD, and RMSE is still satisfactory, the SOC underestimation already appears for an NBR2 of 0.075 and becomes striking for a NBR2 of 0.1 (Figure 10b,c). Demattê et al. [14] remarked that the NBR2 index is not particularly effective for pixels with a residue cover between 80% and 20% and that a stricter NBR2 threshold (<0.05) does not improve detection of pure soil pixels. On the contrary, a less severe threshold (>0.05) will result in the selection of pixels that are almost completely covered by dry vegetation or residues.
The spectrum of a bare soil pixel can also be affected by soil moisture, which entails a nonlinear reduction of the reflectance over the entire spectrum [32]. This challenge in predicting soil moisture effects made many researchers refrain from the use of laboratory spectral libraries for the calibration of remote sensing data or incited them to find alternative solutions to the problem by spectral transfer methods [33] or by a bottom-up approach [11]. We faced the issue from another perspective, detecting and selecting only the bare soil pixels mimicking the conditions of a dry soil sample. Several methods exist that could be applied to mask moist soil pixels based on hyperspectral imagery [34,35]. The adaptation of these methods to multispectral satellites has yet to be investigated. Some authors investigated the suitability of normalized soil moisture indices to be applied for remote sensing [13,34]. Both the normalized indices calibrated by Haubrock et al. [13] and Castaldi et al. [34] exploited the combination of two SWIR wavelengths: 1800 nm and 2119 nm in Haubrock et al. [13] (NSMI = normalized soil moisture index) and between 1700 nm and 2100 nm in Castaldi et al. [34] (SMIR_A = soil moisture index for remote sensing). The SWIR bands of the MSI/Sentinel-2 sensor used for the NBR2 are very close to those employed for the NSMI and SWIR_A index. B11 is centered at 1610 nm and it has the upper boundary at 1682 nm, while B12 is centered at 2185 nm and has the lower boundary at 2078 nm. Thus, the NBR2 can be used as proxy for soil moisture characterization.
Soil roughness is another disturbing factor for SOC estimation by satellite data. Generally, its influence is limited when the satellite images are acquired for smooth soils in a seedbed condition. However, this particular condition occurs only twice yearly and lasts several days to a week in temperate regions. Therefore, the time window for acquiring Sentinel-2 images is rather restricted. Usually, the periods of the year when it is more probable to find soil in a seedbed condition is spring when summer crops such as sugar beet, potatoes, or maize are sown or the end of the summer when winter cereals are sown.
Taking into account, the above-mentioned restrictions, on the one the hand, entails an improvement of the SOC estimation accuracy and, on the other hand, decreases the extent of the mapped surface (Figures 7 and 8). Nevertheless, the pixel selection process can be tuned according to the goal of the survey, which gives a warning on the reliability of the SOC map for each NBR2 threshold.

Conclusions
We tested the effect of the threshold for a spectral index linked to soil moisture and crop residues on the performance of SOC prediction models in croplands using the European LUCAS database and Sentinel-2 data. The resampling procedure of the LUCAS spectra according the MSI/Sentintel-2 bands allowed pointing out the strong correlation between the reflectance in the SWIR region and SOC content. However, the SOC estimation models using SWIR bands showed a consistent decrease of their performances when they were applied to real Sentinel-2 spectra due to the poorer spatial and spectral resolution and lower SNR as compared to the other bands, combined with the confounding effect of crop residues and soil moisture on the SWIR bands. The red edge MSI/Sentinel-2 bands and a combination of red and red edge bands provided the best estimation accuracies for the validation dataset. These results showed that: (i) the application of an SOC prediction model using visible and red edge spectra acquired in the laboratory to real Sentinel-2 data, and (ii) the high signal quality of the three Sentinel-2 bands involved in the RE-CI and RRE-CI indices (B4, B5, and B6). The pure pixel selection based on an NRB2 threshold has proven to be necessary to obtain reliable SOC maps. In particular, an NBR2 threshold of 0.05 allowed selecting (nearly) bare soil pixels from the Sentinel-2 image. It was observed that higher NBR2 threshold values allowed mapping SOC over larger areas at the expense of a decreasing estimation accuracy, especially for NBR2 values above 0.175.