Georeferenced tractor wheel slip data for prediction of spatial variability in soil physical properties

The upcoming technological breakthrough in the cropping system will offer a more detailed insight into soil-to-plant, man-to-soil, and man-to-plant impacts, thus improving the forecasting and ensuring more efficient in-field management. This study presents various on-the-go sensing procedures which were conducted in order to evaluate the quality of spatial estimations of soil physical properties such as soil compaction, soil moisture content, bulk density and texture. Standard statistical tools showed high positive correlations between soil specific resistance and soil compaction (R2 = .75), soil electromagnetic conductivity and moisture content (R2 = .72) and tractor wheel slip and soil compaction (R2 = .64). Variogram modeling of spatial autocorrelation gave the highest prediction error for tillage resistance (9.85%), followed by cone index (4.49%), moisture content (3.7%), bulk density (1.39%), clay + silt content (.98%), soil electromagnetic conductivity (.95%) and the least error was obtained for tractor wheel slip (.58%). The Central Composite Design (CCD) analysis confirmed significant contribution of soil compaction in the modeling of the specific soil resistance and tractor wheel slip, while soil moisture content and fine particle (clay + silt) content had a major impact on soil electromagnetic conductivity measurement. Soil bulk density had considerable importance in CCD modeling of tractor wheel slip.


Introduction
The efficiency of soil and plant management is strongly dependent on the method of data collection and data validity. Large scientific resources are being allocated for the development of farmer-friendly devices (with low price and high efficiency, easy to handle), algorithms and PC software that could be equally efficient in various agro-ecological environments. Site-specific crop production is designed to enable maximum productivity per unit of input and land area used. In order to succeed in all these challenges, field activities must be spatio-temporally adjustable considering the local soil characteristics of the field and specific requirements of the crop produced. Soil heterogeneity has been observed and confirmed in research (Hanquet et al., 2004;Hemmat et al., 2008;Lapen et al., 2002), but short-time soil dynamics complicate long-time spatial forecasting. Factors such as type of soil, climate, topography and applied in-field management represent spatial variability in the field. When the field size is enlarged by merging a couple of adjacent fields in to one entity, the overall field heterogeneity can be obtained (Oliver, 2010). In-row fertilizer placement for increased efficiency results in more variations in the nutrient content of soil (Hu et al., 2015).
Current site-specific management lacks in robust data acquisition techniques which must be a priori cheap, farmer friendly, quick and facilitate instant information availability. Proximal sensing has become very popular for its fast collection of geodata. In addition, the use of proximal sensors for field scouting has become increasingly popular as they are costaffordable to almost any farmer. Although the proximal sensor may seem a simple device, it can offer high-resolution reading with great potential for recognition of field hotspots. Numerous authors have used the systems for measurement of soil tillage resistance in order to determine the variability in soil physical conditions (Agüera et al., 2013;Fountas et al., 2013;Hemmat et al., 2013). Van Bergeijk et al. (2001) stated that soil tillage resistance could provide information about the locations of different types of soil. Lapen et al. (2002) discovered that if there is a relation between tillage resistance and soil compaction in the field, then the map of soil resistance can be used for delineation of compacted zones. Modeling variability in soil properties is very important for precision agriculture if the variability has certain correlation with the key parameters such as soil moisture, soil type and soil fertility. Numerous studies have been conducted with the aim of proving the usefulness of soil mechanical resistance measurement especially with data geo tracking (Fountas et al., 2013;Hemmat et al., 2008Hemmat et al., , 2013. Although the process of soil tillage resistance measurement implies disruption of natural soil condition, the values obtained are mostly related to subsoil layers, hence subsequently generated maps can provide certain insight into subsoil physical condition. Nowadays, with the precision agriculture (PA, "ISPA Official Definition of Precision Agriculture", 2019) concept, on-the-go soil electromagnetic conductivity sensors are most frequently used (Corwin & Lesch, 2005;Sudduth et al., 2001;Triantafilis & Lesch, 2005) as the indirect measurement of specific soil properties over large land areas can be done quickly. The application of electrical conductivity (ECa) probes on production fields could be very efficient in the determination of the variability of certain soil properties if the contributions of the other affecting soil properties to the ECa measurement are known or can be estimated (Sudduth et al., 2001). Also, measurements of apparent ECa using electromagnetic induction (EMI) give information related to soil physical properties and are broadly used to delineate management zones of yield potential (Rodrigues, Bramley, et al., 2015;Rodrigues, Ortiz-Monasterio, et al., 2015). Concerning the previous statement, Singh et al. (2016) found out in their research that the topographic position of the field should be considered while making the correlations of ECa with soil properties. Ma et al. (2011) gave temperature correction models for soil electrical conductivity measurement and they observed the best performance of the exponential model of Sheets and Hendrickx (1995) as corrected by Corwin and Lesch (2005). ECa readings were utilized for delineation of field zones with respect to chemical content in the soil. Some of the nutrients (Na + 2, Mg + 2, Mn + 2, Cu + 2, Ca + 2, Zn + 2, Fe + 2) were the key to the explanation of the ECa (R 2 > .90). The group of authors (Brevik, 2012;Lawes & Bramley, 2012;Rudolph et al. 2015) discovered that ECa is influenced by soil properties such as salinity, temperature, moisture content and bulk density, and the strength of these relations is closely linked to the local field conditions. Sensing of a single parameter of any environmental feature is not sufficient for the modeling due to the complex spatial structure and interactions. For example, Siqueira et al. (2014) used two different soil EM probes (Veris P3000, Veris Technologies, Salina, KS, USA and EMD38 Dual Dipole Version, Geonics, Ontario, Canada, 2005) to generate optimal soil sampling designs and improve the spatial estimation of soil resistance to penetration (PR) on the field studied.
The complexity of agriculture and limited potential for adoption of sophisticated technology by farmers impose the need for the development of robust but simple operating solutions such as global navigation satellite systems (GNSS) steering or yield monitoring which are well accepted since precision agriculture was established. Considering that the tractor is a machine engaged in the field mostly during the growing season, the idea of possible tracking of part of the tractor where soil acts in a predictive manner appeared as meaningful. Modern tractors are already equipped with much electronics mainly utilized for the engine, hydraulics or transmission control, but some of the generated data could be potentially useful for soil monitoring. Similar to continuous yield recording and subsequent field variability evaluation using georeferenced data from a grain elevator, the hypothesis was that tractor wheel slip might be related to soil layers which tractor tires run over, and/ or deeper layers that are disturbed by the primary tillage implements. Therefore, the aim of this research was to observe the spatial relationship between soil physical properties determined in-situ and on-the-go sensed parameters. The main contribution of this study is the introduction of the tractor drive wheel slip (DWS) in the spatial prediction of soil physical properties, which could be of interest to farmers in terms of low-cost sensing equipment and the simplicity of measurement procedure without any need for calibration.

Trial setup
The field experiment was conducted in the northern part of the Republic of Serbia (geographic location: 45°26′09″N, 19°37′47″E) on the calcic chernozem type of soil, which is characterized as highly fertile. The total porosity of this type of soil was about 50% with about 20% of macropores and 30% of micropores. The field surface was classified as flat without undulations or inclinations, covering an area of .8 ha and commonly used in the crop rotation of wheat-maize-soybean under conventional soil tillage management. A similar design was used by Ünal et al. (2020). Site-specific crop production is the concept relating mostly to fields where soil variability unambiguously exists, mostly on large fields. An extensive scientific research on such fields is quite difficult considering the overall resources (time, workforce, cost) and requirements. Therefore, soil physical conditions 1 3 were artificially induced over a small size field to compensate for a large field. Small field trials are more convenient for managing, especially if the observations require spatially dense sampling. For example, soil electromagnetic conductivity measurement with EM38-MK2 probe (which was used in this study) on large areas included side effects like instrumental drift (up to 2 mS m −1 ) as mentioned in very detailed research by Sudduth et al. (2001). Similar disadvantages concerning the effects of temperature on soil EC reading were investigated by Ma et al. (2011). They suggested the exponential model equation as the best for the correction of ECa readings, obtained within 3-50 °C range, to provide measurements referenced to 25 °C.
According to the various types of agronomic expertise related to the local area, the hypothesis about possible soil uniformity that could cause bad modeling of soil variability was accepted. To increase the scale of soil heterogeneity over the field, some soil treatments were conducted during the fall of 2015 and 2016 (Table 1). In September 2015, the field was divided into three equal plots (Fig. 1a) where different soil tillage methods were applied afterwards (deep plowing, disc harrowing and shallow cultivation). The same treatment was repeated in the following year to exclude potential effects of soil reconsolidation due to argillic pedoturbation during the vegetation period. From the very beginning, the objective was to make soil variability more controllable and spatially predictable (and not erratic).
Depth of disking ranged from 100 to 150 mm (D plot), depth of plowing was set to 250 mm (P plot), and soil chiseling was conducted at 180-200 mm (C plot)- Fig. 1a. Prior to soil core sampling and penetration resistance measurement as standard soil condition assessment procedures (Fig. 1b), the field was scanned with EM38-MK2 (Geonics Ltd., Mississauga, Ont, Canada) soil probe (Fig. 1e) to give preliminary insight in to existing soil variability and determine a range of grid sampling interval which would give a tolerable prediction error. The principle by Kerry et al. (2010) by which half the variogram range of an ancillary data such as elevation, areal images or ECa is potentially indicative of variables of interest and could be used as a satisfactory sampling distance. The configuration of the measurement equipment is shown in Fig. 2b. Next step included simultaneous measurement of the soil tillage resistance and tractor drive wheel slip (Figs. 1c,d), both with data georeferencing.
In 2017, physical conditions of soil were determined using standard georeferencing methods for each observed location. Trimble Juno SB with the modular system PathFinder ProXT (Trimble Inc., Sunnyvale, CA, USA, accuracy .5-1 m) was used for in-field navigation during 1 3 soil compaction measurement and soil sampling. Considering the limited resources, the aim was to optimize the efficiency of potentially available people-time-finance per square meter and also to collect enough data for unbiased field representation. A rectangular grid sampling design was applied (10 m × 10 m), which is generally widely recommended when little information on field variability is available as well as for easier data processing and mapping (Rossel & McBratney, 1998).

Methods of determination of soil properties
Soil compaction as a cone index (CI) was measured directly while soil particle content, water content (WC) and bulk density (BD) were determined from the soil cores. Soil CI measurement was carried out prior to soil sampling. Nine penetrations were carried out at each observed location, at points where soil samples were taken subsequently. CI measurement followed standard procedure ASAE S313.3 (2011), using the hand-held device Penetrologger 6.00 (Eijkelkamp Agrisearch Equipment, Giesbeek, The Netherlands) which had a 100 mm 2 cone tip angled at 60°. Steel cylinders of .0001 m 3 volume were used for sampling of undisturbed soil cores in accordance with ISO 11272. Soil profiles were observed in two layers, 0-150 mm and 150-300 mm, from which the cores were taken. The sampling rings were inserted in the middle of the layer with three replications. Bulk density was calculated on the dry basis. Gravimetric soil moisture content was determined by drying the samples at 100 °C in a dryer until a constant mass was reached. The pipette method was used to determine soil texture. Soil texture was classified according to the ISSS-International Society of Soil Science. The percentage of mechanical components of soil specified by the Department of Agriculture Natural Resources Conservation Service (USDA-NRCS, 2014) classified this soil as clay loam. In further analysis, the soil fractions were divided into two groups, FPC (clay + silt) and CPC (sand), with the aim of emphasizing the differences of selected soil indicators and obtain a more consistent relationship with the sensor-based parameters.

Specific soil resistance and wheel slip measurement
Specific soil resistance (SSR) and tractor drive wheel slip (DWS) data were obtained during the 30 passes using a standard agricultural tractor and a non-reversible moldboard plow (Fig. 2a). Working depth (≈ .25 m), working width (≈ 1.05 m) and working speed (1.8 ms −1 ) were kept constant. Godwin et al. (2007) stated that if the tillage tool with constant geometry and working parameters passes through the soil, the reaction force should be in correlation with the physical properties of soil. Working depth and width were checked at every pass in three replications. Some authors have concluded that uncontrolled variations in the working width of the first plowing body do not affect significantly the draft resistance, especially if it is a multi-furrow plow (van Bergeijk et al., 2001). Depth of plowing was limited by a land wheel for the rear plowing body and by the lower hitches on the tractor, where the mechanical fixing device was embedded to avoid uncontrolled movement of the hitches. Specific soil resistance was measured with an original measurement system which was introduced by Kostić et al. (2014). An incremental encoder was used for the measurement of theoretical speed of the left rear wheel which constantly ran out of furrow. In such a tractor set up, the left side wheels of the tractor are unballasted and more disposed to slipping. Hypothetically, variations in the furrow that is being tilled will affect wheel slip directly, so this measurement concept could be, indirectly, reasonably sensitive to changes in soil properties.
The theoretical speed of the tractor was calculated based on the number of counted pulses in the time sequence and the dynamic wheel diameter (1). The sensor-counter used gave 60 pulses per revolution.\text where S is the route length in m, t is time in s, N imp is number of impulses and D w dynamic wheel diameter in m.
Wheel slip was calculated based on the known value of the tractor theoretical speed (v t ) and the true ground speed (v a ) gathered from GNSS.
Data acquisition was undertaken with a high-resolution universal measurement module QuantumX device (HBM, Inc., model MX440A, Darmstadt, Germany). The speed and geographical position were determined by GNSS device Trimble TMX 2050 (RTK, Trimble Inc., Sunnyvale, CA, USA). Data from the sensors and GNSS were read at 10 Hz, which proved to be an adequate sampling frequency in the research by Kostić et al. (2014) and Hayhoe et al. (2002).

Comparison of SSR, WC and ECa with other soil physical properties
During the field measurement, signals from force transducers, incremental encoder and GNSS were recorded simultaneously on PC with data collection software Catman Easy (HBM, Inc., Version 3.4.1, Darmstadt, Germany). First, the data collected at the headland were eliminated because the results of the signal processing were extracted recordings for individual, adjacent passes. Considering the uncontrolled tractor speed deviations which could be influential, the measured soil resistance was recalculated with respect to reference speed according to Summers et al. (1986).
The first step in data processing was the reduction of amplitude fluctuations of the raw signal caused by the shearing of soil flow under the plowing body (Onwualu, 1998). In addition, the signal of SSR was down sampled to 1 Hz (in Catman Easy software). Newly created data recordings represent average tillage resistance at a length of 1/f int ≈ 1.9 m and a width of 1.05 m (average plow working width). In further steps, the SSR, WC and ECa data were extracted to make them comparable with low spatial density data of soil physical properties. Some SSR, WC and ECa data were selected on the map using the GIS program (ArcMap v10.5, ESRI, Redlands, CA, USA) by choosing points-data within a 5 m distance around each observed location. The mean values calculated from selected data points were later used for comparison with soil parameters taking into account the common geographical location.

Statistics and geostatistics
The effects of applied treatments on the plot differences concerning the soil physical condition were assessed by comparing the mean values of the observed soil properties. Basic statistical indicators were used in the following analysis. Duncan's test with a 95% level of confidence was applied to determine the statistically significant differences between data. Pearson's coefficient (R) was used to determine linearity between the observed parameters.
The influence of sample-based data on sensor-based data was analyzed by CCD (Central Composite Design) to obtain the corresponding mathematical models. The CCD was generated by "Design-Expert" software (version 11, Stat-Ease, Inc., Minneapolis, MN) and it was used to find the best inter-relationship of four sample-based data (CI, WC, BD and FPC) with sensor-based data (ECa, SSR and DWS). Thus, soil compaction (CI, MPa), moisture content (WC, %), bulk density (BD, Mg m −3 ) and soil fractions (FPC, %) were chosen as the independent variables shown in Table 2. The electrical conductivity (ECa, mS m −1 ), specific soil resistance (SSR, kN m −2 ) and tractor drive wheel slip (DWS, %) were dependent variables. The inter-relationship of the variables was determined by fitting a polynomial equation to data obtained from 57 experiments. Statistical analysis of the model was undertaken to evaluate the analysis of variance. The minimum and maximum range of investigated variables and the full experimental plan with respect to their actual and coded values are given in Table 2.
As in the papers by Cui (2010), Zhang et al. (2010) and Cui et al. (2016), a multiple regression analysis of the data was carried out and the polynomial equation defining the predicted response (Y) in terms of the independent variables (X 1 , X 2 , X 3 and X 4 ) was obtained.
where X 1 , X 2 , X 3 and X 4 are input variables, B 0 is a constant, B 1 , B 2 , B 3 and B 4 are linear coefficients, B 12 , B 13 , B 14 , B 23 , B 24 and B 34 are cross-product coefficients, and B 11 , and B 44 are squared coefficients. Combination of factors (such as X 1 X 2 ) represents an interaction between the individual factors in that term. The variogram was estimated using a geostatistical tool for spatial data modeling and recognition of their values of dependence. The process of variogram modeling for an individual variable was iterated until the best possible diagnostic outcomes were reached. Ordinary kriging and subsequently generated maps followed immediately after the variogram modeling. Data distribution was examined prior to spatial modeling for all observed parameters. A preliminary survey showed that the anisotropic variogram in the N-S direction fitted spatial data the best. The estimated semivariance was calculated following Goovaerts (1997).
where z(x i ) and z(x i + h) are true values of the variable Z at locations x i and x i + h, and N(h) is the number of compared pairs at lag distance h. The parameters used for describing the variogram prediction models are: C 0 -nugget variance, C + C 0 -sill variance, and a-range.
Suitable variogram models for the kriging interpolation were selected considering the cross-validation results of accuracy of the predicted values (Oliver & Webster, 2014).
Comparison of prediction quality was done by calculating the relative root mean square error (RRMSE, Fortes et al., 2015): where n is the number of observations, P i is the predicted value, Õ i is the calculated value and Õ i is the calculated mean value. The best method was the one which had the lowest RRMSE (Jamieson et al., 1991).
Geostatistical analysis and visualization of the results were done in ArcGIS software.  Strudley et al. (2008). Soil moisture content (WC) differed based on the statistical values (Table 3) where the highest average value (19.38%) was recorded on C plot with minimal data dispersion

Overall statistics of the field trial
(5.67%). On D plot, 18.43% of average moisture content was obtained, while the lowest WC (16.51%) and, at the same time, the highest data variability (7.91%) was recorded for P plot. Soil texture analysis showed small differences in data between the plots. Although data grouping was done (fine fraction = Clay + silt and coarse fraction = Sand), the results obtained suggest insignificant differences with respect to the applied treatment. Obviously, short-time tillage management could not induce rapid soil alterations, since the soil particle content refers to specific soil morphology which is influenced by the factors that have been intertwined in the evolution of soil (parent material, climate, topography, time, organisms, etc.). Naturally, different soil textures can be found in larger fields with distinctive topographic structures, but that was not the case in this study. Thus, it was expected that significant variability in soil texture would not be noted on a small and flat field. Table 3 also shows the results of the sensor-based data analysis relative to the plots. Highest mean ECa was obtained for C plot (136.12 mS m −1 ), followed by D plot (129.90 1 3 mS m −1 ), and the lowest value was recorded on P plot (126.05 mS m −1 ). Although P plot had the lowest ECa, its CV value was slightly higher than those of C and D plots. A brief look at ECa statistics (Table 3) might lead to a conclusion that there is a small difference between the mean values for the plots, but if all data ranges (97.68-155.01 mS m −1 ) were considered, there is obviously a much wider interval between all ECa data taken from the whole field, which is a sufficient internal variability for satisfactory modeling. Soil resistance to plowing (SSR) on the plots was moderately influenced by the applied soil tillage management so the values associated with separate plots have small inter-plot differences (Table 3). Duncan's test confirmed the differences between P plot (60.42 kN m −2 ) and the other two, while C (54.67 kN m −2 ) and D (55.25 kN m −2 ) plots were categorized into the same data group. Wheel slip data (DWS) had similar results to SSR regarding the plots, which was generally expected due to a direct relationship between these two phenomena, as shown by Kumar et al. (2016). P plot had a significantly higher average DWS (13.21%) compared to D and C plots whose average DWS values are quite close (10.58% and 10.85%, respectively). The skewness and kurtosis are within the range (± 2) proposed by Curran et al. (1996) which indicates normally distributed data. Furthermore, the indicators of soil physical condition, like sensed features (ECa, SSR and DWS), show quite similar distributions.
The following text presents an analysis of the relationship between soil physical properties and sensor-based values (Fig. 3). In this part, CPC was skipped due to its redundant contribution regarding the FPC parameter. These two variables have direct-inverse Fig. 3 The correlation matrix chart between soil properties and sensor-based variables proportions over space and provide direct-inverse results with no additional informative value.

Response surface analysis
The results of CCD experiments, where the effects of four independent variables (CI, WC, BD and FPC) on three dependent variables (ECa, SSR and DWS) were studied, are presented in Table A1 in the appendix. The ANOVA (Tables 4, 5 and 6) indicated significant model terms. Having excluded the insignificant terms (not counting those required to support the hierarchy, Table 4), multiple regression analysis of the experimental data gave the following polynomial Eq. (7): where ECa is electrical conductivity, while X 1 , X 2 , X 3 and X 4 are the coded values of the test variables representing soil compaction, moisture content, bulk density and soil fractions, respectively. The regression equation obtained from ANOVA indicates that the multiple regression coefficient R 2 is .83. The model can explain 83.53% variation of the response. The adjusted R 2 = .81 is also very high, to advocate for a high significance of the model (7). The model F value 30.42 implies that the model is significant. Furthermore, the interactive effects of X 1 X 3 , X 1 X 4 , X 2 X 3 and X 2 X 4 are significant.
The 3D response surfaces are presented in Figs. 4a and d where two of four independent variables are fixed to corresponding average coded values.
Again, after the omission of the insignificant terms (not counting those required to support the hierarchy, Table 5), multiple regression analysis of the experimental data gave the following polynomial Eq. (8): where SSR is specific soil resistance, while X 1 , X 2 , X 3 and X 4 are the coded values of the test variables representing soil compaction, moisture content, bulk density and soil fractions, respectively. The regression equation obtained from ANOVA indicates that the multiple regression coefficient R 2 is .79 and adjusted R 2 = .77. The model F value of 31.93 implies that the model is significant. Also, the interactive effects of X 2 X 3 and X 2 X 4 are significant.
The 3D response surfaces are presented in Fig. 4b and e where two out of four independent variables are fixed to corresponding average values.
Similarly, after the omission of insignificant terms (again, not counting those required to support the hierarchy, Table 6), multiple regression analysis of the experimental data gave the following polynomial Eq. (9): where DWS is tractor drive wheel slip, while X 1 , X 2 , X 3 and X 4 are the coded values of the test variables representing soil compaction, moisture content, bulk density and soil fractions, respectively. The regression equation obtained from ANOVA indicates that the multiple regression coefficient R 2 is .74 and adjusted R 2 = .71. The model F value of 29.17 implies that the model is significant. Also, the interactive effect of X 2 X 4 is significant.

Variogram analysis
Variogram analysis was conducted to observe the spatial dependence between the geo-referenced data of each variable. Experimental variograms obtained from the field data are shown in Fig. 5. The mathematical functions were selected (spherical, stable and Gaussian) for the best-fitting of mean values of semivariance on defined lags. The variogram parameters are listed in Table 7 where indicators of quality of the following kriging interpolations are provided in Table 8. The inputs such as lag distance, number of lags and type of model were chosen based on the results of cross-validation. This procedure is well described by Oliver and Webster (2014) who suggested cross-validation as the main rule of selection of the appropriate model. In the process of cross-validation, the measured values are compared with the predicted ones. Figure 6 shows raster images where the predicted values are associated with a color ramp in the form of spatial maps. The cross-validation results for all parameters considered N-E search direction. The CI, WC, BD, ECa and SSR maps show overlapping in spatial distribution of delineated zones with respect to the field trial setup. Figure 6a shows that the highly compacted areas are mostly concentrated in the middle part of the field in which M plot was established. In the case of C and D plot, the zones with "high" compaction are less present. Tebrügge and Diring (1999) observed that soil aggregates that emerge after plowing are less resistant to compaction in comparison with aggregates which form after vertical tillage. Bielders and Michels (2006) stated that the variations in soil compaction are caused by the variations of internal soil characteristics such as bulk density, porosity, cohesion, etc. Low WC zone concentrated along the middle part of the field and its distribution largely overlapped the M plot, while higher WC values were recorded for DH and CP plots (Fig. 6b). Although soil BD had relatively low variability, the central plot as well as the northern headland area showed delineation (Fig. 6c). The generated map for FPC shows that the zones with similar values were randomly placed lacking certain regularity (Fig. 6d). Soil ECa measurement showed lowest values on M plot, which also continued over to the D plot but with lower percentage. C plot had the highest ECa with relatively clear distinction in regards to the rest of the field (Fig. 6e). The SSR map presents distinctive areas where the highest predicted values are assigned to a central part of the field (Fig. 6f). Spatial variability of DWS depicted in Fig. 6g, could be characterized as low compared to CI, BD, ECa, and SSR. A good spatial characterization of DWS variability can be noted with respect to the plot disposition and soil properties determined, which suggests that the DWS measurement had significant influence on soil condition assessment.

Discussion
A great challenge in this research was to set up an experiment on a small field where soil physical condition would be sufficiently variable and aligned with the conditions typical of large fields. This (small) field area proved to be convenient in terms of time consumption  . 6 The maps of observed parameters 1 3 for data collection, avoiding possible interruption of measurement quality, respecting the technical requirements of the equipment used (temperature impact, battery capacity, sensor calibration, etc.), as well as the subsequent data processing. On the other hand, precision cropping mostly refers to the large fields and all scientific efforts must comply with real agricultural circumstances, otherwise the usability of the achievements could be diminished. According to Table 3, it can be noted that there is certain matching with regard to CV values of some parameters such as CI, SSR and DWS. Lower CV values are associated with higher mean value and vice versa. For instance, the P plot has higher mean values (CI, SSR and DWS) than D and C plots, which might be due to intensive soil layer inversion and structure interruption which have reduced the soil particle resistance to compaction during a growing season. The explanation could be that conventional tillage, with respect to the other two tillage systems, forms a consistent soil layer up to the tillage depth. Nevertheless, tillage systems and their effects were not the subject of this research so no further causation analysis was done. Figure 3 shows that the sensors were able to detect the majority of variations of certain soil properties. The best positive correlation (R 2 = .75) for all plots was achieved between CI and SSR which generally coincides with the results of Desbiolles et al. (1999), Arvidsson and Keller (2011) and Khalilian et al. (2014). Also, weaker but still confident positive correlation (R 2 = .72) was found between the ECa and WC (Bai et al., 2013). Electromagnetic soil conductivity was undoubtedly highly impacted by soil moisture content variation (Brevik, 2012) which was confirmed by the results of the test. In previous research, the authors mostly examined the influence of wheel slip on soil compaction (Davies et al., 1973;Mazeiro et al., 1986). This might be the first time to set up a trial where soil compaction was set as the influencing parameter on wheel slip during soil tillage and where a positive correlation (R 2 = .64) was obtained. Although BD parameter had small scale deviation (CV up to 2.49%, Table 3), moderate correlation with soil ECa (R 2 = .58) was obtained. In the research of Krajčo (2007), EM38 was operated in the vertical mode and it was not sensitive enough to measure the differences in soil bulk density. In this case, 49% of ECa measurements were influenced by fine soil particle content. Domsch and Giebel (2004) included weighted clay and silt contents and showed a closer connection with the EC25 (R 2 = .67). The results from Fig. 3 correspond to the results of Arvidsson and Keller (2011) who determined negative correlation between the SSR and WC. Increased soil moisture reduces the cohesive forces between clay particles; consequently, less energy is required for tillage as there is less resistance.
The overall CCD analysis gave multiple R 2 , which implies that the soil variables determined from the samples in their joint interaction have a relatively predictive influence on the sensor-based variables, and the fitted models have good performance ( adjusted R 2 DWS = .71; R SSR 2 = .77; R 2 ECa = .81 ). The values given in the Coefficient column (Tables 4-6) indicate the degree and type of correlation that exists between the explanatory and dependent variables. The CCD analysis provides slightly different insight into the relationship between the variables considering standard correlation analysis given in Fig. 4. A comparison of statistical weights of the sample-based parameters for the ECa model (Table 4) and correlations given in Fig. 4 shows certain compatibility. The WC parameter is depicted as most influential on ECa model (p < .0001), then FPC (p = .0011), both with positive response. Likewise, CCD model (7) as statistically significant (p = < .0001) predicted positive relationship between ECa and CI, as well as ECa and BD, while opposite relationships can be seen in Fig. 4. However, the coefficients of CI and BD in model (7) are not significant (Table 4); they are required to support the hierarchy, only. The real influence of CI and BD on ECa can be seen based on the coefficient of interactive effects of X 1 X 3 . CCD model (7) showed interactive contributions of the parameters and found high significance for all of them (Table 4), which suggests that complex interactions of soil parameters on ECa readings must be considered in the prediction task. CCD modeling of ECa nominates soil water content (WC) and fine particle content (FPC) as dominant factors. As for the SSR, the model (8) obtained emphasized CI parameter with respect to the others (p < .0001) with positive growth response, which coincides with the results in Fig. 3 (R 2 = .75). Other soil parameters were not defined as statistically significant. CCD model of SSR given in (8) shows high level of confidence (p < .0001). Interactive relationships between BD and WC, also WC and FPC were recognized as highly influential, which was generally expected regarding their multicollinearity. The CCD model (8) clearly showed that SSR is a reliable estimator of soil compaction. According to Table 6, CCD model (9) of DWS has high statistical significance (p < .0001). Positive and significant impacts were recorded for CI (p < .0001) and BD (.0409) which coincides with the results from Fig. 4. Interactive impact of explanatory variables FPC and WC on DWS was recognized. According to Table 6, DWS model (9) could be used as a reliable predictor of BD and CI. Table 7 shows almost zero value of nugget variance of the models for CI, BD, WC and FPC (sample-based soil parameters), similar to ECa and DWS (sensor-based parameters), which indicates that the variances between spatially close measurements are minimal, thus the process is continuous in space (Oliver, 2010). Zero nugget of the estimated variograms does not imply the quality of the fitted models (Oliver & Webster, 2014). The cause of SSR nugget variance might lie in the nature of soil-to-tool interaction. In fact, cracks and shear plains formed in the soil blocks during soil layer disturbance by the plow body induce high amplitude fluctuation of soil reaction force in a short time range (Lapen et al., 2001). Although signal filtering was done in order to reduce the "pseudo" variations of SSR and negative impact on the data modeling, small-scale heterogeneity of soil contributed to short-time changes of the sensor signal. Recorded SSR data varied significantly over small spatial distance and this is often described as the "noise effect" (Kerry et al., 2010). Also, Oliver (2010) argued that high nugget variance is a consequence of unsuitable sampling intervals. Spatial dependence ranged from 7.57 m distance for ECa up to 66 m for WC. Dependency range varied with respect to spatial heterogeneity of the observed parameters, measurement frequency and inputs of the variogram model such as lag interval and bin width. Oliver and Webster (2014) stated that if the correlation ranges are higher than the minimal sampling distance, the fitted variograms could be used for spatial prediction o F value s at the unsampled locations. The selection of a suitable model for individual variables was done according to the RRMSE diagnostic parameter of cross-validation (Table 8). The prediction quality was evaluated by using the RRMSE parameter which expresses relative error of the model. As Jamieson et al. (1991) stated (as cited in Fortes et al., 2015), the validation is considered to be excellent when the RRMSE is < 10%, good if the RRMSE is between 10 and 20%, acceptable if the RRMSE is between 20 and 30%, and poor if it is > 30%. Following abovementioned, excellent prediction was made for models DWS (.58%), ECa (.95%), FPC (.98%), BD (1.39%), WC (3.70%), CI (4.49%), followed by the SSR (9.85%) with a good quality of prediction. This led to the conclusion that ECa and DWS had the best fit, which was important for subsequent analysis of kriged maps. Table 8 shows that the parameter root-mean-square standardized error (RMSStE) has positive value (close to 1) for almost all soil properties except CI (1.68) which confirms the objectivity of the models (Johnson et al., 2005). According to RMSStE, the most objective model was obtained for ECa (.99) followed by the SSR parameter (1.005) and DWS (1.07).
In the study, all observed parameters were more or less spatially correlated and had spatial continuity. However, particular interest was given to the study of a tractor wheel slip in soil physical condition modeling since it has never been researched before in the spatial context. By using variogram modeling, DWS provided an evident advance in the precision of data fitting comparing to other sensor-based entities (Fig. 5). DWS experimental data are gradually dispersed around the variogram line which defines DWS as convenient for soil condition assessment.
With regards to spatial distribution of soil properties (Fig. 6), the highest complementarity with SSR could be noted in the case of CI. The spatial structure of predicted maps suggests certain matching of SSR with respect to the parameter WC. The maps for SSR and WC indicated an opposite constellation. These findings have been provided by many other authors such as Arvidsson and Keller (2011). The color distribution on the ECa map confirmed the matching with WC map colors, which supports the confident relationship obtained through regression analysis (R 2 = .72). DWS map depicts a wide range of predicted values, from 6 to 16.8%, where M plot covered the majority with the highest DWS values, obviously very similar to spatial maps of CI and SSR. This fact points to the conclusion that DWS could be a valuable estimator of soil physical properties (like SSR or ECa) yet with a significant difference in the price of measurement equipment and much more farmer-friendly in terms of data collection and processing. DWS measurement could be performed with a low-cost sensor and GNSS assistance for measurement of true speed which has been already used by many farmers even in Serbia with the tendency to be fully accepted. Renowned tractor brands have an option to measure wheel slip during field operation, thus, they do not require any additional costs relating to DWS mapping. Data can be logged in real time domain using special software and subsequently mapped. The range of DWS data depends on soil condition, type of tillage, and quality of implement adjustment. In theory, if a tractor pulls an implement at constant speed and depth with no additional adjustment, wheel slip would be in correlation to soil condition. Also, DWS can be influenced by tire wear. Bad tires, especially on saturated soil, will block wheel grooves and restrict attraction force between soil and tractor wheels. In this situation, wheel slip would be mostly caused by the surface condition and the obtained wheel slip data would not be valuable for the prediction of soil condition up to the tillage depth. The greatest disadvantage of wheel slip measurement for the assessment of soil condition is that it must be conducted during soil tillage, which means that the initial soil condition determined at the moment of measuring will be changed.
The given relationships could not be characterized as temporally or spatially constant due to the complexity of soil structure and interaction with other factors such as farmer field management, climate, plants behavior, etc. Further examinations should be encouraged under the principles of precision agriculture field management, especially regarding the utilization of wheel slip as an alternative parameter for soil physical condition assessment.

Conclusions
The sensor parameters obtained correlated with some but not all soil properties. Statistical analysis showed the strongest association between soil resistance and cone penetration resistance (75%), coinciding with the principles of soil deformation. Soil electromagnetic conduction was highly sensitive to soil moisture content with 72% level of explanation. Soil compaction had similar effects on the wheel slip (R 2 = .64) as it did on soil tillage resistance.
Geostatistics enabled visual perspective of in-field variability and linear modeling of the influence of soil parameters on the sensed parameters (SSR, ECa and DWS).
CCD analysis offered objective and reliable models with an explanatory level of adjusted R 2 DWS = .71; R 2 SSR = .77; R 2 ECa = .81 and confirmed that CI dominated in SSR and DWS reading, while WC and FPC mostly impacted the measurement of soil ECa. Furthermore, according to ordinary least squares analysis, soil bulk density had a significant contribution to the DWS parameter.
The overall conclusion would be that tractor wheel slip data could be used as an alternative indicator for soil physical condition if detected during soil tillage under normal soil conditions It was assumed that the wheel slip variations during tillage were related to the heterogeneity of the shallow soil layer (which influenced wheel traction efficiency) and variations in deeper soil layers (which contributed mainly to the soil resistance to tillage). The results suggest that wheel slip data are more convenient for variogram modeling of subsequent map predictions than soil electromagnetic conductivity and soil tillage resistance. Low price of the sensor device and simple setup procedure make it widely available to almost every farmer regardless of the farm size. Wheel slip does not change rapidly like soil resistance so data processing does not require a complicated procedure of signal filtering, which makes the information necessary for mapping almost instantly available.