Spatial variability of soil physical properties in longitudinal profiles

Abstract

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 profiles. 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, five longitudinal profiles 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 fitted by ordinary least squares (OLS), being detected in higher superficial 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.

Key words
Precision agriculture; geostatistics; soil physics; agricultural management

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), affirming 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 influence. It can be noted by soil moisture with a strong influence on possible changes in the soil physical structure, influencing 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 influence 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.

MATERIALS AND METHODS

The experiment was performed at Brejão farm during a rainy season in December 2010. It is located in the municipality of Três Pontas, southern Minas Gerais, Brazil, geographical coordinates 21º26’08”S and 45º24’53”W, and 915 m high from sea level. The studied area has 10.5 ha cultivated with coffee crop (Coffea arabica L.) cultivar Mundo Novo 376/4, planted in 2005 with 3.8 spacing between lines and 0.8 m between plants. The region climate is Cwa according to Köppen classification, characterized as mild, high altitude tropical climate, hot and rainy summer (Sá Júnior et al. 2012), with a mean annual precipitation of 1.440 mm and mean annual temperature of 18 °C (Souza et al. 2016). 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).

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.

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.

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 decision-making 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.

Once the theoretical models of 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.

Table I
Results of the descriptive analysis for BD (g cm-3), PR (MPa) and GM (%), for the five coffee streets.

Property	Street	Min.	Max.	Median	Average	Variance	Standard deviation	CV (%)
BD	Street 1	0.930	1.620	1.340	1.307	0.036	0.189	14.461
BD	Street 2	0.930	1.173	1.385	1.336	0.030	0.172	12.874
BD	Street 3	0.920	1.680	1.310	1.288	0.040	0.199	15.450
BD	Street 4	1.240	1.780	1.535	1.513	0.024	0.155	10.245
BD	Street 5	0.990	1.630	1.475	1.409	0.038	0.194	13.769
PR	Street 1	0.175	1.611	0.696	0.702	0.151	0.389	55.413
PR	Street 2	0.183	2.120	0.805	0.952	0.300	0.548	57.563
PR	Street 3	0.225	2.181	0.750	0.896	0.248	0.498	55.580
PR	Street 4	0.228	1.717	0.793	0.790	0.136	0.369	46.709
PR	Street 5	0.179	1.805	0.616	0.724	0.185	0.431	59.530
GM	Street 1	16.340	36.990	23.980	24.583	15.062	3.811	15.503
GM	Street 2	20.170	34.030	23.550	24.209	7.798	2.793	11.537
GM	Street 3	16.800	27.240	22.075	21.915	5.463	2.337	10.664
GM	Street 4	17.500	27.990	20.580	20.570	4.260	2.064	10.034
GM	Street 5	12.630	28.920	21.445	21.569	10.456	3.234	14.994

• BD: bulk density; PR: penetration resistance; GM: gravimetric moisture.

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.

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
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.

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, 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.

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.

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.

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

Property	Streets	C0	C1	Threshold	Range	SPD	ME
BD	Street 1	0.002	0.038	0.04	31.16	Strong	-0.00441
BD	Street 2	0.00366746	0.024	0.02766746	50	Strong	0.00530
BD	Street 3	0.01	0.025	0.035	95	Moderate	-0.00870
BD	Street 4	0.015	0.018	0.033	28	Moderate	-0.00319
BD	Street 5	0.02	0.016	0.036	55.01	Moderate	0.00567
PR	Street 1	0.02207626	0.17661007	0.19868633	20.76	Strong	-0.00527
PR	Street 2	0.01605987	0.137	0.15305987	17.12	Strong	0.00000
PR	Street 3	0.0175515	0.15	0.1675515	31.06	Strong	-0.00054
PR	Street 4	0.03	0.056	0.086	19.65	Moderate	0.00369
PR	Street 5	0.047	0.03	0.077	21.86	Moderate	-0.00067
GM	Street 1	0.72964076	4.2	4.92964076	30.8	Strong	-0.03965
GM	Street 2	1	4.8	5.8	68	Strong	-0.01540
GM	Street 3	0.5	10	10.5	103	Strong	0.05846
GM	Street 4	2	4.9	6.9	49.75	Moderate	-0.00128
GM	Street 5	3.7	4.5	8.2	74.02	Moderate	0.09724

• BD: bulk density; PR: penetration resistance; GM: gravimetric moisture.

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).

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 (Figures 2 and 3).

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 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 (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).

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 moisture (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 effect, 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 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 by Silva et al. (2016). In this context, the authors verified by visual comparison of variability maps (Figures 4, 5, 6 and 7), a determined decrease of RP in sites with higher moisture (M) values.

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 the property GM, perhaps during a certain water status, with a representation of mean values from the classes interpolated by color intelligible to the total amplitude of interpolated values. Furthermore, in cases where the spherical model tends to the pure nugget effect, it may be recommended to use other geostatistical methods for the inference of parameters that define the empirical semivariogram, assuming likely data normality, such as in the REML method.

ACKNOWLEDGMENTS

The authors would like to thank for the support of UFLA in the research of geostatistics applied on precision agriculture.

REFERENCES

• ADDIS HK, KLIK A, & STROHMEIER S. 2015. Spatial variability of selected soil attributes under agricultural land use system in a mountainous watershed, Ethiopia. Int J Geos 6(06): 605-6013.
• ALVES MC, SUZUKI LGAS & SUZUKI LEAS. 2007. Soil bulk density and water infiltration rate as indicators of physical quality recovery of an oxisol. Rev Bras de Cienc Solo 31(4): 617-625.
• AMERICAN SOCIETY OF AGRICULTURAL ENGINEERS. 1999. ASAE S313.1: Soil cone penetrometer. In: ASAE Standards 1999: standards engineering practices data, 46th ed., Saint Joseph, 859-860 p.
• ARAÚJO DCDS, MONTENEGRO SM, MONTENEGRO AADA, SILVA JUNIOR, VDP & SANTOS SMD. 2018. Spatial variability of soil attributes in an experimental basin in the semi-arid region of Pernambuco, Brazil. Rev Bras Eng Agric Ambient 22(1): 38-44.
• ARAÚJO J, DIAS J, GUIMARÃES PTG & ALCÂNTARA EM. 2011. Load bearing capacity and critical water content of a Latossol induced by different managements. Rev Bras de Cienc Solo 35(1): 115-131.
• ARSHAD MA, LOWERY B & GROSSMAN B. 1996. Physical tests for monitoring soil quality. In: J.W. Doran and A.J. Jones. (Eds), Methods for assessing soil quality. SSSA Special Publication, Madison: Soil Science Society of America, 123-141 p.
• BARIK K, AKSAKAL EL, ISLAM KR, SARI S & ANGIN I. 2014. Spatial variability in soil compaction properties associated with field traffic operations. Catena 120: 122-133.
• BIONDI F, MYERS DE & AVERY CC. 1994. Geostatistically modeling stem size and increment in an old-growth forest. Can J Forest Res 24(7): 1354-1368.
• BOCKHEIM JG. 2014. Soil Geography of USA: A Diagnostic-Horizon Approach. University of Wisconsin: Soil Science Departament, 13 p.
• CAMBARDELLA CA, MOORMAN TB, PARKIN TB, KARLEN DL, NOVAK JM, TURCO RF & KONOPKA AE. 1994. Field-scale variability of soil properties in central Iowa soils. Soil Sci Soc Am J 58(5): 1501-1511.
• CAMPOS MCC, AQUINO RE, OLIVEIRA IA & BERGAMIM AC. 2013. Variabilidade espacial da resistência mecânica do solo à penetração e umidade do solo em área cultivada com cana-de-açúcar na região de Humaitá, Amazonas, Brasil. Rev Bras Cienc Agr 8(2): 305-310.
• CARVALHO APV, DIAS HCT, PAIVA HN & TONELLO KC. 2012. Resistance mechanics of soil to penetration in the Riacho Fundo watershade, Felixlândia-MG. Rev Árvore 36(6): 1091-1097.
• CARVALHO LCC, SILVA FM, FERRAZ GAS, SILVA FC & STRACIERI J. 2013. Spatial variability of soil physical attributes and agronomic characteristics of coffee crop. Coffee Sci 8(3): 265-275.
• COLLARES GL, REINERT DJ, REICHERT JM & KAISER DR. 2008. Traffic-induced compaction of an oxisol and its relation with growth and yield of common bean and wheat. Rev Bras de Cienc Solo 32: 933-942.
• CORTÉS-D DL, CAMACHO-TAMAYO JH & GIRALDO R. 2016. Spatial prediction of soil penetration resistance using functional geostatistics. Sci Agr 73(5): 455-461.
• CRESSIE NAC. 1993. Statistics for Spatial Data. John Wiley & Sons. 900 p
• CRUZ JS, ASSIS JÚNIOR RN, MATIAS SSR & CAMACHO-TAMAYO JH. 2011. Spatial variability of an Alfisol cultivated with sugarcane. Cienc Investig Agrar 38(1): 155-164.
• ELLEITHY D, HASSAN A, HAGRAS M & RIAD P. 2015. Best applicable geostatistical model for interpolating groundwater-levels in El-Obour city, Egypt. WIT Transactions on The Built Environment, 168: 435-446.
• FERRAZ GAS, SILVA FM, OLIVEIRA MS, CUSTÓDIO AAP & FERRAZ PFP. 2017. Spatial variability of plant attributes in a coffee plantation. Rev Ciênc Agron 48(1): 91-91.
• FERRAZ GAS, SILVA FC, NUNES RA & PONCIANO PF. 2013. Spatial variability of the noise generated by a portable harvester in a coffee field. Coffee Sci 8(3): 276-283.
• FREDDI OS, CENTURION JF, DUARTE AP & PERES FSC. 2009. I-Compactação do solo e produção de cultivares de milho em Latossolo Vermelho. II-Intervalo hídrico ótimo e sistema radicular. Rev Bras de Cienc Solo 33(4): 793-803.
• GOMES FP. 1990. Curso de estatística experimental, 12 ed., São Paulo: Nobel, 467 p.
• GREGO CR, OLIVEIRA RP & VIEIRA SR. 2014. Geostatistics applied to precision agriculture. In: BERNARDI ACC, NAIME JM, RESENDE AV, BASSOI LH & INAMASU RY. (Eds), Agricultura de Precisão: Resultados de um Novo Olhar, Brasília: Embrapa Instrumentação, 194-208 p.
• GUIMARÃES DP, REIS RJ & LANDAU EC. 2010. Índices Pluviométricos em Minas Gerais, Boletim de Pesquisa e Desenvolvimento., Embrapa Milho e Sorgo, 88 p.
• HAPCA S, BAVEYE PC, WILSON C, LARK RM & OTTEN W. 2015. Three-Dimensional Mapping of Soil Chemical Characteristics at Micrometric Scale by Combining 2D SEM-EDX Data and 3D X-Ray CT Images. Plos One 10(9): e0137205.
• HENGL T, HEUVELINK GBM, KEMPEN B, LEENAARS JGB, WALSH MG & TONDOH JE. 2015. Mapping Soil Properties of Africa at 250 m Resolution: Random Forests Significantly Improve Current Predictions. Plos One 10(6): e0125814.
• HEUVELINK GBM, KROS J, REINDS GJ & VRIES WD. 2016. Geostatistical prediction and simulation of European soil property maps. Geoderma Reg 7(2): 201-215.
• IMHOFF S, SILVA AP, GHIBERTO PJ, TORMENA CA, PILATTI MA & LIBARDI PL. 2016. Physical Quality Indicators and Mechanical Behavior of Agricultural Soils of Argentina. Plos One 11(4): e0153827.
• INMET - INSTITUTO NACIONAL DE METEOROLOGIA. 2010. Normais Climatológicas do Brasil, 1981-2010. http://www.inmet.gov.br/portal/index.php?r=clima/normaisClimatologicas Accessed 18 May 2020.
  » http://www.inmet.gov.br/portal/index.php?r=clima/normaisClimatologicas
• ISAAKS EH, SRIVASTAVA RM. 1989. An introduction to applied geoestatistics. New York: Oxford University Press, 561 p.
• JACINTHO JL, FERRAZ, GAS, SILVA, FM & SANTOS AS. 2017. Management zones in coffee cultivation. Rev Bras Eng Agric Ambient 21: 94-99.
• JOMAR FILHO AC. 2003. Nota Metodológica sobre Modelos Lineares Mistos. Universidade Federal do Paraná. Departamento de Estatística.
• JOURNEL AG & HUIJBREGTS CHJ. 1991. Mining geostatistics. London: Academic Press.
• KALIL EB. 1977. Princípios de técnica experimental com animais. Piracicaba: ESALQ/USP, 210 p.
• KLEIN VA & LIBARDI PL. 2002. Bulk density and pore size distribution of an oxisol under different use and management systems. Rev Bras de Cienc Solo 26(4): 857-867.
• LI HY, WEBSTER R & SHI Z. 2015. Mapping soil salinity in the Yangtze delta: REML and universal kriging (E-BLUP) revisited. Geoderma 237: 71-77.
• LIU WF, Y.F. LEUNG YF & LO MK. 2016. Integrated framework for characterization of spatialvariability of geological profiles. Can Geotech J 54(1): 47-58.
• MALI SS, NAIK SK & BHATT BP. 2016. Spatial Variability in Soil Properties of Mango Orchards in Eastern Plateau and Hill Region of India. Veg An Int J Plant Res 29(3): 1-6.
• MARQUES JÚNIOR J & CORA JE. 1998. Atributos do solo para agricultura de precisão. In: SILVA FM, BORGES PHM (Eds), Mecanização e agricultura de precisão. Lavras: UFLA/SBEA, 31-70 p.
• MELLO JM, BATISTA JLF, RIBEIRO JÚNIOR PJ & OLIVEIRA MS. 2005. Adjustment and selection of spatial models of semivariogram envisaging Eucalyptus grandis volumetric estimates. Sci For 69: 25-37.
• MORAES MT, DEBIASI H, FRANCHINI JC & SILVA VR. 2013. Soil penetration resistance in a rhodic eutrudox affected by machinery traffic and soil water content. Eng Agríc 33(4): 748-757.
• MORAES MT, DEBIASI H, FRANCHINI JC & SILVA VR. 2014. Critical limits of soil penetration resistance in a rhodic Eutrudox. Rev Bras de Cienc Sol. 38(1): 288-298.
• OLIVER MA & WEBSTER R. 2014. A tutorial guide to geostatistics: Computing and modelling variograms and kriging. Catena 113: 56-69.
• OLIVEIRA IA, CAMPOS MCC, SOARES MDR, AQUINO RE, MARQUES JÚNIOR J & NASCIMENTO EP. 2013. Spatial variability of physical properties in a cambisol under different land uses in the southern Amazon region. Rev Bras de Cienc Solo 37(4): 1103-1112.
• PALMA MAZ, VOLPATO CES, SILVA FC, SOUZA P & SILVA JA. 2013. Soil penetration resistance in coffee plantations cultivated with mechanized and manual systems. Coffee Sci 8(3): 364-370.
• REICHARDT K. 1985. Transfer process in the soil-plant-atmosphere system, 4th ed., Campinas: Fundação Cargill, 445 p.
• REICHERT JM, ALBUQUERQUE JA, GUBIANI PI, KAISER DR, MINELLA JPG & REINERT, DJ. 2011. Soil hydrology, water availability to plants and agroclimatic zoning. Tópicos Ci. Solo 7:1-54.
• RESENDE AV, HURTADO SMC, VILELA MF, CORAZZA EJ & SHIRATSUCHI LS. 2014. Precision agriculture applications for grain crops production systems in Brazil. In: BERNARDI ACC, NAIME JM, RESENDE AV, BASSOI LH & INAMASU RY (Ed), Agricultura de Precisão: Resultados de um Novo Olhar, Brasília: Embrapa Instrumentação, 14: 194-208.
• REZA SK, NAYAK DC, CHATTOPADHYAY T, MUKHOPADHYAY S, SINGH SK & SRINIVASAN R. 2016. Spatial distribution of soil physical properties of alluvial soils: a geostatistical approach. Arch Agron Soil Sci 62(7): 972-981.
• RIBEIRO JÚNIOR PJ & DIGGLE PJ. 2001. GeoR a package for geostatistical analysis. R-News 1(2): 15-18.
• RINCON DAP, CAMARGO SA, REY D, VELEZ JJ & ZAMBRANO J. 2017. Change of land use and soil coverage in the area of influence of natural national park Los Nevados. Eur Water 58(2): 293-300.
• RONCHI CP, SOUSA JÚNIOR JMD, AMEIDA WLD, SOUZA DS, SILVA NO, OLIVEIRA LBD & FERREIRA PA. 2015. Root morphology of Arabica coffee cultivars subjected to different spatial arrangements. Pesq Agropec Bras 50(3): 187-195.
• SÁ JÚNIOR A, CARVALHO LG, SILVA FF & ALVES MC. 2012. Application of the Köppen classification for climatic zoning in the state of Minas Gerais, Brazil. Theor Appl Climatol 108(1-2): 1-7.
• SANTOS EOJ, GONTIJO I, SILVA MB & PARTELLI FL. 2017. Sampling Design of Soil Physical Properties in a Conilon Coffee Field. Rev Bras de Cienc Solo 41: 1-13 e0160426.
• SANTOS HG, JACOMINE PKT, ANJOS LHC, OLIVEIRA VA & LUBRERAS JF. 2013. Sistema brasileiro de classificação de solos, 3rd ed., Brasília: Embrapa, 353 p.
• SILVA WM, BIANCHINI A & CUNHA CA. 2016. Modeling and correction of soil penetration resistance for variations in soil moisture and soil bulk density. Eng Agríc 36(3): 449-459.
• SOIL SURVEY STAFF. 1993. Soil Survey Manual: Examination and Description of Soil Profiles. USDA-SCS. U.S. Gov. Print. Office. 437 p.
• SOUZA CG, CARVALHO L, AGUIAR P & ARANTES TB. 2016. Machine learning algorithms and variable of remote sensing for coffee cropping mapping. Bol Ciênc Geod 22(4): 751-733.
• STEFANOSKI DC, SANTOS GG, MARCHÃO RL, PETTER FA & PACHECO LP. 2013. Soil use and management and its impact on physical quality. Rev Bras Eng Agríc Amb 17(12): 1301-1309.
• TAYLOR HM, ROBERTSON GM & PARKER JJ. 1966. Soil strength root penetration relations for medium to coarse textured soil materials. Soil Sci 102(1): 18-22
• TORMENA CA, BARBOSA MC, COSTA ACS & GONÇALVES ACA. 2002. Soil bulk density, porosity and resistance to root penetration in an oxisol managed by different soil tillage systems. Sci Agr 59(4): 795-801.
• VAZ CMP, PRIMAVESI O, PRATIZZI VC & LOSSI MF. 2002. Influência da umidade na resistência do solo medida com penetrômetro de impacto. São Carlos: Embrapa Instrumentação e Agropecuária, 5 p.
• YAO X, FU B, LÜ Y, SUN F, WANG S & LIU M. 2013. Comparison of Four Spatial Interpolation Methods for Estimating Soil Moisture in a Complex Terrain Catchment. Plos One 8(1): 1-13 e54660.
• ZUCCO G, BROCCA L, MORAMARCO T & MORBIDELLI R. 2014. Influence of land use on soil moisture spatial–temporal variability and monitoring. J Hydrol 516: 193-199.
• WHELAN B & TAYLOR J. 2013. Precision Agriculture for Grain Production Systems. Csiro Publishing. 208 p.

Publication Dates

• Publication in this collection
  25 Mar 2022
• Date of issue
  2022

History

• Received
  27 Mar 2020
• Accepted
  07 June 2020