Geostatistical Characterization of Soil Clay Patterns in a Typical Watershed of Cultivated Black Soil

Soil clay inuences a range of processes, including vertical and lateral redistribution of water, runoff, and erosion. Knowledge of the spatial pattern of soil clay content (SCC) in cultivated soil region is important. The objective of this study was to determine the degree of spatial variability of SCC across depths in a cultivated black soil region in northeastern China. Data collected from each of 58 sites on a regular grid of 50 by 50 m were analyzed both statistically and geostatistically to describe the spatial distribution. The SCCs between different segments in soil proles were signicantly different. The 0- to 65-cm soil prole could be clustered into three segments: 0- to 30-cm, 30- to 35-cm, and 35- to 65-cm layers.The standard deviation (SD) and coecients of variation (CV) for the SCC of the 30- to 35-cm segment were the largest. Experimental semivariograms of SCC were best tted by spherical models. Nugget-to-sill ratios indicated a strong spatial dependence for SCC at other depths, except the 40- to 45-cm and 60- to 65-cm segments. The 35- to 40-cm layer had the largest spatial dependency compared with the other layers. Cross-validation of the kriging map indicated that prediction of SCC using semivariogram parameters was better than assuming the mean of the observed value for any unsampled location. The mapping displayed heterogeneity of SCC across the experimental site and revealed higher SCC close to the tail of an eroded gully and lower SCC neighbouring eroded gully margins. The fragmentation degree and clay-enriched patch amount increased from the near-surface down to 65-cm of depth, suggesting the higher evenness of SCC in the cultivated layer than that in the tillage pan.


Introduction
Soil properties and their spatial distribution on a large scale are available in soil survey reports. When soil data are required at a much ner scale, on-site detailed sampling is often introduced (Ferre et al. 2015).
Soil variability is usually associated with spatial, temporal, or management-related factors, and each of these sources can partially or fully contribute to the variability of a soil property under investigation (Korres et al. 2015). For cultivated soils, tillage practices can disturb spatial variability of soil clay content (SCC; Rejman et al. 2014). Clay has a major in uence on a range of processes, including vertical and lateral redistribution of water, runoff, and erosion (He et al. 2014). Understanding the spatial characteristics of SCC across the landscape is important for evaluating the in uence of agricultural activities on landscape patterns (Mao et al. 2015).
Geostatistics is a useful tool for analyzing spatial variability (Penizek et al. 2004). Variograms can re ect the structural characteristics of a regional variable (Liang et al. 2014). Spatial dependence of soil properties is reported from meter to kilometer scales (He et al. 2015; Wang et al. 2013). SCC variability focuses on the surface or near-surface soil layer (Brye et al. 2006;Takagi et al. 2012) and region scale. Vauclin et al. (1983) used a classical and geostatistical technique to study the spatial variability of sand, silt, and clay contents. Iqbal et al. (2005) sampled 1-m-deep soil pro les to determine the degree of spatial variability of soil physical properties and the variance structure of surface, subsurface, and deep horizons in soils developed in alluvium. Shukla et al. (2007) indicated the variability of sand, silt, and clay contents in the 0-to 15-cm depth of mine soils was mostly low for clay and low to moderate for silt and sand content. Scienti c information on the spatial variability and distribution of SCC is critical for designing sustainable soil-crop and environmental management decisions. However, little is known on spatial distribution and variability of SCC at the watershed scale in many cultivated regions, including the black soil region of China (Kravchenko et al. 2006). The objectives of this study were to (i) describe and quantify the variability of SCC in cultivated black soil and to (ii) assess the strengths of spatial correlation for SCC at different soil horizons.
The soil began to form at the beginning of the Quaternary Period from loess deposits and has a rich-dark organic horizon. The black soil is classi ed as silty clay loam Mollisols (USDA, 1999) and Phacozems (FAO, 1998). There is an obviously eroded gully approximately 383 m long and 2-9 m wide. The width of the eroded gully was bigger in the upper slope than the attened foot slope. On the southern edge, the eroded gully is irregular, and many branching eroded channels exist.
The watershed has been used in a rotational system for growing corn, soy, and potatoes since 1948.

Soil Sampling and Laboratory Analysis
Soil was sampled using the soil boring method from each of 58 sites on an regular grid of 50 by 50 m in October 2008 at the top 10 segments (0-to 10-cm, 10-to 20-cm, 20-to 30-cm, 30-to 35-cm, 35-to 40-cm, 40-to 45-cm, 45-to 50-cm, 50-to 55-cm, 55-to 60-cm, and 60-to 65-cm depth) (Fig. 1). A soil sample from each segment was placed into its designated bucket, mixed, and packed into boxes. Soil samples were air dried, crushed, sieved (2 mm), and analyzed for soil clay content using the hydrometer method (Gee et al. 1986). The coordinates of each sampling locations were recorded using No. NTS-302 electronic total station (South Surveying & Mapping Instrument CO., LTD.)

Statistical Methods
Soil clay data were analyzed using classical statistical methods to obtain the mean, stand deviation (SD), coe cients of variation (CV), median, minimum, maximum, skewness, and kurtosis for each horizon (n = 58). Wilding (1985) proposed a ranking of variability in soil properties that occur within landscape units of a few hectares or less. He ranked a CV < 15% as the least, 15% < CV < 35% as moderate, and CV > 35% as the most variable.
A one-way ANOVA was also performed using STATISTICA (Version 7.0, Statsoft, Inc., Tulsa, USA) to compare SCC across the soil pro les using a protected least signi cant (P < 0.05) difference test (Table 1). Means within a column that share the same letter do not differ at P ≤ 0.05.
Cluster analysis for the SCC in each horizon was utilized to acquire the relationship among different soil horizons using STATISTICA (Version 7.0, Statsoft, Inc., Tulsa, USA).
The degree of spatial variability for SCC was determined by geostatistical methods using semivariogram analysis and kriging (Iqbal et al., 2005). The experimental semivariogram function (Goovaerts, 1997) was calculated as follows: where γ (h) is the semivariance for lag interval h, z(x i ) and z(x i + h) are variables at locations x i and x i + h, and N(h) is the number of pairs separated by distance h.
Directional (north-south, northeast-southwest, west-east and northwest-southeast with an angular tolerance of 45°) and omnidirectional variograms were used in this analysis. Theoretical variograms were obtained by tting spherical, exponential, and Gaussian models using GS + ( nally selected as the best-tting model. The nondimensional number obtained as the ratio of the nugget to the sill, i.e., the 'relative nugget effect', was used as an indicator of spatial dependence (Shukla et al., 2007). We adopted spatial class ratios to de ne distinctive classes of spatial dependency. If the ratio of spatial class was < 25%, the variable was considered strongly spatially dependent; if the ratio was > 25% and < 75%, the variable was considered moderately spatially dependent; and if the ratio was > 75%, the variable was considered weakly spatially dependent (Iqbal et al., 2005).
We created contour maps of SCC for each horizon through ordinary kriging using their respective semivariogram models in ArcMap (Version 9.3, ESRI Inc., USA).
Within the 30-to 65-cm depth range, SCC increased with increasing depth. The SCC of the 30-to 35-cm horizon was signi cantly lower than the 35-to 40-cm, 40-to 45-cm, 45-to 50-cm, 50-to 55-cm, 55-to 60-cm, and 60-to 65-cm horizons. No signi cant difference was found in SCC among the 10-to 20-cm, 20to 30-cm and 30-to 35-cm horizons. According to the variability of SCC in the 0-to 65 cm soil pro le above, the 0-to 10-cm and 30-to 35-cm horizons were obviously different from other neighboring soil horizons.
The Kolmogorov-Smirnov test indicated that observed distribution of SCC at all depths were signi cantly different from normal distribution (Fig. 2, P < 0.01). The variability of SCC for each horizon was moderate (CV = 16.8-28.6%), and the SCC variability of the 50-to 65-cm depth range was less than the 0-to 50-cm ( Table 1). The CV of SCC was greatest at the 30-to 35-cm horizon (Table 1).

Spatial Structure Analysis
The experimental semivariograms for the SCC of different soil horizons were obtained with a lag of 20 m and a cut distance of 250 m (maximum separation distance was 522 m) to avoid errors from using greater lag distances (Journel et al. Figure 4). Directional sample variograms were plotted for each horizon of 0-to 65-cm soil pro le to examine the anisotropy of the spatial correlation structure. Figure 4 shows north-south, northeast-southwest, west-east, and northwest-southeast variograms for different horizons. The spatial correlation structure was anisotropic with a shorter range in the north-south and northeast-southwest directions than the west-east and northwest-southeast directions. The anisotropy was related to the shape of the watershed. The hillslope in the north-south direction was shorter than the west-east direction. Drainage patterns that tended to run northwest-southeast were longer than those that were northeast-southwest; the slope of the southwest hillslope (10.8°) was greater than the northwest hillslope (9.7°).
For the remainder of the paper, only omnidirectional variograms were considered for two reasons. First, the anisotropy observed was likely to be a function of the orientation for this particular watershed rather than a characteristic of the landscape as a whole. Second, topographic control was not the main focus of this paper. SCC exhibited differences in spatial pattern for each horizon (Fig. 5). Spatial patterns for SCC also varied among different horizons in both magnitude and space (Figs. 4 and 5). The SCC in the foot slope was higher than the upper slope in each horizon, except for the 40-to 45-cm (Fig. 5). The irregularity of SCC increased with horizon depth (Fig. 5). In the 60-to 65-cm horizon, the greatest irregularity was observed (Fig. 5). The spherical model provided good estimates of isotropic semivariogram parameters ( Fig. 4; R 2 = 0.325 ~ 0.696), as indicated by the small sum of squared deviations (SS) between experimental and theoretical semivariograms (Table 2, Fig. 4).
The semivariance increased with distance between sample locations, or lag distance, to a constant value or sill (total variance) at a given distance (Fig. 4), which is known as the range of spatial dependence (Iqbal et al., 2005). All the omnidirectional variograms indicated that the soil clay content was stationary because they exhibited clear sills (Fig. 4), which were close to the sample variance (Sill = 1.0561×variance, R 2 = 0.9648, p < 0.001). Theoretically, the semivariance at h = 0 is equal to zero, but the experimental semivariogram frequently exhibited a discontinuity known as the nugget variance, which is the local variation occurring at scales ner than the sampling interval, such as sampling error, ne-scale spatial variability, and measurement error (Iqbal et al., 2005). All variogram models of SCC in different soil horizons showed a positive nugget effect (Fig. 4).
Nugget variances ranged from low (or 0.1, 0.2% of the sill) to high (or 40.7, 42.0% of the sill) and can be considered moderate (< 50%). Thus, little variation was present at distances shorter than the rst lag (20 m) of the semivariogram.
Sill values of the theoretical semivariograms of SCC in each horizon were similar to the sample variance (Fig. 6), indicating a general absence of trends (Shukla et al., 2007). The sill varied between 36.87 and 96.90. Among the 10 horizons of the 0-to 65-cm soil pro le, the sill value was highest in the 30-to 35-cm horizon (Table 2, Fig. 6). The higher dispersion variance for SCC in the 30-to 35-cm horizon than other horizons was largely due to the variability of tillage depth. For SCC in the 0-to 35-cm soil depth, however, the dispersion variance was smaller in the 0-to 20-cm horizons. The smaller dispersion variance in the 0to 20-cm horizons was consistent with the evenly high level of soil disturbance (i.e., tillage practices). In the 30-to 65-cm soil horizons, the sills of the 50-to 65-cm were lower than the 30-to 50-cm.
The range of in uence is considered to be the distance beyond which observations are not spatially correlated (Shukla et al., 2007). Therefore, a de nite and positive range and no pure nugget (range = 0) for experimental variograms (Table 2; Fig. 6) showed that SCC were not completely random at the scale of measurement. The range varied from 41.9 to 146.8 m ( Table 2; Fig. 6). Range values, which varied between 138.3-146.8 m, were similar in the 0-to 40-cm horizons ( Table 2; Fig. 6) and obviously higher than the 40-to 65-cm horizons. In 40-to 65-cm soil horizons, the range decreased with increasing soil depth and reached the lowest (41.9 m) at the 60-to 65-cm horizon.
When the spatial dependence ranges are similar, the relative nugget effect can be used to classify the spatial dependence of soil properties (Shukla et al. 2007). Lower relative nugget effects correspond to stronger spatial dependence (Shukla et al. 2007). In the 0-to 40-cm horizons, the ranges were similar (Fig. 6, 138.3-146.8 m). The relative nugget effect was augmented as the distance to the soil surface increased. This indicated that the spatial dependence of the 0-to 10-cm horizon was strongest. A high coe cient of determination (R 2 = 0.5193) and a signi cant binomial relationship (p < 0.001) were obtained between the soil horizon depth and the relative nugget effect for soil clay contents in 0-to 65-cm horizons. Because range values for all horizons of the 0-to 65-cm soil pro le were widely different (Fig. 6), no other attempts could be made to characterize the spatial dependence using the relative nugget effect across all 10 horizons.

Discussion
In general, vertical variation of SCC within an undisturbed natural soil pro le is simple. SCC normally increased with increasing soil depth, but in agricultural regions, this may be disturbed by human management (i.e.,continuous soil management practices). Shukla et al. (2007) reported that the SCC of reclaimed minesoil and unmined grass soil at the 0-to 15-cm depth varied between 20% and 21%, and Georgia Piedmont pasture at 0-to 20-cm depth varied between 10% and 20% (Endale, et al., 2006). However, our experimental site at 0-to 10-cm depth was 34.6%, which is higher than that of Shukla et al. (2007) and Enadle et al. (2006). In general because it is subjected to rain water in ltration and wind erosion, SCC within the upper range of the soil pro le increases as depth increases (Enadle, et al., 2006). When disturbed by management practices, this tendency could be altered (Table 1). Currently, on the typically agricultural regions of black soil in northeastern China, the tillage depth is set to 30-to 35-cm by modern tilling machines, which we veri ed with our clustering results of the pro le-depth variables of SCC (Fig. 3). The deep soil tilling plough could take abundant clay in the tillage pan to the surface horizon. This is one of the key reasons for high SCC at the 0-to 10-cm depth (Indorante et al. 2014). The higher SCC could induce a higher probability for wind erosion (e.g., sand storms). When tilled by agricultural machines, the upper and lower soil in tillage depth are often inverted, so the disturbed degree of the upper and lower soil horizons are more intense than the middle soil horizon. Although the 10-to 35-cm soil depth is in the range of the tillage depth, the SCC is lower than other neighboring soil horizons (Table 1). This trend could be related to clay depositing processes with water in ltration from the top to the bottom and tilling activity (Colazo et al. 2010;Sun et al. 2015) Although the SCC in the 30-to 35-cm horizon is lowest (Table 1), the CV (Table 1) and sill (Fig. 6) are highest. This could result from the uneven tillage depth, although the tillage depth was often the 30-to 35-cm soil depth. The CV (Table 1) and sill (Fig. 6) of the SCC in the 0-to 10-cm soil depth are lower than other segments, but the SCC of top soil in the experimental site also presents a tendency along the eroded gully (Fig. 5). The SCC on the upper portion of the eroded gully is lower than the footslope (Fig. 5). The existing variability is related to SCC lateral redistribution by water and runoff, terrain elements, and possibly depositional process variances due to different slopes (Kravchenko et al. 2006).
Under such in uences as air, water, solar radiation, and soil preparation, the evenness and the horizontal spatial dependence of the surface soil are higher than those of the subsoil. This nding was presented in the pattern plots (Fig. 5).The texture of the SCC pattern in different soil segments becomes less dispersed and presents less clustered regions as the segment-to-surface decreases (Fig. 5). This could also be attributed to a tillage effect (Gal et al. 2007).
The spatial structure of SCC presents an obviously directional tendency (Fig. 4). This nding could also be attributed to the terrain change (Kravchenko et al. 2006). The correlation lengths of the 0-to 40-cm soil segments are similar and higher than the 40-to 65-cm soil depth. Because the tillage activity of soil preparation using modern machines could nearly till the 0-to 35-cm soil, this could be one of the reasons for the similar correlation lengths of the 0-to 35-cm soil depth.

Conclusions
The vertical variation of SCC in the soil pro le re ects tillage depth. Tillage action was re ected in the tendency of SCC spatial correlation structure change. The natural patterns of SCC in soil pro les are changed by agricultural management practices. The SCC in the top 0-to 10-cm segment was not signi cantly different from the 40-to 65-cm segments and was signi cantly higher than the 30-to 35-cm segments. The SCC was lowest in the 30-to 35-cm segment. The spatial dependence of the 0-to 10-cm horizon was strongest. Figure 1 Elevation contour lines in meters, eroded channel boundary, and location of clay measurement sites    The spatial patterns of clay content in 0-to 65-cm soil pro le.

Figure 6
The change of the geostatistical structure of the soil clay pattern with soil depth increasing in the 0-to 65cm soil pro le, showing (a) the relation between sample variance and sill; (b) the sill; (c) the correlation length; and (d) the nugget.