Spatial Variability of Soil Aggregate Stability in a Disturbed River Watershed

1. Introduction

Spatial distribution of soil properties, techniques such as classical statistics and Geostatitics have been widely applied (Webster and .Oliver, 1990, Wang et al 2009). Geostatitics provides the basis for the interpolation and interpretation of the spatial variability of soil properties (Webster, 1985, Pohlmann, 1993, Cambardella et al, 1994). Information on the spatial variability of soil properties leads to better management decisions aimed at correcting problems and at least maintaining productivity and sustainability of the soils and thus increasing the precision of farming practices (Schimel et al, 2000, Ozg¨oz, 2009)).A better understanding of the spatial variability of soil properties would enable refining agricultural management practices by identifying sites where remediation and management are needed. This promotes sustainable soil and land use and also provides a valuable base against which subsequent future studies can be proposed (Cambardella et al, 1994)

Aggregate stability has been studied extensively since the early 19th century (Yoder, 1936). This sustained interest shows that aggregate stability plays a central role in the behavior of soils and that its role is still not clearly understood. However opportunities arising from spatial dependence of soil properties can be exploited to explain the role of aggregate stability in the behavior of soils. So far a few studies have shown that spatial dependence of soil properties influence weed distribution (Gaston et al., 2001), crop yield (Cassel et al., 2000) and can explain the spatially variable safe sites concept (Zanin et al., 1998).

Several authors have studied the spatial distribution of soil properties like particle size, organic matter, pH and soil microbial activity at different scales ranging from a few meters (Gajem et al., 1981) to several kilometers (Ovalles and Collins, 1988). Among the various soil properties hydraulic conductivity (Ks) is reported to have the highest statistical variability (Bigger and Nielsen, 1976). It is evident now that spatial variability of soil properties are scale dependent especially those related to water transport in the soils (Iqbal et al., 2005).

Indeed within the last two decades spatial variability has been studied using technologies such as remote sensing (Viscarra- Rossell and McBratney, 1997) and geographical information systems (Mulla and Schepers, 1997). These technologies have been developed to help land users better manage land resources utilizing spatially varying prescription intensities (Stafford, 1997; Cassel et al., 2000).Soil management practices associated with various land uses may modify the soil behavior for example through additions or withdrawal of organic matter and compaction thereby profoundly altering spatial variability of aggregate stability. Stutter et al., (2004) studied the spatial variability and variance structure of cation exchange chemistry in a granitic, hither moorland site in Northwest Scotland. Their results showed strongly significant short-range vertical and lateral variability. Soil properties vary spatially from a field to a larger regional scale affected by both intrinsic (soil forming factors) and extrinsic factors (soil management practices, fertilization, and crop rotation) (Cambardella and Karlen, 1999). The variation is a gradual change in soil properties as a function of landforms, geomorphic elements, soil forming factors and soil management (Buol et al., 1997). The variation of soil properties should be monitored and quantified to understand the effects of land use and management systems on soils. Geostatistical methods like kriging have been used successfully for predicting spatial variability of soil properties (Trangmar et al.1985; Gaston et al., 2001; Cambardella et al., 1994; Saldana et al., 1998; Zebarth et al., 2002; Lark, 2002; Dercon et al., 2003).

Kriging uses a semivariogram which is a function of the distance and direction separating two locations—to quantify the spatial dependence in the data. A semivariogram is constructed by calculating half the average squared difference of the values of all the pairs of measurements at locations separated by a given distance h. The semivariogram is plotted on the y axis against the separation distance h. This semivariogram model is then used to define the weights that determine the contribution of each observed data point to the prediction of new values at unsampled locations. There are some statistical assumptions behind kriging. The main assumption is stationarity (spatial homogeneity). If data is stationary, the data mean and the semivariogram are the same at all locations in the data extent. If this assumption is held, just a few kriging model parameters have to be estimated from the data. Indeed it can be concluded from the literature that most soil properties have short-range spatial structure extending only tens of meters and that long-range spatial structure is doubtful (Campbell, 1978; Trangmar et al., 1986; Gaston et al., 2001; Iqbal et al., 2005).

The study sought to establish that land use and Landscape interactions had influence on the spatial variability of aggregate stability within the watershed.

2.0 Materials and Methods

2.1 Location and Environmental Conditions of the Study Area

The Middle river Njoro catchment is located in Nakuru County within the Kenyan part of the East African Rift System. It lies between longitudes 35˚ 05’E and 36˚ 05’E, and latitudes 0˚ 15’S and 0˚ 25’S.The catchment covers about 8,200 ha and is located about 200km North West of Nairobi in Nakuru County. The catchments lie between altitudes 1720m above sea level (asl) and 3000m asl (Figure 1). The sub watershed has slopes in the range of less than 2% to more than 30 %.( Figure 2). The river has its source from Mau escarpment and drains into Lake Nakuru (Mainuri and Owino, 2014).

Figure 1: Study Area (Source: Mainuri and Owino, 2014)

Figure 2: Slope Map

2.2 Methods

2.2.1: Soil Mapping

The study area was classified into six (6) different landscape units according to physiographic features namely: mountains, plateaus, uplands, valleys, peneplains, alluvial plains, lacustrine plains and hills (Figure 3).

Figure 3: Physiographic Map of the Middle River Njoro Watershed

The selection of profile pits for sampling was based on the physiographic units and positions for the pits selected at random. Soil classification was done according to FAO/UNESCO (1997).

2.2.2 Determination of Aggregate Stability

Soil aggregate analysis was done following the wet sieving technique of Yoder 1936 with the help of a nest of six sieves (5, 2, 1, 0.5, 0.25, and 0.1 mm). Soil aggregates (>5.0–8.0 mm) were spread on the top sieve and soaked with salt-free water twice at the interval of 5 min. The nest of sieves were then submerged under water and oscillated at 35 cycles/minutes. The aggregates on each sieve was collected, oven-dried, and weighed. The aggregates were treated with H2O2 and HCl and passed through the same sieve and the weights of non aggregated primary particles were taken after oven drying. Soil aggregate stability was expressed in terms of mean weight diameter (MWD) in millimeters, defined as the sum of the product of the mean diameter, Xi, and the total sample weight, Wi, of each size fraction (MWD = Xi Wi).

2.2.3 Mapping Aggregate Stability

Analysis of spatial variability of soil aggregate stability within the Middle River Njoro watershed was done using two geostatistical methods: (i) investigating spatial autocorrelation in the datasets using a semivariogram/covariance cloud (Cressie, 1985; Cressie, 1988), and (ii) prediction of attribute values at unsampled locations using a generalized linear regression technique called Empirical Bayesian kriging. The structure of spatial variance between aggregate stability observations was derived from the sample semivariogram

2.2.4 The Semivariogram

A spatial layer of the soil sampling sites with the respective soil aggregate stability was composed using ArcGIS software. The spatial sampling site layer was input as the source of data and the soil aggregate stability as the data field. To explore and quantify the spatial autocorrelation of the soil sampling sites, geostatistical wizard functionality within ArcGIS was used (Isaaks et. al. 1989; Journel and Huijbregts, 1978; Stein, 1999). Empirical Bayesian Kriging linear regression technique was used to generate a soil aggregate stability prediction map

3.0 Results and Discussions

3.1: Soil Mapping

The soils in the catchments were distinguished on the basis of Physiographic, parent material/ geology and soil characteristics. Six major physiographic units, mountains, hills, plateaus, uplands, plains and valleys (Figure 3) were identified in the area. Mountains (Mo) and major scarps have slopes greater than 30% (Forests).The soils of this geomorphic unit are developed on volcanic ashes and other pyroclastic rocks of recent volcanoes with stable to very stable structures. Hills (Hie) and minor scarps (Forests) have slopes which are greater than 16% with stable structures. The soils of this geomorphic unit are developed on ashes and other pyroclastic rocks of recent volcanoes. The plateaus (PU) (Agriculture and grasslands) are undulating with slopes of between 5 – 8%. The soils are developed on ashes and other pyroclastic rocks of recent volcanoes with moderately stable structure. The uplands (UP) (agriculture, grasslands and minor forests) are undulating with slopes of between 5 – 8%. The soils which are moderately stable are developed on volcanic ashes and other pyroclastic rocks of recent volcanoes.

Figure 3: Soil Mapping Units

The plains (PL) (grasslands and agriculture) are flat to very gently undulating in relief with moderately stable structures and slopes of between 0 – 2%. The soils of this unit are well drained, deep to very deep. Soils of Valleys (VA) (grasslands, agriculture and wetlands) are well drained to imperfectly drained, moderately deep to deep and unstable. As observed from the results, the landscapes where most agricultural activities take place are dominated by soils that have structures that range from moderately stable to unstable making them quite susceptible to soil erosion.

3.2 Mapping Aggregate Stability Using Empirical Bayesian Kriging (EBK)

The selection of sampling sites was based on the physiographic units and sampling points for aggregate stability determination selected at random Figure 4.

Figure 4: Sampling Sites for Aggregate Stability showing Aggregate stability values

The Empirical Bayesian Kriging (EBK) was used to build a valid kriging model due to its ability to interpolate and automate the most difficult aspects of building a valid model (Figure 5). Other kriging methods in Geostatitical Analysis require somebody to manually adjust parameters to receive accurate results, but EBK automatically calculates these parameters through a process of sub setting and simulations. Empirical Bayesian kriging also differs from other kriging methods because it accounts for the error introduced by estimating the underlying semivariogram. By not taking the uncertainty of semivariogram estimation into account, other kriging methods underestimate the standard errors of prediction. EBK differs from other kriging methods in Geostatitical Analyst by using an intrinsic random function as the kriging model. The variability of aggregate stability exhibited spatial dependence that was described using semivariogram models. The degree of spatial dependence for each variable was determined with geostatistical methods using semivariogram analysis and kriging.

Figure 5: The spectrum of the semivariogram models produced by EBK.

The semivariogram for aggregate stability analysis consisted of three basic parameters which described the spatial structures γ(h) = Co + C. Co represented the nugget effect, which was the local variation occurring at scales finer than the sampling interval, such as sampling error, fine-scale spatial variability, and measurement error; Co + C is the sill (total variance); and the distance at which semivariogram levels off at the sill which was called the range (beyond that distance the sampling variables were not correlated). The variability of aggregate stability exhibited spatial dependence that was described using acceptable or licit semivariogram models which included Spherical, Exponential, Gaussian and Power models. Each semivariogram had advantages and disadvantages. When choosing a semivariogram therefore, the calculation time and the flexibility of the model (the ability to accurately accommodate a broad range of datasets) was taken into account. The spherical model (Figure 6A) is very fast but the least flexible and actually reaches the specified sill value, c, at the specified range, a. The exponential model (Figure 6B) offers a flexible transformation which is faster but the shape of the semivariogram is not flexible. It is slow compared to Power and spherical models. The Gaussian model (Figure 6C) approaches the sill asymptotically, with a representing the practical range which is the distance at which the semi variance reaches 95% of the sill value. The Gaussian model, with its parabolic behavior at the origin, represented the very smoothly varying properties of aggregate stability (However, using the Gaussian model alone without a nugget effect can lead to numerical instabilities in the kriging process.) The spherical and exponential models exhibited linear behavior at origin, which was appropriate for representing these properties which had a higher level of short-range variability. Models with a finite sill, like the Gaussian, exponential, and spherical, are referred to as transition models and have corresponding covariance functions. The power model (Figure 6D), however, does not reach a finite sill and does not have a corresponding covariance function.

The Gaussian model (Figure 6C) was selected for fitting the data onto the semivariogram as it had a Standard Mean which was most nearest to zero, the smallest root-mean-squared prediction error, the smallest difference between the Average Standard Error (ASE) and Root Mean Squared (RMS), and a Standardized Root Mean Squared (SRMS) prediction nearest to one as recommended by ESRI, 2007. The calculations leading to the choice of Gaussian Model is shown in Table 1.

Table 1 Parametric Comparison of the various Semivariogram model

Mean Root- Mean- Squared [RMS]

Average Standard Error [ASE]

ASE- RMS

Mean Standardized

Root-Mean- Square Standardized

Circular -0.6075 13.52 13.71 0.19 -0.03936 0.9925

Spherical -0.4487 13.57 13.69 0.12 -0.02833 0.9976

Tetraspherical -0.4333 13.62 13.74 0.12 -0.02732 0.9969

Pentaspherical -0.434 13.89 13.62 -0.27 -0.02677 1.026

Gaussian -0.4375 13.35 13.71 0.36 -0.02763 0.9812

Exponential -0.434 13.89 13.62 -0.27 -0.02677 1.026

Rational Quadratic -0.4303 13.87 13.81 -0.06 -0.02683 1.008

Hole Effect -0.4674 13.45 14.05 0.6 -0.03032 0.9615

K-bessel -0.1188 13.75 13.93 0.18 -0.004366 0.9917

J-Bessel -0.152 13.42 14.18 0.76 -0.008349 0.9514

The spatial dependence in this Gaussian process can also be used in showing the patterns of ridges and valleys. The notion of spatial heterogeneity in current spatial statistics is only used to characterize local variance of spatial dependence or regression. The heterogeneity is power law like rather than Gaussian distribution like. With this broad perspective, both spatial dependence and heterogeneity depicted the true picture of the Earth's surface. There are far more small things than large ones across all scales or globally, but things are more or less similar at one scale or locally. Empirical Bayesian kriging (EBK) offered the multiplicative skewing normal score transformation.

The new data was transformed and a new semivariogram model was simultaneously estimated from the simulated data. Predictions and prediction standard errors were made using weights and then back transformed with bias correction.When the data distribution is Gaussian, the best predictor is one that uses a linear combination of the nearby data values. For other distributions, however, the best predictor is often nonlinear and, therefore, more complex. The data can be transformed to follow a Gaussian distribution. Then it is possible to accurately back transform kriging predictions to the original data scale, which can be done in ArcGIS Geostatistical Analyst. Classical kriging also assumes that the estimated semivariogram is the true semivariogram of the observed data. This means the data was generated from Gaussian distribution with the correlation structure defined by the estimated semivariogram (Figure 6C). This is a very strong assumption, and it rarely holds true in practice. Hence, action should be taken to make the statistical model more realistic.

Figure 6: Several Semivariogram models being compared to the Gaussian model.

Spatial modeling to generate empirical semivariogram was computed using the formula below:

   – + = ∑

  

() () ( ) () [ ]

2 2 1 h N

) (

i i x Z h x Z h N h γ

1

=

i

Where Z (xi + h) and Z (xi) are observation positions at positions xi + h and xi, respectively, h is the distance between observations and N (h) denotes the number of pairs of observations separated by the same lag distance h. Semivariogram (maximum distance h) = 0.5 MD* average [(value at location i– value at location j)2]

The semivariogram for aggregate stability analysis describes the spatial structures:

γ(h) = Co + C.

Co = nugget effect (NE (%))

Co + C = the sill (total variance);

a = practical range

𝐶𝑜

NE (%) = (

𝐶𝑜+𝐶1)x100 In which C0 is the nugget effect and C1 is the contribution, both semivariogram parameters. According to Cambardella et al. (1994), NE has the following classification: strong spatial dependence (NE (%) ≤ 25 %), moderate spatial dependence (25 % < NE (%) ≤ 75 %), and weak spatial dependence (NE (%) > 75 %).

The second traditional index of spatial dependence (SPD), presented in Biondi et al. (1994), is given by the expression:

SPD (%) =(

𝐶𝑜+𝐶1)x100, in which C0 is the nugget effect and C1 is the contribution, both semivariogram parameters like the NE index. Adjusting the classification given by Cambardella et al. (1994), we have the following induced SPD classification: weak spatial dependence (SPD (%) ≤ 25 %), moderate spatial dependence (25 % < SPD (%) ≤ 75 %), and strong spatial dependence (SPD (%) > 75 %). We can observe that NE (%) = 100 % - SPD (%), that is, NE and SPD essentially provide the same information. Therefore, because of equivalence between the two indexes, in this article we will not regard the NE index as essentially different from the SPD index.

𝐶𝑜

The spatial dependence index (SDI) is given by the following expression for Gaussian model.

𝐶𝑜

𝑎

SDI (Gaussian%) = 0.504𝑥(

0.5𝑀𝐷) x100 in which C0 is the nugget effect, C1 is the contribution, and a is the practical range, both semivariogram parameters like the NE and 0.5MD is half of the maximum distance (MD) between sampled points, given that in the cases in which the ratio a 0.5MD results in a value superior to 1, this ratio is then truncated into 1, in order to assume values only between zero and 1. In addition, the 0.5MD factor is inspired on practical recommendations for using pairs of points up to half of the largest sampling distance in order to estimate semi variances (Journel and Huijbregts, 2003;Olea, 2006). The constant 0.504 or 50.4% (for Gaussian semivariogram) is the value of the model factor (MF) of the model that expresses the strength of the spatial dependence that the model can reach.

𝐶𝑜+𝐶1)x (

A run of the data using EBK model derived the empirical semivariogram/covariance spatial dependence (SDP) model:

SDP = 75.153*Gaussian (8711.3, 4525.7, 51.8) +158.48*Nugget,

Where: Nugget = microstructure = Measurement error

Microstructure = 172.6 (100%) and

Measurement error =0 (0%)

Given that the SPD for the Gaussian exhibit positive asymmetric distribution, it must be categorized having the median and the 3rd quartile as limits (Table 2), in order to have a classification that is more coherent and in fitting with the distribution form. Thus, based on the results above, we propose the following classification for the SPD and SDI Where SDI = 100 – SDP or SDP = 100 – SDI giving SDI=25% and SDP= 75%

SPD Gaussian(%) ≤ 25% will therefore, give weak spatial dependence, moderate spatial dependence 25% < SDP ( % ) ≤ 75% and strong spatial dependence SDP (%) > 75%, while SDI Gaussian (%) ≤ 25% gives a strong spatial dependence index, moderate spatial dependence index 25% < SDI (%) ≤ 75%, and weak spatial dependence index SDI (% ) > 75%. The Gaussian semivariogram model has high spatial dependence, with a high MF which led to a wider scope of values. The semivariogram model was used to derive a soil aggregate stability prediction graph illustrated by Figure 7 and a semivariogram model Figure 8.

Figure 7: Soil aggregate stability prediction graph

Figure 8: The semivariogram values for the pairs of points (red) and the estimated semivariogram model (yellow line)

The semivariogram models were then used to undertake Empirical kriging to create a prediction map of aggregate stability. A prediction map of the spatial dependence derived from the EBK modeling process is shown in Figure 9 below.

Figure 9: Empirical Bayesian Kriging prediction map with sample values.

In the figure above, the mean weight diameters (MWD) of 0.25 – 0.45 represents unstable soils mostly found in wetlands occurring in valleys, mountains, plains, and depressions in hills, 0 55 –0.62 represents moderately stable soils mostly in agricultural and grassland areas that include plateaus, uplands, and plains, while 0.62 – 0.92 represents stable and very stable soils being found in forested areas, mountains and hills.

3.3: Exploring Data Characteristics

3.3.1: The Histogram

Data was statistically analyzed in two phases: (1) data were described using the classical statistics (mean, mode, median, and standard deviation, coefficient of variation, skewness, and kurtosis); (2) frequency distribution was examined and the test for normality was conducted. Geostatitical methods are optimal when data are normally distributed and stationary that is when mean and variance do not vary significantly in Space as in this case. Significant deviations from normality and stationarity can cause problems, so it is always best to begin by looking at a histogram or similar plot to check for normality and a posting of the data values in space to check for significant trends. Looking at the histogram Figure 10 (with a normal density superimposed) and a normal quantile-quantile plot shows that the aggregate stability distribution does not deviate too severely from normality: Judging from Table 2, the distribution of data is close to normal with an almost zero skew and the mean being equal to the median. Aggregate stability values were collected from55 sample points distributed throughout the field, which is approximately 8, 000 hectares.

Table 2. Distribution of Data

Figure 10: Distribution of Aggregate stability data in Mean Weight Diameter (MWD)

The determination of aggregate stability through the set of sieves involved the estimation of the amount of intact aggregates against the forces of water entry into aggregates, such as forces related to water entry in soil aggregates. It was generally considered that retention of large aggregates on the top sieve against forces of water entry was indicative of good aggregate stability Therefore, the aggregate stability of the soil was evaluated as mean weight diameter (MWD) determined by the wet sieving. The MWD of the aggregates was calculated where the number of size fraction = the mean diameter of any particular size range of aggregates separated by sieving = the weight of aggregates in that size range as a fraction of the total dry weight of the sample analyzed. The values ranged from 0.25 to 0.92 MWD. Here are the data values posted at the Sample locations (Figure 11).

Figure 11: Plotted Aggregate stability values according to their location

3.3.2: Normal Quantile - Quantile Plot

The data was subjected to some normality checks with which Quantile - Quantile plot was a major test of normal distribution in the data. The data in this Q-Q plot appear to be normally distributed, because even though there is a slight, possibly curved trend in the plot, the dots are still pretty close enough to the line to not disqualify these data from being normal. The Q-Q plot (Figure 12) looks alright except for one outlier depicted by the point furthest to the right. All geostatistical methods assume spatial dependence that closer things are more similar than things that are farther away. Spatial dependence is "the propensity for nearby locations to influence each other and to possess similar attributes"(Tobler, 1989). In other words, to paraphrase this famous geographer, “while everything is related to everything else, things that are close together tend to be more related than things that are far apart”. Terrain elevations and soil types for instance, are more likely to be similar at points few meters apart than at points several kilometers apart.

Figure 12: Normal Q-Q plot

Given that geographic data are expensive to create; spatial dependence turns out to be a very useful property. We can sample attributes at a limited number of locations, then estimate the attributes of intermediate locations. The process of estimating unknown values from nearby known values is called interpolation. Interpolated values are reliable only to the extent that the spatial dependence of the phenomenon can be assumed. If we were unable to assume some degree of spatial dependence, it would be impossible to represent continuous geographic phenomena in digital form.

4.0 Conclusion

The soils of the geomorphic units of mountains, hills and minor scarps were developed on volcanic ashes and other pyroclastic rocks of recent volcanoes. They had slopes ranging from 8-30 % and had stable to very stable soil aggregates. The plateaus (PU) (Agriculture and grasslands and the uplands (UP) (agriculture, grasslands and minor forests) were undulating with slopes of between 5 – 8%. The soils which were moderately stable were developed on volcanic ashes and other pyroclastic rocks of recent volcanoes. The plains (PL) (grasslands and agriculture) and Valleys (VA) (grasslands, agriculture and wetlands) were flat to very gently undulating in relief with moderately stable to unstable soil aggregates. They had slopes of between 0 – 2% with well drained, deep, and moderately deep to very deep soils.

The variability of aggregate stability exhibited spatial dependence that was described using semivariogram models. The Gaussian model was selected for fitting the data onto the semivariogram model The spatial dependence (SPD) for the Gaussian exhibited positive asymmetric distribution giving a strong spatial dependence (SPD) and a strong spatial dependence index (SDI) indicated by SPD Gaussian (%) > 75% and SDI Gaussian (%) ≤ 25%. The semivariogram models were used to undertake Empirical kriging to create spatial dependence prediction map of aggregate stability. The mean weight diameters (MWD) of 0.25 – 0.45 represented unstable soil aggregates, 0 55 –0.62 moderately stable soil aggregates and 0.62 – 0.92 representing stable and very stable soil aggregates.

References