Spatial variability of soil physical properties in longitudinal proﬁ les

: The spatial variability of physical properties, such as bulk density, penetration resistance and gravimetric moisture, obtained by applying geostatistics in precision agriculture, can effectively indicate the physical behavior of agricultural soils in longitudinal profi les. In this way, the spatial dependence of physical properties in streets of coffee plantations with different lengths was evaluated in the southern Minas Gerais, Brazil. For this purpose, fi ve longitudinal profi les were measured in streets, each one with depths ranging from 0 to 0.60 m, in six layers of 0.10 m, being the database composed of 432 property, 144 by property, submitted to the ordinary kriging geostatistical method in order to obtain spatial variability maps using the R software. They were evaluated by the lower mean cross-validation error of theoretical models fi tted by ordinary least squares (OLS), being detected in higher superfi cial layers, from 0 to 0.30 m, lower bulk density and lower penetration resistance, with variable gravimetric moisture in the length direction of some streets of coffee plantations, being that these properties presented different structures of spatial dependence for each street.


INTRODUCTION
Soil water status is a primordial parameter for several agricultural operations, either for the moisture control in irrigation, to data of condition penetration resistance composing indicators to evaluate the physical quality of soils for crop production, as calculated by the indicator named "least limiting water range" used by Imhoff et al. (2016), affi rming that the type and magnitude of soil deformation depend on external factors that determine the applied stress, as well as on soil physical and mechanical properties, of which water content exerts great infl uence.It can be noted by soil moisture with a strong infl uence on possible changes in the soil physical structure, infl uencing the productive potential of some crops due to a higher or lower compaction state, in this case, as seen by the variability mapping of properties capable of characterizing a given soil state, proposed to evaluate the physical quality of soils, such as density and the penetration resistance values in an inverse relationship with the amount of porous space, related to root development capacity, considering soil moisture as dependent variables in resistance sampling through penetrometer (Resende et al. 2014).
The use of agricultural equipment without observing the infl uence of their considerable weights can contribute to specific soil compaction, negatively influencing crops perceived by the reduced production (Araújo et al. 2011), with consequent economic impact.Thus, the incoherent agricultural management can contribute to the increase of soil density and soil penetration resistance, as well as management activities in soil tillage, identified as the most important activity in the physical quality of soils (Stefanoski et al. 2013).
The interventions in soil use promote changes in their physical attributes, such as increased soil density, decreased total porosity and organic matter content, confirming the importance of knowledge of spatial variability in order to minimize management errors, according to Oliveira et al. (2013).
Higher investments in agricultural mechanization without adequate management may contribute to soil physical degradation, since there is a considerable increase in traffic of machines and implements, making the soil more compact, perceived through soil density analysis and penetration resistance.Thus, the higher the soil density, the lower the water content, resulting in higher penetration resistance (Palma et al. 2013).In this context, several forms of management seek to improve the soil physical qualities, promoting the production quality of crops suitable for agribusiness.
Soils are naturally heterogeneous in depth and length, due to factors as climate, organic matter, relief and formation, significantly disturbed through human activities, up to 25 cm depending on texture for the anthropic epipedon, according to Bockheim (2014).Therefore, soils are continuously variable in the space or, as mentioned by Cruz et al. (2011), soil variability occurs due to natural and anthropogenic factors acting at various spatial and temporal scales.Also, the study and mapping of the spatial variability of soil properties are important to guide sampling and interpretation of exact results, making it possible to verify and evaluate the spatial dependence of these properties and parameters aiming to the estimation of values at unmeasured locations, in accordance with Elleithy et al. (2015), citing that kriging (geostatistic method) provides the best linear unbiased estimation for unmeasured locations, with a minimum mean interpolation error, using the model semivariogram.
The soil density is an important physical property to evaluate the quality of agricultural soils, directly related to the structure and porosity, often used as an indicator parameter of soil degradation or conservation processes due to management or, as mentioned by Klein & Libardi (2002), by a structural physical condition, soil density can be affected by mechanical forces originating from the pressure caused by wheels of the agricultural machinery and by the very action of implements, reflecting important physical and hydraulic properties, such as aeration porosity, retention and water availability in the soil-plant system, particularly correlated with the mechanical resistance to penetration.
The soil penetration resistance may be indicative of compaction, in direct relationship to the development of crops, influenced by a particular management type, being reflected in the root growth by pressure between particles or aggregates, favoring or limiting the root growth in length and diameter, sometimes precluding shoot growth of plants, being considered crucial in the evaluation of effects from soil management (Tormena et al. 2002), According to Moraes et al. (2013), bulk density and soil penetration resistance have been widely used as indicators of soil compaction, where penetration resistance describes the mechanical resistance provided by the soil against something moving through it, as roots or a tillage tool, showing high correlations with bulk density.
Moisture is related to soil and cropping processes, such as compaction to determine the effects of field traffic operations (Barik et al. 2014).Therefore, its spatial variability is important to determine soil quality and crop development, being that the amount of water retained in the soil is unstable and variable according to the surface recharge that occurs through irrigation (Reichert 2011).
Bulk density, penetration resistance and moisture are variable physical properties in depths according to the characteristics of each soil, type of mechanization, management, climate and, in this way, the spatial variability of these properties can be investigated through the foundations of precision agriculture, in a localized and coherent way, through geostatistical tools.
To this end, according to Marques Júnior & Corá (1998), detailed knowledge of causes and factors controlling the soil-plant system is important to implement precision agriculture.
The present study used the spherical model to determine the variability of soil physical properties, which is justified because it evidences an increasing spatial correlation structure with the distance up to a certain point, then the semivariance becomes constant, limiting the sample space influence, widely used for studies related to soil and coffee crop (Carvalho et al. 2013, Ferraz et al. 2013, Reza et al. 2016, Mali et al. 2016, Rincon et al. 2017, Santos et al. 2017).
With an approach directed to the representation of the vertical physical complexity of the soil by spatial variability profiles, this paper aimed to develop a geostatistical study by ordinary kriging with Ordinary Least Squares (OLS) fitting method aggregated to the visual recognition of the spherical model of semivariograms of the bulk density, penetration resistance and gravimetric moisture with the purpose of verifying and indicating the quality of agricultural soils for the coffee crop.The annual evaporation rate of the region is on average, 1.000 mm per year (INMET 2010).The soil of the area was classified as clayey dystrophic red latosol (Ferraz et al. 2017, Jacintho et al. 2017).

MATERIALS AND METHODS
This grid was determined by the standardized distribution of alignments between sampling points in the transitional range between plants in the length direction of five coffee streets, dimensioned as longitudinal profiles with a fixed depth of 0.60 m and length in meters varying according to the size of the five coffee streets (144.87 m, 182.64 m, 244.26 m, 310.74 m and 353.79 m).The sizing of the five rectangular grids (profiles) was possible by assigning horizontal scales H 1:1 and vertical scales V 100:1 in the statistical software R. The midpoint of the vertical spacing of 0.10 m per sample was the location of the sampling point of soil sublayers at 0.05 m, 0.15 m, 0.25 m, 0.35 m, 0.45 m and 0.55 m, being the total depth interpolated from zero to 0.60 m.In this way, the study area consisted by the distribution of sample points in five alignments or coffee streets, sized for five rectangles (longitudinal profiles) according to the length of each street at a fixed depth of 0.60 m, according to Figure 1.
In order to determine the bulk density (BD), the volumetric ring method was used with the aid of an Uhland sampler, also used to determine the gravimetric moisture (GM), according to Reichardt (1985).The sampling of the property soil mechanical penetration resistance (PR) was determined by the use of PenetroLOG digital penetrometer from Falker, according to ASAE S.313.3 standards (American Society Of Agricultural Engineers 1999).For the moisture property, the gravimetric method was used with systematic measurement of sample masses.
Bulk density "BD", soil mechanical penetration resistance "PR" and gravimetric moisture "GM" were convened, being such properties grouped by streets and their respective depths.For this purpose, digital spreadsheets capable of providing the horizontal distance between the points sampled geographically by alignments in six layers were used, configuring longitudinal profile given by the length of streets versus sampling depth of 0.60 m, represented in Figure 1.
For the geostatistical processing and analyses, the R software was used, The Comprehensive R Archive Network (CRAN), using the geoR package developed by Ribeiro Júnior & Diggle (2001).Together with the mentioned package, other algorithms were developed and ordered in five scripts referring to the five coffee streets, according to the cited methodologies.Similarly, due to the large volume of geographical information, the maps interpolated in R software were recorded and resampled through the geographic information system (GIS) QGIS by algorithms specifically developed for georeferencing and export of krigged results in raster matrix format for this latter.This enables to integrate future information in the GIS that contributes to greater clarification of studied variabilities in order to consolidate geographical information and analyses by consciously use of interoperability provided by such processing platforms in a secure geostatistical computing environment integrated with the dynamics from scientific research processes or incoherent agricultural management, responsibly, in time and space, as in cases where the time factor is essential in agricultural productivity estimates.
The exploratory analysis of the data is important for the development of geostatistical processes in order to identify outlier data, to verify the distribution normality, frequency and variation of the studied properties, being that outlier data are verified in some cases by sampling errors (Grego et al. 2014).
During the preliminary study phase, correlation and trend analyses can be used for an initial understanding of the spatial variation of the properties.Firstly, the graphical representation of the sample frequency distributions allows evaluating a likely theoretical distribution model, to calculate the values of central tendency, dispersion analysis as well as to characterize some zoning type, besides the possible identification of sparse or spurious data as well, as in the descriptive analysis the coefficient of variation "CV", used by Addis et al. (2015) in selected soil attributes under agricultural land use system.Moreover, in the descriptive analysis, the coefficient of variation (CV) expressed in a standardized form CV, in percentage, can indicate the variable's dispersion, being that coefficients of variation above 30% are classified as very high (Gomes 1990), indicating a higher degree of dispersion.
The use of geostatistics to analyze the spatial variability of soil properties investigates the spatial dependence of variables through spatial variability maps, using routines, as suggested by Oliver & Webster (2015), for computing and modelling variograms and kriging.These maps are obtained by application of the kriging technique, by interpolation of variables with non-tendentious estimator and minimum variance, with semivariogram as the main component for the geostatistical validation, which is explained by the difficulty of conventional statistics in considering the spatial aspects of the phenomena.Similarly, the regionalized variable is considered as a basic element of geostatistics, as demonstrated by Journel & Huijbregts (1991), since in order to make parameters applicable in mutual dependence, the spatial stationarity hypothesis is assumed.Thus, the regionalized variables are considered by the same phenomenon that distinguishes them geographically, through the semivariogram interpretation for each value "z(x i )", as a particular realization of a random variable "Z(x i )" at the point "x i ", as mentioned by Journel & Huijbregts (1991), with an estimator "2γ*(h)" as arithmetic mean of the squared differences between the experimental measures "[(z(x i ), z(x i + h)]", as the equation used by Araújo el al. (2018): Where, N(h) is the number of pairs of two experimental measures "[(Z(x i ), Z(x i + h)]", separated by a vector "h", which generates the graph of "γ*(h)" as a function of a vector with three-dimensional coordinates (h u , h υ , h ѡ ) denoted "h" (Journel & Huijbregts 1991).
Once the experimental semivariograms were determined as a function of several "h" vectors, an initial visual fitting of the theoretical model to these semivariances was performed by "eye fit", a visual inference of the theoretical spherical model applied to semivariograms, with the definition of the parameters nugget effect, contribution and range, pursuing conformation that best defines the theoretical model on the regionalized variables of the semivariogram, in order to then perform the fitting of these models by the ordinary least squares (OLS) method, when the parameters that define the entire semivariogram are observed primarily.The models were defined and fitted through three indicators: the Spatial Dependence Degree (SPD), the mean cross-validation error (ME) and the visual recognition of the spherical model structure in its relationship with semivariances, preferring an optimum sill in order to make the decision more assertive in the choice of models, directly influencing the positional distribution of interpolated variabilities.
Furthermore, the spherical theoretical model was used to detect the spatial dependence for all streets and properties, based on semivariogram analysis, which according to Isaaks & Srivastava (1989) is a model of high spatial continuity and less erratic over short distances.
In regression kriging, the deterministic part is fitted first using ordinary least squares (OLS), and then the residuals are estimated by a method of interpolation such as Ordinary Kriging (Hapca et al. 2015).The OLS fitting method consists in obtaining a theoretical model that best represents the phenomenon or physical property with minimal variance, pursuing the parameter values of this model designed on the semivariances, in order to minimize the sum of squares of the difference between the observed and estimated values, by minimizing the vector expression of the parameters that completely defines the semivariogram (Mello et al. 2005), maximum kriging distance, type of estimator, trend degree, nugget effect, sill, range and distribution form of weights, being that the development of such method favors the attenuation of subjectivities on the estimation of parameters, visually, as stated by Cressie (1993).
In order to evaluate the best theoretical model and its parameters for later kriging, the cross-validation technique might be used (Yao et al. 2013), indicating a quality in the ordinary kriging situation of data, also observed by the theoretical parameters C0 (nugget effect), contribution (C1) and "a" (range) by determining the spatial dependence (SPD), measure proposed by Biondi et al (1994), in which contribution and sill are related, following the classification used by Cambardella et al. (1994), called here as spatial dependence degree.The relationship between parameters of semivariogram theoretical model for SPD determination was given by: Where is the sill.
In the search for assertiveness of processes that involve geostatistical modeling with a view to kriging, the elaboration of variability maps and its relationship with quality indicators, such as SPD and ME by cross-validation, were investigated, followed by visual recognition of the spherical model after OLS fitting.The SPD is perceived only by parameters of the theoretical model, while the latter is based on the positional simulation of possible consistencies between the sampled values and those interpolated by simple difference.Such comparison indicates a possible quality of the final theoretical models, which would corroborate a better estimation situation, resulting in a variability map with a certain degree of accuracy for each profile and properties.The visual decisionmaking is important in the choice of estimation windows by visual recognition of the empirical semivariogram, as mentioned by Journel & Huijbregts (1991), after OLS fitting in this case.
O n ce t h e t h e o re t i cal m o d e l s o f semivariograms were determined, the OLS fitting was performed to determine which of these fitted models would imply the lowest mean cross-validation error, being the latter applied by the temporary withdrawal of variables in the sampled sites Z(s i ) of the predicted kriging data set in the site Ẑ(s i ), suggesting an error by comparing the differences between measured and predicted values, returning, as response, the mean error "ME" and others, such as the standard deviation of mean errors, the reduced mean error and the standard deviation of reduced mean error, being the ME given by the equation 3, also used by Hengl et al. (2015) as an evaluation performance criteria of the predictive models: Where, "n" is the number of data, "Z(s i )" is the value observed in the point "s i " and "Ẑ(s i )" is the value predicted by ordinary kriging in the point "s i ".
Regarding the map development through longitudinal profiles, it was observed after performed kriging by software R that the variability of properties showed in the maps is as great as the number of classes that define it, i.e., during map generation stage, the number of interpolation classes to be represented by the variability map must be defined through captions, in a proportionable to represent the total amplitude of variables interpolated by kriging, when a smaller number of classes would result in discrepant maps, as states Whelan & Taylor (2013) showing that the number of classes affects the way of variability is categorized.Therefore, nine interpolation classes were defined for all streets and properties in search for a coherent representation of the variability of these latter.It was observed in R a limitation in relation to the presentation of values from the classes of estimated variables resulting from kriging, when it was possible to observe a greater number of color classes representing variabilities in the generated maps in relation to their respective estimated values, sometimes omitting values of interpolated classes that could be determinant in the application of agricultural practices accurately, as by the geographical application of possible corrections aiming at better agricultural development, by specific spot elevations in the field.
After kriging, the estimated results of each property and street were exported to the QGIS software, so that it was possible to represent and quantify the interpolated total amplitude, in order to indicate compaction levels for each of the properties observed in coffee streets through longitudinal profiles of the soil at maximum depth of 0.60 m.

RESULTS AND DISCUSSION
The relation between the standard deviation and the arithmetic mean gives the coefficient of variation (CV), which expresses the degree of dispersion of the variables.According to Kalil (1977), the main quality of CV as a measure of dispersion is the ability to compare different results that involve the same response variable, used by Cortés-D et al. (2016), allowing quantifying an initial behavior, as seen in Table I.
According to Gomes (1990), CVs below 10% are classified as low, between 10 and 20% as medium, between 20 and 30% as high, and very high above 30%, i.e., lower values represent more homogeneous properties and higher values suggest more heterogeneous variables.For this classification, the soil density "BD" for the five streets shows average values of dispersion with CV from 10 to 20%.For data regarding soil mechanical resistance to penetration "RP", it is verified very high CV (higher than 30%), indicating a high degree of dispersion around their averages, different from that observed for the values of gravimetric moisture "M" with CV from 10 to 20%, i.e., medium degree of dispersion, according to Table I, nearly of moisture CV values evaluated by Zucco et al. (2014) investigating the influence of land use on soil moisture variability, with average value equal to 27%.The medium and high dispersion values observed by CV of the properties SD, RP and M are presented in Table I, making it difficult to estimate parameters in the exact consistency of spherical models to the semivariance of these properties (Figures 2  and 3), evidencing models tending to the pure nugget effect.
For the property gravimetric moisture "GM", perhaps due to the high variance when in relation to other properties (BD and PR, Table I), it was noted a lower amplitude variation of the parameters of spherical models from these properties, for all the streets, when minimal variations of model parameters from this property resulted in minimal differences in empirical semivariograms after the OLS fitting.This fact was not verified for the properties BD and PR in the modeling phase, when minimal variations of semivariogram parameters resulted in a wide variety of empirical models, some tending to pure nugget effect and others more consistent with the expectation of spherical models to the regionalized variables, making it difficult to reach an optimum sill, often tending to pure nugget effect (Figures 2 and 3).
The results of Table I were compared to the penetration resistance classes adapted from Soil Survey Staff (1993), qualifying the penetration resistance lower than 0.01 MPa as extremely low, from 0.01 to 0.1 MPa as very low, 0.1 to 1.0 MPa as low, 1.0 to 2.0 MPa as moderate, 2.0 to 4.0 MPa as high, 4.0 to 8.0 MPa as very high, and above 8.0 MPa as extremely high.According to this classification, the mean values of the observations indicate very high penetration resistance (4.0 to 8.0 MPa) for streets one, four and five, and extremely high RP (above 8.0 Mpa) for streets two and three.According to the value for soil penetration resistance of 2.0 MPa significantly restricts the root growth of crops under conventional tillage systems (Taylor et al. 1966, Arshad et al. 1996), thus indicating possible adaptation of the agricultural management, Table I. Results of the descriptive analysis for BD (g cm-3), PR (MPa) and GM (%), for the five coffee streets.Thus, it can be suggested that the maximum values of soil penetration resistance (Table I), found in the five streets under study, should not cause marked losses to the growth of coffee plantation in this farm.

Property
The soil density in the variability maps, as an indicator of the compaction degree, was higher among the last layers from 0.30 to 0.60 m, similar to that found by Alves et al. (2007), at a maximum depth of 40 cm in a dystrophic red latosol of sandy clay loam texture, when observed higher density values of in these last layers from 0.20 to 0.40 m from the half depth studied by those authors, which can be associated with natural processes of soil or arising from the traffic of machines, or even by the increased work intensity of agricultural implements.Therefore, the changes of these physical parameters due to the soil and crop management system may imply physical alterations that determine ideal or limiting conditions for crop development (Collares et al. 2008).Moreover, high values of soil density up to the 30 cm depth can affect a certain coffee cultivation in full development, since the majority of the root system under these conditions is in the upper layers from 0 to 30 cm, which are important in the extraction of water in a coffee plant (Ronchi et al. 2015), when is observed medium to low BD in the studied coffee streets favorable to root development, for all the streets up to 30 cm depth.Moreover, according to Freddi et al. (2009), densities from 1.18 g.cm -3 in typical dystrophic red latosol of medium texture, indicate limited development of some crops, but only in cases where the soil is below an optimal water range.
In soil compaction studies, it should be noted that penetration resistance is dependent on soil moisture, thus small changes in gravimetric moisture would reflect great variations for penetration resistance (Campos et al. 2013).
Table II presents the results of the geostatistical analysis by the semivariogram parameters, evidencing the mean error (ME) and the spatial dependence degree (SPD) used by Cambardella et al. (1994), ranging from 0 to 1, in which values lower than 0.25 suggest a strong spatial dependence, from 0.25 to 0.75 a moderate spatial dependence, and values greater than 0.75 indicate a weak spatial dependence.According to the parameters of the theoretical models of experimental semivariograms of the properties BD, PR and GM, in order to determine the spatial dependence degree, a strong SPD was verified for all the properties from the street 2 and moderate for all the properties from streets 4 and 5, occurring moderate SPD peculiarly in street 3, when showed such degree only for the property BD, being strong for the other two properties PR and GM.
In accordance with Cambardella et al. (1994) classification, for all the properties from studied coffee streets, a moderate to strong variance of SPD was observed.
The sill was considerably important in determining thresholds of spatial dependence, in this case from moderate to strong, being that a difficulty was noted during modeling phase to adequate the spherical models (Figures 2 and 3) to the experimental semivariograms due to the high degree of dispersion indicated by CV (Table I) for BD and PR.When the application process of the OLS fitting methodology for streets 4 and 5 was performed, minimum variations of the models caused pure nugget effect, being necessary to repeat the definition process of parameters from the theoretical models of these streets until it was observed visually and aided by cross validation by minimum errors adequacy to the spherical model characteristic of kriging.For the property gravimetric moisture (GM), by observing the high values of the parameters of spherical models concerning to others from Table II, the low coefficient of variation of this property is observed, confirmed in characterization phase (Table I), providing easiness in estimation by kriging.
In cases where a theoretical model is established tending to pure nugget effect (Figures 2 and 3), it may be recommended to apply other model of fitting and estimation methods, such as the restricted maximum likelihood (REML) with inferences related to parameters of random effects of the model, based on likelihood-based statistics (Jomar Filho 2003), model used by Liu (2016) in characterization of spatial variability of geological profiles and also by Li et al. (2015) mapping soil salinity.With the fitted semivariogram models, interpolation by ordinary kriging was performed to elaborate variability profiles of the soil properties, by longitudinal profiles, with the maximum values of properties BD and RP represented by the red color and the maximum values of the property GM in blue color (Figures 4,5,6 and 7).
The variability maps would tend to compact in the direction of soil depth, when deeper layers would be more compact, being the upper layers influenced by agricultural management.In this context, in the maps resulting from the study of properties BD and PR, there are more horizontal variabilities, different from the property (GM), with a strong indication of their physical characteristics by a more homogeneous distribution (medium CV), perhaps influenced by determined water status during rainy season in the data collection, with index of soil water supply for such period of approximately 0.72 of the field capacity from 40 mm, for 1 m of depth, considering the evapotranspiration loss of 3 mm/day and water replenishment in the soil by long-term rainfall time series, with average precipitation for the study region of 1359 mm, according to Guimarães et al. (2010), probably reflecting vertical and horizontal variances of the variable M in the analyzed profiles.From the resulting maps (Figures 4, 5, 6 and 7), the RP for the five streets is observed gradually, representing the continuous nature attributed to determined soil compaction in depth.
Through the software R, package geoR, when generating variability maps by nine classes of predicted variables, it was possible to verify the influence of the dispersion of properties for the cases RP street 4, BD street 5, and PR street 5 (Figures 6 and 7).It was evident the importance of the correct distribution of distances among sample points in order to minimize the coefficient of variation, since greater variability of the properties was observed for smaller streets (streets 1, 2 and 3), different from those observed for larger streets 4 and 5, with greater spacing among samples, showing lower spatial variability, with spherical models tending to the pure nugget as verified by semivariogram models (Figures 2 and 3).Moreover, regarding the definition of interpolation classes, the lower the number of classes, the lower the variability of the properties represented on maps, suggesting standardization in the presentation of such maps by sufficient variability classes (Figures 4,  5, 6 and 7), with views to the total interpolation range indicated by captions, promoting the precision agriculture, with management practices directed to the soil profile in Brazilian crops, using software R and QGIS.
For the analyzed street profiles there was no significant correlation between the properties, furthermore, agricultural systems would exhibit lower spatial variability in soil properties than others, since management activities as tillage and irrigation aim to homogenize properties (Loescher et al. 2014), but with no doubt that geostatistical prediction is a useful technique to quantify uncertainties in the soil properties maps, according to Heuvelink et al. (2016).
In general it is recommended to perform RP measurements in the soil medium moisture condition.However, Vaz et al. (2002) suggest that the ideal would be to perform moisture measurements simultaneously with the probing with penetrometer and then make corrections to a constant moisture value.In order to correctly interpret such data, it would be necessary to  It was possible to verify the influence of great distances among samples of the street 4 (Figure 6) and street 5 (Figure 7) for BD, PR and GM, with empirical models tending to the pure nugget effect (Figures 2 and 3), for these streets of greater length, different from the observed for smaller streets, with greater evidence of spherical model in the modeling phase.
It is essential to know the soil properties in longitudinal profiles because it is associated with the use and proper management of the soil.So, it could help the producer in the irrigation management, drainage, soil preparation and soil and water conservation.It is still important to note that farmers agricultural management and decision-making depend on a combination of knowledge that came from different areas of agricultural interest.Therefore, the results of this study can contribute to more efficient and sustainable coffee growing.

CONCLUSIONS
The physical properties bulk density, penetration resistance and gravimetric moisture showed different spatial dependence structures for each street.
The spherical models of the semivariograms revealed modeling difficulties for streets of greater length (street 4 and 5), as well as, it was possible to verify a clear relationship between the coefficient of variation in the exploratory phase and the minimum variances of parameters from theoretical models, influencing the result, sometimes tending to the pure nugget effect.
By applying the kriging method with OLS fitting and visual identification of the spherical theoretical model, the soil compaction was observed by physical properties BD and PR, with variabilities in the direction of soil depth.For property GM, it was found great influence of high values, consistent with the study period, during the rainy season of the year 2010, with an average annual rainfall of 1359 mm.It was verified vertical and horizontal variabilities of

Figure 1 .
Figure 1.Definition of the study area with a distribution of sampling points by coffee streets, configuring respective longitudinal profiles and sampling points at a maximum depth of 0.60 m.
promoting better development of plants.However, in soils non-revolved annually, such as coffee plantation, resistance values up to 4.0 MPa are tolerated due to the permanence and continuity of pores, more active biological activity and greater stability of aggregates, as stated by Carvalho et al. (2012).Moreover, Moraes et al. (2014) for a Rhodic Eutrudox soil (in accordance with Santos et al. 2013) with very clayey texture under continuous no-tillage, analyzing critital limits of soil penetration resistance, suggests the adoption of a threshold of 3.5 MPa instead of 2 MPa normally used and, for minimum tillage with chiseling, regardless of the cropping systems, the critical RP limit should be raised to 3 MPa.

Figure 2 .
Figure 2. Semivariograms resulting from the OLS fitting aggregated with visual recognition of the spherical model structures in its relationship with the semivariances for BD and PR, per coffee street.

Figure 3 .
Figure 3. Semivariograms resulting from the OLS fitting aggregated with visual recognition of the spherical model structures in its relationship with the semivariances for GM, per coffee street.

Figure 4 .
Figure 4. Spatial variability of soil physical properties BD, PR and GM, in longitudinal profiles, interpolated for streets 1 and 2, using ordinary kriging method whit OLS fitting, aggregated with visual recognition of the spherical models (Figures2 and 3).
firstly understand the dependence of RP for a wide range of moistures, as well as the influence of soil formation factors, being the measure of soil moisture along with penetration resistance important to reduce misinterpretation of results obtained in different field conditions and soil management systems, as mentioned bySilva et  al. (2016).In this context, the authors verified by visual comparison of variability maps (Figures4, 5, 6 and 7), a determined decrease of RP in sites with higher moisture (M) values.

Figure
Figure 5. Spatial variability of soil physical properties BD, PR and GM, in longitudinal profiles, interpolated for street 3, using ordinary kriging method whit OLS fitting, aggregated with visual recognition of the spherical models (Figures 2 and 3).

Figure
Figure 7. Spatial variability of soil physical properties BD, PR and GM, in longitudinal profiles, interpolated for street 5, using ordinary kriging method whit OLS fitting, aggregated with visual recognition of the spherical models (Figures 2 and 3).
Figure 7. Spatial variability of soil physical properties BD, PR and GM, in longitudinal profiles, interpolated for street 5, using ordinary kriging method whit OLS fitting, aggregated with visual recognition of the spherical models (Figures 2 and 3).

Figure
Figure 6.Spatial variability of soil physical property BD, PR and GM, in longitudinal profiles, interpolated for street 4, using ordinary kriging method whit OLS fitting, aggregated with visual recognition of the spherical models (Figures2 and 3).

Table II .
Parameters of semivariograms, spatial dependence degree (SPD) and mean error (ME).