Some Peculiarities of Arable Soil Organic Matter Detection Using Optical Remote Sensing Data

Abstract

Optical remote sensing only provides information about the very thin surface layer of soil. Rainfall splash alters soil surface properties and its spectral reflectance. We analyzed the impact of rainfall on the success of soil organic matter (SOM) content (% by mass) detection and mapping based on optical remote sensing data. The subject of the study was the arable soils of a test field located in the Tula region (Russia), their spectral reflectance, and Sentinel-2 data. Our research demonstrated that rainfall negatively affects the accuracy of SOM predictions based on Sentinel-2 data. Depending on the average precipitation per day, the R²cv of models varied from 0.67 to 0.72, RMSEcv from 0.64 to 1.1% and RPIQ from 1.4 to 2.3. The incorporation of information on the soil surface state in the model resulted in an increase in accuracy of SOM content detection based on Sentinel-2 data: the R²cv of the models increased up to 0.78 to 0.84, the RMSEcv decreased to 0.61 to 0.71%, and the RPIQ increased to 2.1 to 2.4. Further studies are necessary to identify how the SOM content and composition of the soil surface change under the influence of rainfall for other soils, and to determine the relationships between rainfall-induced SOM changes and soil surface spectral reflectance.

1. Introduction

Organic matter is an important indicator of soil state as it forms and maintains the main regimes, properties, and functions of soil [1]. Evaluation of the soil state and quality is based mainly on the content and quality of soil organic matter (SOM) and its fractions: humus, total organic carbon, water soluble carbon, microbial biomass carbon, etc. [2,3,4,5,6]. Information on organic matter content is included in soil databases and monitoring systems from regional to global scales [7,8,9,10]. According to a previous review [11], SOM is the most frequently used indicator of soil quality. A decrease in SOM content, caused by human-induced activities, is acknowledged to be one of the five global soil threats as it results in the disruption of soil functions, including soil fertility, and regulation of CO₂ in the atmosphere as a main driver of global warming [12,13].

Optical remote sensing data have been widely used to map the organic matter content in agricultural soils. Numerous studies conducted in various regions and with different types of remote sensing data have demonstrated that SOM content can be detected with different degrees of accuracy [14,15,16,17,18,19,20,21,22,23]. Due to short wavelengths, the optical range only provides information about the very thin soil surface layer [24]. Therefore, the measured spectral reflectance of the surface of arable soils in the optical range only contains information about the properties of this surface layer. In a previous study, Liang [24] determined that the sensible optical depth is about 3, which corresponds to a geometric depth of four to five times the particle effective radius. For example, for clay particles, the geometric depth was 2.2 μm at 0.5 μm and 3.15 μm at 1.105 μm.

The moisture content and surface roughness are considered the main soil-related limiting factors when explaining the low performance of remote sensing-based SOM detection models [25,26,27,28,29]. However, one more important interfering factor that is rarely considered or discussed is the change in bare soil surface properties such as the soil mineralogical composition, soil texture, and SOM content and chemical composition under the impact of rainfall. Rainfall splash causes the breakdown of soil surface aggregates and the redistribution of formed soil material which is partially transported as a lateral surface flow and partially washed in soil surface pores, causing the sealing of the soil surface and the formation of the soil crust [30,31,32,33,34]. Studies on the impact of rainfall on the soil surface show that it results in changes in the soil surface organic matter content, mineral composition, and texture as well as soil surface reflectance [35,36,37,38,39,40,41]. The longer the surface of non-irrigated arable soils is left untreated under the influence of atmospheric precipitation, the more its spectral properties and material composition change and the greater the difference between its properties and the properties of the underlying part of the arable horizon becomes [40,42]. Most studies using optical satellite data or field spectral reflectance data to map SOM content are based only on a single date, and they do not register the soil surface state at the time of spectral data acquisition. However, the extent to which the soil surface transformed by rainfall provides information on the properties of the arable horizon is not yet clear. This discrepancy is one of the main reasons why models obtained in the laboratory fail in usage when applied to field measurements. It also limits the reproducibility of models predicting the SOM content based on satellite data.

The aim of the research was to study the impact of rainfall-induced soil surface spectral changes, in particular, the formation of soil crust, on the accuracy of SOM content mapping based on optical remote sensing data.

2. Materials and Methods

2.1. Study Area and Field Data Collection

The study area is a test field located in the Yasnogorsky district of the Tula region (Russia) with grey forest arable soil (Albic Luvisols [43]) (Figure 1). The area of the test field is 13.5 ha. The test plot is located on a flat part of the slope, 2 to 3 degrees of the western aspect. In 2019, the field was complete fallow, which is a part of crop rotation. Such a field is free of crops during the growing season. It is treated throughout the warm period of the year, maintaining the soil in a loose state, carrying out extermination measures against weeds, pests and pathogens of crops. Periodical tilling as a part of complete fallow cultivation agrotechnology is performed to exterminate emerging weeds and to break down the soil crust to accumulate soil moisture for the next crop and reduce the risk of soil degradation development. While this practice has some disadvantages, it is still widely used in Russia, especially in dry regions, to prepare soil for the sowing of winter crops [44].

During field work conducted on the test field on 15 August 2019, we collected 30 mixed soil samples from the upper part of the arable horizon (0 to 5 cm). Each mixed sample consisted of five single sub-samples of the upper part of the arable horizon collected from the area around the sample point (approximately 1 m). Sub-samples were collected using the diagonal method at five points. We used the simple random sampling pattern. The sampling scheme is presented in Figure 1.

SOM content was determined in the collected soil samples as a percentage of mass by the Tyurin method [45], which is the analogue of the Walkley–Black method [46]. The spectral reflectance of the dry surface and wet subsurface (2 to 5 cm depth) layer was measured in the field at the sample points with a HandHeld-2 field spectroradiometer (proximal data) operating in the range of 325 to 1075 nm with steps of 1 nm. To measure the spectral reflectance of the soil subsurface layer, we removed 2 cm of the soil surface after measuring the soil surface spectral reflectance. At each sample point, the spectral reflectance was measured with fivefold repetition. Reflectance measurements were made using radiance reflected from a sample and radiance incidence reflected from a 100% reflectance reference target (white reference panel) [47]. The white reference panel was positioned horizontally near the sample and was oriented to receive all of the available incident illumination while at the same time avoiding shadowing and light reflected from surrounding objects. Since the reflected radiance from the panel is very close to the total energy falling on the sample, this reading can be used in the denominator of the reflectance measurement. Sample reflectance was calculated by HandHeld-2 instrument software. The location of each sample point was registered by GPS attached to the spectroradiometer (GlobalSat BU-353 GLONASS: https://www.globalsat.ru/catalog/bu-353-glonass (accessed on 25 May 2021)).

In this research, we also used additional spectral data collected from this test field in our previous study: spectral reflectance of the soil crust (Figure 2a) and spectral reflectance of the sample of the ploughed soil layer. Spectral reflectance of these two additional samples was also measured with a HandHeld-2 field spectroradiometer, the same way as spectral reflectance of dry surface and wet subsurface layer as described above. The soil sample of the ploughed soil layer represents a soil surface that forms in the field after the ploughing and has not experienced the impact of rainfall yet.

2.2. Soil Crusting Process

There exist three main types of soil crust: physical, biological and chemical [48]. In our research, we are dealing with the physical soil crust, which is a very thin layer formed as a result of the breakdown of soil surface aggregates under the impact of atmospheric precipitation and the redistribution of the formed soil material. According to our previous research on the formation of soil crust and data acquired in micromorphological and microtomographic analysis, the thickness of the well-developed soil crust formed on soils of the European part of Russia is about 1 to 2 mm, and it varies spatially (Figure 2b,c) [41]. We think that it would be more correct to call such a thin layer a microlayer. The degree of soil crust development depends on how long the soil surface is bare, undisturbed by agricultural practices, and on the intensity of rainfall events as well as on a number of soil properties (soil texture, salinity, SOM content, percentage of calcium carbonate, percentage of exchangeable sodium) [48,49]. Moreover, the soil crust is comprised of heterogenous soil material (residues of soil aggregates and washed soil material, as well as, in some cases, water-resistant soil aggregates) (Figure 2a).

The studied soil did not contain calcium carbonate or water-soluble salts, and the percentage of exchangeable sodium is very low (Table 1). Therefore, the main soil-related factors affecting soil crust formation for the studied soil are SOM content and soil texture. As SOM acts as an aggregate stabilizing agent, its content affects the speed and degree of soil crust development. Soil texture predetermines the mechanism of soil crust formation.

The analyzed soil was formed on loess. The composite sample collected from the ploughed horizon (0 to 30 cm) of the studied field using the diagonal method contained 6.5% of sand (0.05 to 1 mm), 75.6% of silt (0.001 to 0.05 mm) and 17.9% of clay (<0.001 mm) [50]. This is silt loam in the USDA classification. The conversion from the Russian classification to the USDA classification was made based on [51].

The impact of atmospheric precipitation on soils with such a texture causes the breakdown of soil surface aggregates, and the formed material is washed and moved into the pores of the soil surface. In Figure 2b, this washed material has a light brown color, while the dark brown color corresponds to the remains of the soil aggregates. In the field, when we look at the surface of the soil crust, the washed soil material has a very light color, close to white (Figure 2a). The formation of the soil crust in the field goes through certain stages. They have been defined micromorphologically through the analysis of changes in the microstructure [36].

In 2017, we conducted a model experiment, wherein we exposed the samples of ploughed horizons of three soil types, including the studied soil from the test field, to the impact of atmospheric precipitation [40,41]. We put soil samples in rectangular plastic boxes and left them outside, periodically taking photos of the soil surface (Figure 3).

According to our results, which are in accordance with [33], the longer bare soil surface was subjected to rainfall, the less soil aggregates were left there and the more washed soil material formed as a result of the breakdown of soil surface aggregates (Figure 3 and Figure 4). Therefore, the amount of washed soil material on the soil surface can be used as an indicator of soil crust development. As the soil surface on 11 July 2017 was wet (Figure 3f), the area of washed soil material was underestimated due to the interference of soil moisture content (Figure 4). The cracks did not appear after the first rainfall. We registered the first cracks after the cumulative rainfall reached 100 mm (Figure 4).

Rainfall-induced changes in the soil surface microlayer affect the soil spectral reflectance. For example, the soil crust of the studied soil (Figure 2a), collected from the watershed of the test field in our previous studies, had higher reflectance intensity compared to the sample of the ploughed soil layer (Figure 5). The difference in spectral reflectance increased with the increase in wavelength. Moreover, the average spectral reflectance of dry surface measured in the field on 15 August 2019 is very close to the spectral reflectance of the sample of the ploughed soil layer (Figure 5). It shows that tillage occurred close to the date of our field trip, meaning that the upper layer was rather uniform. Thus, the main difference between the soil surface microlayer and subsurface layer on 15 August 2019 was due to the moisture content.

Soil crust is the initial stage of soil degradation. This microlayer is comprised of the remains of soil aggregates, so its original structure is lost. As it is compacted, the water penetration decreases while the surface runoff increases. The latter adds up to soil erosion development. The loss of aggregate structure, compaction, and increased surface runoff contribute to the proliferation of soil degradation.

2.3. General Strategy of the Research

First of all, we started with the development of the model for SOM content prediction from the data collected in the field on 15 August 2019 (Figure 6). For that, we used spectral reflectance of the dry soil surface microlayer. As a result, we determined informative spectral parameters that can be used for the modelling of SOM in the upper layer of the test field from optical remote sensing data. Detailed description is given in the Section 2.6 and Section 2.7.

As soil moisture can negatively affect the accuracy of SOM content prediction as well as result in the underestimation of the soil crust area, we used spectral reflectance of the dry soil surface and wet subsurface layer measured in the field on 15 August 2019 to select Sentinel-2 images with the dry soil surface. The selection was based on the intensity of spectral reflectance and its closeness to the spectral reflectance of the dry surface and wet subsurface layer. Detailed description is given in the Section 2.8.

At the next step, we modeled SOM content from Sentinel-2 images for the dates with the dry soil surface and estimated the accuracy of the developed models. For all the dates, we used the same predictors—the informative spectral parameter (parameters) determined during the analysis of spectral reflectance of the dry soil surface measured in the field on 15 August 2019. We assumed that difference in the accuracy of models at different dates from the model based on spectral reflectance of the dry soil surface was caused by the difference in the state of the soil surface microlayer on these dates and on the date of the field trip in 2019. We looked at the average spectral reflectance for each analyzed date, spectral reflectance of spectral mixtures with different percentages of soil crust and the rainfall data to support our assumption. A detailed description is given in the Section 2.9.

Finally, we accounted for the changes in the soil surface microlayer in models for SOM content estimation by adding as predictors the area of the soil crust and error of spectral unmixing, indicating the area of shadows (including cracks). These parameters were determined for each date with a dry soil surface as a result of spectral unmixing. A detailed description is given in the Section 2.10.

As at the time of field trip in 2019, there was not a soil crust on the soil surface and we did not catch the soil surface immediately after the ploughing; we used additional spectral data, collected during our previous studies and described in Section 2.1, representing these two soil surface states (spectral reflectance of soil crust and of sample of ploughed soil layer) to model spectral mixtures with different percentages of soil crust and to perform spectral unmixing.

2.4. Sentinel-2 Data Collection and Preprocessing

Seven Sentinel-2 scenes were chosen for the studied region: 2 April 2019, 17 April 2019, 20 April 2019, 5 May 2019, 6 June 2019, 19 June 2019, and 28 August 2019 (https://earthexplorer.usgs.gov/ (accessed on 25 May 2021)) (Table 2). Atmospheric correction of the chosen scenes was performed with Sen2Cor in SNAP. Afterwards, we imported all of the scenes in ILWIS 3.3 Academic [52] and extracted the reflectance values for pixels where we had collected spectral data and soil samples in the field.

2.5. Meteorological Data Collection and Analysis

Meteorological data for 2019 for the test field were downloaded from the web-service Vega (http://pro-vega.ru/ (accessed on 25 May 2021)) (reanalysis data). The data were used to study rainfall patterns in between dates of selected Sentinel-2 data acquisition. We also calculated the average precipitation per day (APD) for the periods between image acquisition dates as the ratio of atmospheric precipitation accumulated during the period to the length of the period in days. The APD was further used as a proxy of the degree of impact of rainfall on the bare soil surface.

2.6. The Preprocessing and Analysis of Field Spectral Data

All spectral reflectance curves used in this study were averaged and smoothed with the Savitzky–Golay function of the prospectr package of R software [53]. Due to the high level of noise, regions shorter than 350 nm and longer than 900 nm were excluded from further analysis. After that, smoothed field spectral reflectance was convolved into Sentinel-2 bands assuming Gaussian distribution with the spectralResampling function of the hsdar package of R software [54] by analogy with [55]. Considering the spectral range of the HandHeld-2 spectroradiomener, as a result of spectral generalization of the field spectral data, we had spectral reflectance values only for nine Sentinel-2 bands (from Band 1 to Band 8a) (Table 2). All further analysis and calculations were made for spectrally generalized data.

As the soil moisture content has a masking effect on SOM spectral features [26,29], we started the analysis with spectral reflectance of the dry surface. Using these spectrally generalized spectral reflectance data, we calculated a number of spectral indices and ratios: spectral ratios of analyzed Sentinel-2 bands, BI, SI, HI, CI, RI, AVI, GSI, CRUST, BSCI, NDVI, SAVI, EVI2, and NDWI [56,57,58,59,60,61,62,63,64] (Table A1). Then, we checked the correlations between the calculated parameters to reduce the number of predictors for further regression modelling and to avoid multicollinearity. Only the indices and ratios with correlations between each other lower than 0.6 were further used in linear regression modelling. The selection of such parameters was made in R software using the function vifcor of the usdm package with a threshold of 0.6 [65].

2.7. Modelling SOM Content with Transformed Field Spectral Data

We used the parameters selected in the previous step to model SOM content with the linear regression method. First, all of the selected parameters were included in the model. Then, step by step, we excluded the parameters with a p-value above 0.05 until the model included only the parameters with a p-value below 0.05. Regression modelling was performed with the lm function of the stats package [66].

To estimate the accuracy and quality of the obtained models, we repeated 10-fold cross-validation 100 times and calculated the cross-validated coefficient of determination (R²cv), the cross-validated root-mean-squared error (RMSEcv), and the ratio of the interquartile range to the error of prediction (RPIQ). The model developed from spectral reflectance of the dry soil surface with the highest R²cv and RPIQ and the lowest RMSEcv was considered the best, and the parameter (or parameters) included in such a model was considered the most informative parameter for the prediction of SOM content. Cross-validation was performed with the package caret [67] by analogy with the example presented in (http://www.sthda.com/english/articles/38-regression-model-validation/157-cross-validation-essentials-in-r/ (accessed on 25 May 2021)). The same informative parameter (parameters) was used to model SOM content from spectral reflectance of the wet subsurface layer. This model was also cross-validated.

We also checked for spatial autocorrelation of the residuals for models from both spectral reflectance of the dry surface and wet subsurface layer to ensure that it did not affect our results. For this, we used a semivariogram by analogy with [68]. First, we calculated the residuals for the best regression model with the stats package [66]. Then, we added this information to the attribute table with coordinates of sample points in ILWIS Academic 3.3 and calculated the semi-variance with a lag spacing of 10 m. The semivariogram was created in Microsoft Excel.

2.8. Assessment of Soil Surface Status at the Time of Satellite Data Acquisition

As soil moisture has a masking effect on the spectral features of other soil properties (including SOM content and mineralogy), its interference can result in underestimation of the soil crust area. In our model experiment in 2017, described in Section 2.2, for the date with a wet soil surface (Figure 3f), we registered a significant drop in the area of washed soil material based on the analysis of soil surface photos (Figure 4). As the soil surface of the soil sample was kept undisturbed during the whole experiment, such a one-time decrease in the detected area of washed soil material was caused by the masking effect of soil moisture. Therefore, we decided to exclude Sentinel-2 images from further analysis where soil surface moisture content may interfere significantly. We suggested that the soil surface microlayer (upper 1 to 2 mm) is most often dry at the times when satellite images in an optical domain are acquired, except on the days when images were acquired after rainfall and if soils are hydromorphic (with the capillary rise from ground water). As recent research has shown, the soil surface microlayer dries out quicker than the underlying soil horizon [69]. Therefore, to correctly estimate the soil moisture content in the soil surface microlayer, it is necessary to collect samples that include the material only from that microlayer. Otherwise, the results will be unreliable. This is a considerably challenging task. In [70], the authors collected soil samples from the very first mm of the soil surface for the analysis of soil moisture content, but they pointed out that it is not always possible to collect only this thin layer without including the material of the underlying soil horizon, which can result in the overestimation of the soil moisture content. Unfortunately, they did not describe how they did this.

Considering the above, for the current research, we decided to choose Sentinel-2 images with a dry soil surface microlayer based on the intensity of spectral reflectance. We based our assessment on spectral reflectance of the dry surface and wet subsurface layer collected in the field, which we subsequently named “dry” and “wet”. After a rain event, the soil surface is wet. Then, it starts drying, and the spectral reflectance curve correspondently changes. At the moment when the spectral reflectance curve stops changing, the soil surface becomes dry. We created a graph wherein we positioned sample points according to their spectral reflectance values measured in the field (spectral reflectance of dry surface and wet subsurface layer) as well as extracted from all Sentinel-2 images. Axes were chosen from informative Sentinel-2 bands included in the regression model to predict SOM content, which was developed in the previous step (Section 2.6).

In the studied spectral range, soil moisture mainly affects the intensity of spectral reflectance in the optical range, meaning that wetter soils have lower spectral reflectance while drier soils have higher spectral reflectance. Thus, on the graph, they occupy two different regions. Dates when the spectral reflectance of sample points extracted from Sentinel data was close to the spectral reflectance of wet subsurface layer were excluded from further analysis. Otherwise, the image was included for further study. We also looked at meteorological data for the test field for 2019. For each date, we mainly checked when the closest rainfall was and its amount. Nocita et al. 2012 [71] showed that the application of the model developed based on spectral reflectance of dry soil to predict organic carbon content in wet soil with a moisture content from 0.05 to more than 0.15 gW gS⁻¹ results in a significant drop in model accuracy. We also checked this in our data. For this, SOM content for each selected Sentinel-2 image was predicted from models for both the dry surface and wet subsurface layer developed from field spectral data. The RMSE calculated for each date for each set of models was used as an indicator of model accuracy.

2.9. Assessment of the Influence Rainfall-Induced Soil Surface Changes on Its Spectral Reflectance and on the Modelling of SOM Content (Only for the Dates with Dry Soil Surface)

The influence of rainfall on a bare soil surface results in the development of a soil crust. Therefore, the longer the soil surface experiences the impact of rainfall, the more washed soil material will be found on its surface and the more developed the soil crust will be. Spectral reflectance and other properties will also change (Figure 5). We considered the spectral reflectance of the sample of the ploughed soil layer as a starting point and the spectral reflectance of the soil crust as the finishing point of soil surface rainfall-induced changes (Figure 5). The spectral reflectance between these two points is a spectral mixture of spectral reflectance of the soil surface with and without a soil crust.

We calculated the spectral reflectance of the spectral mixtures with different fractions of soil crust, using a model for the linear spectral mixture:

SRss = a ∗ SR_pl + (1 − a) ∗ SR_sc,(1)

where SRss is the modeled spectral reflectance of spectral mixture in Sentinel-2 bands, SR_pl is the spectral reflectance of the sample of the ploughed soil layer convolved into Sentinel-2 bands, SR_sc is the spectral reflectance of the soil crust convolved into Sentinel-2 bands, a is the fraction of the soil surface without soil crust, and (1 − a) is the fraction of the soil crust.

Then, for each sample point, we created spectral curves from the spectral reflectance data extracted from Sentinel-2 images. Afterwards, for each point, we calculated the deviation of its spectral reflectance from the spectral reflectance of the dry surface microlayer measured in the field on 15 August 2019 and averaged the results by the date of Sentinel-2 data acquisition. We also calculated the deviation between the spectral reflectance of the dry soil surface microlayer measured in the field and spectral reflectance of the spectral mixtures with various fractions of the soil crust to relate the spectral deviation for each Sentinel-2 date to the degree of rainfall-induced soil surface changes. The higher the deviation of spectral reflectance for a certain date was from the field spectral reflectance of dry surface microlayer, the more soil surface on this date differed from the soil surface at the time of field data collection.

The spectral reflectance of spectral mixtures with different fractions of the soil crust as well as the APD was used to relate the observed differences to the soil surface transformation under the impact of rainfall. Then, we modeled the SOM content for each Sentinel-2 image using informative parameters obtained from the field spectral data analysis (Section 2.6). We further analyzed how the deviations in the spectral reflectance of the soil surface at different dates of Sentinel-2 image acquisition from the spectral reflectance of dry soil surface microlayer measured in the field on 15 August 2019 impacted the parameters, accuracy, and quality of the developed models. The comparison of the models’ accuracy and quality was based on parameters such as the R²cv, RMSEPv, and RPIQ (see Section 2.7). We also checked for the spatial correlation of the residuals in the same way as described in Section 2.7.

2.10. Accounting for Soil Surface State When Mapping SOM Content Based on Sentinel-2 Data

As rainfall splash results in the formation of the soil crust, we need to know its area to account for the soil surface state on a particular date. As we excluded dates with a wet soil surface microlayer (see Section 2.8), we assumed the soil surface microlayer to be dry on the rest of the images and the impact of soil moisture content on spectral reflectance, if there was any, to be minimal.

To determine the area of the soil crust, we performed linear spectral unmixing (function unmix of package hsdar of R software) for each analyzed Sentinel-2 image. As endmembers, we used spectral reflectance of the soil crust and soil sample of the ploughed horizon (Figure 5) convolved into Sentinel-2 bands. A least square solution was performed during spectral unmixing. The error representing the difference between the image spectral reflectance and spectral reflectance predicted by the linear spectral mixture model for each pixel was also calculated. The error value gives the Euclidean norm of the error vector after least square error minimization. This error is caused by the components that were not included as endmembers in the linear unmixing procedure. In our model, we did not account for shadows and soil surface cracks. Therefore, we suggest that the error can be used as an indicator of the area of shadows and/or cracks on the soil surface. As a result, for each date, we obtained the area of soil crust and error for the test field. Then, we extracted the values of these parameters for pixels where we had collected spectral data and soil samples in the field.

Next, we included the error and the area of soil crust as predictors into the models for SOM content prediction from Sentinel-2 data together with the informative parameter (parameters) determined during the analysis of spectral reflectance of the dry surface microlayer (see Section 2.7). Then, we removed the parameters with p-values above 0.05 from the models. We also checked for the spatial correlation of the residuals. The cross-validation, estimation of the new models’ accuracy and the check of residuals for spatial correlation were performed the same way as described in Section 2.6.

Periodical tillage is a part of complete fallow cultivation agrotechnology. Thus, the soil crust microlayer forming on the surface of complete fallow field is periodically disrupted and mixed with the rest of the ploughed soil layer. As a result, there are several cycles of soil crust formation during the season. If the tillage happened close to the date of Sentinel-2 image acquisition, we would observe low values of crust area and high values of error due to the increase in soil surface roughness.

3. Results and Discussion

3.1. Modelling of SOM Content Based on Field Spectral Data

The SOM content of soil samples of the upper soil layer of the test field varied from 1.5% to 7.3% with an average value of 3.2% (Table 1). The best model for SOM content prediction developed based on the field spectral reflectance of the dry surface microlayer included only one parameter, the coloration index (CI), which is calculated using spectral reflectance in the green and red spectral bands:

CI = (Band 4 − Band 3)/(Band 4 + Band 3)(2)

where Band 3 and Band 4 refer to the field spectral reflectance recalculated into corresponding Sentinel-2 bands.

When we used CI to model SOM content from the field spectral reflectance of the wet subsurface layer, we found that two models had a similar accuracy and prediction ability (Figure 7). Natural log was used in the model because the relationship between SOM content and CI was better described by this function.

Originally, CI was related to soil color [72]. Later, Galvao and Vitorello [73] showed that organic matter has the highest negative correlation with the soil spectral reflectance in the region between 550 and 700 nm. Green and red spectral regions involved in CI calculation were also found to be informative for SOM prediction in a number of studies [16,74,75,76].

Residuals for models developed for both dry and wet subsurface layers did not show spatial autocorrelation (Figure 8).

3.2. Assessment of Soil Surface Status at the Time of Satellite Data Acquisition

The dry soil surface layer has a higher spectral reflectance than that of the wet subsurface layer (Figure 9). The deviation of spectral reflectance from the average in both cases increases when moving from the visible to near-infrared spectral range as this range is more sensitive to the variations in soil moisture content. Therefore, the visible range is suitable to distinguish between the dry and wet soil surface independent of differences in moisture content.

When reflectance values for all Sentinel-2 data and proximal sensing data were positioned in the Band 3 to Band 4 spectral space, we found that two distinct clouds of points were formed (Figure 10). The first cloud was concentrated around the points corresponding to the wet subsurface layer; the second cloud was located around the points corresponding to dry surface layer. Therefore, we considered that the first cloud was formed from the points where the soil surface was wet at the time of Sentinel-2 data acquisition and the second cloud was formed from the points where the soil surface was dry.

The point cloud for the wet soil surface showed less scattering than the point cloud for the dry soil surface. This can be explained by the masking effect that the soil moisture has on spatial variations of the soil spectral reflectance.

Scattering within the point cloud for the dry soil surface could have been caused by the variation in the state of the soil surface between different dates of Sentinel-2 data acquisition and, for each date, by the differences in the dry soil surface state between sample points.

The analysis of meteorological data showed that our assignment of the soil surface to dry and wet classes based on its reflectance in Band 3 and Band 4 was correct (Figure 11). The dates 17 April 2019, 6 June 2019, 19 June 2019, and 28 August 2019 fell in the middle of the dry period in between rainfall. The dates 2 April 2019 and 17 April 2019 were in the middle or immediately after the long rain period; thus, the soil surface did not have enough time to dry out. As for 5 May 2019, it fell at the end of a short period with low precipitation. Rainfall at that time was registered in the early morning or in the second half of the day, meaning that the soil surface was dry during the majority of daylight hours.

The highest accumulated atmospheric precipitation was observed in the period from 19 June 2019 to 28 August 2019, meaning that the soil surface had undergone maximal rainfall-induced transformation by 28 August 2019 compared to the previous dates.

The application of the regression model developed for the dry surface layer for the dates that fell in the point cloud for the wet surface resulted in an RMSE almost two times higher than when the model developed for the wet subsurface layer was used for the same dates (Table 3). The application of the regression model developed for the wet subsurface layer for the dates that fell into the point cloud for the dry surface led to greater RMSE values than when we used the model developed for the dry surface layer for these dates.

The image acquired on 19 June 2019 was the exception. The soil surface on that date was classified as dry but the application of the regression model developed for the dry surface resulted in a higher RMSE than when we used the model developed for the wet subsurface layer for this date.

3.3. Influence of Rainfall-Induced Soil Surface Changes on Its Spectral Reflectance and on the Modelling of SOM Content

The spectral reflectance of the soil surface microlayer on 5 May 2019 and 19 June 2019 in bands used for SOM prediction was very close to the spectral reflectance of the dry soil surface microlayer measured in the field (Figure 12). The models of SOM content for these dates had the highest accuracy and the lowest standard errors of the model parameters (Table 4). The APD was the lowest: 1.15 for 5 May 2019 and 2.09 for 19 June 2019.

The average deviation of spectral reflectance of the soil surface microlayer from the average spectral curve derived from field data was the highest on 28 August 2019 (Figure 12). At the same time, it was between the spectral reflectance of the spectral mix with 20% and 10% of the soil crust in Band 3 and Band 4 used for SOM content modeling. The model parameters for this date were similar to the parameters of the model developed from field spectral data for the dry surface microlayer, while the standard errors for model parameters were twice as high and the model accuracy was 12% lower (Table 4). The APD was the highest for this date (3.7).

The spectral reflectance for 20 April 2019 and 6 June 2019 was close to the spectral reflectance of the spectral mixes with 10% of the soil crust (Figure 12). The model parameters for these dates differed the most from the parameters of the model based on the field spectral reflectance of the dry soil surface microlayer, and the standard errors of the model parameters were rather high (Table 4). The regression model for 20 April 2019 had the lowest accuracy among all models developed from Sentinel-2 data. The APD for 20 April 2019 was 2.43. As for 6 June 2019, despite having an R²cv of 0.75, the RPIQ was rather low, demonstrating the poor performance of the model for this date. The APD for 6 June 2019 was 2.7.

Spatial autocorrelation of the residuals was observed for the models developed for 20 April 2019 and 28 August 2019 (Figure 13a,e). Residuals of the models for the rest of the dates did not show spatial autocorrelation. (Figure 13b–d).

3.4. Accounting for Soil Surface State When Modelling SOM Content Based on Sentinel-2 Data

The highest mean and maximum area of the soil crust at sample points were registered on 28 August 2019 followed by 6 June 2019 (Table 5). The APD was also the highest for these dates. The lowest mean area of the soil crust was estimated on 5 May 2019 followed by 19 June 2019 (Figure 14a). The APD for these dates was also the lowest. The variation coefficient for the dates with a low crust area was rather high (400 to 100%), demonstrating high spatial heterogeneity of soil crust development (Figure 14c). As for the error, it was the highest for the dates with a low crust area and APD value (Figure 14b).

At the field level, we observed a similar tendency. A lower mean area of the soil crust was registered for the dates with lower APD (Table 5, Figure 14a), while higher values were determined for the dates with the higher APD (Table 5, Figure 14c). As for the error, the opposite pattern was determined. Maximum values of the parameters were higher when calculated at the field level because more pixels were included in the analysis. Higher maximum values of the soil crust area for the dates with lower APD mainly occurred in the area along northern field boundary (Figure 14a), which is close to the road (Figure 1).

When we included the information on the soil surface state in the models, their quality and accuracy increased (Table 6). For the scene from 20 April 2019, the R²cv increased from 0.67 to 0.80, the RPIQ increased from 1.4 to 2.1, and the RMSEcv decreased from 1.1 to 0.71 (Table 4 and Table 6). For scenes from 5 May 2019 and 6 June 2019, the R²cv increased from 0.77 and 0.75 to 0.79 and 0.81, respectively; the RPIQ increased from 2.2 and 1.7 to 2.4 and 2.3, respectively; and the RMSEcv decreased from 0.68 and 0.85 to 0.63 and 0.61, respectively (Table 4 and Table 6). For the scene from 28 August 2019, the R²cv and RPIQ increased from 0.71 and 1.4 to 0.84 and 2.4, respectively, while the RMSEcv decreased from 1.04 to 0.62 (Table 4 and Table 6).

The error related to the area of shadows (including soil surface cracks) was informative for most analyzed dates. The area of soil crust became a significant parameter only when its mean value reached 10% (Table 5). The APD for these dates was the highest. Moreover, in the case of 28 August 2019, the area of soil crust was the only parameter in the model. As for 19 June 2019, the quality of the model did not improve after adding new predictors; it stayed at the level of the model based on the field spectral reflectance of the dry soil surface microlayer. Incorporation of information on the soil surface state in the model resulted in an increase in model accuracy up to the level of the model based on the field spectral reflectance of the dry soil surface microlayer. Model accuracy stayed at the same level independent of the date of Sentinel-2 image acquisition.

We did not observe spatial autocorrelation of the residuals for most models developed with information on the soil surface state (Figure 15b–e). Only the residuals of the model for 20 April 2019 showed an autocorrelation pattern (Figure 15a). However, it was weaker when compared to the model without information on the soil surface state for this date (Figure 13a). Residuals of the model with an area of soil crust as a single predictor developed for 28 August 2019 did not show spatial autocorrelation in contrast to the model including only CI as a predictor of SOM content for this date (Figure 13e).

4. Discussion

The soil crust formed under the impact of rainfall is a microlayer, which differs in microstructure, texture, and mineralogy from the underlying soil horizon as well as in soil spectral reflectance [36,38,41,77]. Several studies tried to connect the changes in soil surface properties due to the formation of soil crust to the changes in soil spectral reflectance. For example, [38] related the difference in spectral reflectance between the soil crust and underlying layer in the range of 1.2 to 2.5 µm to changes in mineralogy and soil texture during crust formation. They based their conclusions on the analysis of the absorption peaks. They also demonstrated that changes in soil spectral reflectance due to soil crust formation can be used to model the infiltration rate. The authors of [77] also related observed changes in soil spectral reflectance during crust formation to changes in soil texture and mineralogy. For their research, they choose soils of different texture and carbonate content and low SOM content in the range of 0.2 to 0.8%.

In contrast, the authors of [78] in their study found no difference in soil texture and organic matter content between soil crust and control sample. Thus, they attributed the observed changes in soil spectral reflectance solely to the impact of soil surface roughness. However, they did not specify the depth of the sampling of the soil crust for the analysis of SOM content and soil texture, only mentioning that the depth of the layer subjected to the artificial rainfall was 2 cm. Therefore, the absence of a difference in SOM content as well as in soil texture between the soil crust and control sample could be due to the fact that the analyzed samples included both a soil crust microlayer and the underlining layer. When the process is happening at a microlayer, the collection of soil samples for the estimation of soil properties of that microlayer should be very precise. Otherwise, the results can be misleading. Previous studies showed that both SOM content and its composition affect soil’s spectral reflectance [14,72,79,80]. At the same time, atmospheric precipitation alters soil’s surface composition, including soil organic matter [39,81]. Usually, when developing models for the prediction of SOM content based on optical remote sensing data, the relationship is established between the SOM content of the mixed sample of the upper layer (0 to 5 cm, 0 to 15 cm, or 0 to 20 cm) and spectral reflectance of the soil surface measured once [15,27,82]. Moreover, in some cases, satellite data used in modelling were acquired several years before the collection of the field data. This makes the developed models not only site-specific but also not reproducible. The reproducibility of the models for SOM content estimation is very important for soil monitoring and in precision agriculture.

In the current study, the CI obtained from proximal and satellite remote sensing data was related to the content and composition of SOM in the upper soil layer of the test field. Raindrop action (splash) causes the breakdown of soil surface aggregates, containing about 90% of soil surface organic carbon [83,84]. As different conditions are formed inside and outside soil aggregates, the organic matter of the surface of the aggregates is different from the organic matter of their inner part [85]. Therefore, we suggest that the registered deterioration of the accuracy and quality of models developed from Sentinel-2 data was caused by the differences in SOM content and composition between soil surface at the time of satellite data acquisition and the upper soil layer. Moreover, in those cases, we observed a shift in the average satellite spectral curves from the average spectral curve derived from the field spectral reflectance of dry soil surface microlayer to spectral mixes with different fractions of soil crust, which generally demonstrated higher spectral reflectance. The degree of deviation increased with the APD. The increase in spectral reflectance of the soil surface microlayer as a result of soil crust formation and development was observed by [76] in their experiment with four soil types differing in soil texture.

While field spectra were measured on 15 August 2019, almost at the end of the study, it did not represent the endpoint with maximum soil crusting as can be expected (Figure 12). This is also true for Sentinel-2 data. For example, the average spectral curves for 20 April 2019 and 6 June 2019 are close to the spectral mix with 10% of soil crust, while the spectral curves for 5 May 2019 and 19 June 2019 are close to the average spectral curve derived from the field spectral reflectance of the dry soil surface microlayer (Figure 12). This was probably caused by periodical tillage operations occurring several days before the Sentinel-2 image acquisition on 5 May 2019 and 19 June 2019 as well as prior to field data collection. The accuracy of the models for these dates is the closest to the accuracy of the model derived from the field spectral reflectance of dry soil surface microlayer (Table 4). According to the results of spectral unmixing, the average area of soil crust was the lowest for these dates, while the average error, related to the area of shadows and/or cracks, was the highest (Table 5). The soil crust microlayer can be easily destroyed during tillage operations, and mixed with the rest of ploughed horizon. After tillage, when the soil surface experiences the impact of rainfall, the soil crust microlayer starts to form again. There can be several cycles of soil crust microlayer formation during one season. Depending on how long the soil surface is left bare, and the intensity, and the amount of rainfall, the degree of soil crust microlayer development in each cycle can vary. This makes it even more important to estimate the soil surface state at the time of spectral data acquisition to ensure the stability and reproducibility of the developed models.

The results presented in Table 3 are in accordance with the results obtained by [71]. They show that the application of the model developed for the dry surface to predict SOM content for the wet surface resulted in a significant drop in the model’s accuracy. This supports our decision to exclude certain images prior to the analysis to minimize the interference of soil moisture content. The results of spectral unmixing showed that variations of the soil surface state both in time and within the field can be rather high. We attribute the observed temporal variations to rainfall patterns. The area of soil crust on the studied field increased with the APD. We observed similar results in our model experiment (Figure 3 and Figure 4), in which the area of the soil crust increased with cumulative rainfall. At the same time, the APD had the opposite effect on the error, which is related to the area of shadows (including surface cracks). As rainfall causes the breakdown of soil surface aggregates, the amount and size of large aggregates that can produce shadows decreases under its impact, while the soil crust becomes more pronounced [36,40]. Additionally, cracks in the soil surface can disappear after heavy rainfall, and the cracking process will begin again when the soil surface starts to dry. Vindeker et al. [40] demonstrated that the position of cracks changes in time, and the soil surface usually undergoes several cycles of transformation during the season. Within-field heterogeneity of the soil surface state of the test field is mainly connected to the SOM content. The SOM content affects the speed of soil surface transformation and soil crust development as it determines the stability of soil aggregates. At the beginning of soil crust development, the difference between soils with various SOM content will be low. While the soil crust develops, the difference between areas with higher and lower SOM content will increase due to the variation in the aggregate stability. Aggregates of soils with higher SOM content are more resistant to the impact of raindrops, and soil crust formation will be slower compared to soils with lower SOM content. Therefore, the longer the soil surface is exposed to the rainfall, the higher the spatial heterogeneity of soil crust development will be. As this heterogeneity is caused by the difference in SOM content, at a certain stage of soil crust development, the soil crust area will be directly related to SOM content. We observed such a situation for 28 August 2019, when the soil crust area became the only predictor in the model for SOM content estimation (Table 6). The lowest area of the soil crust in the test field on this date was observed as expected in the areas with high SOM content (Figure 14c).

Our results on the informativeness of the area occupied by the shadows at certain dates when modelling SOM content are in accordance with the results obtained by [27]. They showed that correction of the effect of shadows allowed for the improvement of models for soil organic carbon by 27% for field spectral data and by 25% for airborne spectral data. However, they used only the data collected during the period of 5 to 6 October and did not study how the impact of the shadows on the modelling of soil organic carbon changes in time when the soil surface experiences a more prolonged impact of rainfall. According to our results (Table 4 and Table 6), the incorporation of the information on the area occupied by the shadows in the models leads to a decrease in the RMSEcv by 7 to 35% depending on the date of image acquisition.

Here, we used a simple approach to account for the rainfall-induced changes in the soil surface when mapping SOM content from Sentinel-2 data. The approach needs further development and testing in different environments and on soils with different textures and mineralogy before we suggest its wide application. However, it is valuable for loamy soils formed on loess in similar climatic conditions (the temperate climatic zone), as soil texture in non-saline soils predetermines the mechanism of soil crust formation [33]. In our model experiment in 2017, we studied two different loamy soils (grey forest (from the test field) and leached chernozem (Luvic Chernozems according to [43]) and found that the transformations that occurred under the impact of rainfall were common [40,41]. The difference was only in the degree of soil aggregate disruption and the amount of washed material. This depends mainly on the organic matter content, because the higher the organic matter content is, the more stable the soil aggregates are.

Additional studies are necessary to identify how the SOM content and composition of the soil surface microlayer changes under the influence of rainfall. It is also important to determine the relationships among rainfall-induced SOM, mineralogy changes, and soil surface spectral reflectance.

5. Conclusions

The impact of rainfall leads to a discrepancy in the spectral reflectance and properties between the soil surface microlayer and the rest of the ploughed horizon. Our study showed that this rainfall-driven discrepancy negatively affects the accuracy and quality of models developed to map SOM content from Sentinel-2 data. The higher the rainfall impact was, the more the soil surface differed from the rest of the ploughed layer. The incorporation of information on the soil surface state in the model allowed the compensation of negative rainfall impact and the development of stable models with higher quality and accuracy.

Accounting for the rainfall-induced changes in the soil surface is necessary to ensure the reproducibility of the developed models and the correct assessment of organic matter in the soil ploughed layer. This is a crucial issue in SOM detection and mapping based on remote sensing data.

Author Contributions

Conceptualization, methodology, writing—review and editing, supervision, project administration, resources, data curation, and funding acquisition, I.S.; writing—original draft preparation, methodology, visualization, validation, formal analysis, and investigation, E.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Ministry of Science and Higher Education of Russia (agreement No 075-15-2020-909).

Acknowledgments

This paper has been supported by the RUDN University Strategic Academic Leadership Program. We would like to thank the Center of collective usage of equipment of the V.V. Dokuchaev Soil Science Institute for their help with laboratory analysis of soil samples and field data collection.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A

Table A1

Spectral indices used in the study.

Index	Equation	Equation Adopted to Sentinel-2 bands	Reference
BI (Brightness Index)	$\sqrt{\frac{\left(B^2\hspace{0.17em}+\hspace{0.17em}{G}^2\hspace{0.17em}+\hspace{0.17em}{R}^2\right)}{3}}$	$\sqrt{\frac{\left(Band\hspace{0.17em}{2}^2\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}{3}^2\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}{4}^2\right)}{3}}$	[58]
CI (Colouration Index)	$\frac{\left(R\hspace{0.17em}-\hspace{0.17em}G\right)}{\left(R\hspace{0.17em}+\hspace{0.17em}G\right)}$	$\frac{\left(Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}3\right)}{\left(Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}3\right)}$	[58]
HI (Hue Index)	$\frac{\left(2\hspace{0.17em}\ast \hspace{0.17em}R\hspace{0.17em}-\hspace{0.17em}G\hspace{0.17em}-\hspace{0.17em}B\right)}{\left(G\hspace{0.17em}-\hspace{0.17em}B\right)}$	$\frac{\left(2\hspace{0.17em}\ast \hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}3\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}2\right)}{\left(Band\hspace{0.17em}3\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}2\right)}$	[58]
SI (Saturation Index)	$\frac{\left(R\hspace{0.17em}-\hspace{0.17em}B\right)}{\left(R\hspace{0.17em}+\hspace{0.17em}B\right)}$	$\frac{\left(Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}2\right)}{\left(Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}2\right)}$	[58]
RI (Redness Index)	$\frac{R^2}{\left(B\hspace{0.17em}\ast \hspace{0.17em}{G}^3\right)}$	$\frac{Band\hspace{0.17em}{4}^2}{\left(Band\hspace{0.17em}2\hspace{0.17em}\ast \hspace{0.17em}Band\hspace{0.17em}{3}^3\right)}$	[58]
BSCI (biological soil crust index)	$\frac{1\hspace{0.17em}-\hspace{0.17em}2\hspace{0.17em}\ast \hspace{0.17em}\left|R\hspace{0.17em}-\hspace{0.17em}G\right|}{(\frac{R\hspace{0.17em}+\hspace{0.17em}G\hspace{0.17em}+\hspace{0.17em}NIR}{3})}$	$$\begin{gathered} \frac{1\hspace{0.17em}-\hspace{0.17em}2\hspace{0.17em}\ast \hspace{0.17em}\left|Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}3\right|}{(\frac{Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}3\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}5}{3})}; \\ \frac{1\hspace{0.17em}-\hspace{0.17em}2\hspace{0.17em}\ast \hspace{0.17em}\left|Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}3\right|}{(\frac{Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}3\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}6}{3})}; \\ \frac{1\hspace{0.17em}-\hspace{0.17em}2\hspace{0.17em}\ast \hspace{0.17em}\left|Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}3\right|}{(\frac{Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}3\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}7}{3})}; \\ \frac{1\hspace{0.17em}-\hspace{0.17em}2\hspace{0.17em}\ast \hspace{0.17em}\left|Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}3\right|}{(\frac{Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}3\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}8}{3})}; \\ \frac{1\hspace{0.17em}-\hspace{0.17em}2\hspace{0.17em}\ast \hspace{0.17em}\left|Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}3\right|}{(\frac{Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}3\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}8a}{3})} \end{gathered}$$	[59]
GSI (Grain Size Index)	$\frac{\left(R\hspace{0.17em}-\hspace{0.17em}B\right)}{\left(R\hspace{0.17em}+\hspace{0.17em}B\hspace{0.17em}+\hspace{0.17em}G\right)}$	$\frac{\left(Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}2\right)}{\left(Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}3\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}2\right)}$	[60]
CRUST	$1\hspace{0.17em}-\hspace{0.17em}\frac{\left(R\hspace{0.17em}-\hspace{0.17em}B\right)}{\left(R\hspace{0.17em}+\hspace{0.17em}B\right)}$	$1\hspace{0.17em}-\hspace{0.17em}\frac{\left(Band\hspace{0.17em}4\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}2\right)}{\left(Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}2\right)}$	[57]
AVI (Advanced vegetation index)	$\sqrt[3]{\left(NIR\hspace{0.17em}+\hspace{0.17em}1\right)\hspace{0.17em}\ast \hspace{0.17em}\left(1\hspace{0.17em}-\hspace{0.17em}R\right)\hspace{0.17em}\ast \hspace{0.17em}\left(NIR\hspace{0.17em}-\hspace{0.17em}R\right)}$	$\sqrt[3]{\left(Band\hspace{0.17em}8\hspace{0.17em}+\hspace{0.17em}1\right)\hspace{0.17em}\ast \hspace{0.17em}\left(1\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)\hspace{0.17em}\ast \hspace{0.17em}\left(Band\hspace{0.17em}8\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}$	[63]
NDVI (Normalized difference vegetation index)	$\frac{\left(NIR\hspace{0.17em}-\hspace{0.17em}R\right)}{\left(NIR\hspace{0.17em}+\hspace{0.17em}R\right)}$	$$\begin{gathered} \frac{\left(Band\hspace{0.17em}8\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}8\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}4\right)}; \\ \frac{\left(Band\hspace{0.17em}8\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}5\right)}{\left(Band\hspace{0.17em}8\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}5\right)}; \\ \frac{\left(Band\hspace{0.17em}8\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}6\right)}{\left(Band\hspace{0.17em}8\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}6\right)}; \\ \frac{\left(Band\hspace{0.17em}8\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}7\right)}{\left(Band\hspace{0.17em}8\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}7\right)}; \\ \frac{\left(Band\hspace{0.17em}8\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}8a\right)}{\left(Band\hspace{0.17em}8\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}8a\right)}; \end{gathered}$$	[64]
SAVI (Soil-adjusted vegetation index)	$1.5\ast \frac{\left(NIR\hspace{0.17em}-\hspace{0.17em}R\right)}{\left(NIR\hspace{0.17em}+\hspace{0.17em}R\hspace{0.17em}+\hspace{0.17em}0.5\right)}$	$$\begin{gathered} 1.5\ast \frac{\left(Band\hspace{0.17em}5\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}5\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}0.5\right)}; \\ 1.5\ast \frac{\left(Band\hspace{0.17em}6\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}6\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}0.5\right)}; \\ 1.5\ast \frac{\left(Band\hspace{0.17em}7\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}7\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}0.5\right)}; \\ 1.5\ast \frac{\left(Band\hspace{0.17em}8\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}8\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}0.5\right)}; \\ 1.5\ast \frac{\left(Band\hspace{0.17em}8a\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}8a\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}0.5\right)}; \end{gathered}$$	[56]
EVI2 (Enhanced vegetation index)	$2.5\ast \frac{\left(NIR\hspace{0.17em}-\hspace{0.17em}R\right)}{\left(NIR\hspace{0.17em}+\hspace{0.17em}2.4\hspace{0.17em}\ast \hspace{0.17em}R\hspace{0.17em}+\hspace{0.17em}1\right)}$	$$\begin{gathered} 2.5\ast \frac{\left(Band\hspace{0.17em}5\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}5\hspace{0.17em}+\hspace{0.17em}2.4\hspace{0.17em}\ast \hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}1\right)}; \\ 2.5\ast \frac{\left(Band\hspace{0.17em}6\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}6\hspace{0.17em}+\hspace{0.17em}2.4\hspace{0.17em}\ast \hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}1\right)}; \\ 2.5\ast \frac{\left(Band\hspace{0.17em}7\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}7\hspace{0.17em}+\hspace{0.17em}2.4\hspace{0.17em}\ast \hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}1\right)}; \\ 2.5\ast \frac{\left(Band\hspace{0.17em}8\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}8\hspace{0.17em}+\hspace{0.17em}2.4\hspace{0.17em}\ast \hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}1\right)}; \\ 2.5\ast \frac{\left(Band\hspace{0.17em}8a\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}4\right)}{\left(Band\hspace{0.17em}8a\hspace{0.17em}+\hspace{0.17em}2.4\hspace{0.17em}\ast \hspace{0.17em}Band\hspace{0.17em}4\hspace{0.17em}+\hspace{0.17em}1\right)}; \end{gathered}$$	[61]
NDWI (Normalized difference water index)	$\frac{\left(G\hspace{0.17em}-\hspace{0.17em}NIR\right)}{\left(G\hspace{0.17em}+\hspace{0.17em}NIR\right)}$	$\frac{\left(Band\hspace{0.17em}3\hspace{0.17em}-\hspace{0.17em}Band\hspace{0.17em}8\right)}{\left(Band\hspace{0.17em}3\hspace{0.17em}+\hspace{0.17em}Band\hspace{0.17em}8\right)}$	[62]

B—spectral reflectance in the blue part of the spectrum; G—spectral reflectance in the green part of the spectrum; R—spectral reflectance in the red part of the spectrum; NIR—spectral reflectance in the near-infrared part of the spectrum.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Geographical location of the test plot and sampling scheme: the boundaries of the test field are delineated in red color.

Figure 2. Soil crust on the surface of the studied soil: (a) photo of soil crust of grey forest soil collected on watershed of the test field. The photo was made on an iPhone XR located above the soil surface; (b) microstructure of vertical thin section of upper 5 mm by an Olympus BX51 polarizing microscope with an Olympus DP26 digital camera; (c) model of micromonolith of upper 4 cm soil layer showing inner structure and pore space. The model was produced by X-ray microtomography using a SkyScan 1172G. Red lines on (b,c) delineate the boundary between soil crust and the underlying layer.

Figure 3. Photos of soil surface of grey forest soil in model experiment [40]: (a) 7.06.2017 (the start of the experiment); (b) 16.06.2017; (c) 26.06.2017; (d) 30.06.2017; (e) 07.07.2017; (f) 11.07.2017; (g) 25.07.2017.

Figure 4. Change in soil surface status of grey forest soil in model experiment; R1—the amount of rainfall from the last measurement, R2 —accumulated rainfall from the start of the experiment. The projected area, occupied by washed soil material, was determined from photos of the soil surface taken during the experiment [40]. The information on rainfall was collected from the meteorological station near the experiment site (http://www.pogodaiklimat.ru/weather.php?id=27605 (accessed on 25 May 2021)).

Figure 5. Spectral reflectance of soil samples with different state of soil surface. SR_ds and SR_ws—average spectral reflectance of dry soil surface and wet subsurface layer of 30 sample points, respectively, measured in the field on 15 August 2019; SR_sc—spectral reflectance of soil crust; SR_pl—spectral reflectance of sample of ploughed soil layer.

Figure 6. Process flow chart: (1) modelling of SOM content from the field data and selection of Sentinel-2 images with dry soil surface; (2) modelling SOM content from Sentinel-2 images with dry soil surface; (3) accounting for soil surface state when modelling SOM content from Sentinel-2 data for the dates with dry soil surface. SR_ds and SR_ws—spectral reflectance of dry soil surface and wet subsurface layer, respectively, measured in the field on 15 August 2019; SR_sc—spectral reflectance of soil crust, SR_pl—spectral reflectance of sample of ploughed soil layer.

Figure 7. Relationship between CI and SOM content for the dry surface and wet subsurface layer based on the field data collected on 15 August 2019.

Figure 8. Spatial autocorrelation within model residuals for: (a) dry surface layer; (b) wet subsurface layer.

Figure 9. Spectral reflectance of dry surface and wet subsurface layers measured in the field: dry, wet—average spectral reflectance of dry surface and wet subsurface layers correspondingly for 30 sample points; dry+sd, wet+sd—average spectral reflectance of dry surface and wet subsurface layers plus standard deviation of soil surface spectral reflectance of the corresponding layer; dry-sd, wet-sd—average spectral reflectance of dry surface and wet subsurface layers minus standard deviation of soil spectral reflectance of the corresponding layer.

Figure 10. Spectral reflectance for sample points measured with a HandHeld-2 field spectroradiometer and registered by Sentiniel-2 data in the Band 3 to Band 4 spectral space: dry soil surface layer—red points; wet subsurface layer—black points; 2 April 2019—green points; 17 April 2019—yellow points; 20 April 2019—orange points; 5 May 2019—purple points; 6 June 2019—brown points; 19 June 2019—dark green points; 28 August 2019—dark blue points; the red oval delineates points with a dry soil surface; black oval—with a wet soil surface.

Figure 11. Rainfall plot for the year 2019 for the test field: red arrows indicate Sentinal-2 image acquisition on 2 April 2019—day 92; 17 April 2019—day 107; 20 April 2019—day 110; 5 May 2019—day 125; 6 June 2019—day 157; 19 June 2019—day 170; 28 August 2019—day 240. Values above the colored blocks show atmospheric precipitation accumulated between the studied dates. Rainfall amount (mm) is presented on the OY axis and the number of a day in the year on the OX axis.

Figure 12. Average deviation of spectral reflectance of the bare soil surface microlayer of sample points at the dates of Sentinel-2 data acquisition from the spectral reflectance of dry soil surface microlayer on 15 August 2019, and spectral mixtures with different soil crust fractions. Central wavelength for each analyzed Sentinel-2 band is given in brackets.

Figure 13. Spatial autocorrelation within model residuals for: (a) 20 April 2019; (b) 5 May 2019; (c) 6 June 2019; (d) 19 June 2019; (e) 28 August 2019.

Figure 14. Spatial distribution of error (b,d) and the area of soil crust (a,c) determined as a result of spectral unmixing on the test field for 5 May 2019 (a,b) and 28 August 2019.

Figure 15. Spatial autocorrelation within model residuals for models incorporating data on soil surface state for: (a) 20 April 2019; (b) 5 May 2019; (c) 6 June 2019; (d) 19 June 2019; (e) 28 August 2019.

Figure 15. Spatial autocorrelation within model residuals for models incorporating data on soil surface state for: (a) 20 April 2019; (b) 5 May 2019; (c) 6 June 2019; (d) 19 June 2019; (e) 28 August 2019.

Table 1

Chemical properties of grey forest soils of the test field (the data for 30 samples collected on 15 August 2015).

Descriptive Statistics	pH of Salt Extract	Mobile P (mg/kg)	Mobile K (mg/kg)	SOM (%)	Exchangeable Cations
					Ca++ (mmol/100 g)	Mg++ (mmol/100 g)	Na+ (mg/kg)	K+ (mg/kg)
Minimum	4.9	85.2	141.6	1.5	9.7	1.9	0.0	0.3
Maximum	6.9	3541.4	395.0	7.3	47.6	3.3	0.1	0.8
Average	5.3	536.0	227.0	3.2	19.6	2.7	0.1	0.4
Standard deviation	0.6	943.3	61.8	1.5	10.2	0.4	0.0	0.1
Coefficient of variation (%)	11.9	176.0	27.2	46.1	52.4	15.0	26.2	29.6

Table 2

Sentinel-2 band characteristics.

Sentinel-2 Bands	Central Wavelength (nm)	Band Width (nm)	Resolution (m)
Band 1—coastal aerosol	443	20	60
Band 2—blue	490	65	10
Band 3—green	560	35	10
Band 4—red	665	30	10
Band 5—vegetation red edge	705	15	20
Band 6—vegetation red edge	740	15	20
Band 7—vegetation red edge	783	20	20
Band 8—near infrared	842	115	10
Band 8a—vegetation red edge	865	20	20
Band 9—water vapor	945	20	60
Band 10—shortwave infrared—Cirrus	1375	30	60
Band 11—shortwave infrared	1610	90	20
Band 12—shortwave infrared	2190	180	20

Table 3

Variations in the root-mean-square error of prediction (RMSE, %) of the SOM content modelling due to differences in the soil surface wetness.

Acquisition Date of Sentinel-2 Images	RMSE When Using Regression Model Developed from Spectral Reflectance of Dry Soil Surface Microlayer	RMSE When Using Regression Model Developed from Spectral Reflectance of Wet Subsurface Layer	Surface Class from Figure 10
2/04/2019	6.78	6.69	wet
17/04/2019	1.95	1.64	wet
20/04/2019	1.97	3.98	dry
5/05/2019	1.08	1.43	dry
6/06/2019	1.21	1.53	dry
19/06/2019	1.61	0.87	dry
28/08/2019	1.10	2.03	dry

Table 4

Variation in model parameters and accuracy when modelling SOM content of the soil surface layer from Sentinel-2 data (only for the dates with a dry soil surface).

Acquisition Date of Sentinel-2 Images	R2cv	RMSEcv	RPIQ	Intercept		Ln (CI)
				Value	St. Error	Value	St. Error
model based on spectral reflectance of dry soil surface microlayer measured in the field on 15 August 2019	0.81	0.64	2.3	−13.35	1.65	−7.77	0.77
20.04.2019	0.67	1.1	1.4	−10.43	3.23	−5.84	1.38
05.05.2019	0.77	0.68	2.2	−12.25	1.73	−7.62	0.85
06.06.2019	0.75	0.85	1.7	−17.44	3.00	−10.18	1.48
19.06.2019	0.78	0.64	2.3	−11.45	1.49	−7.54	0.77
28.08.2019	0.71	1.04	1.4	−13.07	3.32	−7.80	1.59

Table 5

Temporal variations in the area of soil crust and error estimated from Sentinel-2 data for sample points and for the test field (1400 pixels): results for sample points/results for the test field.

Acquisition Date of Sentinel-2 Images	Crust (%)				Error				APD
	Maximum	Mean	Standard Deviation	Variation Coefficient	Maximum	Mean	Standard Deviation	Variation Coefficient
20.04.2019	17/30	5.0/7.0	5/6	100/86	0.006/0.007	0.003/0.002	0.001/0.0008	33/39	2.43
5.05.2019	12/30	0.5/9	2/9	400/100	0.015/0.015	0.005/0.005	0.004/0.003	80/60	1.15
6.06.2019	22/60	10.0/13.0	6/8	60/62	0.009/0.032	0.003/0.003	0.002/0.004	67/133	2.70
19.06.2019	11/60	3.0/6.0	3/8	100/133	0.014/0.033	0.005/0.005	0.003/0.004	60/80	2.09
28.08.2019	34/40	20.0/22.0	11/9	55/41	0.007/0.023	0.003/0.003	0.001/0.002	33/67	3.70

Table 6

Variation in model parameters and accuracy when modelling SOM content with information on soil surface state (only for the dates with a dry soil surface).

Acquisition Date of Sentinel-2 Images	R2cv	RMSEcv	RPIQ	Intercept		Ln (CI)		Error		Crust
				Value	St. Error	Value	St. Error	Value	St. Error	Value	St. Error
20.04.2019	0.80	0.71	2.1	−3.7	2.2	−2.0	1.0	891.3	130.4	-	-
05.05.2019	0.78	0.63	2.4	−5.8	3.2	−4.0	1.8	175.1	73.9	-	-
06.06.2019	0.81	0.61	2.4	−3.0	3.0	−3.1	1.4	293.1	60.2	−9.0	2.3
19.06.2019	0.81	0.64	2.3	−6.6	3.0	−4.6	1.7	157.9	82.3	-	-
28.08.2019	0.84	0.62	2.4	5.4	0.3	-	-	-	-	−11.2	1.1