Evaluating the tillage management direction effects on soil attributes by space series analysis (case study: a semiarid region in Iran)

The assessment of soil attributes affected by land use changes or different cultivation management strategies is commonly based on a comparison between agricultural fields, neglecting the natural soil spatial variability. This study aimed to develop a methodology based on improved space series to differentiate between spatial variability of soil attributes and the effect of tillage direction when the evaluation is based on comparison between adjacent fields. The study area consists of two adjacent fields of different tillage directions, i.e. up-down tillage (UDT) and contour tillage (COT). Soil sampling was performed at 40 points in each filed at 5 m intervals along a contour line at the mid-slope position. All measured soil attributes, i.e. sand, silt, clay, MWD, GMD, bulk density; SP, CCE, OM, of UDT were significantly (p < 0.05) different from those of COT compared by independent sample t test. This analysis could not differentiate between the spatial variability of the soil and the changes induced by tillage. To determine the net effect of UDT on soil attributes, we (i) performed space series analysis on COT data, (ii) used autoregressive, moving average and autoregressive-moving average models to model the space series data on COT field, and (iii) used the best model obtained for each soil attributes on COT to forecast the value of the property in ten adjacent points in the UDT field. Comparison between the forecasted and measured data in UDT showed that the evaluation of tillage direction effect on soil attributes based on comparison between adjacent fields can be over or under estimated when the sampling coordinates and the spatial correlation among adjacent observations of data are ignored. The methodology used was able to differentiate between natural and management induced differences of soil attributes. Overall, the use of this methodology will improve the prediction and understanding of the effects of different cultivation practices on soil quality.


Introduction
Tillage types and management impact both soil properties (e.g. bulk density, organic carbon, structure and water storage capacity) and soil behaviour including infiltration, runoff and erosion (Bogunovica et al. 2018;Carretta et al. 2021). Combined tillage and water erosion are dominant contributors of soil loss and degradation in steeply sloping drylands (Wang et al. 2021). Liu et al. (2015) estimated a 10.6% loss of topsoil during 60 years of cultivation. Tillage erosion rate and redistribution are greatest when carried out up-and down-slope in steep hillslopes (Novara et al. 2019;Wang et al. 2021). In such situations, tillage influences soil physical properties, soil loss and crop yield (Bogunovica et al. 2018). Tillage can have an impact on the size and stability of soil aggregates and reduce soil organic carbon and nutrients as well (Catania et al. 2018). Mukherjee and Lal (2015) evaluated the impact of tillage type and duration on an organic soil compared to a mineral soil, while the longterm conventional tillage did not affect the organic soil, short-term conventional tillage impacted soil bulk density, water-stable aggregates, and available water capacity in the mineral soil. Catania et al. (2018) compared the effects of a spading machine, chisel plough, rotary tiller and no tillage (NT) on soil organic carbon (SOC), water stable aggregates and soil penetration resistance (PR). No tillage showed the highest PR and tillage with spading machine showed the lowest PR. SOC as well as water content and water stable aggregates were higher under spading machine compared with those under rotary tillage and chisel plough. In the study of Bogunovic et al. (2020), soil bulk density and PR were the lowest under reduced tillage compared with conventional tillage and minimum tillage, while minimum tillage created a loose topsoil and dense subsoil.
Tillage erosion also has an important role in the spatial and temporal variability of soil attributes due to soil redistribution (Malvezi et al. 2019). Accurate methods to assess the effects of tillage erosion and their causing factors in a given field, are still to be better evaluated. Most of these methods are expensive to be used, time-consuming and require private research stations. Several studies (e.g. Asadi et al. 2012;Mukherjee and Lal 2015;Bogunovic et al. 2020) analysed adjacent fields/plots to compare tillage types, land uses, or managements. The main problem with such comparative studies is the assumption of initial uniformity in space of the soil attribute in adjacent areas, disregarding possible differences among the spatial variability structure of each soil attribute. Therefore, methods should focus on saving time and money and provide an accurate measure or estimate a soil attribute at the same spatial point along the time or at the same time at different spatial points for different tillage types, managements and land uses.
The analysis of time series usually is applied to model the phenomenon under study based on observations of that phenomenon, to obtain statistical conclusions and to assess the ability of the model to forecast values of the study phenomenon (Chatfield 2003). While, time (or spatial) series analysis also has been used to model the variability of soil attributes using autoregressive, autoregressive state-space, and autoregressive moving average models (Heuvelink and Webster 2001;Aquino et al. 2015;Centeno et al. 2020), very few studies (Timm et al. 2006) have been carried out to forecast soil attributes using space series.
Taking into account that the spatial variability of soil attributes has been ignored in many comparative studies of tillage/management/land use changes, we present here a methodology based on time series analysis to evaluate the effects of different tillage changes on the magnitude of soil attributes. In this way, the objectives of this were to: (i) assess and compare the gross effect of tillage direction (up-down vs. contour) on the magnitude of soil attributes in two adjacent fields using the independent sample t test; (ii) analyze and model the trend and spatial variability of soil attributes using geostatistical tools and space series analysis in the field with contour tillage (COT); (iii) assess the potential transferability of those obtained COT models in forecasting soil attributes of the adjacent field of up-down tillage (UPT); and (iv) compare forecasted and observed UPT soil attributes to determine the net effect of tillage type.

Study area
The study was carried out on a sloping land in the Kouhin region, Qazvin province, Iran (36° 22′ N, 49° 35′ E 1500 m asl) ( Fig. 1a and b). Although presenting a very low mean annual precipitation and high potential evapotranspiration of 325 mm and 1200 mm, respectively, rain-fed agriculture is dominant in the region. The mean annual, minimum and maximum temperatures are 12.2, 6.4 and 18.1 °C, respectively (Iran Meteorological Organization 2014). Soil water content and temperature regimes of the study area are Xeric and Mesic (Keshavarzi et al. 2017). The soil of the experimental area was classified as an Inceptisol (silty loam textural class) according to the Soil Survey Staff (2014). The major cultivated crops are wheat, lentil and barley with crop rotation wheat-lentil-barley usually or wheat-fallow-barley depends on annual rainfall of the region. In dry years, a fallow is included in rotation.
The study area was a south face hillslope with average slope of 10%, and was divided in two sub-areas of about 5 ha each (Fig. 1c). In the eastern part of the study area is located a Research Station in which the contour tillage system (COT) (Fig. 1c) has been performed over more than 30 years. The western part of the study area belongs to local farmers, in which the plowing operations using moldboards have been performed up-down-slope (UDT) (Fig. 1c). Tillage was conventional and to the depth of 30 cm. Signs of sheet erosion by subsoil exposure in the upslope part and deposition in the foot-slope part are visible in the UDT field. The development of rills at lower parts of the slope together with poor plant cover indicates the occurrence of severe water and tillage erosion in the UDT field.

Soil sampling and analysis
To conduct the study, the middle position of the hillslope (10% of slope) was selected for the establishment of a perpendicular 400 m spatial transect to the main slope. Soil sampling was performed in July 2016 at 80 points (40 soil samples from each field), equidistantly 5 m along the transect ( Fig. 1c and d) from the 0-15 cm soil layer.
Soil sampling point, the following attributes were determined: particle size distributions (sand, silt and clay contents) were determined by hydrometer method (Gee and Or 2002); mean weight (MWD) and geometric weight (GMD) diameters of soil aggregates were determined by the wet sieving method (Nimmo and Perkins 2002); soil bulk density (ρ b ) was determined by the volumetric ring method (VR diameter 5 cm, length 5 cm) with three replications (Grossman and Reinsch 2002), while the saturated soil water content (SP) was determined gravimetrically (Topp and Ferre 2002); organic carbon (OC) content was determined by the Walkley and Black method (Bremmer and Mulvaney 1982); and calcium carbonate equivalent (CCE) was determined by the neutralization with hydrochloric acid method (Nelson 1982).

Exploratory data and statistical analyses
Exploratory data analyses were run for the data set of each field (COT and UDT) separately to calculate descriptive statistics including mean, minimum, maximum, standard deviation and coefficient of variation (CV) values. Dispersion diagrams were used to describe the spatial behavior of soil attributes under COT and UDT along the spatial transect. Data normality was evaluated using the Shapiro-Wilk test. The effect of tillage direction between COT and UDT on soil attributes was evaluated by the t test for independent samples using the SPSS 16.0 package (SPSS 2007). From this analysis, the gross difference between COT and UDT tillage directions was determined which includes eventual natural soil spatial variability.

Spatial analysis
To find the difference between random spatial (in this study tillage direction in COT and UDT fields) and structural spatial variations, experimental semi-variograms were calculated and adjusted to theoretical ones using the GS + software, version 5.1.1 (Gamma Design Software, Plainwell, MI, USA). Because of the choice of soil sampling through a spatial transect, data sets were considered to have an isotropic behavior. The spatial dependence degree (SD) of each data set was calculated as the ratio between the nugget variance and the sill [(C 0 /C 0 + C) × 100] and The geographical location of the study area in Qazvin province, c Aerial view of the experimental area; d The study area, sampling design as well as field with up-down tillage (UDT) and adjacent field with contour till-age (COT) sub-areas. X refers to the 10 sampling points in the UDT field used to compare forecasted and observed soil attributes based on models developed from COT field data classified as: SD ≤ 25%, strong spatial dependence degree; 25 < SD ≤ 75%, moderate spatial dependence degree; and SD > 75%, weak spatial dependence degree (Cambardella et al. 1994).

Space series analysis
The method is based on the assumption that the magnitude of soil attributes cultivated under COT would have the same magnitude if they were cultivated under UDT. Briefly, soil attribute measurements of 40 sampling points demarcated in the COT field were used to model and to forecast soil attributes at 10 adjacent locations of the UDT field. Then, forecasted and measured data in the UDT field were compared and the net effect of tillage direction was determined. The flowchart of the proposed method is presented in Fig. 2.
Below the proposed method is presented step by step in more details: Step one. Spatial distributions of all COT data sets were plotted using dispersion graphs. Static regression models of time series provide the basic framework for handling behavioral relationships, and are appropriate when the system is in equilibrium (Harvey 1990). We used the static regression models in this study. To fit a static model to the data, it is necessary to check if data series present or not a trend. The Mann-Kendall (MK) non-parametric test was used to check any monotonic trend in data (Drápela and Drápelová 2011). Mann-Kendall statistic S is defined as: where x j and x k are spatially sorted values and n is the number of samples (here, n=40). The sgn (x j -x k ) is equal to − 1 for (x j -x k ) < 0, is equal to zero for (x j -x k ) = 0, and is equal to 1 for (x j -x k ) > 0. Then the variance of S is computed, and MK test statistic (Z) is calculated as follows: The positive values of Z represent additive trend and negative values represent decreasing trend. The null hypothesis states that there is no trend is rejected if |Z|≥ Z 1-a/2 , where a is the pre-assigned significance level (Yue et al. 2002).
The Mann-Kendall test showed that there is a trend in OM and CCE data. The nonparametric Sens's method was used to estimate the slope of linear trend in OM and CCE data. The trend was removed by subtracting the Sens's estimate of each point from the measured data. The remaining data were considered as residual values (Chatfield 2003). The normality test was run for the new data of OM and CCE by the Shapiro-Wilk test again.
Step two. The data normality was checked for all data sets using the Shapiro-Wilk test at 5% significance level (p value = 0.05). Autocorrelation (ACF) and partial autocorrelation (PACF) functions were computed for each data set to determine the spatial correlation range and to identify the order of each model, respectively.
Based on ACF and PACF results, an Autoregressive (AR), Moving Average (MA) or Autoregressive Moving Average (ARMA) model was selected to describe the spatial variability behavior of each variable along the transect (Shumway and Stoffer 2016). The model parameters of each selected model were obtained using the SPSS 16.0 package (SPSS 2007). Forecasting accuracy of each model was evaluated by calculating the mean absolute percentage error (MAPE) and classified as: MAPE < 10%, highly accurate forecasting; 10 ≤ MAPE ≤ 20%, good forecasting; 20 ≤ MAPE ≤ 50%, reasonable forecasting; and MAPE > 50%, inaccurate forecasting (Lewis 1982).
Step three. The best model obtained for each soil attributes (Step two) under COT was used to forecast the value of that soil attributes in 10 adjacent points in UDT (Figs. 2d and 3) following the methodology proposed by Salas et al. (1980). The difference between forecasted and measured Comparison of the forecasted and measured data in UDT field using the Mann-Whitney U test

Contour tillage (COT) Up-down tillage (UDT)
Comparison of data between adjacent fields using independent T test and spatial variations analysis Space series analysis of the data, using the best obtained models to forecast the soil property values for the next ten points in the field with up and down tillage Comparison the result of the two tests (T test vs. U test) and finding the net effect of tillage type values of the soil attribute under study in the same 10 sampling points of the UDT field was considered the net effect of tillage direction. The estimation of normal distribution parameters, such as mean and especially variance, was a limitation due to the use of small sample size. Furthermore, it was a hard task to evaluate if the variable under study tended to the normal distribution or not. To overcome both limitations, the Mann-Whitney U test was used to compare mean values between forecasted and measured values of the soil property under evaluation in the UDT field. Mean values and results from the application of the independent samples t test for all UDT and COT measured soil data sets are presented in Table 1. Mean values of all soil attributes presented difference between the two tillage directions at p < 0.01, except for the CCE data set which presented difference at p < 0.05. But these differences are the gross difference between COT and UDT which include eventual natural soil spatial variability. In other words, the difference between COT and UDT fields presented in Table 1 may in part be due to spatial variability rather than the pure effect of tillage direction.

Exploratory data analysis and comparison between mean values of COT and UDT soil attributes
The clay content of the soil in the filed under UDT was about 24% higher than that under COT. Up-down tillage induces water and tillage erosion which has resulted in removing soil surface layer and the exposure of subsoil of high clay content.

Spatial variability of soil attributes
The spatial distribution of OM and CCE data sets under COT showed a trend using the Mann-Kendall test at 5% and 1%, respectively (Table 2). Both spatial trends were removed by calculating the residuals (see, Spatial analysis Statistical analysis of data, Step one). Figure 4a and b shows that in the mid-slope along the contour line spatial transect (from the right to the left sense, Fig. 1), the CCE percentage increases and the OM percentage decreases in the COT field, demonstrating a clear opposite trend behaviour through the transect. Therefore, the trend of each spatial data series was removed by subtracting each measured value at each point from each Sens's estimated value. CCE and OM residual data series are presented in Fig. 4.
The spherical model better described the spatial variability structure for seven out of nine soil attributes under COT, while the spherical model better described the spatial variability structure, for four out of nine soil attributes, under UDT (Table 3). UDT reduced both nugget effect (C 0 ) and sill (C 0 + C) parameters for sand, silt and ρ b theoretical semi-variograms as compared to those ones for the same soil attributes under COT. The opposite was found for those semi-variogram parameters for clay data. In the case of soil aggregate stability indices (MWD and GMD), while the nugget effect was increased under UDT as compared to that one under COT, the sill remained unaffected by tillage direction. For SP theoretical semi-variogram, UDT increased C 0 and decreased C 0 + C, while the opposite behaviour was found for those parameters under COT. Theoretical semivariograms for CCE and OM (Table 3) showed a pure nugget effect under UDT which means that up-down tillage adjacent observations of each soil property are not spatially correlated at our sampling scale. On the other hand, theoretical semi-variograms for CCE and OM residuals under COT manifested spatial variability structures which were better described by exponential and spherical models, respectively.
Soil textural fractions (sand, silt and clay contents), MWD, GMD and ρ b data sets showed moderate spatial dependency degree (25% < SD ≤ 75%, Cambardella et al. 1994) at both tillage directions. The SD for remained soil attributes was affected by the tillage direction. UDT reduced spatial dependency of SP, CCE and OM as compared with COT. The spatial dependency reduced form strong (SD = 4.6) to moderate (SD = 50) for SP, from moderate (SD = 33) to a pure nugget effect (SD = 100, i.e. without any degree of spatial dependency) for CCE, and from strong (SD = 19.5) to a pure nugget effect for OM (Table 3). The spatial dependency and structure of the soil attributes (Table 3) supports the hypothesis that all differences presented in Table 1 would not be just due to the effect of tillage direction.

Space series models of soil attributes
Autocorrelation (ACF) and partial autocorrelation (PACF) functions were computed for COT data to assess if  observations of each data set were spatially correlated and to select the model type and its respective order (Table 4).
As an example, Fig. 5 shows ACF and PACF plots for sand and SP percentages up to 16 lags. Sand observations showed significant spatial correlation up to 1 lag or 5 m in our study using a t test at 5% probability level (Fig. 5a). According to Nielsen and Wendroth (2003), the ACF of sand indicates a possibility of obtaining a spatial interpretation for sand data. ACF and PACF (Fig. 5b) behaviours for sand data set suggested an AR(1) model to describe its spatial variation along the transect. SP adjacent observations (Fig. 5c) showed significant spatial correlation up to 3 lags (15 m here) using the same t test at 5% probability level. From ACF (Fig. 5c) and PACF (Fig. 5d) results for SP data, an ARMA (3,2) model was used for describing SP data along the spatial transect. Similarly, AR(1) model was suggested for MWD and OMresidue series, ARMA(1,1) for silt series, and ARMA (2,2) for ρ b series. The ACF and PACF did not show significant autocorrelation at any lags for clay, GMD and CCE-residue series, therefore no models were used for these soil attributes. The models presented in Table 4 were used to forecast soil attributes in the UDT field for 10 first points of the transect. Based on the MAPE values (Table 4), forecasting accuracy was high for silt, ρ b , SP and OM, good for sand and reasonable for MWD. Figure 6 presents the spatial distribution behaviour for all soil property data sets measured under COT (40 points) and forecasted for UDT field (10 first points) using the best fitted time series model presented in Table 4. Mean values and the results from the application of the Mann-Whitney U test for all ten measured and forecasted data of each soil property in UDT field are presented in Table 5. In contrast with the results of the independent sample t test run for measured data of both two fields, in which all soil attributes were significantly different between those fields (Exploratory data analysis and comparison between mean values of COT and UDT soil attributes, Table 1), mean measured and forecasted values for sand and ρ b did not present significant difference at p < 0.05 (Table 5 and Fig. 5a and d). The forecasted mean value for silt content data was significantly higher than the measured one for UDT field (Table 5 and Fig. 5b). The reverse behaviour was found for mean values of silt content data when using all 40 observations from each field (Table 1 and Fig. 3b). Also, mean magnitude differences between forecasted and measured data were quite different for MWD, SP and OM (residuals). The MWD forecasted mean value (0.84 mm, Table 5) using the AR(1) model (Table 4) was about 56% higher than that one obtained using all 40 MWD observations under UDT (0.58 mm, Table 1). The results were quite similar for SP data, i.e., the mean difference was very low (45.1 vs. 46.4%) based on measured SP data as compared to that one (44.2 vs. 49.4%) based on SP forecasted data. Opposite results were obtained for OM data, in which the mean difference decreased from 45% (0.65 vs. 0.94%) to 29% (0.66 vs. 0.85%) when compared the modelled and measured data (Table 5) instead of comparing adjacent fields ( Table 1).

Comparison of different methods for OM
To give more insights into the results, further analysis was performed for the OM data. Figure 7 presents the way of predicting the net effects of tillage direction on OM by linear regression equations or the nonparametric Sen's slopes of linear trends. The content of SOM was estimated for each tillage direction at the fields' border by extrapolating using the equation fitted to the data of that tillage type. The results are presented in Table 6 alongside with the mean values estimated from various sample sizes. The mean value of OM under UDT was not affected by sample size, while its opposite behaviour was observed under COT. The results are vice versa in the case of modelling (Table 6), the forecasted OM values are almost the same for COT, but they are different for UDT. According to this results, if we compare the adjacent fields to find out the effect of tillage type on OM for example, there will be an overestimation of the magnitude effect of the tillage type on OM values which is affected also by the sample size. But the method presented in this paper (Space series models of soil attributes) is supported by the data presented in Fig. 7 and Table 6.

Discussion
We studied the effects of long-term tillage direction on soil attributes in two adjacent fields of a south-facing hillslope in a rain-fed area. The comparison was based on 40 soil samples taken along a contour line in each field, located in the middle of the slope. Based on the analysis by classical statistics (independent t test), significant differences were observed in all measured soil attributes under the two tillage directions (Table 1). Up-down tillage generally results in surface soil displacement and induces water erosion, both processes resulting in changes of soil attributes, which are responsible for soil quality and productivity reduction Statistical comparisons between adjacent fields/plots have been carried out in the literature neglecting eventual effects of soil spatial variability not allowing to distinguish the effect induced by the tillage type. Soil attributes may show different spatial variability structures, which depend also on the soil property and the scale. Strong spatial correlation and variability of soil attributes exists even at small scales (Malvezi et al. 2019). Kisic et al. (2018) also pointed out the existence of small-scale soil variability before establishing their experiments and expressed its effect on the results.
The results of the MK test (Table 2, Fig. 3) indicated a spatial trend for several soil attributes under both tillage directions (UDT and COT). Also, empirical semi-variogram modelling (Table 3) showed different spatial variability structure for different soil attributes under UDT and COT fields. The moderate spatial dependency for SP, and pure nugget effect for CCE and OM under UDT as compared to their strong, moderate and strong spatial dependencies, respectively, under COT, showed that up-down tillage may increase spatial uniformity of the soil. This can probably be related to the effect of scale and factors influencing the random component (such as tillage type) (Marzvan et al. 2015). Usually strong and weak SD of soil attributes can be attributed to inherent factors of soil (such as soil formation) and external factors (such as soil management methods, tillage, rotation, crop rotation, etc.), respectively (Yemefack et al. 2005).
According to the result of space series analysis for COT data (Table 4), the AR(1) and ARMA models were the most fitted model to forecast soil attributes in this study. Although time series models have not frequently been used in the literature to deal with space series to forecast data beyond the measured spatial domain. They have commonly been used to model time/spatial data series inside the spatial domain in which data were sampled (e.g. Heuvelink and Webster 2001;Timm et al. 2004). The best model obtained for each soil property was used to forecast its value for 10 consequent points in the UDT field. The forecasted values were considered as the net effect of tillage direction neglecting the spatial variability (Table 5). Analyzing Tables 1 and 5 together, it can be seen that as compared to COT, the gross effect of UDT (comparison between adjacent fields) was significant on all soil attributes (Table 1), while the net effect of UDT was not significant for sand content and ρ b forecasted observations (Table 1). The results were reverse for silt content in these two cases (Tables 1 and 5).
Mean clay content in the UDT field was significantly higher than that one in the COT adjacent field (Table 1) as well as higher than that one for those forecasted values ( Table 5). The combined effects of tillage and water erosion in the UDT field resulted in soil surface removal and the exposure of the clay enriched sub soil. An extensive study of the area (Asadi et al. 2017) showed that the clay content in the 15-30 cm soil layer was in average 5% higher than that one in the 0-15 cm layer.
Mean soil OM content in the UDT was lower than that ones for either the adjacent COT field or the modeled SOM values for COT at the same location. Soil OM content and nutrients usually increased by adopting COT in comparison with those ones measured in the UDT because of increased water infiltration (Pengfei et al. 2018), decreased soil erosion   (Kisic et al. 2017(Kisic et al. , 2018 and increased crop production (Quinton and Catt 2004). Nie et al. (2019) carried out a study emphasizing the effect of soil displacement and redistribution by tillage on the SOC dynamic in hillslopes. Soil aggregate stability indices (represented by MWD and GMD in our study) were higher in COT fields than those ones in the adjacent UDT field (Table 1). Values of both indices were forecasted to be higher under COT in the UDT field (Table 5). Despite the direct destructive effect of UDT by moldboard action on soil aggregates (Zheng et al. 2018), the difference in soil aggregate stability was also interpretable by the difference in soil degradation and erosion rates and SOM contents on the study soil layer in the two tillage types.

Conclusion
The comparison is usually made between adjacent lands/ fields/plots to study the effects of management changes on soil characteristics and is based on the assumption of initial uniformity of the soil in adjacent areas. We hypothesized that this assumption is not in agreement with the spatial variability of soil characteristics. The effect of UDT on soil attributes was compared with COT by forecasting soil attributes in the same field. The generated soil data of the UDT field were forecasted by the ARMA and MA models  developed by space series analysis of soil data of the COT field. The assessment of tillage type effect on soil attributes by space series analysis illustrates how we can realize and differentiate the natural soil spatial variability from soil management induced changes when comparing the adjacent fields. Two methods of analysis were used to compare the effect of tillage direction on soil attributes including (i) the classical t test and regression model in which the sampling coordinates were ignored and (ii) the space series analysis in which the sampling coordinates were not disregarded. The gross effect of tillage direction (up-down vs. contour) evaluated by t test showed that there were significant differences between the magnitude of mean values of all studied soil attributes at two adjacent fields of different tillage direction. However, this effect included natural soil spatial variability and may not be the pure effect of tillage direction.
The behaviour of the spatial variability pattern of each studied soil property was affected by the tillage direction; however, it was not only affected by the effect of tillage    Fig. 7 Estimation of the net difference in organic matter of two tillage directions by the Sen's slope of linear trend or by linear regression soil attributes showed that sand content and ρ b data were not significantly affected by tillage direction, and reverse result was observed for silt content.
Comparing the results of the two analysis methods for other soil attributes showed that the effect of tillage direction on MWD and SP was diminished, and on OM was exaggerated when the classical method was used as compared to when the sampling coordinates and the spatial correlation among adjacent observations were accounted in analysis. This means that the evaluation of tillage direction effect on soil attributes based on comparison between adjacent fields can be over or under estimated when the sampling coordinates and the spatial correlation among adjacent observations of data are ignored. Therefore, it can be concluded that some differences observed between the adjacent fields are due to spatial variability of soil attributes not due to tillage type effects.
The use of space series analysis showed that future studies at both adjacent fields should take into account the non-uniformity of soil attributes. Inaccurate estimates of the effects of tillage operations on soil attributes can mislead the engineers and farmers to choose appropriate management or planning.
Our study can improve the evaluation of experimental field data and provide a framework for precise prediction of the effects of different tillage types, cultivation managements, land use changes, and soil erosion on soil attributes.