Terrain Effects on the Spatial Variability of Soil Physical and Chemical Properties

Understanding topography effects on soil properties is vital to modelling landscape hydrology and establishing sustainable on-field management practices. This research focuses on an arable area (117 km2) in Southwestern Ethiopia where agricultural fields and bush cover are the dominant land uses. We postulate that adapting either of the soil data resources, coarse resolution FAO-UNESCO (Food and Agriculture Organization of the United Nations Educational, Scientific and Cultural Organization) soil data or pedo-transfer functions (PTFs) is not reliable to indicate future watershed management directions. The FAO-UNESCO data does not account for scale issues and assigns the same soil property at different landscape gradients. The PTFs, on the other hand, do not account for environmental effects and fail to provide all the required data. In this regard, mapping soil property spatial dynamics can help understand landscape physicochemical processes and corresponding land use changes. For this purpose, soil samples were collected across the watershed following a gridded sampling scheme. In areas with heterogeneous topography, soil is spatially variable as influenced by land use and slope. To understand the spatial variation, this research develops indicators, such as topographic index, soil topographic wetness index, elevation, aspect, and slope. Pearson correlation (r), among others, was used to investigate terrain effects on selected soil properties: organic matter (OM), available water content (AWC), sand content (%), clay content (%), silt content (%), electrical conductivity (EC), moist bulk density (MBD), and saturated hydraulic conductivity (Ksat). The results show that there were statistically significant correlations between elevation-based variables and soil physical properties. Among the variables considered, the ‘r’ value between topographic index and soil attributes (i.e., OM, EC, AWC, sand, clay, silt, and Ksat) were 0.66, 0.5, 0.7, 0.55, 0.62, 0.4, and 0.66, respectively. In conclusion, while understanding topography effects on soil properties is vital, implementing either FAO-UNESCO or PTFs soil data do not provide appropriate information pertaining to scale issues.


Introduction
Soils vary widely as a function of their position on the landscape [1][2][3] and agricultural management, land use, and cultivation intensity [4,5]. As a result, identification of landscape features is an important unreliable PTfs. Instead, we developed soil-landscape relation for a tropical catchment in the Rift Valley basin where data are a scarce resource.
Here, two points need to be balanced. Firstly, the indispensability of knowledge of soil properties for planning territorial development and environmental conservation. Secondly, implementation of methodologies that rapidly and effectively capture information about the spatial variability of soils in a way that reduces the need for intensive and expensive sampling [30]. In this regard, the geostatistical approach, which uses probabilistic methods, has been successfully used in soil science for quantitative description of the spatial variability of soil physical and chemical properties [31]. Since topography parameters are known to control water and sediments distribution over the landscape, features such as altitude, slope, and shape are correlated with soil physical properties [3].
To understand the spatial variability of soil and land use in a rugged, mountainous environment of the Rift Valley region, this research extracted indicators from digital elevation models, such as a topographic index, soil topographic wetness index, elevation, aspect, and slope. These indicators were then used to identify the spatial characteristics of soils and analyze the spatial variability as it is affected by climate, soil property, and topography. The use of the digital elevation model increased the predictive accuracy of soil parameters from terrain attributes.
In this paper, we measured the spatial variability of soils across the watershed and evaluated the use of topographic indices to develop explanatory variables. Specific objectives were (1) developing a fine scale soil map of the watershed, (2) spatially correlate soil properties and topographic attributes, and (3) assess the interactive effect of soil, land use type, and slope position on soil physiochemical properties. To improve soil property mapping, we used auxiliary variables that help understand the processes occurring over the landscape and improve the precision of soil mapping.

Site Description
The study area is a 117 km 2  The climate in the area is mainly controlled by seasonal migration of the Inter-Tropical Convergence Zone and its associated atmospheric circulation [32]. The mean annual rainfall lies between 1280 and 1339 mm and the mean annual temperature in Hosanna station varies from 11 to 22 • C in the watershed.
Based on the 2018 Landsat based land cover classification, the dominant land use in the area were agricultural (48.4%) and bushland units (45.2%). Grazing land, mixed forest and water body covered the remaining 6.4%. Landscape datasets employed in the study included spatial soil attributes and elevation data with 30 m horizontal and 1 m vertical resolution [33].

Soil Survey and Mapping
A combination of field and laboratory work intended to identify the basic physicochemical properties of soils, establish the distribution of those soils at specific map scales, and interpret the information for a variety of uses. Soil samples were collected in a grid scheme field observation (approximately 1 sample for each 1 km 2 grid) that uses soil sample pits and description augers (see section 3.1). The sampling sites were selected on an arable area (117 km 2 ) in Southwestern Ethiopia, where agricultural fields and bush cover are the dominant land uses. Projected coordinates (UTM) for each soil sample were recorded to help soil property spatial characterization.
The field soil survey followed the latest United States Department of Agriculture (USDA)-natural resource conservation service norms and standards for field mapping [34]. Soil sampling allowed the estimation and characterization of soil physical and chemical attributes at unsampled sites through existing models. Because of the high complexity relationships among soil properties and the unpredictable spatial pattern [35], the deterministic method of soil characterization does not result in accurate estimation. Probabilistic methods, on the other hand, admit some uncertainty about how soil properties could vary in space, and soil properties at the unmeasured sites are considered to be outcomes of some random process [36].
The sub-horizontal structure [37] indicated the existence of varying soil landscape relationships and that the vertical position and horizontal location of each layer could be derived from the elevation

Soil Survey and Mapping
A combination of field and laboratory work intended to identify the basic physicochemical properties of soils, establish the distribution of those soils at specific map scales, and interpret the information for a variety of uses. Soil samples were collected in a grid scheme field observation (approximately 1 sample for each 1 km 2 grid) that uses soil sample pits and description augers (see Section 3.1). The sampling sites were selected on an arable area (117 km 2 ) in Southwestern Ethiopia, where agricultural fields and bush cover are the dominant land uses. Projected coordinates (UTM) for each soil sample were recorded to help soil property spatial characterization.
The field soil survey followed the latest United States Department of Agriculture (USDA)-natural resource conservation service norms and standards for field mapping [34]. Soil sampling allowed the estimation and characterization of soil physical and chemical attributes at unsampled sites through existing models. Because of the high complexity relationships among soil properties and the unpredictable spatial pattern [35], the deterministic method of soil characterization does not result in accurate estimation. Probabilistic methods, on the other hand, admit some uncertainty about how soil properties could vary in space, and soil properties at the unmeasured sites are considered to be outcomes of some random process [36].
The sub-horizontal structure [37] indicated the existence of varying soil landscape relationships and that the vertical position and horizontal location of each layer could be derived from the elevation above sea level using the digital elevation model. The surface was built using Kriging point iteration, Soil Syst. 2020, 4, 1 5 of 21 a model suitable to analyze datasets with a spatial autocorrelation. Kriging estimates the statistical relationships among the sample points with the assumption that the distance between sampling points reflects a spatial correlation that can be used to explain data variations [38]. The algorithm implemented for the Kriging model first creates an empirical semivariogram representing the variance for each pair of observation points followed by adjustment of the model to fit the semivariogram (see Section 2.3). The methodology adopted in this study is presented in the following flowchart ( Figure 2). Soil Syst. 2020, 4, 1 5 of 20 above sea level using the digital elevation model. The surface was built using Kriging point iteration, a model suitable to analyze datasets with a spatial autocorrelation. Kriging estimates the statistical relationships among the sample points with the assumption that the distance between sampling points reflects a spatial correlation that can be used to explain data variations [38]. The algorithm implemented for the Kriging model first creates an empirical semivariogram representing the variance for each pair of observation points followed by adjustment of the model to fit the semivariogram (see Section 2.3).
The methodology adopted in this study is presented in the following flowchart ( Figure 2).

Spatial Mapping of Topography Index
The standard topographic index (TI) within the yopographic hydrologic model (TOPMODEL) (Equation (1)) uses the catchment digital elevation model. The soil-topographic wetness index analysis (Equation (2)), on the other hand, extends the purely topography-based TOPMODEL analysis by accounting for spatial variation in hydrologically relevant soil properties [39]. TOPMODEL concepts are generalized in a semi-physical representation of rainfall-runoff dynamics, which is parameterized with catchment topographic data and initial conditions of soil moisture. The topographic index components in the model indicate that as the subsurface contributing area increases, there is a likelihood increase in the relative saturation. In addition, the TI incorporates Darcy's law to estimate subsurface flow proportionality with the hydraulic gradient approximated by the tangent of the ground-surface slope (tan β). A steep topographic slope makes a greater hydraulic gradient and increases drainage potential, which, in turn, reduces the level of saturation. The topographic index (TI) within TOPMODEL is expressed as

Spatial Mapping of Topography Index
The standard topographic index (TI) within the yopographic hydrologic model (TOPMODEL) (Equation (1)) uses the catchment digital elevation model. The soil-topographic wetness index analysis (Equation (2)), on the other hand, extends the purely topography-based TOPMODEL analysis by accounting for spatial variation in hydrologically relevant soil properties [39]. TOPMODEL concepts are generalized in a semi-physical representation of rainfall-runoff dynamics, which is parameterized with catchment topographic data and initial conditions of soil moisture. The topographic index components in the model indicate that as the subsurface contributing area increases, there is a likelihood increase in the relative saturation. In addition, the TI incorporates Darcy's law to estimate subsurface flow proportionality with the hydraulic gradient approximated by the tangent of the ground-surface slope (tan β). A steep topographic slope makes a greater hydraulic gradient and increases drainage potential, which, in turn, reduces the level of saturation. The topographic index (TI) within TOPMODEL is expressed as where α is the upslope contributing area for the cell per unit of contour line (m), tan β is the topographic slope of the cell, and T is the transmissivity (soil depth × saturated soil hydraulic conductivity) of the uppermost layer of soil (m 2 d −1 ) [40,41].

Geostatistical Analysis
A geostatistical analysis was conducted to elucidate the spatial patterns of soil variability within the study area and to allow the interpolation of soil variables to unmeasured locations. Unsampled soil properties across the entire spatial domain were produced using the variogram model in kriging. Variogram modeling is an advanced geostatistical procedure that generates an estimated surface from a scattered set of points with z-values. Prediction maps were created using a set of training data to employ Kriging point iteration. Kriging works with the regionalized variable theory [42,43], assuming that the spatial variation in the phenomenon represented by the z-values is statistically homogeneous. The spatial variation is quantified using semivariogram [44] computed from soil samples collected from field and laboratory datasets across the watershed following a gridded sampling scheme ( Figure 3). Omnidirectional experimental semi-variograms that quantify the dissimilarity γ(h) between observations as a function of the separation distance h were computed in ArcGIS 10.6.1 as half the average squared difference between the components of every data pair. The sample semivariogram for a separation distance of h, lag, is the average squared difference in z-value between pairs of input sample points separated by h. Given n measurements of a spatial attribute (as in Equation (3)), the method of moments estimate of the semi-variogram is calculated as

Soil Wetness and Dryness with Topography
In order to assess the spatial variability of terrain elevation with soil wetness/dryness, the soil topographic wetness index, a map, Figure S1E, was developed from the coincidence of TI ( Figure S1B) and soil type ( Figure S1D). Reconceptualization of the TI map into the soil topographic wetness index and topographic wetness index provides accurate information on how land use changes with terrain and soil wetness/dryness at the lowest functional units. The watershed topographic wetness index was Sill variance depicting the total variance of the process is the maximum value reached by the variogram after an initial increase [26]. The range in geostatistical analysis is the distance at which the variogram reaches the sill, beyond which the process is no longer spatially correlated. The combination of measurement error and variation over distances less than the shortest sampling interval is explained by the nugget variance at lag distance 0 [45]. High variations in nugget indicate high, unquantified, fine-scale variability and measurement error. On the other hand, the spatial heterogeneity of soil property dynamics is explained by the sill value [46]. During data analysis, the distribution of soil properties was assessed using an exploratory data analysis employing the Kolmogorov-Smirnov test. Table 1 shows the ranking of qualifiers [34] based on the Food and agricultural organization reference soil group guideline [37]. Following the field soil survey, a laboratory analysis and descriptive augers, seven soil taxa at the family level of USDA soil taxonomy hierarchy were identified and mapped. These soil types were chromic luvisol, chromic vertisol, haplic gleysol, humic nitisol, lithic leptosol, ochric regosol, and pellic vertisols ( Figure 3). Nearly 92% of the catchment was covered with three soil types: chromic vertisol (33.3%), chromic luvisol (30.2%) and pellic vertisol (28.2%). The remaining soil classes (humic nitisol, haplic gleysol, lithic leptosol, and ochric regosol) covered 8% of the area.

Soil Wetness and Dryness with Topography
In order to assess the spatial variability of terrain elevation with soil wetness/dryness, the soil topographic wetness index, a map, Figure S1E, was developed from the coincidence of TI ( Figure S1B) and soil type ( Figure S1D). Reconceptualization of the TI map into the soil topographic wetness index and topographic wetness index provides accurate information on how land use changes with terrain and soil wetness/dryness at the lowest functional units. The watershed topographic wetness index was lumped into 10 equal area intervals ranging from 1 to 10, with index class 1 covering 10% of the watershed area with the lowest topographic wetness index (i.e., lowest propensity to saturate) and index class 10 containing 10% of the watershed with the highest topographic wetness index (i.e., highest propensity to saturate). Pertinent soil information was extracted with an aerially weighted soil topographic wetness index and included in the index input table by soil wetness class. The results show that the saturated hydraulic conductivity (Ksat) of soils has an unpredictable pattern with elevation. Note: 1*, 1, 2, 3, 4, 5 suggestions for ranking qualifiers in soil unit names for the most common qualifiers. 6 names of reference soil groups can further be explained using prefixes added with formative element qualifiers. a influential and shaping qualifier elements for naming reference soil groups units. b referring to diagnostic horizons, properties or soil materials. c directly related diagnostic horizon, property or soil material. d secondary characteristics qualifiers referring to soil characteristics that indicates a particular water regime, soil solution or groundwater chemistry, exchange complex specification, weak soil development or morphological characteristics than diagnostic horizons, properties or soil materials. e having the typical expression of the reference soil group in the sense that there is no further or meaningful characterisation. f indicates depth of occurrence or degree of expression of soil characteristics. RSGs indicate reference soil groups.

Soil Topography Relationships
For each soil sample, texture was defined using sand and clay content obtained from a laboratory analysis and using soil-plant-air-water (SPAW) soil water characteristics' pedo-transfer functions. During the textural classification, SPAW soil water characteristic equations [47,48] are valid within a range of approximately 0%-60% clay content and 0%-95% sand content. In addition, the spatial relationships between soil properties and topographic attributes for the data points were assessed for the first and second soil layers. In each sampling point, average topographic height, slope and aspect were calculated from the digital elevation model in ArcGIS environment. The spatial dependence of each individual variable as well as the relationship between them was evaluated with Pearson correlation and descriptive statistics of soil topography variables.
The descriptive statistical results for altitude and soil physical attributes are presented in Tables 2  and 3. Samples were collected along the slope with altitude values ranging from 2065 to 2950 m above sea level. Frequency distribution was evaluated for these attributes and the results are summarized in Tables 2 and 3. According to Webster [35], skewness values up to 0.5 suggest a specific attribute with normal distribution. Besides altimetry, the only attributes that showed both features, a skewness smaller than 0.5 and a normal distribution were sand, clay, OM, and soil bulk density. Negative skewness values ( Table 3) indicate data that are skewed to the left, which means that the left tail is long relative to the right tail. On the other hand, positive skewness values indicate data that are skewed right. The kurtosis shows how the tails of a distribution differ from the normal distribution. Negative kurtosis (Tables 2 and 3) indicates a distribution flatter than a normal curve and tails that are heavier than a normal distribution are indicated by positive kurtosis values. The small standard deviation values (<30%) were found for most of the variables, other than elevation and aspect. The high values of coefficient of variation for elevation can be explained by the great amplitude of variation in the area (minimum and maximum values).
The summary statistics for the second soil layer indicate that topographic and soil properties are highly variable between sampling points in the watershsed.
The spatial dependence of each individual variable was evaluated with semivariogram and cross semivariograms (Figure 4). The unsampled estimates of variables across the entire spatial domain were produced using the variogram model in kriging. The spatial variation was quantified by the semivariogram estimated by using the sample semivariogram computed from point input dataset. Figure 4 summarizes the experimental semi-variogram plots and how to relate the practical semi-variogram to the ideal models. It can be seen that whilst some variables may follow simple behaviour (e.g., elevation), many others, such as sand, clay, silt, moist bulk density, available water content, organic matter, electrical conductivity, and saturated hydraulic conductivity, may require a non-linear or relatively complex models to describe the experimental semi-variogram.   There are several important features to note from the sample semivariogram plots. For some variables, at relatively short lag distances of h, the semivariance is small but increases with the distance between the pairs of sample points. At a distance referred to as the range, the semivariance levels off to a relatively constant value referred to as the sill. This implies that beyond this range distance, the variation in z-values is no longer spatially correlated. Within the range, the z-value variation is smaller when the pairs of sample points are closer together. The extent of the x-axis of the semivariogram is determined by the distance between the most widely separated pair of points in the input sample data.
The results show that there was a spatial correlation between slope, sand, silt, clay, AWC, EC, OM, and MBD. A positive correlation was characterized between soil depth and AWC, EC and sand content, OM and EC, AWC and OM, and Ksat and soil depth ( Table 4). For the second soil layer (Table 5), a positive correlation was characterized between TI and silt content, soil depth and EC, sand content and depth, clay content and TI, sand content and OM, EC and soil depth, and clay content and MBD. In addition, AWC and MBD, slope and elevation, soil depth and terrain elevation, and AWC and clay content were positively correlated. A low positive and negative relation was reported between the other variables.
In summary, Pearson product moment correlation (Tables 4 and 5) better estimates the linear relationship between two variables. In order to assess the linearity and stationarity of the dataset, Pearson correlation was also calculated for the log-transformed data (Tables 4 and 5). We also believe that any two variables may not linearly relate to each other to implement Pearson correlation. As a result, the monotonic relationship between two variables was assessed using Spearman's rank correlation coefficient (or Spearman's rho) and Kendall's rank correlation coefficient (or Kendall's tau) [49]. For a detailed presentation of a statistical test comparison between Pearson correlation, Spearman's rho, and Kendall's tau, please see the Excel Supplementary File (Tables S1 and S2).
Pearson correlation being the most used statistical tool to evaluate the relationship between two variables, the coefficient of correlation is also used in geostatistics to improve data interpolation. The correlation between the soil-topography index and selected soil physicochemical properties for top and second layers is presented in Figure 5.
The graph relating topographic index (TI) with soil physicochemical properties indicated that soil organic matter (OM), moist bulk density (MBD), slope, depth, and saturated hydraulic conductivity (K sat ) shows a decrease in TI for both soil layers. However, the sand content in the upper and second soil layers shows a relatively different trend of a tendency for a slight increase for the first soil layer, while a noticeable decrease is observed for the second soil layer. On the other hand, silt and clay contents, OM, OC, and available water content (AWC) show an increase with TI.
In the second soil layer, the increase in the TI resulted in an increase in clay and AWC with a relatively stronger correlation coefficient of 0.6 and 0.66, respectively ( Figure 5, lower panel). On the contrary, the influence of TI on the silt content was relatively weaker (r = 0.36). Linear regression also indicated that the correlation between Ksat and MBD with TI is 0.66 and 0.37, respectively.     Note: TI, OM, EC, MBD, AWC, and Ksat represent topographic index, soil organic matter, electrical conductivity, moist bulk density, and saturated hydraulic conductivity, respectively. Pearson correlation being the most used statistical tool to evaluate the relationship between two variables, the coefficient of correlation is also used in geostatistics to improve data interpolation. The correlation between the soil-topography index and selected soil physicochemical properties for top and second layers is presented in Figure 5. The graph relating topographic index (TI) with soil physicochemical properties indicated that soil organic matter (OM), moist bulk density (MBD), slope, depth, and saturated hydraulic conductivity (Ksat) shows a decrease in TI for both soil layers. However, the sand content in the upper and second soil layers shows a relatively different trend of a tendency for a slight increase for the first soil layer, while a noticeable decrease is observed for the second soil layer. On the other hand, silt and clay contents, OM, OC, and available water content (AWC) show an increase with TI.
In the second soil layer, the increase in the TI resulted in an increase in clay and AWC with a relatively stronger correlation coefficient of 0.6 and 0.66, respectively ( Figure 5, lower panel). On the contrary, the influence of TI on the silt content was relatively weaker (r = 0.36). Linear regression also indicated that the correlation between Ksat and MBD with TI is 0.66 and 0.37, respectively.

Soil, Topography and Land Use Relationships
Soil-land use-slope relationships and the spatial variability of land use units can be collated based on six snapshots (1973-2018) for five slope classes. For example, in 1973 ( Figure 6, Table S3), 141 ha of the agricultural area was topographically located in slope class 2 (slope ranging from 2 to 10 • ) and the specific soil type was chromic vertisol.

Soil, Topography and Land Use Relationships
Soil-land use-slope relationships and the spatial variability of land use units can be collated based on six snapshots (1973-2018) for five slope classes. For example, in 1973 ( Figure 6, Table S3), 141 ha of the agricultural area was topographically located in slope class 2 (slope ranging from 2 to 10°) and the specific soil type was chromic vertisol.  Table S3 and Figure S2.
As an important indicator of topography, terrain slope has a substantial effect on land use, thereby affecting soil fertility and the density and pattern of land use. It was also indicated that time series land use change pattern changed in agriculturally accessible and arable soils (Table S3), while showing a relative stability in mountainous and poorly fertile soils. According to World Food and Agricultural Organization study [37], soils such as vertisols, luvisols, nitisols, and gleysols are fertile soils. On the other hand, leptosols and regosols are weakly developed mineral soils that are extensive in eroding lands, arid and semi-arid areas and mountain regions. Intensification of agricultural practices leading to a decrease in soil fertility could also be related to the acute shortage of land on the gentle slope.
Studies [50,51] highlighted that in line with agricultural unsuitability, properties of dryness/wetness and magnitude and seasonal pattern of rainfall alter the pattern land use change and soil fertility. This suggests that management for any one of these soil properties may yield unintended cascading effects throughout the soil subsystem.

Slope Effects on Soil Physicochemical Properties and Land Cover Distribution
Slope has a significant effect on the intensity of soil erosion, which can affect soil fertility and patterns of land use. Balba [52] indicated that steep slopes are typified by rapid runoff, less moisture entering the soil, and lower crop productivity. To evaluate the impact of slope position on soil texture (i.e., sand, clay, and silt), OM, and AWC, the spatial variation of each variable is presented in Table  S4 (Excel Supplementary File). Figure 7 illustrates the spatial relation between catchment slope and selected soil physicochemical properties. The change in OM and AWC is highly related to the type of soil and land cover density. For instance, a highly leached leptosol is characterized by low OM.  Table S3 and Figure S2.
As an important indicator of topography, terrain slope has a substantial effect on land use, thereby affecting soil fertility and the density and pattern of land use. It was also indicated that time series land use change pattern changed in agriculturally accessible and arable soils (Table S3), while showing a relative stability in mountainous and poorly fertile soils. According to World Food and Agricultural Organization study [37], soils such as vertisols, luvisols, nitisols, and gleysols are fertile soils. On the other hand, leptosols and regosols are weakly developed mineral soils that are extensive in eroding lands, arid and semi-arid areas and mountain regions. Intensification of agricultural practices leading to a decrease in soil fertility could also be related to the acute shortage of land on the gentle slope.
Studies [50,51] highlighted that in line with agricultural unsuitability, properties of dryness/wetness and magnitude and seasonal pattern of rainfall alter the pattern land use change and soil fertility. This suggests that management for any one of these soil properties may yield unintended cascading effects throughout the soil subsystem.

Slope Effects on Soil Physicochemical Properties and Land Cover Distribution
Slope has a significant effect on the intensity of soil erosion, which can affect soil fertility and patterns of land use. Balba [52] indicated that steep slopes are typified by rapid runoff, less moisture entering the soil, and lower crop productivity. To evaluate the impact of slope position on soil texture (i.e., sand, clay, and silt), OM, and AWC, the spatial variation of each variable is presented in Table S4 (Excel Supplementary File). Figure 7 illustrates the spatial relation between catchment slope and selected soil physicochemical properties. The change in OM and AWC is highly related to the type of soil and land cover density. For instance, a highly leached leptosol is characterized by low OM. Other similar studies conducted in the Ethiopian highlands [53] indicated that variations in soil physicochemical properties are related to slope position and land use. According to their findings, low Other similar studies conducted in the Ethiopian highlands [53] indicated that variations in soil physicochemical properties are related to slope position and land use. According to their findings, low soil pH was observed at forested areas in the upslope position, which was characterized by high OM. They also noted that high pH values were recorded on lower slope positions where low OM cultivation land was the dominant land use. Griffiths et al. [54], highlighted that topography could influence soil physicochemical properties, such as OM, soil depth, and texture. Our findings agree with Hu et al. [55], who showed that a decrease in OM and the associated changes in soil nutrients were related to land use conversions into agriculture and increased incidence of tillage.
Details of topography effects on soil physicochemical properties and land cover distribution are presented in Table S4 (excel supplementary file). It is shown that land use was significantly affected by slope position, soil type, and soil fertility. For illustration, in 2018, 35.5 ha of the agricultural area was located in slope class 2 (slope ranging from 2 to 10 • ) in a clayey soil type 2 (i.e., abundance between 34.7% and 38.1%).

Conclusions
In areas with heterogeneous topography, soil is spatially variable and the way we manage catchment soil and land use alters the fluxes of water and sediment. This effect is more pronounced in highly rugged agrarian landscapes where there is a significant land disturbance and topographic influences are dominant. In this regard, knowledge of soil properties associated with topographical attributes is vital for modelling soil-landscape relationships and establishing sustainable on-field management practices. To this end, previous studies on the spatial variability of soil properties have either been site-specific, coarse in scale, or implemented pedo-transfer-functions that neither accounted for environmental effects nor provided all the required data inputs.
In this study, we assessed soil-landscape relation and mapped soil property spatial dynamics to understand landscape physicochemical processes and used the result for hydrological modelling and future soil and watershed management applications. The significance of the work therefore stems from the fact that we did not rely on our research using either coarse resolution FAO-soils that assigns the same soil property at different landscape gradient or site-specific and unreliable PTfs. Instead, we developed a soil-landscape relation for a tropical catchment in the Rift Valley basin where data are a scarce resource.
The relationship between topography and the spatial variability of soil physicochemical properties was evaluated using Kriging predictions and it is shown that the use of Kriging is a promising tool in soil physics. Descriptive statistics were used to evaluate central tendency and to explore the correlation between spatially explicit soil and terrain attributes. In an study aiming to relate practical semivariogram to the ideal models, it is seen that whilst some variables may follow simple behaviour (e.g., elevation), many others, such as sand, clay, silt, moist bulk density, available water content, organic matter, electrical conductivity, and saturated hydraulic conductivity may require relatively complex models to describe the experimental semi-variogram.
Supplementary Materials: The following are available online at http://www.mdpi.com/2571-8789/4/1/1/s1, Table S1 (excel supplementary file). Pearson, Spearman's rho, and Kendall's rank correlation coefficient for top soil layer. Table S2 (excel supplementary file). Statistical test comparison between Pearson correlation, Spearman's rho, and Kendall's tau for second soil layer. Table S3. Land use distribution (area in ha) with elevation and soils in the Batena Watershed. Note: the acronym for the soil types, slope range and land use are explained. Table S4 (excel  supplementary file). Terrain effects on soil physicochemical properties and land cover distribution. A note for the value range of each of the five soil physicochemical properties, slope range ( • ), and time series land cover change is included. Figure S1. Terrain elevation (A), topographic index classes (B), small enlarged view of topographic index (C), catchment soil map (D), soil topography wetness index (E), and comparative small enlarged view of soil topography wetness index (F) with topographic index (C). Figure