Spatial variability of soil penetration resistance influenced by season of sampling

Variabilidade espacial da resistência do solo a penetração influenciada pela época da coleta de amostras

Spatial variability of soil penetration resistance influenced by season of sampling

Variabilidade espacial da resistência do solo a penetração influenciada pela época da coleta de amostras

Juan José BonninI, ^*; José Manuel Mirás-Avalos^II; Kléber Pereira Lanças^I; Antonio Paz González^II; Sidney Rosa Vieira^III

^IFaculdade de Ciências Agronômicas. UNESP. Fazenda Experimental Lageado. 18603-970, Botucatu, (SP), Brasil. E-mail: [e-mail]

^IIFacultad de Ciencias. Universidade da Coruña, Campus A Zapateira, C.P. 15071, A Coruña, (Spain). E-mail: [e-mail]; [e-mail] ^IIIInstituto Agronômico. Avenida Barão de Itapura, 1481, 13020-902 Campinas (SP), Brasil. E-mail: [e-mail]

ABSTRACT

The aim of this work was to analyze the spatial distribution of soil compaction and the influence of soil water content on the resistance to penetration. The latter variable was described by the cone index. The soil at the study site was a Nitisol and the cone index data were obtained using a penetrometer. Soil resistance was assessed at 5 different depths, i.e. 0-10 cm, 10-20 cm, 20-30 cm, 30-40 cm and deeper than 40 cm, whereas soil water content was measured at 0-20 cm and 20-40 cm. Soil water conditions varied during the different samplings. Coefficients of variation for cone index ranged from 16.5% to 45.8% while those for soil water content varied from 8.96% to 21.38%. Results suggested a high correlation between soil resistance, as assessed by the cone index, and soil depth. However, the expected relation with soil water content was not observed. Spatial dependence was observed in 31 out of 35 data series, both cone index and soil water content. This structure was fitted to exponential models with nugget effect varying from 0 to 90% of the sill value. Four of the data series showed a random behaviour. Inverse distance technique was used in order to map the distribution of the variables when no spatial structure was observed. Ordinary kriging showed a smoothing of the maps compared to those from inverse distance weighing. Indicator kriging was used to map the cone index spatial distribution for recommendation of further soil management.

Key words: cone index, indicator kriging, soil compaction, spatial distribution.

RESUMO

O objetivo deste trabalho foi analizar a distribuição espacial da compactação do solo e a influência da umidade do solo na resistência à penetração. Esta última variável foi descrita pelo índice de cone. O solo estudado foi Nitossolo e os dados de índice de cone foram obtidos usando um penetrômetro. A resistência do solo foi avaliada a 5 profundidades diferentes, 0-10 cm, 10-20 cm, 20-30 cm, 30-40 cm e mais de 40 cm, porém o conteúdo de umidade do solo foi medido a 0-20 cm e 20-40 cm. As condições hídricas do solo variaram nas diferentes amostragems. Os coeficientes de variação para o índice de cone foram 16,5% a 45,8% e os do conteúdo de umidade do solo variaram entre 8,96% e 21,38%. Os resultados sugeriram elevada correlação entre a resistência do solo, estimada pelo índice de cone e a profundidade do solo. Sem embargo, a relação esperada com a umidade do solo não foi apreciada. Observou-se dependência espacial em 31 de 35 séries de dados de índice de cone e umidade do solo. Esta dependência foi ajustada por modelos exponenciais com efeito pepita variável de 0 a 90% o valor do patamar. Em séries de dados o comportamento foi aleatório. Portanto, a técnica das distâncias inversas foi utilizada para cartografar a distribuição das variáveis que não tiveram estrutura espacial. Na krigagem constatou-se uma suavização dos mapas comparados com esses das distâncias inversas. A krigagem indicadora foi utilizada para cartografar a variabilidade espacial do índice de cone e recomendar melhor manejo do solo.

Palavras-chave: índice de cone, krigagem indicadora, compactação do solo, distribuição espacial.

INTRODUCTION

Soil variability is a consequence of complex interactions between its forming processes. Apart from these factors, management practices are additional causes of variability. Management can alter chemical, physical, and biological attributes of the soil (BLEVINS et al., 1983; SETA et al., 1993, CAMBARDELLA et al., 1994; PAZ GONZÁLEZ et al., 2000). Knowledge of the spatial and temporal variability of soil properties influencing crop yields is the first step for the successful application of precision-agriculture techniques (CORÁ et al., 2004).

Usually, soil-attributes variability has been described using traditional statistical techniques. However, it has been proved that agricultural surfaces, considered homogeneous from pedological and management viewpoints, show differences in the spatial distribution of their attributes due to the type of tillage and grown crop (PAZ GONZÁLEZ et al., 2000; VIEIRA and PAZ GONZÁLEZ, 2003).

Most practices in precision agriculture, such as the delineation of areas requiring further management techniques, require a prior mapping of the target attribute over the study area (CORÁ et al., 2004; SIQUEIRA et al., 2008). Until recently, soil properties at unsampled grid nodes have been mainly determined using minimum error variance (kriging) interpolation algorithms (GOOVAERTS, 1999), which have been proved to be reasonably accurate for certain purposes.

Nevertheless, in many agricultural problems, it is of interest to map the zones where the variable under study shows the probability of being higher than a threshold value instead of obtaining maps with specific values. For instance, soil compaction caused by tractors increased soil resistance to penetration (RAPER, 2005; GABRIEL FILHO et al., 2007), which spatial variability has been previously observed (USOWICZ and LIPIEC, 2009). The level of soil compaction caused by traffic increases with soil water content (RICHARD et al., 1999), thus it is important to consider this relation in order to diminish this problem. Soil penetration resistances higher than 2 MPa might cause difficulty for the establishment of crops (SILVA et al., 2000). Therefore, further management or tillage techniques should be performed in order to solve this problem. In such cases, indicator kriging is considered an appropriate tool (MOTOYIMA et al., 2006) since it provides a binary transformation of the results, namely between 0 and 1, which is the estimated-values probability of being greater or lower than the threshold. Thus, the spatial of soil penetration resistance variability can be used to determine minimum number of penetrations and sampling positions for accurate evaluation of the effectiveness of agricultural practices and mechanical impedance for root growth (USOWICZ and LIPIEC, 2009). The maps of spatial variability of soil penetration resistance obtained by kriging are valuable means in precision agriculture to establish adequate management measures (SIRJACOBS et al., 2002).

Moreover, the measurement of soil resistance is expensive due to the specific equipment required for its assessment, thus, relating this variable to soil water content might provide an inexpensive approach for the farmers. The use of secondary information for enhancing the estimation of an attribute is well described in the literature (GOOVAERTS, 1997; BURROUGH and MCDONNELL, 1998). Briefly, when two variables are correlated and one of them is easier to measure or exhaustively known in the whole study area; cokriging techniques can be applied in order to obtain more precise results for the variable of interest. Cokriging is a multivariate extension of kriging (GOOVAERTS, 1997).

By using geostatistical techniques, spatial variability of soil attributes can be assessed, facilitating a regionalised management, following the premises of precision agriculture, namely, performing tillage operations according to this spatial variability (PAZ GONZÁLEZ et al., 2000). Thus, it is possible to improve crop yields, optimize the use of tractors, and reduce the cost of these operations whereas providing a better control of agriculture impact on the environment (CORÁ et al., 2004; MOTOMYIA et al., 2006).

Therefore, the aim of this study was to identify zones with different management requirements using the cone index as an indicator of soil resistance. The use of indicator kriging may allow to describe zones where the use of further management practices should be recommended. In addition, the relation between cone index and soil water content was evaluated in order to enhance the estimations by using the latter as a secondary information in the interpolations.

2. MATERIAL AND METHODS

Experimental site, sampling, and statistical characterization of data sets

The experiment was conducted at Lageado Research Station, which belongs to the Faculty of Agronomical Sciences of the São Paulo State University (Brazil), located in Botucatu at 22º 50' 24'' south and 48º 25' 23''west, 791 m above sea level. According to the classification of Köppen, the climate of the region is mesotermic (Cwa), with a dry season from April to August and a rainy season from September to March, having January as the wettest month.

Soil in this area was classified as a Nitossolo vermelho distroférrico according to EMBRAPA (1999) and Nitisol according to FAO-ISRIC-ISSS (1998), with clay content of approximately 43%. The physical characterization of this soil is presented in Table 1. The study field was located on a rolling relief, covered with chives and has been in a fallow stage during four years. Five samplings were conducted in order to assess the variability of cone index and soil water content during the dry and rainy seasons.

As described in GABRIEL FILHO et al. (2007), soil resistance was measured using the SSMU (Soil Sampling Mobile Unit) equipped with a hydraulic-electronic penetrometer, actioned and mobilized by a piston and hydraulic valves. It has an electronic data-acquisition system which registers the strength values obtained by a load cell, according to ASAE S313.2 (1995) and their corresponding depth data are generated by a powermeter.

Cone index was evaluated at five different depths, i.e. 0-10 cm, 10-20 cm, 20-30 cm, 30-40 cm, and deeper than 40 cm, whereas soil water content was described at 0-20 cm and 20-40 cm. In the end, 102 points were surveyed. Soil water conditions varied during the five different samplings as it was assessed by a gravimetric method.

The main statistical moments, which are generally accepted as indicators of the central tendency and of the data spread, were analyzed. This includes examination of the mean, variance, coefficient of variation and extreme minimum and maximum values. The examination of mean and extreme values may be useful in determining possible outliers. For deciding whether or not data follow the normal frequency distribution, it may be enough to examine the coefficients of skewness and kurtosis. For a population that follows the normal frequency distribution, these coefficients should have values of 0 and three, respectively. The linear correlation between cone index and soil water content was also examined using Pearson's r coefficient, to guide the use of cokriging techniques (WEBSTER and OLIVER, 2001).

Cone index spatial characterization

The characterization of cone index spatial variability is of great interest for planning agricultural practices. The goal in this section is to use the data set to spatially characterize cone index on a grid using a traditional geostatistical tool (kriging). In order to achieve this aim, the analysis is divided into two stages: (1) semivariogram analysis and (2) mapping.

Semivariogram analysis

Spatial variability was assessed through the analysis of scaled semivariograms (VIEIRA et al., 1997) of cone index and soil water content. Semivariogram analysis is an essential tool in the classic geostatistical analysis. The purpose of this analysis is to (1) identify the spatial structure of a variable by computing an "experimental" semivariogram, and (2) fit the "best" semivariogram model for subsequent kriging mapping. Since the aim of this paper is not to describe thoroughly this methodology, readers are referred to ISAAKS and SRIVASTAVA (1989) and GOOVAERTS (1997) for more details on semivariogram analysis.

Briefly, the experimental semivariogram γ(h) is defined as the expected squared difference in value between pairs of samples with a given relative orientation:

where n is the number of observations in the sample of the attribute z separated by the distance h, known as the lag distance. In general, the range in lag distances should be no greater than one-half the width of the area being sampled and no less than the smallest sampling distance (BURROUGH and MCDONNELL, 1998). A lag distance approximately equal to the average sample spacing is often acceptable. The semivariogram is a function of both distance and direction, and so it can account for direction-dependent variability (anisotropic pattern) (GOOVAERTS, 2000).

Based on the mathematical model used to fit the data, the following semivariogram parameters were defined: a) nugget effect (C₀), which is the value of the semivariogram when distance is 0; b) range of the spatial dependence (d), which is the range or distance at which semivariogram remains approximately constant, after increasing as distance increases; c) treshold (C₀+C₁) which is the sill value approaching the data variance.

Initially, the method used to select the parameters (nugget effect, range, sill) of the different semivariogram models was the weighted least-squares method, where the weights were the number of data pairs that provided the information needed to calculate the experimental semivariograms. Then, the cross-validation technique (CHILÉS and DELFINER, 1999) was used to check the model performance. According to KARNIELI (1990), two criteria were considered for the goodness-of-fit: (i) sample correlation coefficient between the estimation values and the known values (r², this term should approach 1), and (ii) sample correlation coefficient between the estimation values and the standardized estimation values (MSE, this term should approach 0).

Because of the limited number of measured data, only the omnidirectional semivariogram was computed, and hence the spatial variability is assumed to be identical in all directions. Classic criteria for calculating semivariograms were considered (GOOVAERTS, 1997).

Both soil physical attributes involved in this study were fitted to the exponential model (equation 2). This model describes processes with lower oscillations than those described by the spherical one, which is widely used for describing natural phenomena.

where a is the maximum distance till which the semivariogram is defined. As described in equation 2, the exponential model reaches the sill quickly and asymptotically.

Both experimental semivariograms and their corresponding fitted models were calculated as described in VIEIRA AND PAZ GONZÁLEZ (2003).

Moreover, after fitting the mathematical structure to the experimental semivariogram, the dependence ratio (DR) was determined according to CAMBARDELLA et al. (1994). This ratio represents the percentage of the nugget effect (C₀) in relation to the sill (C₀+C₁). The values of DR can be grouped as follows: strong dependence (< 25%), moderate dependence (25% - 75%), and weak dependence (> 75%).

Soil resistance mapping

Inverse distances method was used to interpolate values for mapping soil resistance data when they did not have spatial dependence. In this method the variable is estimated as a linear combination of several neighbouring observations, with the weights being inversely proportional to the square distance between observations and the point to be estimated (BURROUGH and MCDONNELL, 1998).

Geostatistical interpolation was performed employing two methods: ordinary kriging (OK), and i ndicator kriging (IK) for characterizing the probability of measuring soil resistance values greater than 2 MPa.

The complete mathematical description of these methodologies is beyond the scope of this paper and can be found elsewhere (ISAAKS and SRIVASTAVA, 1989; BURROUGH and MCDONNELL, 1998; GOOVAERTS, 1999), thus, they are described only briefly.

OK is by far the most common type of kriging in practice (WEBSTER and OLIVER, 2001). Kriging interpolations provide each cell with a local, optimal prediction and an estimation of the error that depends on the semivariogram and the spatial configuration of the data (BURROUGH and MCDONNELL, 1998). Kriging is a generalized technique that allows one to account for the spatial dependence between observations, as revealed by the semivariogram, into spatial predictions. The OK weights are determined such as to minimize the estimation variance (GOOVAERTS, 2000).

Instead of assuming a normal distribution, when using Indicator Kriging (IK), at each estimate location, the cumulative distribution function will be built up at each point based on the behaviour and correlation structure of indicator transformed data points in the neighbourhood (WEBSTER and OLIVER, 2001). To achieve this, IK needs a series of threshold values (IK cutoffs) between the smallest and largest data values in the set. For each IK cutoff, data in the neighbourhood are transformed into 0 (zero) and 1 (one): Zeros if the data are greater than the threshold, and ones if they are less. IK then estimates the probability that the estimation point is less than the threshold value, given this neighbourhood of transformed data and a model of the IK cutoff correlation structure (BURROUGH and MCDONNELL, 1998). From this, probability maps are obtained. For interpreting these maps, the zones with greater probability values are those which are more likely to not surpass the limit value of the threshold.

Mapping, both by inverse distances and kriging (both ordinary and indicator) techniques, was performed using GSTAT (PEBESMA, 2004).

3. RESULTS AND DISCUSSION

A descriptive statistical summary for all the studied variables is presented in table 2. A high variability on the mean values of cone index can be observed among the different samplings and depths considered. These values were usually lower in the surface soil layers and increased with depth due to a more compacted soil in the deeper layers. The effect of sampling on the increasing cone index values can be explained by the traffic of tractors on the area (GABRIEL FILHO et al., 2007). Coefficients of variation (defined as the ratio of the standard deviation to the mean) for these data were highly variable, ranging from 16.5 to 45.8 %, proving that cone index and soil water content are heterogeneous variables. Regarding soil water content, these data were least dependent on depth and their coefficients of variation were smaller than those for the cone index.

Cone index data in the two top layers, i.e. 0-10 and 10-20 cm, the most interesting from the agricultural viewpoint, since these are the layers affected by tillage, differed widely among samplings, varying from 1.58 MPa, in the case of the 0-10 cm layer of the first sampling, to 6.91 MPa, for the 10-20 cm layer during the fifth sampling (Table 2). This latter value is very far from the threshold of 2 MPa which prevents the establishment of the crop. This increase in the cone index is caused by the traffic of tractors in the area, as previously observed by GABRIEL FILHO et al. (2007). Conversely, the soil water content was lower in the case of the fifth sampling than in the case of the first one.

Skewness and kurtosis coefficients were approximately those considered for a standard Gaussian distribution in most cases (Table 2). However, no data transformation was performed on those variables which presented a non-Gaussian distribution.

By linear regression analysis it was found that correlation between cone index and soil water content ranged from 0.02 to 0.39 (data not shown). Since these values proved to be non significant, the use of soil water content as secondary information into the mapping of soil resistance was not worth accounting. Therefore, crokiging techniques were not used for any further data analysis.

In table 3 the parameters of the models which have been fitted to the experimental semivariogram data are shown. A model was fitted to 31 out of 35 data sets, the rest of the data showed a random behaviour expressed as a pure nugget effect. The spatial structure fitted was always an exponential one (Table 3).

Mean squared error (MSE) was close to the value considered a good fit in all the cases (Table 3). In addition, regression coefficients (r²) oscillated between 3% and 91%, approximately, depending on which variable, depth and sampling were considered (Table 3). Higher values were observed in deeper layers suggesting a greater homogeneity in those layers.

Dependence ratio varied from 0% to 89.69%. According to the criteria established by CAMBARDELLA et al. (1994), a weak dependence relation was observed in 8 out of the 35 data saries, the rest of the data sets varied from moderate (17) to strong (10) spatial dependences (Table 3). This indicates that soil resistance to penetration possesses a high spatial variability at short distances compared to other soil physical properties such as clay content or pH (VIEIRA and PAZ GONZÁLEZ, 2003). However, the magnitude of the nugget effect might also be caused by measurement errors or by the size of the penetrometer compared to the grid size. However, it also may be a consequence of variability at a smaller scale than the one considered in this study suggesting the need of a detailed study of the soil resistance to penetration at lower scale. It must also be noted that semivariograms represent but do not explain spatial variation (WEBSTER and OLIVER, 2001).

Range values fluctuated between 16 and 120 m, approximately. This indicates a highly variable spatial dependence. Apart from the case of cone index at more than 40 cm depth in the fourth sampling, the highest range values for this variable were observed in the first sampling (Table 3) when the soil water content was higher (Table 2) suggesting a correlation between both variables. However, as explained before, no significant relation between these variables was found; which may suggest that the two analysed variables are correlated to a certain extent within a specific range of values.

An example of experimental semivariogram with its best fitted model; is shown in figure 1. In this case, the represented variable is cone index at 30-40 cm depth for the third sampling.

The depicted model showed a reasonable value of range (80 m), thus it is considered to be appropriate for describing the spatial variability of soil resistance at this sampling and depth, as shown by the parameters in table 3. Figure 1 shows clearly the moderate magnitude of the nugget effect (37.45%), according to CAMBARDELLA et al. (1994) criteria.

Inverse distances method was used for mapping the variables when they did not show any defined spatial structure. Figure 2 shows as an example the isoline maps obtained using this technique for some of the studied data sets.

Maps obtained by this technique showed a high dependence on the sample grid, as shown in the example of Figure 2, which resulted in a discontinuous appereance of the maps. In addition, they tended to not show a general pattern of behaviour of the variable within the whole map as those obtained by OK. However, they reproduced maximum and minimum values for cone index and soil water content, in contrast to isoline maps obtained by OK.

Two examples of isoline map for cone index corresponding to the third sampling and different depths are shown in figure 3, one of them corresponding to the semivariogram depicted on figure 1.

OK maps seemed to be too smooth for reproducing measured maximum and minimum values. Moreover, associated error maps tended to show a rather high and uniform pattern (Figure 3). Overall, this technique tended to underestimate both cone index and soil water content values. High error values can be explained by the low range values which limit the estimations to the areas close to the measured points. Nevertheless, this approach is worth accounting for the continuity of the maps that it provides which is more useful for planification purposes than the inverse distances approach.

A 2 MPa threshold value was considered in the IK approaches to these cone index data since it was regarded as the limit which can impede the establishment and growth of crops. Figure 4 shows the probability maps obtained by this technique for all the cone index data sets corresponding to the third sampling.

The resulting maps showed that high probability of the occurrence of high resistance values is highly variable with depth. In the case presented as example (Figure 4), zones located at the West of the plot are more likely to surpass the considered threshold value of 2 MPa. This suggests that further management or tillage practices should be performed on this zone in order to facilitate the establishment of the crops. However, probability gradient was different depending on the the accounted depth, as figure 4 describes. This gradient also varied with sampling. It should be noted that the threshold values accounted during the IK process suppose a critical aspect of this analysis, as pointed out by MOTOMIYA et al. (2006).

4. CONCLUSION

A high correlation between soil resistance to penetration, as assessed by the cone index, and soil depth was evidenced. However, the expected relation with soil water content was not observed. A single model for explaining the spatial distribution of the soil resistance to penetration was not identified. Higher values of autocorrelation were found in deeper soil layers suggesting that soil resistance to penetration is more homogeneous in these layers than in the surface ones, due mostly to the significant effect of traffic in the surface soil layers. The existence of spatial dependence allowed the construction of maps using kriging interpolation. This way, visual comparison of cone index map for different depths was possible. In general, IK proved to be an efficient tool for characterizing the soil resistance under this soil type conditions and, thus, provide useful information for advicing further management. The other two techniques used in this study for mapping soil resistance offered valuable results but more difficult to implement with further soil management. Further research should investigate whether other environmental descriptors are able or not to improve the mapping of soil resistance on an interdisciplinary approach involving soil scientists, agronomists, farmers and decision makers.

Received for publication in September 9, 2008 and accepted in December 22, 2009.

• ASAE STANDARD. Soil Cone Penetrometer St. Joseph: American Society of Agricultural Engineering, 1988. 683p. (ASAE STANDARD S313.2)
• BLEVINS, R.L.; SMITH, M.S.; FRYE, W.W. Changes in soil properties after 10 years of no-tillage and conservation tilled corn. Soil and Tillage Research, v.3, p.135-146, 1983.
• BURROUGH , P.A.; MCDONNELL, R.A. Principles of Geographical Information Systems: Spatial Information Systems and Geostatistics. 2.ed. Oxford: Oxford University Press, 1988. 333p.
• CAMBARDELLA, C.A.; MOORMAN, T.B.; NOVAK, J.M., PARKIN; T.B.; KARLEN, D.L.; TURCO, R.F.; KONOPKA, A.E. Field-scale variability of soil properties in central Iowa soils. Soil Science Society of America Journal, v.58, p.1501-1511, 1994.
• CHILÉS, J.P.; DELFINER, P. Geostatistics. Modeling Spatial Uncertainty.New York: John Wiley, 1999. 695p.
• CORÁ, J.E.; ARAUJO, A.V.; PEREIRA, G.T.; BERALDO, J.M.G. Variabilidade espacial de atributos do solo para adoção do sistema de agricultura de precisão na cultura de cana-de-açúcar. Revista Brasileira da Ciência do Solo, v.28, p.1013-1021, 2004.
• EMBRAPA. Sistema brasileiro de classificação de solos Brasilia: EMBRAPA Produção de Informação, 1999. 412p.
• FAO-ISRIC-ISSS. World Reference Base for Soil Resources Rome: FAO (World Soil Resources Reports 84)
• GABRIEL FILHO, A.; LANÇAS, K.P.; BONNIN, J.J.; LEITE, F.; PAULA, C.A.; MONTEIRO, L.A. Variaciones del índice de cono en suelo movilizado posterior al tráfico de un tractor equipado con dos tipos de neumáticos y a cuatro velocidades de desplazamiento. Cadernos do Laboratorio Xeolóxico de Laxe, v.32, p.35-46, 2007.
• GOOVAERTS, P. Geostatistics for natural resources evaluation New York: Oxford University, 1997. 483p.
• GOOVAERTS, P. Geostatistics in soil science: state-of-the-art and perspectives. Geoderma, v.89, p.1-45, 1999.
• GOOVAERTS, P. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. Journal of Hydrology, v.228, p.113-129, 2000.
• ISAAKS, E.H.; SRIVASTAVA, R.M. Applied Geostatistics New York: Oxford University, 1989. 561p.
• KARNIELI, A. Application of kriging technique to areal precipitation mapping in Arizona. GeoJournal, v.22, p.391-398, 1990.
• MOTOMIYA, A.V.D.A.; CORÁ, J.E.; PEREIRA, G.T. Uso da krigagem indicatriz na avaliação de indicadores de fertilidade do solo. Revista Brasileira da Ciencia do Solo, v.30, p.485-496, 2006.
• PAZ GONZÁLEZ, A.; VIEIRA, S.R.; TABOADA CASTRO, M.T. The effect of cultivation on the spatial variability of selected properties of an umbric horizon. Geoderma, v.97, p.273-292, 2000.
• PEBESMA, E.J. Multivariable geostatistics in S: the gstat package. Computers & Geosciences, v.30, p.683-691, 2004.
• RAPER, R.L. Agricultural traffic impacts on soil. Journal of Terramechanics, v.42, p.259-280, 2005.
• RICHARD, G.; BOIZARDA, H.; ROGER-ESTRADEB, J.; BOIFFIN, J.; GUERIF, J. Field study of soil compaction due to traffic in northern France: pore space and morphological analysis of the compacted zones. Soil and Tillage Research, v.51, p.151-160, 1999.
• SETA, A.K.; BLEVINS, R.L.; FRYE; W.W.; BARFIELD, B.J. Reducing soil erosion and agricultural chemical losses with conservation tillage. Journal of Environmental Quality, v.22, p.661-665, 1993.
• SILVA, V.R.; REINERT, D.J.; REICHERT, J.M. Soil strength as affected by combine wheel traffic and two soil tillage systems. Ciência Rural, v.30, p.795-801, 2000.
• SIRJACOBS, D.; HANQUET, B.; LEBEAU, F.; DESTAIN, M.F. On-line soil mechanical resistance mapping and correlation with soil physical properties for precision agriculture. Soil and Tillage Research, v.64, p.231-242, 2002.
• SIQUEIRA, G.M.; VIEIRA, S.R.; CEDDIA, M.B. Variabilidade de atributos físicos do solo determinados por métodos diversos. Bragantia, v.67, p.203-211, 2008.
• USOWICZ, B.; LIPIEC, J. Spatial distribution of soil penetration resistance as affected by soil compaction: the fractal approach. Ecological Complexity, v.6, p.263-271, 2009.
• VIEIRA, S.R. Geoestatística em estudos de variabilidade espacial do solo. In: NOVAIS, R.F., ALVAREZ, V.H., SCHAEFER, G.R. (Ed.). Tópicos em Ciência do Solo. Viçosa: Sociedade Brasileira de Ciência do Solo, v.1, p.1-54, 2000.
• VIEIRA, S.R.; NIELSEN, D.R.; BIGGAR, J.W.; TILLOTSON, P.H. The scaling of semivariograms and the kriging estimation. Revista Brasileira de Ciência do Solo, v.21, p.523-533. 1997.
• VIEIRA, S.R.; PAZ GONZÁLEZ, A. Analysis of the spatial variability of crop yield and soil properties in small agricultural plots. Bragantia, v.62, p.127-138, 2003.
• WEBSTER, R.; OLIVER, M.A. Geostatistics for Environmental Scientists London: John Wiley, 2001. 271p. (Statistics in Practice Series)

Publication Dates

• Publication in this collection
  14 Feb 2011
• Date of issue
  2010

History

• Received
  09 Sept 2008
• Accepted
  22 Dec 2009