On the impact of soil texture on local scale organic carbon quantification: From airborne to spaceborne sensing domains

Soil & Tillage Research 241 (2024) 106125 Available online 26 April 2024 0167-1987/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). On the impact of soil texture on local scale organic carbon quantification: From airborne to spaceborne sensing domains Vahid Khosravi a , * , Asa Gholizadeh a , Daniel ˇ Zíˇ zala a , b , Radka Kodeˇ sov´ a a , Mohammadmehdi Saberioon c , d , Prince Chapman Agyeman e , Petra Vokurkov´ a a , Anna Juˇ ricov´ a b , Marko Spasi´ c a , Luboˇ s Borůvka a a Department of Soil Science and Soil Protection, Faculty of Agrobiology, Food and Natural Resources, Czech University of Life Sciences Prague, Kamýck´ a 129, Suchdol, Prague 16500, Czech Republic b Research Institute for Soil and Water Conservation, ˇ Zabovˇ resk´ a 250, Prague 15600, Czech Republic c Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences, Section 1.4 Remote Sensing and Geoinformatics, Telegrafenberg, Potsdam 14473, Germany d ILVO, Flanders Research Institute for Agriculture, Fisheries and Food, Technology and Food Science-Agricultural Engineering, Merelbeke 9820, Belgium e Division of Plant Sciences and Technology, University of Missouri-Columbia, Columbia, MO 65211, USA ARTICLE INFO Keywords: Soil organic carbon Airborne hyperspectral data Sentinel-2 Soil texture Stratification ABSTRACT Soil organic carbon (SOC) distribution and interaction with light is influenced by soil texture parameters (clay, silt and sand), which makes SOC prediction complicated, especially in samples with considerable pedological variability. Hence, understanding the relationship between SOC and soil texture is important within the context of SOC prediction using remote sensing data. The main objective of this study was to find the impact of soil texture on the performance of local SOC prediction models that were developed on Sentinel-2 (S2) multispectral and CASI/SASI (CS) hyperspectral airborne data as the main predictor variables. One approach to that objective was to lowering the texture variance by stratification of the samples. Therefore, soil samples collected from four agricultural sites in the Czech Republic were segregated based on the i) site-based and ii) texture-based stratification strategies. Random forest (RF) models were then developed on all stratified classes with and without considering the soil texture parameters as predictor variables and results were compared with those obtained by the RF models developed on the non-stratified (NS) samples. Both stratification strategies provided more homogeneous classes, which enhanced the accuracy of SOC prediction, compared to using the NS samples. In addition, the texture-based RF models yielded higher accuracy predictions than the site-based ones. Except sand, adding texture parameters to the main predictors improved accuracy of the models, so that the highest prediction performance was obtained by a texture-based model developed on clay added CS data. Overall, texture-based stratification could significantly enhance the SOC prediction, when the texture parameters were added to the S2 and CS data as the main predictor variables. 1. Introduction Soil organic carbon (SOC) monitoring is a critical step for sustainable land use practices and alleviating the climate change impacts (Gautam et al., 2022). Traditional SOC quantification methods are precise but require tedious soil sampling campaigns, expensive equipment and hazardous chemicals (Schillaci et al., 2017a). There has thus been a constant demand for rapid and inexpensive alternative techniques to compensate for drawbacks of those methods. Within recent decades, remote sensing has been capable of quick and non-destructive monitoring of some specific soil properties, including SOC. Spaceborne multispectral and hyperspectral sensors have shown huge capabilities for broad and continuous coverage of the earth surface and production of soil surface maps even from the inaccessible and hard to reach regions (Biney et al., 2021). As an example, Sentinel-2 (S2) with revisit interval of 5–7 days, provides multispectral optical imagery for high spatial resolution terrestrial observations. Encouraging SOC prediction performances have been obtained using S2 data for temperate soils in the Czech Republic (Gholizadeh et al., 2018; ˇ Zíˇ zala et al., 2019), Belgium and Germany (Castaldi et al., 2019a) and Northern France * Corresponding author. E-mail address: [email protected] (V. Khosravi). Contents lists available at ScienceDirect Soil & Tillage Research journal homepage: www.elsevier.com/locate/still https://doi.org/10.1016/j.still.2024.106125 Received 4 October 2023; Received in revised form 30 March 2024; Accepted 14 April 2024 Soil & Tillage Research 241 (2024) 106125 2 (Vaudour et al., 2019). The hyperspectral data acquired by the classic spaceborne sensors (e.g., Hyperion) suffer from low spatial and temporal resolution and low signal to noise ratio (SNR) (Gomez et al., 2008). Endeavors have been made to address these shortcomings by the current or forthcoming Earth observation missions such as Environmental Mapping and Analysis Program (EnMAP) (Storch et al., 2023) and Copernicus Hyperspectral Imaging Mission for the Environment (CHIME) (Nieke and Rast, 2018). Successful imaging of SOC has also been reported in several studies that implemented multior hyperspectral sensors mounted on aircrafts (Gholizadeh et al., 2018; Hong et al., 2020; Guo et al., 2021; Majeed et al., 2023) or drones (Wehrhan and Sommer, 2021; Biney et al., 2021; 2023). The lower imaging altitude of these aerial platforms has led to higher spatial resolution data and hence more robust prediction models and higher precision and quality distribution maps. The accuracy of remote sensing - based SOC estimation is usually influenced by factors such as soil texture and moisture (Jiang et al., 2016; Castaldi et al., 2019a), wet and dry vegetation (Dvorakova et al., 2020; Castaldi, 2021) and atmospheric conditions (Dvorakova et al., 2020) which must be taken into account cautiously. The relationship between SOC and soil texture parameters (i.e., clay, silt and sand) is complicated and variable (Stenberg et al., 2010), mainly depending on some factors like soil type, aggregates, land use and management practices. However, there is generally a strong correlation between SOC and finer texture content of soil, such as clay. That is because fine clay particles provide a large surface area with negative charge, allowing them to adsorb SOC and protect it from microbial deterioration and decomposition (Barr´ e et al., 2014; Han et al., 2016; Adiyah et al., 2022). Silt particles also have large surfaces which facilitate the aggregation and retention of SOC, although the SOC and silt have a more complex relationship since silt can be easily lost by wind and water erosion. The relationship between SOC and sand depends on various factors including climate as well as soil type, zone, structure and management practices (Yost and Hartemink, 2019). Generally, soils with coarser textures, such as sand, store lower content of SOC (Yost and Hartemink, 2019). That is mainly due to lower surface area (Almajmaie et al., 2017), water-holding capacity (ˇ Simanský et al., 2019) and nutrient content of sandy soils (Almajmaie et al., 2017; ˇ Simanský et al., 2019), which limit their plant growth capacity and hence organic matter inputs. In addition, sandy soils have more coarse pores which cause better aeration and hence, faster decomposition of their SOC contents. To overcome the above-mentioned complexity between soil texture parameters and SOC, especially in large datasets, soil samples can be stratified (i.e., divided) into several classes based on a certain feature. Models built on stratified samples have shown improved accuracies of SOC estimation, which is because of lower variability and similar characteristics of each sample group (Wetterlind and Stenberg, 2010; Castaldi et al., 2018). In a study in the United States, texture-based stratification of nearly 20000 nationwide soil samples into more homogeneous groups led to a significant improvement in prediction of SOC using visible–near infrared−shortwave infrared (Vis–NIR–SWIR; 350–2500 nm) reflectance spectra (Wijewardane et al., 2016). On a nation-wide scale study in Germany, some soil properties like texture were used for stratification of samples, which caused an increase in accuracy of SOC model calibrated on the NIR spectra. Separation of the sandy soil samples from the other samples caused a 14% reduction in root mean square error (RMSE) (Jaconi et al., 2017). In another study by Moura-Bueno et al. (2020), the stratification of a soil spectral library (SSL) resulted in higher accuracy of SOC predictions than when using general models. Stratification can also be applied based on regional (Igne et al., 2010) and geological (Udelhoven et al., 2003) characteristics of the dataset. Region-based division of samples led to superior local models with enhanced SOC prediction performances (Wijewardane et al., 2016). Despite numerous literatures documenting the relationship of texture parameters with SOC (Bosatta and Ågren, 1997; Dexter, 2004; Wiesmeier et al., 2019; Wang et al., 2021; Johannes et al., 2023), only a few studies have delved into their impact on SOC prediction (Schillaci et al., 2017b; Hamzehpour et al., 2019; Swetha and Chakraborty, 2021; Garosi et al., 2022). Above-mentioned lacuna is more acute, when multiand hyperspectral imagery are being used as the main SOC predictor variables. To the best of our knowledge, no study has investigated the effect of texture-based stratification on SOC prediction models, developed solely on the image spectra. The main purpose of this study, hence, was to explore the influence of soil texture on performance of the SOC prediction models developed using multiand hyperspectral remote sensing data. To do so, we stratified the soil samples into different classes based on i) geographical location and ii) soil texture strategies and then examined how well the S2 multispectral and CASI/SASI (CS) airborne hyperspectral imagery, can predict SOC contents within each defined class, either with or without involving the texture parameters as predictor variables. 2. Materials and methods 2.1. Study sites, soil sampling and analysis This study was conducted in four agricultural sites of Nov´ a Ves nad Popelkou (Nov´ a Ves), Jiˇ cín, Kluˇ cov and Pˇ restavlky (Fig. 1), located in rural areas of the Czech Republic. Oilseed rape, spring and winter cereals, potatoes and maize are mostly being cultivated in them. An attempt was made to select homogeneous soilscape representative sites with the same climatic conditions (precipitation and temperature), land management and terrain characteristics (ˇ Zíˇ zala et al., 2017). Dissected relief with tributary valleys, toe and back slopes and plateaus are some of the main landform features of our study sites. Soils of the study sites, as described by the World reference base (WRB) for soil resources (IUSS Working Group WRB, 2014), are mostly composed of Cambisols and Stagnosols on crystalline and sedimentary rocks and Chernozems and Luvisols on loess. Details about the sampling sites and data collection campaigns are presented in Table 1. Overall, 320 soil samples were collected from the topsoil (0–20 cm) of four study sites in June 2021 (80 samples from each). Stratified random strategy conditioned Latin Hypercube Sampling (cLHS) was used for the sampling design. A GeoXM global positioning system (GPS) (Trimble Inc., Sunnyvale, California, USA) with 1 m accuracy was employed to record the position of each sampling point. Soil samples were taken as composite samples, then air-dried, ground, sieved (≤ 2 mm) and thoroughly mixed before analyzing (ISO 11464:2006). SOC was measured using the Walkley–Black method. Soil texture (i.e., clay, silt and sand) content of the samples was determined by the analysis of particle size distribution (less than 0.002 mm, 0.002–0.01 mm, 0.01–0.05 mm, 0.05–0.25 mm and 0.25–2 mm fractions) using the pipette method (ISO 11277:2009). 2.2. Airborne hyperspectral imaging, data pre-processing and coregistration The airborne hyperspectral imaging campaign was carried out on June 3rd, 2021, at least five days after the last rain. Data acquisition took place using Vis–NIR (380–1050 nm) CASI (Itres Ltd., Calgary, Canada) and shortwave infrared (SWIR; 950–2450 nm) SASI (Itres Ltd., Calgary, Canada) sensors (Table 2), aboard Cessna 208B Grand Caravan photogrammetric aircraft. A ground survey was also carried out, on the day of the campaign, to obtain information on soil cover of vegetation and plant residues, moisture and surface roughness. Data pre-processing (radiometric and geometric corrections) were performed by Global Change Research Institute of the Czech Academy of Sciences, via flying laboratory of imaging systems (FLIS). Radiometric corrections were done in the RadCorr Ver. 9.3.6.0 (Itres Ltd., Calgary, Canada) using laboratory determined calibration parameters. The basic procedure of radiometric correction consisted of the dark subtraction V. Khosravi et al. Soil & Tillage Research 241 (2024) 106125 3 and conversion of raw values taken by the sensor (DN: digital numbers) into physically defined radiance units. Atmospheric corrections were implemented using the MODerate resolution atmospheric TRANsmission (MODTRAN) radiative transfer model (RTM) incorporated in the ATCOR-4 Ver. 7.3 (ReSe Applications Schlapfler Inc., Wil, Switzerland). The resulting atmospherically corrected data were expressed as reflectance at the surface level values. Geometric corrections, orthorectification and georeferencing (UTM zone 33 N, ETRS-89) of data were carried out using data collected by the GNSS/IMU sensor and the digital elevation model (DEM) in the GeoCor Ver. 3.7.2 (Itres Ltd., Calgary, Canada). The noisy bands at the two edges of the sensors’ spectra were eliminated, resulted in 64 (from 383.05 to 981.98 nm) and 98 (from 987.5 to 2442.5 nm) remaining bands for CASI and SASI sensors, respectively. Before the main processing step, the CASI data were resampled to the spatial resolution of the SASI sensor (6 m). Data fusion was then performed to merge the CASI and SASI imagery into one hyperspectral full Vis–NIR–SWIR spectral range data cube. Afterwards, reflectance (R) was transformed to absorbance via log (1/R) and the first derivative (FD) transformation was applied on the spectra to remove the baseline offset (Gholizadeh et al., 2015; Khosravi et al., 2020). 2.3. Multispectral satellite data pre-processing and indices retrieval The two cloud-free bottom of atmosphere (BOA) Sentinel-2 level-2A (S2) images were downloaded from the Copernicus open access hub as close as possible to the sampling campaign date (Table 1). Since the downloaded S2 products were geometrically and atmospherically corrected, there was no need for pre-processing of the downloaded images. The six 20 m-pixel bands were resampled to 10 m to obtain a finer pixel size S2 image and hence more precise selection of the bare soil pixels and draw better comparison with the 6 m resolution CS data. S2 related Fig. 1. Distribution of study sites in the a) Czech Republic and sampling fields in b) Kluˇ cov, c) Pˇ restavlky, d) Nov´ a Ves nad Popelkou and e) Jiˇ cín. Table 1 Study sites, soil sampling and imaging campaign specifications. Site Area (ha) Dominant soil unit Parent material Soil sampling date Airborne data acquisition date Sentinel-2 data acquisition date Nov´ a Ves 129 Eutric Cambisol Permian-Carboniferous rocks (sandstone, siltstone) June, 2021 03.06.2021 04.06.2021 Jiˇ cín 153 Luvisol, Albic Luvisol, Luvic Chernozem Pleistocene loess June, 2021 03.06.2021 04.06.2021 Kluˇ cov 175 Calcic Chernozem, Luvic Chernozem, Regosols Pleistocene loess, Cretaceous marl, terrace gravels June, 2021 03.06.2021 19.06.2021 Pˇ restavlky 58 Haplic Stagnosol, Eutric and Stagnic Cambisol, Leptosol Complex of Proterozoic and Paleozoic rocks (schist, granodiorite) June, 2021 03.06.2021 19.06.2021 Table 2 CASI and SASI sensors specifications. Sensor Number of bands Spectral range (nm) Spectral resolution (nm) Field of view (FOV) FWHM (nm) CASI 72 Vis–NIR, 380–1050 3.2 40◦10 SASI 100 SWIR, 950–2450 15 40◦15 V. Khosravi et al. Soil & Tillage Research 241 (2024) 106125 4 variables including ten bands and 19 spectral indices were extracted and calculated at each sampling point and employed as covariates to develop the SOC prediction models. A Pearson correlation analysis was conducted to explore the possible autocorrelations between the input indices and consequent adverse effects such as increasing the calibration complexity and reduction in accuracy of the SOC estimations. According to the results, there were significant correlations between some indices such as NDVI and TVI, WDVI and V, BI2 and TSAVI, and EVI and SAVI. Some combinations of the indices were used as inputs to investigate their effects on the prediction. Considering the negligible differences observed between the results, all indices were incorporated in calibration of the SOC models to employ the maximum number of S2 related covariates. The S2 bands’ details and equations of the indices can be found in Table A1 and Table A2, respectively. 2.4. Selection of bare soil pixels The normalized difference vegetation index (NDVI) and normalized burn ratio 2 (NBR2) (Table A2) indices have been considerably used to delineate green vegetation (Huang et al., 2021) and crop residues (Dvorakova et al., 2023), respectively. Although the sensitivity of NBR2 to detect crop residues on wet soils is very limited, they form a linear relationship for dry soils (Dvorakova et al., 2021), like the case in this study. Accordingly, both NDVI and NBR2 indices were implemented for both S2 and CS to remove the pixels with green and dry vegetation and crop residue cover. The NDVI threshold was determined by visual inspection of the true color composite (TCC) of S2 (R: 664.6 nm, G: 559.8 nm and B: 492.4 nm) and airborne (R: 658.8 nm, G: 554.2 nm and B: 478.1 nm) images. Various NBR2 threshold values ranging from 0.05 to 0.25, with increment of 0.01 were selected and tested to isolate the bare soil areas. Results were validated by comparing the retrieved pixels with field-based ground truth observations. 2.5. Sample stratification and correlation analysis The samples corresponding to the remaining bare soil areas were stratified based on the study sites and the United States Department of Agriculture (USDA) texture triangle. The texture triangle of the USDA is one of the most popular systems of soil classification, evolved from eight classes right-angled triangle (Whitney, 1911) and ten classes equilateral triangle (Davis, 1927) models to a twelve groups final version introduced in 1951, which is currently in effect (USDA, 1951, 2017). The Pearson correlation between SOC and each texture parameter was calculated within each stratified class to provide an insight into their relationship, while investigating the possible impacts of interaction between the two variables on SOC prediction. 2.6. SOC modeling and performance assessment For each stratified class, random forest (RF) algorithm (Breiman, 2001) was used to develop the prediction models using S2 and CS imagery, with and without separate and combined incorporation of the texture parameters as additional predictor variables. Results were then compared with those obtained by the RF models calibrated on non-stratified (NS) samples (i.e., on the whole dataset). Three specific parameters used for the RF models in this study were: i) the number of trees in the forest (n tree ), ii) quantity of variables used for each tree (m try ), and iii) the minimum data per node (nodesize). A grid search was conducted to determine the optimum quantities for the parameters and the value of 500 was obtained for n tree . One third of the total number of predictor variables was the optimal quantity of the m try while the nodesize was set to 5. The ‘randomForest’ package in RStudio environment (R Core Team, 2020) was used to run the RF model. Using variable importance in projection (VIP), the contribution of all input predictor variables, including the texture parameters, were assessed on outcome of each developed class-model. This can help to determine the extent to which the texture can influence the prediction accuracy. To develop the models, samples were divided into training and testing (75:25% ratio) datasets, using the Kennard–Stone (KS) algorithm, which chooses representative samples based on a distance measure (Kennard and Stone, 1969). A 5-fold cross-validation technique was used to calibrate the SOC prediction models on the training dataset. The two most common evaluation criteria, RMSE and coefficient of determination (R 2 ) (Eq. 1 and Eq. 2, respectively) were selected to assess performance of the developed models. RMSE = 1 n∑ n i=1 (oi−pi)2 √(1) R2=1−∑n i=1(oi−pi)2 ∑n i=1(oi− μ o)2(2) where, n is the number of samples, µ o is the mean for the observed values and o i and p i are the observed and predicted values, respectively. A robust model, generally, has high R 2 and low RMSE. A flowchart of the methodology is provided in Fig. 2. 3. Results 3.1. General statistics and correlations After rigorous perusal, pixels with NDVI values below 0.23 were kept for both satellite and airborne images. Different models were then developed on the samples remaining after applying various NBR2 thresholds (from 0.05 to 0.25 with increment of 0.01) and the optimal value of 0.09 was obtained. By applying the optimal NDVI and NBR2 thresholds, 200 bare soil pixels/samples which were fully compatible between the two datasets were finally retrieved and used for the rest of the study. Fig. 3 and Fig. 4 show the statistics, histograms and the remaining soil samples (based on the USDA texture triangles), for the site-based and texture-based stratified samples, respectively. The statistics of NS samples have also been shown in Fig. 3. According to the results, SOC contents of the NS samples ranged from the minimum value of 0.31% to the maximum value of 2.59% with average and coefficient of variation (CV) of 1.19% and 0.27, respectively. Among the study sites, SOC in Kluˇ cov had the highest CV of 0.29 followed by Pˇ restavlky (CV =0.24) and Nov´ a Ves (CV =0.22), respectively. The highest CV of clay and silt were observed in Pˇ restavlky, while Kluˇ cov showed the highest CV of sand. That was while Jiˇ cín indicated the lowest CV values for all three texture parameters. According to the histogram plots and related distribution curves, SOC contents of all classes almost followed the normal distribution. Other variables did not follow a certain distribution in different classes (Fig. 3). According to the USDA texture triangle (Fig. 4), samples were mostly plotted in lower positioned classes of sandy loam (SL), loam (L), silt loam (SiL), silty clay loam (SiCL) and clay loam (CL). For this reason, samples were stratified into four groups of SL with 46 samples, L with 59 samples, SiL with 51 samples and SiCL_CL as combination of SiCL and CL with 44 samples (Fig. 4). Accordingly, the SiCL_CL and SiL classes had higher amounts of clay and silt, class L had relatively equal quantities of sand and silt and the SL class included samples with high sand contents. Considering the SOC level in each texture class, the SiCL_CL had the highest average of SOC (1.23%) with highest CV of 0.31. SiL and SL classes had the lowest SOC content with averages of 1.15% and 1.17%, respectively. According to the distributions of soil textures within the USDA texture triangle, most of Nov´ a Ves, Jiˇ cín, Kluˇ cov and Pˇ restavlky samples were in loam, silt loam, silty clay loam and loam classes, respectively (Fig. 4). The Pearson correlation coefficients between SOC and clay, silt and sand in the NS and stratified classes are presented in Table 3. None of the V. Khosravi et al. Soil & Tillage Research 241 (2024) 106125 5 texture parameters showed significant correlation with SOC in the NS class. Regarding the site-based stratified classes, SOC showed significant correlation with all texture parameters in Nov´ a Ves, Kluˇ cov and Pˇ restavlky sites. In addition, the interrelationship was positive between SOC and clay as well as SOC and silt and negative between SOC and sand within those study sites. The highest positive correlation of SOC was observed with clay in Kluˇ cov (r =0.84, P <0.001), while the highest negative correlation was detected between SOC and sand for Pˇ restavlky (r = − 0.39, P <0.001). No significant correlations were observed between SOC and texture parameters in Jiˇ cín. Considering the texture-based stratified classes, SOC was positively correlated with the clay content, while the correlation between SOC and sand was negative in all classes excluding L. Though not quite strong, the highest negative correlation between SOC and sand content was in SL class (r = − 0.32, p-value <0.05), while SOC and clay showed the highest positive correlation in SiCL_CL (r =0.75, p-value <0.05). 3.2. Modeling results 3.2.1. General comparison Tables A3 and 4 display the SOC prediction results for the models created on the NS, siteand texture-based stratified samples (calibration and validation sets, respectively). Fig. 5 and Fig. A1 also show the scatterplots of the best models obtained for each strategy and stratified class on the S2 +clay and CS +clay datasets, respectively. Prediction results of models developed on NS samples are likewise shown. As can be Fig. 2. Flowchart of the study. V. Khosravi et al. Soil & Tillage Research 241 (2024) 106125 6 Fig. 3. Histograms and summary statistics of SOC (first column), clay (second column), silt (third column) and sand (fourth column) within all stratified and nonstratified classes (NS: non-stratified, SL: Sandy Loam, L: Loam, SiL: Silt Loam, SiCL_CL: combination of SiCL (Silty Clay Loam) and CL (Clay Loam)). V. Khosravi et al. Soil & Tillage Research 241 (2024) 106125 7 observed, the models developed on the CS outperformed those developed on the S2 data. This superiority was more evident (up to 60%) in the models developed on NS samples (Table 4, Fig. 5 and Fig. A1). Even though fewer samples were used for modeling, models developed on the stratified samples indicated better performance than those developed on the NS samples. As an example, the site-based stratification caused a minimum improvement of 7%, when using the CS data as the main predictor variable (+clay, silt and sand). This improvement was more pronounced (at least 29%), when the S2 image was considered as the main predictor variable (+clay, silt and sand). The lowest model performances among all (RMSE >0.17% and R 2 <0.41), linked to those developed on the NS samples, regardless of the employed data type and predictor variables. Comparing the results obtained under the siteand texture-based stratification strategies, the models developed on SiCL_CL and SiL outperformed all the site-based stratified ones, showing that dividing samples based on their texture parameters can lead to more successful SOC predictions. This conclusion holds true, except for some models developed on the classes with higher sand content, especially SL, which performed very poorly (RMSE >0.16% and R2 <0.42). In other words, the models developed for Kluˇ cov and Pˇ restavlky, under the site-based stratification, had higher performances than the L and SL models, developed under the texture-based stratification. The performance of some models developed for Jiˇ cín was also similar to those developed on the L class (like the model developed on the CS data +clay with RMSE = 0.13% and R 2 =0.60). Among all the texture parameters, adding clay as the predictor variable caused the highest enhancement in performance of the models. The best improvement was obtained for Jiˇ cín, where incorporation of clay enhanced the prediction performance up to 30% compared to the Fig. 3. (continued). V. Khosravi et al. Soil & Tillage Research 241 (2024) 106125 8 situations that no texture parameter added as the input variable. Adding all texture parameters together also made substantial improvements compared to when the multior hyperspectral data were used as the only predictors. Adding silt and sand individually, on the other hand, did not introduce any improvement to the models’ prediction performances. It is noteworthy to mention that adding the texture parameters caused more improvement for the models developed on the site-based stratified samples than those developed on the texture-based stratified ones. Furthermore, adding texture parameters caused an average higher performance of 4.73% for the models developed using the S2 data, while the average enhancement was 3.78% considering the CS data as the main predictor. 3.2.2. Site-based stratification Considering Table 4, Pˇ restavlky and Kluˇ cov yielded acceptable SOC prediction performances on both S2 and airborne data. However, the best results under the site-based stratification strategy were obtained on the airborne hyperspectral data +clay with RMSE =0.11% and R 2 = 0.71, followed by the airborne hyperspectral data +clay +silt +sand with RMSE =0.12% and R 2 =0.70, both for the Kluˇ cov site. Likewise, for the Pˇ restavlky site, the models developed on the airborne hyperspectral data +clay (RMSE =0.12% and R 2 =0.65) and the airborne hyperspectral data +clay +silt +sand (RMSE =0.12% and R 2 =0.63) showed the highest performances. The same trend was observed for the models developed using S2 (with and without the texture parameters) dataset. Accordingly, Kluˇ cov Fig. 4. USDA texture triangles for the non-stratified, siteand texture-based stratified soil classes (NS: non-stratified, SL: Sandy Loam, L: Loam, SiL: Silt Loam, SiCL_CL: combination of SiCL (Silty Clay Loam) and CL (Clay Loam)). V. Khosravi et al. Soil & Tillage Research 241 (2024) 106125 9 and Pˇ restavlky had relatively acceptable predictions, especially when texture parameters were incorporated as the predictor variables (S2 + clay with RMSE =0.13% and R 2 =0.64, S2 +clay +silt +sand with RMSE =0.13% and R 2 =0.59, for Kluˇ cov, and S2 +clay with RMSE = 0.15% and R 2 =0.61 and S2 +clay +silt +sand with RMSE =0.14% and R 2 =0.58, for Pˇ restavlky). For other sites (Nov´ a Ves and Jiˇ cín), the performance of site-based stratification models dropped, regardless of using S2 or airborne datasets. For example, the best model for Jiˇ cín and Nov´ a Ves was developed using airborne hyperspectral data +clay with RMSE =0.13% and R 2 =0.60 and RMSE =0.16% and R 2 =0.46, respectively. 3.2.3. Texture-based stratification Among all models developed on the texture-based stratified samples, those developed on the SiCL_CL samples showed the highest prediction performances, followed by the SiL models (Table 4). This was regardless of the type of data and the texture parameters added as input predictor variables. The best prediction accuracy was obtained for the SiCL_CL class model developed on the airborne hyperspectral data +clay (RMSE =0.08% and R 2 =0.81), followed by the same class model developed on the airborne hyperspectral data +clay +silt +sand (RMSE =0.10% and R 2 =0.78). According to Table 4, it is highlighted that the SiCL_CL models also yielded promising results including RMSE =0.10% and R 2 =0.78 for the airborne hyperspectral data +clay and RMSE =0.10% and R 2 =0.73 for the airborne hyperspectral data +clay +silt +sand. Contrarily, performance of models developed on the L (RMSE =0.12% and R 2 =0.62) and SL (RMSE =0.16% and R 2 =0.42) classes was not reasonable. 3.2.4. Variable importance Variable importance plots of the RF models developed on different classes of samples are shown in Fig. 6. Accordingly, the 30 highest influential predictor variables of each model are presented along with their percentage increase in mean squared error (%IncMSE) as a measure of variable importance obtained during the calibration of the models. As evidenced, clay was amongst the top important variables in almost all models (excluding Nov´ a Ves CS). It also made the largest contribution to the accuracy of the model developed on SiCL_CL class with %IncMSE of 3.48 for the CS dataset. Compared to sand, silt better contributed to development of the models, especially when combined with the CS data. The most significant impact of sand was on the models developed for SL and Pˇ restavlky classes (Fig. 6). 4. Discussion 4.1. SOC prediction: effect of correlation with the texture parameters It has generally been proved that SOC content is positively and negatively related to clay and sand contents of soil, respectively (Song Table 3 Correlation matrix between SOC content and soil texture parameters in the nonstratified and stratified soil classes. Correlation coefficient Texture parameter NS Nov´ a Ves Jiˇ cín Kluˇ cov Pˇ restavlky Clay 0.05 0.46* 0.19 0.84** 0.48** Silt -0.10 041** -0.08 0.73* 0.60*** Sand 0.05 -0.22 0.10 -0.23** -0.39*** Texture parameter SL L SiL SiCL_CL Clay 0.60* 0.22* 0.30** 0.75* Silt 0.50*** -0.10 0.36** -0.40** Sand -0.32* 0.05 -0.04 -0.01 NS: non-stratified, SL: sandy loam, L: loam, SiL: silt loam, SiCL_CL: combination of silty clay loam (SiCL) and clay loam (CL). *, ** and *** indicate significancy at P <0.05, P <0.01 and P <0.001, respectively. Table 4 SOC prediction results of the RF models developed on the non-stratified and stratified classes of samples (validation dataset). Data type Texture input variable NS (200 samples) Nov´ a Ves (60 samples) Jiˇ cín (18 samples) Kluˇ cov (59 samples) Pˇ restavlky (63 samples) SL (46 samples) L (59 samples) SiL (51 samples) SiCL_CL (44 samples) RMSE R 2 RMSE R 2 RMSE R 2 RMSE R 2 RMSE R 2 RMSE R 2 RMSE R 2 RMSE R 2 RMSE R 2 Sentinel-2 None 0.27 0.24 0.28 0.32 0.19 0.43 0.16 0.55 0.17 0.55 0.28 0.24 0.17 0.48 0.13 0.56 0.13 0.59 All 0.26 0.26 0.27 0.34 0.14 0.55 0.13 0.59 0.14 0.58 0.25 0.25 0.13 0.55 0.11 0.59 0.11 0.63 Clay 0.26 0.27 0.27 0.35 0.15 0.56 0.13 0.64 0.15 0.61 0.26 0.24 0.12 0.57 0.11 0.61 0.10 0.66 Silt 0.27 0.24 0.27 0.32 0.18 0.42 0.16 0.55 0.17 0.54 0.28 0.23 0.16 0.48 0.13 0.57 0.13 0.60 Sand 0.27 0.23 0.28 0.31 0.21 0.41 0.15 0.57 0.18 0.55 0.27 0.23 0.18 0.44 0.15 0.56 0.14 0.55 CS None 0.19 0.38 0.18 0.43 0.16 0.53 0.14 0.68 0.14 0.60 0.19 0.36 0.16 0.55 0.12 0.69 0.12 0.76 All 0.17 0.41 0.16 0.44 0.13 0.59 0.12 0.70 0.12 0.63 0.17 0.37 0.12 0.62 0.10 0.73 0.10 0.78 Clay 0.18 0.41 0.16 0.46 0.13 0.60 0.11 0.71 0.12 0.65 0.16 0.42 0.13 0.60 0.10 0.78 0.08 0.81 Silt 0.19 0.38 0.20 0.43 0.15 0.53 0.13 0.67 0.15 0.62 0.18 0.36 0.14 0.58 0.13 0.71 0.11 0.76 Sand 0.23 0.37 0.21 0.41 0.15 0.52 0.15 0.63 0.12 0.67 0.18 0.34 0.15 0.55 0.13 0.69 0.12 0.75 NS: non-stratified, SL: sandy loam, L: loam, SiL: silt loam, SiCL_CL: combination of silty clay loam (SiCL) and clay loam (CL). V. Khosravi et al. Soil & Tillage Research 241 (2024) 106125 16 data for the estimation and mapping of topsoil organic carbon: Feature selection coupled with random forest. Soil Tillage Res. 199, 104589. Hong, H., Chen, S., Fang, Q., Algeo, T.J., Zhao, L., 2019. Adsorption of organic matter on clay minerals in the Dajiuhu peat soil chronosequence, South China. Appl. Clay Sci. 178, 105125. Hontoria, C., Saa, A., Rodríguez-Murillo, J.C., 1999. Relationships between soil organic carbon and site characteristics in peninsular Spain. Soil Sci. Soc. Am. J. 63, 614–621. Huang, S., Tang, L., Hupy, J.P., Wang, Y., Shao, G., 2021. A commentary review on the use of normalized difference vegetation index (NDVI) in the era of popular remote sensing. J. Res. 32, 1–6. Igne, B., Reeves, J.B., McCarty, G., Hively, W.D., Lund, E., Hurburgh, C.R., 2010. Evaluation of spectral pretreatments, partial least squares, least squares support vector machines and locally weighted regression for quantitative spectroscopic analysis of soils. J. Infrared Spectrosc. 18, 167–176. IUSS Working Group WRB, 2014. World reference base for soil resources. In: International Soil Classification System for Naming Soils and Creating Legends for Soil Maps; World Soil Resources Reports No. 106. FAO, Rome, Italy Available online. http:// www.fao.org/3/a-i3794e.pdf, Accessed date: 1 October 2015. Jaconi, A., Don, A., Freibauer, A., 2017. Prediction of soil organic carbon at the country scale: stratification strategies for near-infrared data. Eur. J. Soil Sci. 68, 919–929. Jiang, Q., Chen, Y., Guo, L., Fei, T., Qi, K., 2016. Estimating soil organic carbon of cropland soil at different levels of soil moisture using VIS-NIR spectroscopy. Remote Sens. 8, 755. Jindaluang, W., Kheoruenromne, I., Suddhiprakarn, A., Singh, B.P., Singh, B., 2013. Influence of soil texture and mineralogy on organic matter content and composition in physically separated fractions soils of Thailand. Geoderma 195, 207–219. Johannes, A., Sauzet, O., Matter, A., Boivin, P., 2023. Soil organic carbon content and soil structure quality of clayey cropland soils: A large-scale study in the Swiss Jura region. Soil Use 39, 707–716. Juˇ ricov´ a, A., Chuman, T., ˇ Zíˇ zala, D., 2022. Soil organic carbon content and stock change after half a century of intensive cultivation in a chernozem area. Catena 211, 105950. Kennard, R.W., Stone, L.A., 1969. Computer aided design of experiments. Technometrics 11, 137–148. Kennedy, M.J., Pevear, D.R., Hill, R.J., 2002. Mineral surface control of organic carbon in black shale. Science 295, 657–660. Khosravi, V., Ardejani, F.D., Aryafar, A., Yousefi, S., Karami, S., 2020. Prediction of copper content in waste dump of Sarcheshmeh copper mine using visible and nearinfrared reflectance spectroscopy. Environ. Earth Sci. 79, 165. Konen, M.E., Burras, C.L., Sandor, J.A., 2003. Organic Carbon, Texture, and Quantitative Color Measurement Relationships for Cultivated Soils in North Central Iowa. Soil Sci. Soc. Am. J. 67, 1823–1830. Leifeld, J., Bassin, S., Fuhrer, J., 2005. Carbon stocks in Swiss agricultural soils predicted by land-use, soil characteristics, and altitude. Agric. Ecosyst. Environ. 105, 255–266. Liu, F., Rossiter, D.G., Song, X.D., Zhang, G.L., Yang, R.M., Zhao, Y.G., Li, D.C., Ju, B., 2014. A similarity-based method for three-dimensional prediction of soil organic matter concentration. Geoderma 263, 254–263. Majeed, I., Purushothaman, N.K., Chakraborty, P., Panigrahi, N., Vasava, H.B., Das, B.S., 2023. Estimation of soil and crop residue parameters using AVIRIS-NG hyperspectral data. Int. J. Remote Sens. 44 (6), 2005–2038. Maurya, S., Abraham, J.S., Somasundaram, S., Toteja, R., Gupta, R., Makhija, S., 2020. Indicators for assessment of soil quality: a mini-review. Environ. Monit. Assess. 192, 1–22. McLauchlan, K.K., 2006. Effects of soil texture on soil carbon and nitrogen dynamics after cessation of agriculture. Geoderma 136, 289–299. Moura-Bueno, J.M., Dalmolin, R.S.D., Horst-Heinen, T.Z., ten Caten, A., Vasques, G.M., Dotto, A.C., Grunwald, S., 2020. When does stratification of a subtropical soil spectral library improve predictions of soil organic carbon content? Sci. Total Environ. 737, 139895. Nieke, J. and Rast, M., 2018, July. Towards the copernicus hyperspectral imaging mission for the environment (CHIME). In Igarss 2018-2018 ieee international geoscience and remote sensing symposium (pp. 157-159). IEEE. Pan, G., Li, L., Zhang, X., Dai Jingyu, Z.Y., Zhang, P., 2003. Soil organic carbon storage of China and the sequestration dynamics in agricultural lands. Adv. Earth Sci. 18, 609. Schillaci, C., Acutis, M., Lombardo, L., Lipani, A., Fantappi` e, M., M¨ arker, M., Saia, S., 2017b. Spatio-temporal topsoil organic carbon mapping of a semi-arid Mediterranean region: The role of land use, soil texture, topographic indices and the influence of remote sensing data to modelling. Sci. Total Environ. 601–602. Schillaci, C., Lombardo, L., Saia, S., Fantappi` e, M., M¨ arker, M., Acutis, M., 2017a. Modelling the topsoil carbon stock of agricultural lands with the Stochastic Gradient Treeboost in a semi-arid Mediterranean region. Geoderma 286, 35–45. ˇ Simanský, V., Juriga, M., Jonczak, J., Uzarowicz, Ł., Stępie´ n, W., 2019. How relationships between soil organic matter parameters and soil structure characteristics are affected by the long-term fertilization of a sandy soil. Geoderma 342, 75–84. Sollins, P., Homann, P., Caldwell, B.A., 1996. Stabilization and destabilization of soil organic matter: Mechanisms and controls. Geoderma 74, 65–105. Song, J., Gao, J., Zhang, Y., Li, F., Man, W., Liu, M., Wang, J., Li, M., Zheng, H., Yang, X., Li, C., 2022. Estimation of Soil Organic Carbon Content in Coastal Wetlands with Measured VIS-NIR Spectroscopy Using Optimized Support Vector Machines and Random Forests. Remote Sens (Basel) 14, 4372. Stenberg, B., Viscarra Rossel, R.A., Mouazen, A.M., Wetterlind, J., 2010. Visible and Near Infrared Spectroscopy in Soil Science. pp. 163–215. Storch, T., Honold, H.P., Chabrillat, S., Habermeyer, M., Tucker, P., Brell, M., Ohndorf, A., Wirth, K., Betz, M., Kuchler, M., Mühle, H., 2023. The EnMAP imaging spectroscopy mission towards operations. Remote Sens Environ. 294, 113632. Swetha, R.K., Chakraborty, S., 2021. Combination of soil texture with Nix color sensor can improve soil organic carbon prediction. Geoderma 382, 114775. Thuri` es, L., Pansu, M., Feller, C., Herrmann, P., R´ emy, J.C., 2001. Kinetics of added organic matter decomposition in a Mediterranean sandy soil. Soil Biol. Biochem 33, 997–1010. Udelhoven, T., Emmerling, C., Jarmer, T., 2003. Quantitative analysis of soil chemical properties with diffuse reflectance spectrometry and partial least-square regression: A feasibility study. Plant Soil 251, 319–329. USDA, 1951. Soil survey manual. Soil Conservation Service Volume Handbook, U.S. Department of Agriculture 18. USDA, 2017. Soil survey manual. Soil Conservation Service Volume Handbook, U.S. Department of Agriculture 18. Vaudour, E., Gomez, C., Loiseau, T., Baghdadi, N., Loubet, B., Arrouays, D., Ali, L., Lagacherie, P., 2019. The impact of acquisition date on the prediction performance of topsoil organic carbon from Sentinel-2 for croplands. Remote Sens 11 (18), 2143. Wang, Q., Wen, Y., Zhao, B., Hong, H., Liao, R., Li, J., Liu, J., Lu, H., Yan, C., 2021. Coastal soil texture controls soil organic carbon distribution and storage of mangroves in China. Catena 207, 105709. Wang, S.Q., Yu, G.R., Zhao, Q.J., Niu, D., Chen, Q.M., Wu, Z.F., 2005. Spatial characteristics of soil organic carbon storage in China’s croplands. Pedosphere 15, 417–423. Wehrhan, M., Sommer, M., 2021. A parsimonious approach to estimate soil organic carbon applying unmanned aerial system (UAS) multispectral imagery and the topographic position index in a heterogeneous soil landscape. Remote Sens. 13 (18), 3557. Wetterlind, J., Stenberg, B., 2010. Near-infrared spectroscopy for within-field soil characterization: small local calibrations compared with national libraries spiked with local samples. Eur. J. Soil Sci. 61, 823–843. Whitney, M., 1911. Use of soils east of the Great Plains region. Wiesmeier, M., Urbanski, L., Hobley, E., Lang, B., von Lützow, M., Marin-Spiotta, E., van Wesemael, B., Rabot, E., Ließ, M., Garcia-Franco, N., Wollschl¨ ager, U., Vogel, H.J., K¨ ogel-Knabner, I., 2019. Soil organic carbon storage as a key function of soils - A review of drivers and indicators at various scales. Geoderma 333, 149–162. Wijewardane, N.K., Ge, Y., Wills, S., Loecke, T., 2016. Prediction of soil carbon in the conterminous United States: visible and near infrared reflectance spectroscopy analysis of the rapid carbon assessment project. Soil Sci. Soc. Am. J. 80, 973–982. ˇ Zíˇ zala, D., Minaˇ rík, R., Z´ adorov´ a, T., 2019. Soil organic carbon mapping using multispectral remote sensing data: Prediction ability of data with different spatial and spectral resolutions. Remote Sens. 11 (24), 2947. Yost, J.L., Hartemink, A.E., 2019. Soil organic carbon in sandy soils: A review. pp. 217–310. V. Khosravi et al.