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

: Optical remote sensing only provides information about the very thin surface layer of soil. Rainfall splash alters soil surface properties and its spectral reﬂectance. 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 ﬁeld located in the Tula region (Russia), their spectral reﬂectance, 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 2 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 2 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 inﬂuence of rainfall for other soils, and to determine the relationships between rainfall-induced SOM changes and soil surface spectral reﬂectance.


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 2 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 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 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.
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). 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.
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.

Soil Crusting Process
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).

(%)
The analyzed soil was formed on loess. The composite sample colle ploughed horizon (0 to 30 cm) of the studied field using the diagonal meth 6.5% of sand (0.05 to 1 mm), 75.6% of silt (0.001 to 0.05 mm) and 17.9% o mm) [50]. This is silt loam in the USDA classification. The conversion from classification to the USDA classification was made based on [51].
The impact of atmospheric precipitation on soils with such a texture cau down of soil surface aggregates, and the formed material is washed and m pores of the soil surface. In Figure 2b, this washed material has a light brow the dark brown color corresponds to the remains of the soil aggregates. In t we look at the surface of the soil crust, the washed soil material has a ve close to white (Figure 2a). The formation of the soil crust in the field goes th stages. They have been defined micromorphologically through the analysis the microstructure [36].
In 2017, we conducted a model experiment, wherein we exposed t ploughed horizons of three soil types, including the studied soil from the te impact of atmospheric precipitation [40,41]. We put soil samples in recta boxes and left them outside, periodically taking photos of the soil surface (F According to our results, which are in accordance with [33], the longer face was subjected to rainfall, the less soil aggregates were left there and the soil material formed as a result of the breakdown of soil surface aggregates 4). Therefore, the amount of washed soil material on the soil surface can indicator of soil crust development. As the soil surface on 11 July 2017 wa 3f), the area of washed soil material was underestimated due to the inter moisture content (Figure 4). The cracks did not appear after the first rainfall. the first cracks after the cumulative rainfall reached 100 mm ( Figure 4).
Rainfall-induced changes in the soil surface microlayer affect the soil s tance. For example, the soil crust of the studied soil (Figure 2a), collected fr shed of the test field in our previous studies, had higher reflectance intensit the sample of the ploughed soil layer ( Figure 5). The difference in spectral 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 (Figures 3 and 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). the sample of the ploughed soil layer ( Figure 5). It shows that tillage occurred close to th date of our field trip, meaning that the upper layer was rather uniform. Thus, the mai difference between the soil surface microlayer and subsurface layer on 15 August 201 was due to the moisture content.  Soil crust is the initial stage of soil degradation. This microlayer is comprised of th remains of soil aggregates, so its original structure is lost. As it is compacted, the wate penetration decreases while the surface runoff increases. The latter adds up to soil erosio development. The loss of aggregate structure, compaction, and increased surface runof contribute to the proliferation of soil degradation. 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. 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. 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.

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 Sections 2.6 and 2.7.

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 Sections 2.6 and 2.7. (1) modelling of SOM content from the field data and selection of SenTable 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_scspectral reflectance of soil crust, SR_pl-spectral reflectance of sample of ploughed soil layer.
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 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.
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, Remote Sens. 2021, 13, 2313 8 of 26 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.  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. Table 2. Sentinel-2 band characteristics.

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.

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 Remote Sens. 2021, 13, 2313 9 of 26 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].

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 2 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 2 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.

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 −1 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. 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: 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 2 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.

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 crossvalidation, 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.

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: 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. 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.

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: where Band3 and Band4 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). 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).

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. 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.
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 corre-

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.

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. 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.
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. . 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; drysd, wet-sd-average spectral reflectance of dry surface and wet subsurface layers minus standard deviation of soil spectral reflectance of the corresponding layer.
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.
Scattering within the point cloud for the dry soil surface could have been ca the variation in the state of the soil surface between different dates of Sentinel-2 quisition and, for each date, by the differences in the dry soil surface state between points. we middle or immediately after the long rain period; thus, the soil surface did n enough time to dry out. As for 5 May 2019, it fell at the end of a short period w precipitation. Rainfall at that time was registered in the early morning or in the half of the day, meaning that the soil surface was dry during the majority of hours.  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 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 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.

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.  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 2 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).  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 2 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).

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). emporal variations in the area of soil crust and error estimated from Sentinel-2 data for sample points and for eld (1400 pixels): results for sample points/results for the test field.

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 2 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 (Tables 4 and 6  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 2 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 (Tables 4 and 6 (Tables 4 and 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).  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).

Discussion
The soil crust formed under the impact of rainfall is a microlayer, w microstructure, texture, and mineralogy from the underlying soil horizon a spectral reflectance [36,38,41,77]. Several studies tried to connect the chan face properties due to the formation of soil crust to the changes in soil spec For example, [38] related the difference in spectral reflectance between th underlying layer in the range of 1.2 to 2.5 μm to changes in mineralogy during crust formation. They based their conclusions on the analysis of peaks. They also demonstrated that changes in soil spectral reflectance d formation can be used to model the infiltration rate. The authors of [77] served changes in soil spectral reflectance during crust formation to change and mineralogy. For their research, they choose soils of different texture 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 organic matter content between soil crust and control sample. Thus, they observed changes in soil spectral reflectance solely to the impact of soil sur However, they did not specify the depth of the sampling of the soil crust of SOM content and soil texture, only mentioning that the depth of the la the artificial rainfall was 2 cm. Therefore, the absence of a difference in S well as in soil texture between the soil crust and control sample could be that the analyzed samples included both a soil crust microlayer and the un When the process is happening at a microlayer, the collection of soil samp mation of soil properties of that microlayer should be very precise. Otherw can be misleading. Previous studies showed that both SOM content and affect soil's spectral reflectance [14,72,79,80]. At the same time, atmospher alters soil's surface composition, including soil organic matter [39,81]. Us veloping models for the prediction of SOM content based on optical remo the relationship is established between the SOM content of the mixed sam layer (0 to 5 cm, 0 to 15 cm, or 0 to 20 cm) and spectral reflectance of the so ured once [15,27,82]. Moreover, in some cases, satellite data used in mod quired several years before the collection of the field data. This makes the d els not only site-specific but also not reproducible. The reproducibility of SOM content estimation is very important for soil monitoring and in precis In the current study, the CI obtained from proximal and satellite remo was related to the content and composition of SOM in the upper soil layer Raindrop action (splash) causes the breakdown of soil surface aggrega about 90% of soil surface organic carbon [83,84]. As different conditions ar and outside soil aggregates, the organic matter of the surface of the aggreg from the organic matter of their inner part [85]. Therefore, we suggest tha

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 (Figures 3 and 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 (Tables 4 and 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.

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. 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) (Band 2 2 + Band 3 2 + Band 4 2 ) 3 [58] CI (Colouration Index) HI (Hue Index) 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.