Predicted Maps for Soil Organic Matter Evaluation: The Case of Abruzzo Region (Italy)

: Organic matter, an important component of healthy soils, may be used as an indicator in sustainability assessments. Managing soil carbon storage can foster agricultural productivity and environmental quality, reducing the severity and costs of natural phenomena. Thus, accurately estimating the spatial variability of soil organic matter (SOM) is crucial for sustainable soil management when planning agro-environmental measures at the regional level. SOM variability is very large in Italy, and soil organic carbon (SOC) surveys considering such variability are di ﬃ cult and onerous. The study concerns the Abruzzo Region (about 10,800 km 2 ), in Central Italy, where data from 1753 soil proﬁles were available, together with a Digital Elevation Model (DEM) and Landsat images. Some morphometric parameters and spectral indices with a signiﬁcant degree of correlation with measured data were used as predictors for regression-kriging (RK) application. Estimated map of SOC stocks, and of SOM related to USDA (United States Department of Agriculture) texture—an additional indicator of soil quality—were produced with a satisfactory level of accuracy. Results showed that SOC stocks and SOM concentrations in relation to texture were lower in the hilly area along the shoreline, pointing out the need to improve soil management to guarantee agricultural land sustainability.


Introduction
One of the main challenges for the future is to maintain soil functions, but despite many efforts to promote more sustainable land management, soil degradation in the European Union (EU) is increasing [1], with a severe impact on food production and the supply of ecosystem services. Among the soil properties impacting soil quality, soil organic carbon (SOC)-and soil organic matter (SOM) derived from its determination-deserves special attention, representing a key indicator for evaluating soil quality [2], but also impacting the chemical and physical properties of the soil and its overall health. Properties affected by organic matter include soil structure, water holding capacity, diversity and activity of soil organisms, buffering capacity, and nutrient availability. SOM also regulates the efficiency of soil amendments, fertilizers, pesticides, and herbicides. According to the Food and Agriculture Organization of the United Nations (FAO) [3], one of the characteristics of sustainable soil management (SSM) is a stable or increasing storage of SOM, ideally close to the optimal level for the local environment, for all land uses. Thus, the most efficient way to improve soil quality is to stimulate better SOM management.
SOC is also an essential component of the global carbon cycle [4], representing one of the largest reservoirs of terrestrial carbon that may influence global warming [5,6]. The global carbon budget, and the CO 2 emissions associated with its major components (i.e., atmosphere, ocean, fossil Land 2020, 9,349 3 of 14 covariates from remote sensing. Lamichhane et al. [29] reviewed the applications of different DSM techniques reported in the scientific literature from 2013 until 2019 for mapping SOC concentration and stocks. Spatial interpolation by RK is a method that merges a regression of the primary dependent variable on secondary ancillary variables with simple kriging of the regression residuals [30][31][32][33]. Ancillary information is often derived from the digital elevation model (DEM), which provides topographic information that allows for the calculation of terrain parameters, and from satellite imagery, both easily available at relatively low cost. Since maintaining or increasing SOM is one of the main targets of SSM, our research question was whether such a technique, already tested in limited areas in Italy [15,34], could be suitable for SOM assessment at the regional level. We hypothesized that RK could yield accurate estimates for site-specific analyses without further sampling expenses. The Abruzzo region in central Italy was chosen as the study area due to the complexity of its territory, and for the availability of suitable data for RK application.
In this framework, the study aimed to provide a spatial evaluation-at a regional level-of SOC, soil USDA (United States Department of Agriculture) texture, and SOM levels based on the soil texture from point data, estimating values in non-sampled locations by applying RK. Then, by transferring them into a GIS software, we can produce a reliable estimate and a valid evaluation tool for a SSM.

Study Area
The study area, located in central Italy, consists of about 10,800 km 2 corresponding to the territory of the Abruzzo region. Among the 20 administrative regions of Italy, Abruzzo is one of the most mountainous, located in the central peninsular part of the country. Bordered by the Adriatic Sea in the east and by the Apennines in the west, its territory is very complex and heterogeneous. Within a few kilometers, the environment changes from high mountains (Gran Sasso and Maiella) to the seashore, passing through all the intermediate landscapes: mountain grasslands and woods, hills, plains and river basins. While the mountain ranges lie along a NW-SE direction, the rivers cross them toward the sea along a SE-NW direction. About 40% of the total surface is represented by utilized agricultural area. Along the coastline, the climate is Mediterranean, warm and dry, gradually becoming continental moving inlands. Mean annual air temperature goes from 6 • C in the mountains to 15 • C near the sea. Mean annual rainfall ranges from 600-800 mm in the plains and in the river basins to 1000-1200 mm in the hills, reaching 1600 mm in the mountains. Summer is everywhere the dry season.
Soil regions (SR) are the largest units of soil description, depicting areas with similar soil-forming conditions. These are usually defined as typical associations of dominant soils, occurring in areas limited by a specific climate and/or a characteristic association of parent material.  [36].
In Figure 1, the digital elevation model of the region is reported, together with the location of the sampling points and the boundaries of the SRs. A land use map of the Abruzzo Region is reported as Supplementary Material ( Figure S1).

Data Collection
The study dataset, provided by the Regional Agency for Agricultural Extension Services of Abruzzo Region (ARSSA), consists of 1753 georeferenced soil samples collected by an auger (0-25 and 25-50 cm) in accessible agricultural and forest land. The physical and chemical routine analyses included the measurement of particle size distribution according to the pipette method [37] and of SOC content according to the modified Walkley-Black method [38]. The SOC stock in kg m −2 for each soil profile was calculated as follows: where soc is the organic carbon concentration in % for each soil horizon; bd is the bulk density of the soil horizon in g cm −3 estimated by suitable pedofunctions, different for each type of horizon [39]; th is the thickness of the horizon in cm; gr is the gravel content in %; and n is the number of horizons in the soil profile.
The SOM content was then evaluated from SOC [40] by means of the following formula: From an agronomic standpoint, the pedo-climatic context cannot be neglected in the evaluation of SOM, because in different soil types, the same amount of SOM can differently impact soil functions. Thus, the SOM content of the study area was classified into four different levels (i.e., very low, low, medium, and high) based on the USDA textural classes [41], as detailed in Table 1. A land use map of the Abruzzo Region is reported as Supplementary Material ( Figure S1).

Data Collection
The study dataset, provided by the Regional Agency for Agricultural Extension Services of Abruzzo Region (ARSSA), consists of 1753 georeferenced soil samples collected by an auger (0-25 and 25-50 cm) in accessible agricultural and forest land. The physical and chemical routine analyses included the measurement of particle size distribution according to the pipette method [37] and of SOC content according to the modified Walkley-Black method [38]. The SOC stock in kg m −2 for each soil profile was calculated as follows: where soc is the organic carbon concentration in % for each soil horizon; bd is the bulk density of the soil horizon in g cm −3 estimated by suitable pedofunctions, different for each type of horizon [39]; th is the thickness of the horizon in cm; gr is the gravel content in %; and n is the number of horizons in the soil profile. The SOM content was then evaluated from SOC [40] by means of the following formula: From an agronomic standpoint, the pedo-climatic context cannot be neglected in the evaluation of SOM, because in different soil types, the same amount of SOM can differently impact soil functions. Thus, the SOM content of the study area was classified into four different levels (i.e., very low, low, medium, and high) based on the USDA textural classes [41], as detailed in Table 1. Land use, vegetation, climate, and terrain features are the main factors affecting soil properties at the landscape scale, particularly in hills, hence DEM-derived terrain attributes can be used for the prediction of the spatial distribution of soil features. Ancillary data for the area were derived from: (i) a 30 m degraded version of the 20 m DEM provided by the Land Information Service of the Abruzzo Region; (ii) from Landsat 7 TM imagery (three visible bands and four infrared bands), and (iii) the 1:250,000 Soil Subsystems Map of Abruzzo available from ARSSA [42].
From the DEM, the following morphometric attributes were derived: • Elevation (ELEV); • slope gradient (SLOPE); • curvature plan and profile (PLANC and PROFC), obtained from the second derivative of the maximum slope direction and the perpendicular one respectively [43]; • solar radiation (SOLAR); • Topographic Wetness Index (TWI); • flow accumulation (FLOWACC), which represents the contributing area (i.e., the surface over which water from rainfall, snowfall, etc. can be aggregated) [44]; and • Stream Power Index (SPI).
TWI is a parameter correlating topography and the water movement in slopes, used to display the spatial distribution of soil moisture and the shallow saturation degree: where A s is the specific catchment area and β is the slope [45]. SPI is used to describe potential flow erosion and related landscape processes. When specific catchment area and slope steepness increase, both the amount of water contributed by upslope areas and the velocity of water flow also increase, hence stream power and potential erosion increase: where A s is the specific catchment area and β is the slope [46]. From the Landsat 7 TM imagery (July 2016, cloud cover 0%), Clay Index (CI) and Normalized Difference Vegetation Index (NDVI) were calculated. CI is correlated with the clay content of the soil: where MIR is Mid Infra Red (band 6) and MIR2 is Mid Infra Red (band 7) [47]. NDVI gives a quantitative and qualitative estimation of the vegetation: Land 2020, 9, 349 6 of 14 where NIR is Near Infra Red (band 5) and R is Red (band 4) [48]. From the 1:250,000 Soil Subsystems Map of Abruzzo Region by ARSSA [42], an additional variable-SST86-was derived, defined by 29 different soil systems and 87 soil subsystems. Converting the map into a raster, all these soil units were grouped in a single categorical variable.
Finally, a multivariate correlation analysis allowed us to consider as suitable covariates for prediction only those auxiliary data with a relatively stronger spatial correlation with the target soil variables.

Data Processing and Validation of Results
RK, a sort of BLUP (Best Linear Unbiased Prediction) method for spatial data, assumes that the local mean varies continuously into each neighborhood, and can be estimated by combining both directly measured data and correlated ancillary information [28,29]. The technique uses multiple regression to depict the relationship linking the field primary variable and the secondary data. Kriging is then applied to the regression residuals, and the results from both regression and kriging are joined to obtain the estimation [49,50].
The available measured dataset was randomly divided in a training dataset (75% of total samples), and a test dataset (25% of total samples), using only the training part for prediction and the test one to validate the results.
For estimating the target variables (SOC, sand and clay) by RK, the computational steps reported below were followed [34]: 1.
set up and import predictor data layers (land-surface parameters and soil subsystems map); 2.
match soil samples in the training dataset with land-surface parameters and build the regression matrix. Since using the calculated parameters and indices directly as predictors may cause multicollinearity and redundancy effects, the covariates were transformed in principal components (PCs) [51]. Eleven orthogonal and independent components were defined, and the choice of predictors for each variable was then performed by a stepwise regression, considering only the components significant at p < 0.001; 3.
linear regression analysis and derivation of the regression residuals, resolving the regression coefficients by means of a maximum likelihood algorithm [52]; 4.
analysis of residuals for detecting spatial autocorrelation, and fitting of the theoretical variogram models; 5.
run the interpolation; and 6.
visualization and validation of the results using the test dataset.
A flowchart of the procedure is reported in Figure 2. The software SAGA 7.6.3 [53] and ILWIS 3.8.6 Open [54] were used for the derivation of parameters and indices, and for PCs definition. The statistical software R 3.6.3, with the packages sp and rgdal for spatial data preparation and gstat for geostatistical modeling and prediction, was used to perform the regression analysis [55]. The software ArcGIS 10.2 ® was used for drawing the estimated maps of SOC in kg m −2 , soil texture (USDA classification), and SOM levels based on the USDA texture as specified in Table 1. To assess the precision of prediction, estimated values from the training dataset were compared with the correspondent values from the test dataset and not used in the estimation procedure. Such validation allowed us to evaluate the accuracy of the prediction model by measuring the root mean square prediction error (RMSE): where Ẑ(xi) and Z(xi) are the estimated values and actual observations, respectively, and N is the number of validation points. The RMSE expresses the difference between the model estimations and the observed values, presented in the same unit of measurement. If the value of RMSE is close (lower) to the standard deviation of the data, then the model is a good fit [56]. The Geostatistical Analyst extension of the software ArcGIS 10.2 ® was used to validate the estimation results.

Results and Discussion
Pre-processing of data included basic statistics calculation and frequency distribution analysis. Two variables needed to be transformed to approach a normal distribution: sand (square root) and SOC (cube root). The 11 parameters and indices derived from ancillary data (ELEVATION, SLOPE, PROFC, PLANC, TWI, SOLAR, FLOWACC, SPI, CI, NDVI, SST86) were converted in PCs. Table 2 shows the matrix of transformation coefficients, calculated from the covariance matrix. The PC1 and PC2 components were the most significant ones for the prediction of all the target variables. The last component-PC11-was excluded a priori to avoid any rounding effect in the PCs computation with ILWIS [46]. To assess the precision of prediction, estimated values from the training dataset were compared with the correspondent values from the test dataset and not used in the estimation procedure. Such validation allowed us to evaluate the accuracy of the prediction model by measuring the root mean square prediction error (RMSE): whereẐ(x i ) and Z(x i ) are the estimated values and actual observations, respectively, and N is the number of validation points. The RMSE expresses the difference between the model estimations and the observed values, presented in the same unit of measurement. If the value of RMSE is close (lower) to the standard deviation of the data, then the model is a good fit [56]. The Geostatistical Analyst extension of the software ArcGIS 10.2 ® was used to validate the estimation results.

Results and Discussion
Pre-processing of data included basic statistics calculation and frequency distribution analysis. Two variables needed to be transformed to approach a normal distribution: sand (square root) and SOC (cube root). The 11 parameters and indices derived from ancillary data (ELEVATION, SLOPE, PROFC, PLANC, TWI, SOLAR, FLOWACC, SPI, CI, NDVI, SST86) were converted in PCs. Table 2 shows the matrix of transformation coefficients, calculated from the covariance matrix. The PC1 and PC2 components were the most significant ones for the prediction of all the target variables. The last component-PC11-was excluded a priori to avoid any rounding effect in the PCs computation with ILWIS [46].   The target variables were strongly correlated with the principal components, as can be seen from the R 2 values from the linear regression. As explained in the previous section, the choice of the components to be used as predictors was performed by the stepwise regression, considering a significance level of 0.001. For SOC estimation, the sum of PC1-PC2-PC3-PC5-PC7 was used; for sand estimation, the predictor was the sum of PC1-PC2-PC5-PC6-PC10 components; for clay estimation, the sum of PC1-PC2-PC3-PC4-PC6-PC10 components was chosen. Notably, all the components selected as predictors were related to the soil type and position.
RK was applied to the SOC content of the training dataset as well as to the sand and clay data to estimate soil texture. The prediction accuracy was evaluated comparing the values estimated using the training dataset with the measured values from the test dataset that did not enter in the estimation. The calculated RMSE was 0.31 for SOC, 1.76 for sand, and 11.56 for clay, lower than the standard deviation of the measured data of 0.38, 2.09, and 15.02, respectively, and thus considered a very good result, notwithstanding the irregular distribution of the sampling points.
In Figure 3, the estimated map of SOC is reported, showing that about 88% of the area has a SOC content below 3.0 kg m −2 .
In Figure 4, the map of USDA soil texture, estimated from sand and clay data using the training dataset, is reported, showing that the prevailing textures are loam and sandy loam.
Finally, SOM in g kg −1 was evaluated from SOC (Equation (2)), then SOM values were ranked in four classes (very low, low, medium, high) based on the estimated USDA texture (according to Table 1). The obtained map is reported in Figure 5.
The observed SOM levels were related to soil morphology: in the hilly belt along the coast, the map shows essentially a very low content, while in the interior, where soils are mainly under forest, the SOM content reached a high class in most of the territory (about 53%). In the coastal area, the SOM depletion can be considered both as a cause and effect of the active erosion on hills, usually exacerbated by intensive agricultural practices. This area is indeed intensely cropped, as can be seen from the land use map ( Figure S1). Previous studies of the same type carried out in the northern part of the region, concerning a vineyard district, showed that an intervention to enhance SOM content is necessary [15,16,23,34]. Such a result is in line with the findings by Berhongaray et al. [57], who reported that cultivation caused a reduction of 16% in SOC content at 50 cm depth in the Argentine Pampas.
In Abruzzo, an improvement in agro-environmental planning is required to ensure an appropriate SOC/SOM content to soils, so that they can maintain their ecological and socio-economical functions, and crop yield sustainability in agricultural lands. Stimulating the adoption of SSM practices through the dissemination activity of agricultural extension services is fundamental to preserve soil resources. Such practices can enhance soil quality (water holding capacity, infiltration, soil structure, and soil fertility), improve the SOM cycle, and effectively contribute to climate change mitigation. A rational soil management can also help to reduce GHG emissions, especially carbon dioxide emissions, favoring a decrease in soil organic carbon losses, increasing the organic matter input, or combining both of them [15].
Land 2020, 9, 349 9 of 14 estimate soil texture. The prediction accuracy was evaluated comparing the values estimated using the training dataset with the measured values from the test dataset that did not enter in the estimation. The calculated RMSE was 0.31 for SOC, 1.76 for sand, and 11.56 for clay, lower than the standard deviation of the measured data of 0.38, 2.09, and 15.02, respectively, and thus considered a very good result, notwithstanding the irregular distribution of the sampling points.
In Figure 3, the estimated map of SOC is reported, showing that about 88% of the area has a SOC content below 3.0 kg m −2 . In Figure 4, the map of USDA soil texture, estimated from sand and clay data using the training dataset, is reported, showing that the prevailing textures are loam and sandy loam.  Finally, SOM in g kg −1 was evaluated from SOC (Equation (2)), then SOM values were ranked in four classes (very low, low, medium, high) based on the estimated USDA texture (according to Table  1). The obtained map is reported in Figure 5.  The observed SOM levels were related to soil morphology: in the hilly belt along the coast, the map shows essentially a very low content, while in the interior, where soils are mainly under forest, the SOM content reached a high class in most of the territory (about 53%). In the coastal area, the SOM depletion can be considered both as a cause and effect of the active erosion on hills, usually exacerbated by intensive agricultural practices. This area is indeed intensely cropped, as can be seen from the land use map ( Figure S1). Previous studies of the same type carried out in the northern part of the region, concerning a vineyard district, showed that an intervention to enhance SOM content is necessary [15,16,23,34]. Such a result is in line with the findings by Berhongaray et al. [57], who reported that cultivation caused a reduction of 16% in SOC content at 50 cm depth in the Argentine Pampas.
In Abruzzo, an improvement in agro-environmental planning is required to ensure an appropriate SOC/SOM content to soils, so that they can maintain their ecological and socio-economical functions, and crop yield sustainability in agricultural lands. Stimulating the adoption of SSM practices through the dissemination activity of agricultural extension services is fundamental to preserve soil resources. Such practices can enhance soil quality (water holding capacity, infiltration, soil structure, and soil fertility), improve the SOM cycle, and effectively contribute to climate change mitigation. A rational soil management can also help to reduce GHG emissions, especially carbon dioxide emissions, favoring a decrease in soil organic carbon losses, increasing the organic matter input, or combining both of them [15].
This study, defining the current SOC/SOM status in a region of central Italy, provides evidence that SOM management in its agricultural areas should be improved. The results of our study show that in Abruzzo, most arable lands are susceptible to soil degradation, sped-up by dry conditions and high temperatures causing a rapid mineralization of SOM. These areas need appropriate management to guarantee agricultural land sustainability. Adopting conservative practices such as conservation tillage or no-tillage (e.g., direct seeding), improving rotations with forage crops, This study, defining the current SOC/SOM status in a region of central Italy, provides evidence that SOM management in its agricultural areas should be improved. The results of our study show that in Abruzzo, most arable lands are susceptible to soil degradation, sped-up by dry conditions and high temperatures causing a rapid mineralization of SOM. These areas need appropriate management to guarantee agricultural land sustainability. Adopting conservative practices such as conservation tillage or no-tillage (e.g., direct seeding), improving rotations with forage crops, returning crop residues to soil, growing green manure crops, and supplying the soil with proper exogenous organic matter could lead to an appropriate SOM restoration [16]. Permanent grasslands are effective for soil carbon accumulation in mineral soils and the adoption of agroforestry-the integration of trees and shrubs on agricultural land-and crop diversification might contribute to SSM [1].
An accurate state-of-the-art of SOC/SOM distribution would allow us to foresee future trends, and to evaluate the effectiveness of soil conservation practices stabilizing and increasing carbon stock in soils, which should be adopted for a more sustainable soil management. As in our case, mapping SOC over a large area in Bosnia and Herzegovina by classical geostatistical methods [58] showed that the spatial distribution of SOC concentration was strongly influenced by the intense farming practices. In this case, however, the environmental variables had a reduced capacity to explain the spatial variability of SOC. In Brazil, Bonfatti et al. [59] applied RK for mapping SOC, finding again a similar situation: soils under arable crops and vineyard showed the lowest concentration.
In the northern part of Abruzzo, similar mapping work was carried out comparing ordinary kriging and RK as spatial interpolators, and their performance resulted in being approximately the same [34], depending on the available auxiliary information. Mapping targeted soil parameters such as SOC/SOM applying RK, both in limited areas [15] and at the regional level, can represent a valid tool, providing useful and accurate information to assess soil degradation in a cost-effective and little time-consuming way. Such maps can also improve the process of carbon budgeting and reporting-in line with global initiatives of minimizing greenhouse gas emission-and thus the impacts of climate change [60]. Moreover, they can really help farmers in managing soils, can influence and support local land-use planners, and can also be easily updated whenever new data become available.

Conclusions
A spatial representation of SOM is essential to supply a useful and proper reference tool to decision makers, facilitating and optimizing the regional planning of agro-environmental measures for a SSM. Soil surveys taking into account the extremely large variability of the Italian territory are very difficult and onerous. However, most soil attributes are spatially correlated with ancillary variables derived from DEM, Landsat imagery. and existing soil subsystem maps.
In the Abruzzo region, analyzing the available data and estimating values in non-sampled locations by means of RK-integrating measured data and ancillary variables-allowed us to map soil texture, SOC, and SOM levels based on the USDA texture with an acceptable precision. These maps, obtained at relatively low costs, could also be less accurate than the traditional ones, but providing the associated estimation error may anyway bring added value. The observed SOC/SOM distribution appears to be linked to the soil morphology and to the land use: low or very low in intensively cropped hills near the shore, and high in mountains and forest land.
RK proved to be a rapid and cost-efficient tool for mapping soil properties across large areas, allowing us to easily monitor their changes over time. Maps can be updated every time new information and/or new data become available. In the next future, it is hoped that this technique will be successfully coupled and associated with traditional soil surveying and mapping procedures, and then applied at the national level.