Soil organic carbon mapping and prediction based on depth intervals using kriging technique: A case of study in alluvial soil from Sudan

Magboul M. Sulieman a,b,*, Abdallah M. AlGarni c,d a Department of Soil and Environment Sciences, Faculty of Agriculture, University of Khartoum, Khartoum, Sudan b Department of Soil Science, College of Food and Agriculture Sciences, King Saud University, Riyadh, Saudi Arabia c Institute of Environmental Studies, University of Khartoum, Khartoum, Sudan d Member of Saudi Society for Geosciences, College of Science, King Saud University, Riyadh, Saudi Arabia


Introduction
It is well known that soil is the major pool of organic carbon in the terrestrial biosphere, storing more carbon than that in the plants and atmosphere together (Scharlemann et al., 2014). Soil organic carbon (SOC) plays an important role in agricultural productivity (Acín-Carrera et al., 2013), therefore SOC could be a useful indicator for soil health and quality for soil fertility and climate regulation (Beguería et al., 2015;Chen et al., 2017). Level of SOC mainly revert to soil properties such as aggregate stability, soil structure, porosity or nitrogen content (Cotching et al., 2014;Bienes et al., 2016, Sulieman et al., 2018 and resultant soil processes related to land degradation like soil erosion Jin et al., 2009). Thus, nature-based solutions should be applied in order to conserve the soil carbon stocks in some specific degraded environments such as forests, agricultural fields or alluvial areas (Keesstra et al., 2018, Sulieman et al., 2016. In this way alluvial areas over the world are highly demanded for intensive agricultural production, and this is mostly attributed to their high quality of the soils and water as well (Parker, 1995;Bertalan et al., 2016). In the South-east Anatolia region of Turkey, Dengiz (2010) studied the morphology, physico-chemical properties and classified the soils developed on terraces of the Tigris River. He concluded that, the soils developed on floodplain were characterized by weak pedogenesis, however, development of the B horizon (Bw, Bt, Bss) and carbonate accumulation were the main pedogenic processes in terraces soils, and the main soil order were Entisols, Inceptisols, Alfisols and Vertisols. Specifically, in North African regions, non-suitable land use managements are affecting the soil properties and microbial diversity and only well-planned decisions by farmers and policy makers could be revert this dramatically situation (Bounouara et al., 2017;Camilli et al., 2016, Mlih et al., 2016. More specifically, in Sudan, the most irrigated intensive crop farming areas for vegetables and fruits are largely situated within the alluvium plains of the Blue Nile, White Nile, and River Nile terraces. In these alluvial terraces, research on soil properties and their relationships among them that condition their geographical distribution is not highly conducted Sulieman et al., 2015aSulieman et al., ,b, 2016. The distribution of SOC and other soil properties change across the landscape and it also varies by depth , Sulieman et al., 2018. Dengiz et al. (2015) studied the effects of soil types and land use -land cover on SOC density at Madendere watershed, they concluded that, the two factors had great influence on SOC density spatial variation at the different land use systems. In most natural soils, SOC is higher in the surface horizons and it decreases with depth, however, in agricultural areas the spatial variability can be very high (Lado et al., 2004). Such depth-wise variability is mostly continuous in forest or mountainous soils (Bishop et al., 1999;Hartemink and Minasny, 2014) except in soils with a strong human impact like some soils in the Netherlands (Kempen et al., 2011) or after the abandonment (Rodrigo- Comino et al., 2017). Therefore, soil human impacts can be considered as a driving factor of SOC changes (Gómez et al., 2009;Johannes et al., 2017). Despite the most SOC studies and inventories are confined to 30 cm soil depth, the amount of SOC stored below 30 cm is also relevant in many ecosystems (IPCC, 2003;Smith et al., 2005). Sudan is a poorly mapped country in terms of soil data (Sulieman et al., 2018). Likewise, the spatial variability of SOC is not yet understanding all over the country and, most specifically, in alluvial areas. Geographic Information Systems can contribute as a useful tool for soil mapping (Pereira et al., 2018). The creation of soil maps that show the distribution of soil properties along the landscape is considered a powerful tool to organize the territory suitably (Behrens et al., 2010;Grunwald, 2009). To achieve this goal in semiarid environments, different methods have been applied such as modelling (Keshavarzi et al., 2018;Shiri et al., 2017), pedotranfers (Nasri et al., 2015) or geospatial analysis (Zeraatpisheh et al., 2017). In this current study, it was pretend to fill lack of information by applying digital soil mapping techniques to quantify the SOC content and distribution. Therefore, the main objective of this research is to applied ordinary kriging (OK) geostatistical method to predict and map the spatial distribution of SOC at soil depth intervals of 0-30 cm, 30-60 cm, 60-90 cm, and 90-120 cm in selected alluvium soils along the Blue Nile and River Nile terraces in Sudan.

Study Area
The studied samples representing soils from across the upper, middle, and lower Nile River terraces and the upper and lower Blue Nile terraces (10º00ʺ to 15º 45ʺ N, 30º00ʺ to 35º00ʺ E) in Sudan were collected and analyzed ( Figure 1). According to the updated map of climatic classification (Peel et al., 2007), the climate of the study area varies from arid to semiarid zone, characterized by hot dry in most periods of the year to relatively heavy rains during short period in summers. The mean annual temperature varies from 28.3 to 29.6 °C and mean annual precipitation varies from 163 mm to more than 700 mm (Van der Kevie, 1976). Geologically, the alluvium soils along Nile River belong to the Cretaceous Period which comprises continental clastic sediments including sandstones, siltstones, mudstones and conglomerates and locally known as Nubian Sandstone Formation. Whereas, the soils from Blue Nile terraces belong to the Tertiary Period (Gezira Formation), which comprises unconsolidated fine materials and gravels, and Quaternary Periods comprises unconsolidated sands with some gravels, and shales locally known as Umm-Ruwaba Formation ( Figure 2) (Whiteman, 1971;Ministry of Energy and Mines, 1981). Topographically, the study area is mostly gently undulating. The soil moisture and temperature regimes are aridic to ustic and thermic to hyperthermic, respectively, and the soils are classified as fluvents suborder according to USDA soil Taxonomy (Soil Survey Staff, 2014). The main form of land use in the study area is intensive irrigated agriculture for vegetables and fruits. Ecologically, the area lies within the desert to semi desert ecological zone (Harrison and Jackson, 1958), consists of a wide variety of natural plant species, notably Fagonia cretica, Indigofera oblongifolia, Aerva javonica, Acacia toritllis and Maerua crassifolia.

Soil sampling and analysis
A total of 38 soil profiles were excavated along the study area, fully described following the field books of Schoeneberger (2012), classified according to USDA soil taxonomy (Soil Survey Staff, 2014). Based on the genetic horizons, a total of 152 soil samples were collected from four depth intervals of 0 to 120 cm within each profile as follows; 0-30 cm, 30-60 cm, 60-90 cm, and 90-120 cm. In laboratory, soil samples were first air-dried, ground, and then passed through a 2 mm sieve to obtain the fine fraction. The soil organic carbon was measured by wet oxidation with chromic acid and back titrated with ferrous ammonium sulfate following the standard Walkey and Black method (Nelson and Sommers, 1996). A histogram tool in Arc GIS software version 10.4 was used to identify the outliers of SOC, where a more isolated bar from the main group of bars in the histogram was taken to represent a possible outlier (Chabala et al., 2017). Therefore, in case if outliers observed, the laboratory analyses were repeated and any outliers recorded one more were removed from the samples set.

Statistical analyses
General basic statistics parameters such as minimum and maximum values, mean, median, standard deviations, first quartile, third quartile, Skewness, and Kurtosis were calculated for both tested and validated samples to provide a basic understanding of the characteristics of data using Statistix software version 8.1 (Steel and Torrie, 1980).

Predicting spatial distribution of SOC
In this study, the spatial variability of SOC at four depths intervals was studied using ordinary kriging (OK) method. OK is a geostatistical models that use a set of statistical tools to predict the value of a given soil property (in this case SOC) at a location that was not sampled (Goovaerts, 1998;Regalado and Ritter, 2006;Johnston et al., 2001). Based on OK, the predicted SOC at an unsampled location (pred(x0)), using measured values obs(xi) is given as: Where ( 0 ) is the kriging weight. In order to model the spatial variability of SOC with OK, further data evaluation was performed to check whether the SOC data conformed to the necessary assumptions required for OK include that data are required to be normally distributed with no trend. The histogram and normal quantile-quantile (Q Q) plots tools in ArcGIS were used to evaluate the distribution of the data. The trend analysis tool was also applied to check for trend in the SOC data set. The data check revealed that the data set had no trend. The central tool of geostatistical procedures such as OK is the semivariance, which is half the expected squared difference between the SOC values at two locations. The semivariance quantified the spatial variability of SOC for all possible pairing of SOC data. The plot of the semivariance as a function of distance is referred to as a semivariogram (Chahouki et al., 2011). The semivariogram described how SOC varied with distance among the sampling locations. Positive definite models were applied to fit the semivariogram using ordinary least squares to capture spatial features of the SOC. An automated fitting procedure was followed when fitting the semivariogram of SOC using the exponential model. All data processing and analysis for OK were done in ArcGIS software package version 10.4 with the variogram redrawn in Microsoft excel.

Data validation
The prediction accuracy of SOC at the unsampled locations was evaluated using the leave-out-one cross validation (LOCV) techniques (Reza et al., 2010;Liu et al., 2013;Martín et al., 2016). The indices used during LOCV were root mean square error (RMSE), average standard error (ASE), and standardized RMSE (RMSSE). Thus, for each of the sampled locations, there was a measured value for SOC (obs(xi)) and a predicted value (pred(xi)), with standard values being obs1(xi) and obs2(xi), respectively. The ASE, RMSE, and RMSSE were calculated according to Yang et al. (2009) as follows: Where n is the number of validation points. The RMSSE varies between 0 and 1, with a value of 1 indicating a perfect model and 0 indicating a poor model. In addition, a cross plot of measured and predicted SOC values was drawn in a scatter plot and associated summary statistics of predicted SOC values were compared with those of measured values (Leuangthong et al., 2004) using the Sample t-test (Paired t-test).

Summary of basic statistics
The summary of descriptive statistics for the measured and predicted SOC data are presented in Table 1. Results show that the SOC showed a decreasing trend with increasing soil depth. Mean SOC for both measured and predicted was 0.75, 0.64, and 0.50 g kg −1 for soil layers 0-30 cm, 30-60 cm, and 90-120 cm, respectively, while the mean for measured and predicted SOC at layer 60-90 cm was 0.56 and 0.59 g kg −1 , respectively. This indicates that the model was reliable for SOC prediction in this depth compared to the other three depths, meanwhile, it can be taken into account. The standard deviation (±SD) in the measured SOC was 0.53, 0.52, 0.46, and 0.43 g kg −1 for layers 0-30 cm, 30-60 cm, 60-90 cm, and 90-120 cm, respectively compared to 0.29, 0.29, 0.28, and 0.26 g kg −1 , respectively for the same layers in the predicted SOC. The coefficients of variation (CV) of SOC for both measured data for all layers ranged from 71.56 to 86.72 % with gradually increased from the surface to the bottom layers in the soil profile. This could be due to the complex textural layers. This result is in line with previous studies of He et al. (2009) and Chen et al. (2015). A coefficient of variation (CV) value of 10% indicates low variability and values ranging from 10-90% indicate a moderate variability; CV values > 90% indicate high variability (Fang et al., 2012). Therefore, measured SOC in our study indicated a moderate variability. The coefficients of skewness at four depths were 0.53, 0.61, 71, and 1.23 in the measured SOC, respectively and 0.49 and 0.20, 0.23, and 0.53 in the predicted SOC, respectively. This indicated that the measured and predicted SOC followed similar pattern distribution. Further, the skewness (<0.5) and kurtosis (<3) of SOC indicated the approximate normal distributions of data (Rosemary et al., 2017). The first and third quartiles at the four depths were (0.51, 0.91), (0.34, 0.90), (0.40, 0.78), and (0.30, 0.61) g kg−1, respectively in the predicted values, which were comparable to the recorded values in the measured data. Overall, the summary statistics in the measured and predicted SOC indicated that the model was well fitted for the prediction of SOC at the different soil layers.

Variography and kriging interpolated of SOC at different depth intervals
Figures 3 shows the fitted semivariogram and the associated variograhpic parameters for the different depth intervals. Depending on the criteria suggested by Cambardella et al. (1994), the spatial dependence or autocorrelation was judged to be strong if the nugget to sill ratio was less than 0.25, to be moderate if between 0.25 and 0.75, and to be weak if higher than 0.75. In this study, the spatial autocorrelation (SAC) at depth intervals of 0-30 cm, 30-60 cm, and 60-90 cm has a nugget to sill ratios of 0.42, 0.28, 21 respectively, and at layer 90-120 cm was 0.86 ( Fig. 3a,b,c,d), demonstrating a moderate to weak spatial dependence of SOC for the different soil layers which was controlled by both intrinsic and extrinsic factors. This was evidenced by local, intensive agricultural management, such as tillage and fertilization, widely used in the study area to improve the fruit and vegetables productivity (Sulieman et al., 2016). The maps of SOC spatial distribution generated depending on the measured SOC values and fitted variogram is shown in Figure  4. However the standard error maps of SOC spatial prediction is shown in Figure 5. According to Fu et al. (2014), a strong spatial dependence of soil properties is attributed to soil intrinsic properties, such as soil parent materials, soil texture, topography and vegetation. Whereas, a weak spatial dependence of soil properties indicates that the spatial variability is mainly regulated by extrinsic variations, such as soil fertilization and cultivation practices, and moderate spatial dependence is controlled by both intrinsic and extrinsic factors (Kılıç et al., 2004;. The exponential nature of the fitted semivariogram may indicate that SOC at the site had a gradual transition or that several patterns interfered. This variation in SOC levels can be attributed to the vertical illuviation of mineral and organic material from surface layers downward the soil profiles (Brodsky et al., 2013).

Modeling validation
The indices generated during the leave-out-one cross validation (LOCV) of OK model at four different depths are presented in Table 2. It was observed that goodness of fit between observed and predicted values was higher in surface soil layers than sublayers. The average standard error was 0.58 at depth 0-30 cm (RMSE = 0.47), 0.57 at depth 30-60 cm (RMSE = 0.43), 0.51 at depth 60-90 cm (RMSE = 0.39), and 0.48 at depth 90-120 cm with the RMSE of 0.35. The standardized RMSE was 0.83, 0.79, 0.80, and 0.79 for the same depths grades, respectively. The RMSSE with a value closed to 1 indicating a perfect model (Tang et al., 2017). Overall, the RMSSE values in our study indicating that the model was correctly estimating the variability of SOC at the unsampled locations at the four different layers, and the best model was observed at the surface layer (0-30 cm). Chabala et al. (2017) studied the spatial distribution of SOC at depth 0-20 cm in soils of a selected part of Zambia by using OK technique, they conclude that the model generated by OK was reliable in predict SOC and produce a RMSSE of 1.02. Figure 6 shows the predicted versus measured SOC, and the cross plot clearly indicated a weak to medium positive correlations. Further, there were no significant differences (P >0.05) between the measured and predicted SOC in all soil layers when compared by the Sample t-test (Paired t-test).

Conclusion
In this study, OK was applied to study the spatial interpolation of SOC at 0-30 cm, 30-60 cm, 60-90 cm and 90-120 cm using the measured data from 152 samples in alluvium soils along Blue Nile and River Nile in Sudan. The results indicated that the short-range spatial dependence was moderate to weak with a nugget shifted from zero. Spherical model was selected to describe the spatial pattern of SOC in the study area. A moderate to weak spatial dependence of SOC was observed, indicating that SOC was controlled by both intrinsic factors (e.g. soil parent materials and soil texture) and extrinsic factors (e.g. application of fertilizers and tillage treatment).