Abstract

Soil organic carbon contents are expected to vary from place to place because of variation in soil properties. However, the extent of variability has not been explored in the study area. This study has, therefore, been initiated to assess the spatial variability of soil organic carbon stock in Gurje subwatershed Hadiya Zone, Southern Ethiopia. A total of 40 randomly predefined sampling points were identified for soil sampling using GIS and a total of 80 composite soil samples and 80 core samples were collected from those points at two sampling soil depths (0–20 cm and 20–40 cm). The ordinary kriging (OK) method was used as a geostatistical tool and applied to model the spatial variability of soil organic carbon in this study. With respect to soil depth, the coefficient of variation (CV%) for SOC and SOCS varied from 40.87 to 51.36%, which indicated moderate variability in the study area. For the land use types, the CV% varied from 7.94 to 42.06%, indicating low to moderate variability for the variables in the study area. The exponential semivariogram model described the spatial structure of SOC at 0–20 cm depth while the spherical one was used for SOCS. Moreover, the exponential model was best suited for SOCS at a soil depth of 20–40 cm, while the circular model was appropriate for SOC at this depth. The nugget/sill ratio (C0/C0 + C) of SOC and SOCS varied from nil to 15.58, reflecting a strong spatial dependence, which could be mainly due to the influence of intrinsic factors (e.g., natural variations in soils) in the study area. Overall, the spatial distributions of SOC and SOCS were higher in the northwestern and eastern parts of the subwatershed.

1. Introduction

Soil properties vary from place to place in the landscape even within a pedon. The spatial variability of soil properties is the combined effect of natural and human factors [1]. Specifically, soil-forming parent material and the topographic position of soil on the catena cause variability in the soil properties [2]. Moreover, soil properties change with time and place [3]. As a result, its physical, chemical, and biological characteristics change simultaneously. This change causes soil spatial variability, which plays a critical role in ecological and environmental modeling processes on the landscape [4]. Therefore, assessing the spatial variability of soil properties is used in developing site-specific soil management [5].

Soil acts as a pool of organic carbon in the terrestrial ecosystem. The variability of soil organic carbon (SOC) is associated with parent materials, topography, climate, and biota [6]. The carbon sequestration potential of soil also varies with soil type and physicochemical characteristics, management, and environmental factors [7]. Furthermore, land use causes spatial changes in the SOC contents [1]. Thus, soil properties are significantly affected by changes in land use and management practices [8], with similar effects on SOC. The SOC densities are also affected by elevation and slope [9]. Therefore, SOC varies spatially across the landscape [10].

The spatial variability of SOC, on the other hand, is an indicator of variability in soil quality and carbon pool in the terrestrial ecosystem [11]. Moreover, spatial distribution is useful in developing carbon cycle models and assessing soil quality in the ecosystem because of its vital role in global carbon cycling [12]. Sustainable agricultural land use and management practices also require an understanding of the spatial variability of SOC [13]. Different land use and soil management practices have variable contributions to the nature, quality, and quantity of carbon storage [14]. Moreover, SOC determines many ecosystem services, which are in turn influenced by the nature of the ecosystems.

SOC acts as a determinant of agricultural productivity [15]. It influences soil properties in terms of soil fertility and quality [16]. Thereby, SOC controls soil fertility bases such as soil structure, water-holding capacity, cation exchange capacity (CEC), and carbon to nitrogen ratio (C/N) [17]. Therefore, knowledge of SOC spatial distribution could be crucial for site-specific management strategies.

Spatial variation of SOC resulting from different pedogenic factors could be estimated by geostatistical methods. Geostatistical techniques are used to analyze the spatial variability, distribution of soil properties, and evaluate their spatial dependence [18]. These techniques allow the estimation of SOC stocks using soil information and environmental covariates on landscapes. Spatial interpolation is one of the geostatistical methods used to predict the value of soil properties (e.g., SOC) at unsampled locations by accounting for the spatial correlation between sampled and estimated points [19]. Therefore, SOC is estimated using spatial interpolation techniques considering various affecting factors [17].

The SOC dynamics is affected by land use change and still understudied in sub-Saharan Africa (SSA) [20] because of an incomplete understanding of how SOC changes are affected by various factors and the vital role of ecosystems in global carbon cycling based on the dynamics of SOC accumulation and turnover quantification [21]. Ethiopia is characterized by having a variety of organic carbon contents [22] due to its diverse nature of soil types on different landscapes. However, the study on the spatial variability of SOC is limited and did not touch different arrays of watersheds. Although previous studies on SOC distribution are reported to concentrate in the southern part of the country [14], it did’t cover Gurje subwatershed, Hadiya Zone, Southern Ethiopia. Consequently, the study area was suffered from lack of information on the spatial variability of SOCS to devise a site-specific soil management plan. To tackle this, therefore, the current study was conducted to assess the spatial variability of SOCS in this subwatershed.

2. Materials and Methods

2.1. Description of the Study Area
2.1.1. Location

The study was conducted in Gurje subwatershed in Anabalesa Kebele (the lower administrative unit in Ethiopia), Lemo District, Hadiya Zone, Southern Nations Nationalities, and Peoples (SNNP) Regional State of Ethiopia. The total area of the Gurje subwatershed is 451 hectares. Lemo district is about 230 km southwest of Addis Ababa (the capital city of Ethiopia) and 175 km west of Hawassa (the regional capital). However, Gurje subwatershed is also located in this district about 7 km southeast of Hossana (the zonal town of Hadiya). Geographically, Gurje subwatershed is located between 7°34′30″and 7°36′0″N and 37°55′0″ and 37°56′0″E (Figure 1).

2.1.2. Relief and Climate

The topographic features of the Lemo district are flat, hilly, undulating landscapes with altitude ranges of 1900–2720 meters above sea level (masl) [23]. The district has two different agro-climatic zones, namely: Woyina Dega (midland) and Dega (highland). The mean annual rainfall of the study area is 1200 mm (range from 900–1400 mm) and the mean minimum and mean maximum temperatures are 13 and 22°C, respectively (Figure 2). Moreover, the rainfall distribution of the area is a bimodal type, which occurs in two main rainy seasons, Belg (March–May) and Maher (June–August).

2.1.3. Geology and Soils

The study area is found under the division of Southern Main Ethiopian Rift (SMER) segment which has less evolved rifting (decreased in maturity) compared to the northern MER system. Moreover, the SMER was under the category of Miocene time because this geological date was the onset of the Main Ethiopian Rift (MER) system [25]. However, quaternary volcanic activity in the SMER is limited to the rift margins (e.g., Soddo) except for the basaltic activity that gave rise to the land bridge separating Lake Abaya from Chamo. However, the study area was found close to the Bilate River Basin of the southern MER and it shares the geological units of the Nazareth group (Upper Miocene succession), ignimbrites, tuffs, ash flow, rhyolites, and trachyte, a silicic volcanic complex of Damota and quaternary volcanic eruptive and sedimentation [26].

According to [23], the dominant soil types of the Lemo district are largely derived from igneous volcanic rocks, mainly ignimbrites by weathering processes. The weathered ignimbrite rocks are the main sources of colluvial and alluvial deposits. This is due to the slope and the fluvial action originating in the elevated area and depositing in a low-lying area. Ignimbrite rock outcrops along the escarpment; the hillsides are highly jointed sets following tectonic trends. As a result, the dominant soil types in the study area include Nitisols and Luvisols with loamy clay due to the clayey texture. These soils are darkish brown to reddish in color and characterized by a good agricultural potential for crop production.

2.1.4. Land Use/Land Cover and Farming Practices

Lemo district is characterized by different land use/land cover including forest, bush, cultivated, and grazing land uses. Cultivated land constitutes about 74% of the land use. The vegetation cover (bush, shrubs, and forest) is about 4.5% out of the total land area. The coverage is diminishing from time to time at alarming rates because of the high population pressure and agricultural land expansion, urbanization, fuelwood consumption, etc. The remaining 21.5% of land uses are used for other functions.

Subsistence crop–livestock mixed farming is the most dominant agricultural activity in this area. Growing trees and/or shrubs is commonly practiced in the agricultural landscape, and rain-fed crop production is the major economic activity of the area. The staple food crop of the area is Enset (Ensete venrtricosum) while the main cash crop is khat (Catha edulis) followed by coffee (Coffea arabica). In addition to these, fruits (e.g., avocado—Persea americna), vegetables (e.g., cabbage—Brassica oleracea), cereal crops (e.g., wheat—Triticum spp., maize—Zea mays), and root and tuber crops (e.g., potato—Solanum tuberosum, beetroot—Beta vulgaris) are the most common crops of the study area [23].

2.1.5. Field Survey and Soil Sampling Procedure

A preliminary field survey was conducted using a topographic map of 1 : 50,000 scale. Besides, randomly selected farmers were interviewed to collect information on general overviews on land use types and their management practices. Soil sampling points were randomly predefined on the entire subwatershed using Geographic Information System (GIS) according to the sampling distribution procedure of [27]. Based on the predefined sampling point loaded in the tag, the surveyor navigated to a selected sampling point using a geographical positioning system (GPS) receiver (model Garmin GPSMAP 60Cx), which is always the center of the 10–15 subsampling points. A total of 40 georeferenced sampling points were administered for soil sampling. At each point, soils were sampled at two depths (0–20 and 20–40 cm) using an auger sampler for disturbed and a core sampler of 5 cm diameter for undisturbed samples. The sampling points randomly fell in different land uses/covers such as cultivated (with annual and perennial crops), forest, and grazing lands of the watershed. For each soil sample from each sampling point, about 0.5 kg composite sample was produced from 10–15 subsamples around the georeference point to minimize bias, which is generated from the heterogeneity of soils, while one core sample per depth was taken from each sampling point for bulk density determination. A total of 80 composite soil samples and 80 core samples (40 soil samples per depth of each) were packed and labeled for further laboratory analysis.

2.1.6. Soil Sample Preparation and Laboratory Analyses of the Selected Physical and Chemical Properties

Soil samples were air-dried, crushed, and passed through a 2 mm sieve for the analysis of the selected soil physicochemical parameters except for organic carbon (OC) and total nitrogen (TN), in which case soil samples were sieved with a 0.5 mm mesh. Soil particle size distribution was determined by the Bouyoucos hydrometer method [28]. The bulk density (ρb) was determined from undisturbed soil samples using the core method after drying it in an oven at 105°C to a constant weight [29]. Soil pH-H2O was determined potentiometrically using a pH meter in the supernatant suspension of 1 : 2.5 soil to water ratio [30]. Soil OC was determined by the wet oxidation method [31]. Modified micro-Kjeldahl procedure was employed for the determination of TN as described by [32].

The coarse fractions were determined after the clumps were broken by hand, crushed, grinded, dried, and sieved until the samples can pass through a 2 mm sieve. Thereafter, the coarse fractions (>2 mm) were weighed and their fractions were determined as described in [33]:

2.1.7. Determination of Soil Organic Carbon Stock

The soil organic carbon density SOCD (kg C/m2) was estimated for a single soil depth layer based on the following equation as described in [34]:where SOCDi and SOCi are the SOC density (kg cm−2) and concentration (gkg−1) of the ith layer, ρi is the soil bulk density of the ith layer (gcm−3), di is the depth of the ith layer (cm), and CFi is the fraction (%) of the coarse fragments (>2 mm) in the ith layer; 100 is a conversion factor.

For the individual profiles with a depth, D (cm), SOC densities were calculated by summing up the SOC density of each soil layer of each land use type:

The SOC stock in each soil depth layer of each land use type of the watershed was estimated by multiplying the SOC density in each unit area (kg C/m2) by the total area covered by a particular land use. The summation of SOC stocks in each depth layer gives the total SOC stock in each land use type in the watershed. Therefore, the total SOC storage is calculated as:where SOC stock is the total amount of SOC storage (kg C) at depth, D (cm) of a given land use type, m is the number of layers, i is the ith layer, SOCDi is the average SOC density of the jth land use type (kg cm−2), and Sj is the total area (m2) of a given land use type j in the watershed [34].

2.1.8. Statistical Data Analysis

One-way analysis of variance (ANOVA) was carried out to show the SOCS variability among different fields using SAS software [35] version 9.2. Moreover, the descriptive statistics (mean, standard deviation, maximum, minimum, skewness, and kurtosis) were employed for SOC and SOCS using this software. The Kolmogorov–Smirnov test for normality of the data distribution was employed using a Statistical Package for the Social Science (SPSS) version 20. The coefficient of variation (CV) was obtained for SOC and SOCS in the entire field of the subwatershed.

2.1.9. Geostatistical Analysis

Geostatistical analysis was performed using ArcGIS10.7.1 software to determine the degree of spatial variability and dependency of SOC and SOCS using geostatistical tools for variogram calculation, ordinary kriging (OK), validation, and mapping. The OK technique was performed because it is the most reliable of all predicting techniques based on mean error (ME) [36]. The OK technique is also the best unbiased predicting method in cases in which the soil sample points were selected randomly and sparse to predict the values of the soil properties at the unsampled point [19]. Moreover, the OK method was performed to interpolate the values of the unsampled locations and produce maps of soil properties (e.g., SOC and its stock) [37]. The experimental semivariogram modeling was used to assess the spatial pattern of soil variables (e.g., SOC and SOCS), which were calculated as described in [38]:where γ (h) is the experimental semivariogram value at the distance interval h; N (h) is the number of sample value pairs within the distance interval h; z (xi), z (xi + h) are the sample values at two points separated by the distance interval h.

All pairs of points separated by distance h (lag h) were used to calculate the experimental semivariogram. The semivariogram models provide information about the spatial structure as well as the input parameters for the kriging interpolation [39]. The experimental semivariograms were fitted by three theoretical models such as exponential, spherical, and circular separately, and the best model for each parameter was selected based on the fit [40]. The fitted models were checked based on the error values computed from the entire dataset [41]. In this regard, cross-validation was employed to evaluate the performance of each model based on the prediction errors such as mean error (ME) (equation (6)), root mean square error (RMSE) (equation (7)), average standard error (ASE) (equation (8)), mean standard error (MSE) (equation (9)), and root mean square standardized error (RMSSE) (equation (10)). Accordingly, the model showing RMSE close to the ASE, ME and MSE values close to zero, and RMSSE value close to one is considered as the best-fitting model. After conducting error evaluation, the best-fitted model was selected for prediction. Then, kriging maps indicate the values of a variable at unsampled locations based on sample observations. Therefore, the maps provide a visual representation of the distribution of the soil variables (e.g., SOC and SOCS).where Z (Xi) is the value of the variable Z at location Xi, Ž (Xi) is the predicted value at location i, n is the sample size, and is the kriging variance for location Xi.

As described by [41], ASE is either a result of over or underestimation. If the MSE is more than the RMSE, the variability of the prediction is overestimated, while if the RMSE is more than the ASE, the variability is underestimated. On the other hand, if the ASE is close to the RMSE, the variability prediction is correctly assessed. The RMSSE should be close to one of the prediction standard errors and valid; if the RMSSE is more than one, the variability of the prediction is underestimated, while if it is less than one, the variability is overestimated.

Spatial interpolation was carried out with a point kriging approach using the fitted models [42]. Using the model semivariogram, basic spatial parameters such as nugget variance (C0), structure variance (C), range (A), and sill (C + C0) were calculated. Nugget variance is the variance at zero distance, sill is the lag distance between measurements at which one value of a variable does not influence neighboring values, and range is the distance at which values of one variable become spatially independent of another [43]. As indicated in [44], nugget/sill ratios of <25%, 25–75%, and >75% were classified as strong, moderate, and weak spatial dependence, respectively. More often, for the ratio of 100% or if the slope of the semivariogram is close to zero, the soil variable is considered non-spatially correlated [43].

3. Results and Discussion

3.1. Descriptive Statistics of SOC and SOCS

The soil organic carbon (SOC) concentration and soil organic carbon stock (SOCS) are presented in Table 1. At both sampling depths, the SOC showed a positive skewness and negative kurtosis. The lower mean value of SOC (17.09 g/kg) compared to the median (17.45 g/kg) at 20–40 cm soil layer implies the presence of positively skewed data [16]. However, SOCS (t/ha) showed a positive skewness and kurtosis for the referred soil depths. The Kolmogorov–Smirnov test of normality (K-S) values ranged from 0.003 to 0.155 for the variables at both sampling depths (Table 1). The K–S test values (K-S, ) indicated a normal distribution of data. Moreover, the values of skewness and kurtosis were 0.12 and −0.98, respectively, for SOC at the depth of 0–20 cm and showed nearly normal distribution of data (Table 1). According to [45], the values of skewness (<0.5) and kurtosis (<3) reflect the approximate normal distribution of data.

The coefficient of variation (CV) is the ratio of the standard deviation to the mean expressed as a percentage, and is a useful measure of the overall variability of SOC in the study area. It varied from 40.88 for SOC to 51.49% for SOCS within the upper 40 cm soil depth (Table 1). According to [16], the variability of soil properties could be classified based on the CV as low variability if the CV is less than 10, while it will have moderate variability if the CV is between 10% and 90% and high variability if the CV is greater than 90%. As such, the CV values of SOC and SOCS of the studied soils were under the category of moderate variability for both soil depths (Table 1). However, the CV increased with increasing soil depth for all variables, implying the increasing trend of variability in the deeper soil layers. Therefore, a soil depth of 20–40 cm had a higher CV compared to the surface (0–20 cm) for both variables. Similarly, [34] reported the highest variability of SOC concentrations in the subsurface soil layer. Generally, the results showed a higher CV% for those variables, which could be associated with the heterogeneity of land use types, fertilizer applications, or erosion [13].

3.2. Selected Soil Properties at Different Soil Depths

The mean values of soil particles such as sand, silt, and clay at 0–20 cm were 31.10%, 28.75%, and 40.15%, respectively, while the values at a depth of 20–40 cm were 30.05%, 28.95% and 41.05% for the respective soil particles. The mean values of soil bulk density (BD) at 0–20 cm and 20–40 cm depths were 1.117 and 1.121 g/cm3, respectively (Table 2).

However, the mean soil pH values were 5.67 and 5.90 at 0–20 cm and 20–40 cm soil depth layers, respectively, which showed a statistically significant difference () between the two soil depths (Table 2). Soil pH increased with depth which might be due to the downward leaching of basic cations. The mean values of total nitrogen (TN) at 0–20 cm and 20–40 cm soil depths were 0.163% and 0.128%, respectively (Table 2). The vertical distribution of TN revealed a significant variation () between soil depths. The mean values of carbon to nitrogen ratio (C/N) at the soil depths of 0–20 cm and 20–40 cm were 13.161 and 13.251, respectively. However, the difference was not statistically significant () between soil depths for the C/N ratio (Table 2).

The mean values of SOC and SOCS at 0–20 cm were 21.51 (g/kg) and 34.06 (t/ha), respectively, while their means were 17.09 (g/kg) and 25.85 (t/ha), respectively, at a soil depth of 20–40 cm (Table 2), implying the SOC and SOCS decreased with increasing soil depths. This might be due to the topography, climate, and vegetation that influence the vertical distribution of SOC [46]. Similarly, [47] reported a decreasing trend of OC contents with depth in mineral soils. In this study, the surface soil layer (0–20 cm) had higher SOC and SOCS contents compared to the subsurface ones (20–40 cm) (Table 2). It might be due to the surface soil layer being affected by the application of organic inputs, falling litters, and microbes. The application of organic inputs may have a positive influence on the surface soil [48]. Particularly, a strong accumulation of organic materials from litter, roots, and microbes at the surface [12] will contribute to a higher SOC compared to the subsurface. The CV of the selected soil properties varied from 7.41 to 46.23% at 0–20 cm while it varied from 6.10 to 51.50% at 20–40 cm (Table 2).

3.3. Variation of Soil Organic Carbon in Different Land Use Types

The SOC in cultivated land with perennial crops (25.12 g/kg) was significantly higher () than that in cultivated land with annual crops (15.32 g/Kg) (Table 3). Moreover, the mean of the SOC concentration in the forest land (33.517 g/kg) was significantly higher () than in the grazing land (29.3 g/kg) (Table 3). This is probably due to the presence of more vegetation generating more litter falls which are returned to the soils as organic matter [49]. Besides this, the difference in land use types causes variation in the SOC contents [48]. Therefore, land use exerts a great influence on the SOC variation [50].

However, the lowest mean SOC recorded in the cultivated land with annual cropping systems (Table 4) might be due to the depletion of organic matter, resulting from continuous cultivation, removal of crop residues, and little or no supply of organic inputs. As described in [51], a continuous soil cultivation results in lower soil organic matter content, and in turn reduces the amount of SOC due to low organic matter turnover. Consequently, the lower organic matter content reduces the nutrient and water-holding capacity and soil biological activity, which in turn affects the ecosystems because high biological activity favors the absorption and accumulation of carbon in the ecosystems [52]. On the other hand, the higher SOC in the perennial cropping systems compared to annual ones (Table 3) could be because the perennial systems enhance more inputs of organic residues than annual crop-based systems [51], implying the difference in management. Generally, the concentration of SOC varied significantly () among the land use types. However, no significant variations () were observed between perennial crop and grazing land uses. The grazing land revealed higher SOC contents compared to cultivated land. It could have resulted from the high amount of grass roots and their biomass turnover rate in the soils [53]. The CV of SOC varied from 7.94 to 42.06%, which was highest in cultivated land (annual cropping system) and the lowest in the grazing land use. Thus, the CV of SOC was under the category of low to moderate variability in accordance with [16].

The SOCS in the cultivated land with perennial cropping practices (33.655 t/ha) was significantly higher () than that in cultivated land with annual cropping (23.694 t/ha) (Table 3). However, SOCS in the forest (67.662 t/ha) was significantly higher than that in the grazing land use (41.229 t/ha) (Table 3). The mean of SOCS varied significantly () among the land use types except between perennial crop use and grazing (Table 3). The SOCS in the cultivated land with annual crop-based systems was lower compared to other land uses in this study. This agrees with [54] who reported lower SOCS in agricultural areas, especially in cultivated lands. The CV of SOCS varied from 8.82 to 40.41% (Table 3), implying low to moderate variability at different land uses as previously reported in [16].

3.4. Correlation Analysis

At a soil sampling depth of 0–20 cm, the SOC and SOCS correlated significantly with sand (r = −0.652 and −0.527, respectively) and clay (r = 0.531 and 0.606, respectively) at (). Moreover, the SOC and SOCS significantly correlated with clay (r = 0.341 and 0.401, , respectively) at a depth of 20–40 cm (Table 4). This reflects soil textures such as sand and clay influence the organic matter content of soil by affecting its decomposition rate. For instance, soils with high clay content may have higher organic matter content. This is because a clay particle protects soil organic matter from breakdown by binding strongly and creating physical barriers to microbial access [55]. The SOC and SOCS showed an inverse relation with soil bulk density (Table 4). However, TN had a significant () and positive correlation with SOC and SOCS at different sampling depths (Table 4).

3.5. Geostatistical Analysis of SOC and SOCS

Geostatistical analysis showing the spatial variation, distribution, and dependence of SOC and SOCS at both sampling depths are presented in Table 5 and Figure 3. The exponential semivariogram model provided the best fit for SOC at a soil depth of 0–20 cm, and the circular model for a depth of 20–40 cm . The spherical semivariogram model described the variation of SOCS at 0–20 cm and the exponential model at a depth of 20–40 cm (Table 5 and Figure 3). The nugget/sill ratio (C0/C0 + C) varied from nil to 15.58 at both sampling depths in the study area (Table 5). The results revealed lower nugget/sill ratio for SOC and SOCS. According to [44], the nugget/sill ratio (C0/C0 + C) < 25% reflected a strong spatial dependence of SOC and SOCS. This could be mainly due to the inherent/intrinsic sources of variability (e.g., natural variations in soils, such as soil texture, soil type, parent materials, and topography) [16, 56, 57].

The nugget values (C0) of SOC and SOCS varied from nil to 0.0342 (Table 5). This revealed lower nugget was an implication for micro-variability. As suggested by [58], the lowest nugget values showed they had low spatial variability within small distances. The partial sill (C) of SOC and SOCS varied from 0.245 to 0.270 and 0.182 to 0.278, respectively (Table 5). In this study, the partial sill (C) increased with increasing soil depth, which indicated the total spatial variation, as was reported by [34].

The range of the semivariogram or spatial range varied from 270.9 to 394.5 m for all variables at both sampling depths (Table 5). At a sampling depth of 0–20 cm, the SOC had a high range of semivariogram (381.4 m), while the SOCS (394.5 m) did show a high range at a depth of 20–40 cm (Table 4). Overall, a higher spatial range was observed at the subsurface compared to the surface soil layer (Table 5). This might be due to pedogenic processes like eluvial-illuviation that homogenizes the distribution of SOCS and increases the spatial autocorrelation distance in the subsurface soil layer. Hence, a larger spatial range implies a larger spatial autocorrelation [34] and a broader homogenous scale [59]. Conversely, the smaller spatial range at 0–20 cm indicated the relatively high spatial heterogeneity of the SOC content in this layer [12].

The ME values of SOC and SOCS were 0.0766 and 0.0425, respectively, at 0–20 cm depth, and 0.0027 and −0.0617, respectively, at 20–40 cm (Table 6). The MSE values were −0.0694 and −0.0752 at 0–20 cm depth for SOC and SOCS, respectively, and −0.093 and −0.0616, respectively, at 20–40 cm (Table 6). The cross-validation indices like ME and MSE had values close to zero (Table 6), verifying that the OK interpolation is a reliable technique to predict the spatial distribution of SOC and SOCS at both sampling depths as it was reported similarly by [16]. Moreover, the ME and MSE values are close to zero, which demonstrates that the prediction is relatively unbiased with a small bias [60].

However, the RMSE values of SOC and SOCS were 9.98 and 17.40, respectively, at 0−20 cm and 8.80 and 14.37, respectively, at a depth of 20−40 cm (Table 6). The RMSE of SOCS had higher values at 0−20 cm than SOC due to its larger magnitude. The RMSSE of SOC and SOCS were 0.9509 and 0.9945 at a soil depth of 0−20 cm and 1.0303 and 1.0379, respectively, at 20−40 cm (Table 6). The results showed that the RMSSE values were close to unity (Table 6), reflecting a perfect or best-fitting model [61] and considered accurate [19]. The average standard error (ASE) values of SOC and SOCS were 10.77 and 17.23, respectively, at 0−20 cm depth, and 9.14 and 12.78, respectively at 20−40 cm (Table 6). The results showed that the ASE values were close to the RMSE for the variables at both depths except for SOCS, which did show a higher value of RMSE than ASE at the sampling depth of 20−40 cm. According to [41], if the ASE is close to the RMSE, the variability in prediction is correctly assessed, while if the RMSE is more than the ASE, the variability is underestimated.

3.6. Spatial Distribution of SOC and SOCS

The kriged/interpolation maps of SOC and SOCS produced using the best-fitting models of OK at different depths are presented in Figure 4. These maps signified the spatial distribution of SOC and SOCS in the study area. The lowest SOC concentrations were distributed mainly in the periphery of northwest to southwest and at the center of the study area (Figure 4). The decline in SOC concentrations might be attributed to cultivation practices (e.g., tillage practices) which aggravate soil erosion and oxidation of organic matter. This resulted in water and soil loss because of inadequate protection of tillage management practices [13]. However, the highest SOC concentrations were distributed mainly between the northwest and northeast of the subwatershed (Figure 4). Moreover, the lowest SOCS was distributed in the north and southwestern parts than the north and southeastern parts of the subwatershed, whereas the highest SOCSS were found in the northeastern part of the study area (Figure 4) due to litter falls from forest vegetation. Generally, the differences in SOC and SOCS spatial patterns might be due to the difference in land use patterns in line with topographic variations. A study by [62] verified the influence of topographic patterns on the SOC dynamics.

4. Conclusions

The coefficient of variation (CV) from descriptive statistical analysis revealed a moderate variability for SOC and SOCS in terms of soil depth. The variation of the semivariogram of the SOC was best described by the exponential and circular models, while SOCS by the spherical and exponential models. The semivariogram of SOC and SOCS showed a strong spatial dependence mainly due to the influence of intrinsic factors in the study area. Generally, spatial distribution of SOC and SOCS was relatively higher in the northwestern and eastern parts of the subwatwershed. Finally, the kriged maps of the current study showed the spatial variation of SOC, which improves site-specific soil management schemes.

Data Availability

The data that support the findings of this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

The authors would like to thank the Ministry of Education (MOE) for financial support and Haramaya University for hosting the program. Moreover, a corresponding author is thankful to Wachemo University for providing a chance to pursue PhD study. The authors are also grateful to the soil testing laboratory staff of Agriculture and Natural Resource Development Bureau of Southern Ethiopia for their cooperative work during soil analysis.