Spatial and fractal characterization of soil properties across soil depth in an agricultural field, Northeast Iran

The present study was conducted to explore the fractal behavior and establish fractal dimensions of soil physical and chemical properties (i.e., sand, silt, and clay contents, bulk density, degree of moisture saturation, pH, organic carbon content, total nitrogen, available phosphorus, and available potassium) to characterize their spatial patterns. Soil samples were collected from 0-30 (surface) and 30-60 cm (subsurface) depths from an agricultural field, Mashhad Plain, Northeast Iran. Descriptive statistics and fractal analysis were used to describe the extent and form of variability. Spatial patterns of the soil properties were estimated using GS+ 10.0 software.  Soil properties showed low to high variations in both surface and subsurface layers across the field, where bulk density and pH being the most reliable soil physical and chemical properties in the study area. The variability was high (CV > 35%) for total N, available P, available K and organic carbon in both surface and subsurface soils and it could be attributed to management practices and micro-topographical variations as these are the dynamic properties of soil. The fractal dimension (D) values of soil physical properties ranged from 1.398 to 1.913 at the surface, and from 1.874 to 1.934 at the subsurface indicating both short and long range variations. The D values for the chemical properties ranged from 1.331 to 1.975, and 1.148 to 1.990 in the surface and subsurface layers, respectively. The results showed that fractal analysis could be employed to effectively describe the structure of soil heterogeneity in spatial scale for effective agricultural and environmental management of soil.


Introduction
Soil covers land as a continuum having properties that vary enormously and continuously with depth and horizontal distances (Gessler et al., 2000). This variability is reflected in soil test results and variations in crop yields (Tuffour, 2015). These variations are largely the results of various factors and processes of soil formation, and management, functioning independently or in combination over broad spatial and temporal scales with different intensities (Beckett and Webster 1971;Burrough, 1983a,b;Tuffour, 2015). The spatial variation of soil properties such as organic matter, clay content, pH and water retention capacity is caused by pedogenic processes which are influenced by hydrological and temperature regimes modified by topography (Pilesjő et al., 2005;Florinsky, 2012;Vasu et al., 2017a). The source of these variations, which hitherto recognized as a nuisance, has received widespread recognition among researchers in recent times (Bengough et al., 2000;Linsenmeier et al., 2011;Tuffour, 2015;Miloš and Bensa, 2017;Vasu et al., 2017a). Knowledge on the spatial distribution of soil properties at different scales is a prerequisite for site specific soil management and minimizes the negative effects of management practices on environmental quality (Cambardella et al., 1994;Bogunovic et al., 2014). In recent past, GPS (Global Positioning System), GIS (Geographical Information System) and geostatistics have played a very significant role in the study of spatial variability of soil properties (Vasu et al., 2017b). Geostatistical tools incorporate the spatial coordinates of soil observations in data processing, facilitate description and modelling of spatial patterns, predict values for unsampled locations and validates these predictions (Goovaerts, 1998). Because of these advantages, geostatistical techniques have been widely employed by researchers to address the heterogeneity in soil properties and processes in space under different environmental conditions (Wang and Shao, 2011;Tuffour et al., 2013Tuffour et al., , 2016Kooch et al., 2014). The structure of spatial variation is vital to decide the magnitude of variability and that can be assessed by semivariograms which indicate the spatial correlation in values measured at sample locations (Vasu et al., 2017b). The spatial patterns are described in terms of dissimilarity between samples as a function of separation distance (Goovaerts, 1998). The average dissimilarity between data is separated by a vector and is measured by the semivariogram. However, anomalies in soil data and their deviation from normal (Gaussian) or log-normal distribution need to recognized for identification, accurate delineation, and modeling of soil properties for site specific soil management and precision agriculture (Fu et al., 2010;Bogunovic et al., 2014;Vasu et al., 2017a). Soil particles were assumed to have "self-similar" features, and fractal theories were employed to investigate their characteristics in several studies (e.g., Oleschko et al., 2008;Tuffour, 2015;Li et al., 2016;Miloš and Bensa, 2017). For example, fractal dimensions of particle size distribution (PSD) were found to be influenced by land use and management practices (Wang et al., 2008). For newly formed wetlands in Yellow river delta, the fractal dimension values (D) of PSD varied from 1.82 to 1.90 indicating coarse texture (Yu et al. 2015). Therefore, the objective of this study was to investigate the fractal behavior and establish fractal dimensions of soil physical and chemical properties to characterize their spatial patterns.

Study area
The study site is located in Mashhad Plain, Khorasan-e-Razavi Province, Northeast Iran with an area of 6131 km 2 (Figure 1). The region is located between latitudes 35° 59′ N to 37° 04′ N and longitudes 58° 22′ E to 60° 07′ E. The elevation of the study area ranges from 900 to 1500 m a.s.l., while the elevation of surrounding areas is up to ~1200 m a.s.l. The study area is characterized by semi-arid climate with mean annual precipitation of 222.1 mm and mean annual temperature of 15.8°C (Keshavarzi et al., 2016). The major soil types include Calcaric Cambisols, Gypsic Regosols, Calcaric Regosols and Calcaric Fluvisols occurring in pediment plains, plateau, upper terraces and gravelly colluvial fans, respectively. Soil texture varied from loam, to sandy loam to sandy clay loam. The main land use in the study area is irrigated farming.

Soil sampling and laboratory analyses
A Digital Elevation Model (DEM) with 10×10 m grid size was extracted from a paper-based topographic map (1:25000 scale) using GIS platform and 10-meter contour lines interval. Vasu et al. (2016) established the importance of subsurface soil properties in evaluation of soil quality, spatial variability and their management. Disturbed and undisturbed samples were collected from surface (0-30 cm) and subsurface (30-60 cm) of 48 representative soil profiles (96 samples) by stratified random sampling technique. The sampling points were designed to represent all the major soil and land use types. Large plant materials and pebbles in the samples were separated by hand and discarded. Collected soil samples were air dried, crushed and sieved by using 2 mm sieve size before laboratory analysis. Soil organic carbon (SOC) content was determined by the Walkley and Black method (Walkley and Black, 1934) with dichromate extraction and titrimetric quantization (Nelson and Sommers, 1986). Percentages of clay (< 0.002 mm), silt (0.002-0.05 mm), and sand (0.05-2 mm) particles were measured by hydrometer method (Gee and Bauder, 1986). Soil moisture content based on the degree of saturation (DS) and bulk density (BD) (Blake and Hartge, 1986;Gardner, 1986) were also determined. Soil reaction (pH) was measured in saturated paste extract using a digital pH-meter (Thomas, 1996). Total N, Olsen P (available phosphorus) and available potassium were measured using standard procedures (Sparks et al., 1996).

Descriptive statistics
Measured variables in the data set were analyzed using descriptive statistics including the minimum, maximum, mean, coefficient of variation (CV), skewness and kurtosis. The CV values were classified according to Wilding (1985). Due to the presence of small variations that could result in chance fluctuations in the skewness and kurtosis (Tuffour et al., 2013(Tuffour et al., , 2016, the frequency distribution was validated using the Watson normality test.

Fractal analysis
Fractal dimension ( ) was estimated from the slope of log-log semivariogram. Semivariogram analysis was conducted to ascertain the spatial structure of the soil properties using GS+ (Geostatistics for environmental science; Gamma Design, Plainwell, Mich.) 10.0 software. Isotropy and anisotropy of the semivariograms were determined in four directions (i.e., 0°, 45°, 90° and 135°) with a tolerance of 22.5°, and the structure of spatial variance between observations was calculated from the semivariogram defined by the equation (Gupta et al., 2003;Tuffour, 2015): (1) where, (ℎ) is the semi-variance for separate distance class ℎ, (ℎ) is the number of sample pairs at each distance interval ℎ, ( ) is the value of the variable at sampled location and ( + ℎ) is the value of the variable at a distance ℎ away from (Lark, 2000).
The fractal dimensions ( ) were estimated from the slopes of log-log semivariograms as (Tuffour, 2015): where: = Fractal dimension (Hausdorff-Besicovitch statistic), between 1 and 2 = Slope of the log-log semivariogram plot To measure the smoothness of the measured data sets, the Hurst exponent ( ) was estimated from as follows (Tuffour, 2015): The closer the value of to 0, the more random will be the distribution of the property.

Results and Discussion
Descriptive statistics for soil properties Soils in the study area were predominantly Aridisols, with a few patches of Entisols. The descriptive statistics of the analyzed soil properties are presented in Tables 1.1 and 1.2. The variability was interpreted using the coefficient of variation (CV). The results indicate that the CVs varied from 1.47 to 91.49%. The criteria proposed by Wilding (1985) was used to classify the parameters into most (CV >35%), moderate (CV 15-35%) and least (CV <15%) variable classes. Accordingly, the variability was high (CV >35%) for total N, available P, available K and organic carbon in both surface and subsurface soils and it could be attributed to management practices and micro-topographical variations as these are the dynamic properties of soil . Clay is an intrinsic soil property whose depth distribution is influenced by parent material, topography and hydraulic conductivity. In the present study, clay varied moderately in surface soils (CV=33.43%) but found to be high in variability in subsurface soils (CV=41.02%). Bulk density and pH were low in variability in both surface and subsurface soils (Tables 1.1, 1.2). Similarly, Tuffour et al. (2016) observed low variability for BD in the 0-20 and 20-40 cm depth of an Acrisol.  The moderate variability observed for DS in both layers was due to the influence of both static (topography and soil properties) and dynamic (precipitation and soil moisture content) properties of the study area, which causes variation in soil moisture patterns (Western et al., 2004;Mapfumo et al., 2006;Tuffour et al., 2016). The observed variability in the soil physical properties considered in this study, viz., BD, DS and particle size distribution have important implications on infiltration rate, and runoff/erosion processes. They could have severe impacts on the variability soil-plant-water relationships, and on crop growth (Tuffour et al., 2016). The mean values of sand, clay, BD and DS were higher in the subsurface layers than the surface layers but mean silt content was higher in the surface layer than in the subsurface layer. The high heterogeneity observed for the primary soil particles (i.e., from moderate to high) in both depths could be attributed to land use and soil management in the study area (Tuffour et al., 2016).
Skewness, kurtosis and the Watson normality test results indicated that only sand and silt contents in the surface layer and total N in subsurface layer followed normal distribution. The mean values show high contents of N, P, and K in the surface layer than the subsurface layer. However, OC was higher in the surface than the subsurface layer. The high N, P and K contents observed in the subsurface layer can be attributed to leaching as evidenced by the high degree of moisture saturation in the subsurface layer than the surface layer (Table 1.1), which suggests high water transmission and infiltration rate of the soils studied. The high OC content in the surface layer, however, could have resulted from the accumulation of organic residues from different sources with different compositions and rates of decomposition, and the differential growth behavior of crops in the field (Santra et al., 2008;Tuffour et al., 2013). Similarly, studies by Camacho-Tamayo et al. (2008), Balasundram et al. (2009 and Tuffour et al. (2013) also reported high OC content in the soil surface than the subsurface layer.
Soil pH was not significantly different among the soil layers, however, it was slightly higher in the subsurface than the surface layer. Soil pH in both soil layers showed a more consistent pattern than other properties which was evident from the low CVs (Table 1.2). This low variation in pH could be attributed to the sampling distance adopted in this study, since soil pH has been reported to show higher variations at scales of centimeters than at meters (Göttlein and Stanjek, 1996;Tuffour et al., 2013). However, relatively higher variability of pH in the surface layer than the subsurface layer could be the result of soil management practices, such as tillage and fertilizer application (Huang et al., 1999), seasonal fluctuations in soil moisture, temperature, microbial activity, plant growth and parent materials of the soil in the study area (Tuffour et al., 2013).

Fractal characterization of soil physical and chemical properties
The fractal dimensions (D) determined from the slope obtained by plotting the semivariance against the sample spacing (distance) (Figures 2.1-2.10) were used to describe the spatial variability of the soil properties. The Hurst exponents (H) estimated from Equation (3) are presented for various depths in Table  2. The values of D ranged from 1.398 to 1.934 for the physical properties, and 1.148 to 1.990 for the chemical properties. The range of H for the physical properties was between 0.066 and 0.602, and 0.01 and 0.852 for the chemical properties.   The results obtained for D revealed that the soil properties show fractal behavior, wherein, lower the value of D, longer the range of variation, and vice versa (Vieira et al., 2010;Tuffour, 2015). Based on the scaling properties of fractional Brownian motion reported by Sugihara and May (1990) and Tuffour (2015), H values > 0.5 represent smooth or persistent nature of the soil property, whereas, H < 0.5 represents antipersistent nature. Liao et al. (2017) successfully employed fractal analysis to investigate the spatio-temporal variability of soil moisture content in contrasting land uses. The above studies showed that fractal dimension can be used for assessing the spatial distribution of soil properties at different scales. Mohammadi et al. (2017) analyzed spatial variability of soil textural fractions and fractal parameters derived from Particle Size Distributions (PSD). In this research, fractal features of particle size distribution of soil samples were studied using fractal geometry and then the geostatistics approach was applied to characterize the spatial variability of fractal and soil textural parameters. According to the semivariogram models and validation parameters applied to the models, it was found that the fractal parameters had powerful spatial structure and could better describe the spatial variability of soil texture. The results of the present study ( Table 2) showed that the nature of the distribution of the soil properties ranged from smooth (persistent) to rough (anti-persistent). For instance, DS in the surface layer with D and H values of 1.398 and 0.602, respectively implies a long range variation and a more persistent nature than in the subsurface layer (D = 1.876 and H = 0.124), wherein it would exhibit a short range variation and an anti-persistent nature. Furthermore, with the exception of DS in the surface layer (H = 0.602), and OC in both the surface and subsurface layers (H = 0.669 and 0.852, respectively), all the soil properties were anti-persistent in nature across the study area.

Conclusion
In this study, classical and fractal methods were applied to reveal and describe the spatial variability of soil properties (i.e., sand, silt, and clay contents, bulk density, degree of moisture saturation, pH, organic carbon content, total nitrogen, available phosphorus, and available potassium) in an agricultural field, Northeast Iran using fractal parameters (fractal dimension and Hurst exponent). Comparing the results of the descriptive statistics and fractal analysis, revealed that the selected soil properties in the study area exhibit short to long range variations at the local scale. The results also showed that fractal characterization is an applicable technique for describing the spatial patterns of soil properties.