Sample arrangements and spatial variability characterization of dendometrics parameters of Eucalyptus camaldulensis and physical soil attributes

ABSTRACT

To evaluate the performance of different sampling arrangements, this work aimed to study, in a field of Eucalyptus camaldulensis, the variability and spatial dependence of some attributes of plant and soil (LATOSSOLO VERMELHO Distrófico [Typic hapludox]). To collect the soil (water content [WC] and soil mechanical penetration resistance [PR]) and plant data (tree height [TH] and circumference at breast height [CBH]), a sampling arrangement containing 122 points was used in an area of 1.98 ha and three other arrangements, with 84, 48 and 42 points. Except for the CBH, the reduction of points did not expressively affect the parameters of the classical statistics of all other parameters evaluated. Spatially, for TH and WC, the reduction of collection points did not significantly affect the geostatistical parameters, allowing estimating mapping similar to the initial mapping with high correlation coefficients. The WC at 0.10-0.20 m layer could be estimated from 48 samples of the WC at 0-0.10 m layer and 48 samples in the layer of 0.10-0.20 m with a correlation coefficient (r) equal to 0.94 when compared to initial kriging from an arrangement with 122 sample points.

Keywords: geostatistical planning; sampling; spatial dependence; soil management.

INTRODUCTION

It is currently common sense in the agrarian sciences to consider that the spatial variability of soil attributes is virtually unquestionable. Hence, characterizing this spatial variability can improve soil management actions and strategies in order to reduce the operational costs and the application of inputs (Molin et al., 2015). The study of spatial variability, which is given by the existence of spatial dependence in a sample field, is prioritized as object of study in geoestatistics, which represents a group of mathematical procedures applicable when the data are georeferenced (Yamamoto & Landim, 2013).

The use of this methodology of analysis has been consolidated for some time in the agrarian sciences and has been increasingly applied in the study of the spatial variability of soil attributes (Almeida et al., 2017; Lemos Filho et al., 2017), as well as in the evaluation of the spatial interactions between soil attributes and plant productivity (Lima et al., 2016; Silva et al., 2017) allowing the interpretation of results based on the natural structure of the attribute (Grego et al., 2014; Landim, 2015).

On the other hand, according to Souza et al. (2014) one of the factors limiting the application of geoestatistics is the amount of data necessary for the sampling process to be representative in a way to capture the real spatial variability of a given attribute. Because, due to the high costs of sampling and attribute analysis, there is some difficulty in harmonizing geostatistical rigor with both economic and operational feasibility for the spatial characterization of soil variability on a commercial scale.

According to Cherubin et al. (2015), the spatial variability of soil attributes may be presented at different scales due to formation processes and, also, from anthropic interference. Thus, considering the numerous variables that directly or indirectly influence a certain attribute, the difficulty of proposing a sample model based on universal collection is evident. According to Oliveira et al. (2011), densified samples provide a better assessment of the spatial variability of an attribute, however, depending on the size of the sample area, this condition requires more labor for data collection and generating higher costs (Van Groenigen et al., 1999; Montanari et al., 2005).

This way, in order to efficiently characterize geostatistically samples arrangements, the literature has presented very studies such as the work of Cherubin et al. (2015), which showed that the grid reduction of the sample maintained the accuracy in the characterization of the spatial variability of P and K by Kriging. Souza et al. (2014), in another study, observed that samplings in two regular grids (with 208 and 206 points) showed similarity in the kriging error and inputs estimation for application accordingly to samplings with 105 and 102 points. Gelain (2016), on the other hand, evidenced several sampling points necessary for estimation of P, K and saturation base (187 points); pH and organic matter (95 points), and soil clay (23 points). Thus, several other studies have been carried out with this aim (Gimenez & Zancanaro, 2012; Cherubin et al., 2014; Alves & Reis, 2017).

However, the fact is that when the intention is to initially apply the geostatistical analysis in any area, the sample planning and the decision on the number of collection points precede the possible field results, thus there is no guarantee of success in evidencing the spatial dependence of the attributes, even using as basis some similar studies. Still, once detected the space dependence, and considering the possibility of consolidating the use of geostatistics in the management of this area with the idea of improving productive efficiency, reducing at maximum the amount of sample points without significant losses in the result seems an adequate decision. All things considered, the purpose of this study was to evaluate several sample arrangements to estimate the dependence and spatial variability of tree height and circumference measured at breast height of Eucalyptus camaldulensis, as well as water content and mechanical penetration resistance in a Typic hapludox in the Brazilian savannah.

MATERIAL AND METHODS

This study was developed at the Research and Extension Farm of Engineering School (Unesp), located in Selvíria, state of Mato Grosso do Sul, Brazil (22º 23 'south latitude and 51º 27' west longitude), 363 m above sea level. The climate of the region was classified as C1dAa' by the Thornthwaite system (Rolim et al., 2007), indicating a dry sub-humid region without water surplus, megathermal with evapotranspiration in the summer less than 48% of annual mean, average temperature of 23.5 ° C and yearly average rainfall of 1,400 mm, with period of highest precipitation between the months of September to June, dry winters from June to August.

The geostatistical grid used was installed under a field of Eucalyptus camaldulensis (mounted in mid-1986 with spacing between plants of 4 x 4 m) in a LATOSSOLO VERMELHO Distrófico típico (Santos et al., 2018); Typic hapludox (USDA, 2014).

The attributes of the Eucalyptus evaluated as parameters for all different sample arrangements were: tree height (TH, m) and circumference at breast height at 1.3 m (CBH, cm). For each collection point, the values of TH and CBH were based on the mean of the data of the trees surrounding the sampling point.

The soil attributes evaluated were: water content (WC, %) and the soil mechanical penetration resistance (RP, MPa). To determine the water content, samples of soils with deformed structure were collected, using a bucket auger. The methodology used is described as explained in Donagema et al. (2011):

$WC=\left(\frac{a-b}{b}\right)$ (1)

where: a is the wet sample mass and b is the dry sample mass.

The soil mechanical penetration resistance was obtained in the field using an impact penetrometer (Stolf, 1991) and calculated according to a computer spreadsheet (Stolf, 2014). The soil data were analyzed following the layers: 0 - 0.10 m (WC1 and RP1) and 0.10 - 0.20 m (WC2 and RP2).

The spatial distribution of plant and soil data collection followed the arrangements defined by the tested grids (Figure 1). The initial grid (bigger grid) was composed of 7 parallel lines, with 12 sampling points each, with a total number of 84 points (20 × 15 m), covering an area of 1.98 ha (220 × 90 m). For the purpose of detailing the spatial analysis, two smaller grids (refinement grid) were allocated inside the bigger grid, totalizing another 38 sample points. The dot spacing of the refinement grid was 4 x 4 m. Altogether, the base grid consisted of 122 points, being called Grid ₁₂₂ (Figure 1a). From this base grid, a gradual decrease of the sampling density was realized for the elaboration of three new grids with different sample arrangements. Thus, the final result was a grid with 84 points, Grid ₈₄ (Figure 1b), a grid with 48 points, Grid ₄₈ (Figure 1c) and another with lower density containing 42 points, Grid ₄₂ (Figure 1d).

Figure 1:
Arrangements of sample collection (Grids).

Based on the data extracted and organized in each experimental arrangement, for each attribute studied, the descriptive analysis was carried out with the aid of classical statistics, using the SAS Software (Schlotzhaver & Littell, 1997). To test the hypothesis of the normality of the attributes, the Shapiro-Wilk test was used at the 1% level of probability significance. The data of each arrangement sampling were submitted to variance analysis using the F test at the significance level of 5%.

Spatial analysis was performed using the Gamma Design Software 7.0 - GS⁺ (Robertson, 2004). Thus, separately for each attribute, the spatial dependence was analyzed by calculating a simple semivariogram, based on the stationarity intrinsic hypothesis (Yamamoto & Landim, 2013) estimated by the following expression (Almeida et al., 2017):

$\widehat{\mathrm{\gamma }}\left(h\right)=\frac{1}{2N\left(h\right)}\sum _{i=1}^{N\left(h\right)}[Z\left(X_i\right)-Z(X_i+h){]}^2$ (2)

where: N (h) is the number of observed experimental pairs of Z (xi) and Z (xi + h) separated by a distance h.

Thus, the experimental semivariograms were adjusted by the following theoretical models: exponential (2), spherical (3) and Gaussian (4).

$\widehat{\mathrm{\gamma }}\left(\mathrm{h}\right)={\mathrm{C}}_0+\mathrm{C}\left\{1-\mathrm{e}\mathrm{x}\mathrm{p}\left[-3\left(\frac{\mathrm{h}}{a}\right)\right]\right\};\mathrm{}\mathrm{h}>0$ (3)

$\widehat{\mathrm{\gamma }}\left(\mathrm{h}\right)={\mathrm{C}}_0+\mathrm{C}\left[\frac{3}{2}\left(\frac{\mathrm{h}}{a}\right)-\frac{1}{2}{\left(\frac{\mathrm{h}}{a}\right)}^3\right];\mathrm{}0\underset{\_}{<}\mathrm{h}\underset{\_}{<}{\mathrm{A}}_0\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{ã}\mathrm{o}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\widehat{\mathrm{\gamma }}\left(\mathrm{h}\right)={\mathrm{C}}_0+\mathrm{C};\mathrm{}\mathrm{h}\underset{\_}{>}a$ (4)

$\widehat{\mathrm{\gamma }}\left(\mathrm{h}\right)={\mathrm{C}}_0+\mathrm{C}\left\{1-\mathrm{e}\mathrm{x}\mathrm{p}\left[-3\left({\frac{\mathrm{h}}{a}}^2\right)\right]\right\};\mathrm{}0<\mathrm{h}<\mathrm{d}$ (5)

being the maximum distance in which the variogram was defined, C₀ corresponds to the nugget effect; C₀ + C to the level and a range of the variogram (Vieira, 2000).

The analysis for choosing each semivariographic adjustment was done using as criteria: a) smaller sum of the residuals squares (SSR); b) the highest coefficient of determination (R²), and c) the highest spatial dependence (SDE), which was evaluated as described in Zimback (2001):

$SDE=[C/(C+C_0].100$ (6)

where: SDE is the indication for space dependency; C is the structural variance; C + C₀ is the threshold.

The interpretation for the SDE was, as follows: a) SDE ≤ 25% indicating weak spatial dependence; b) 25% <SDE ≤ 75% moderately dependence, and c) SDE> 75% strongly dependence.

The highest correlation coefficient (r) between observed and estimated values from the cross-validation process was adopted as the final criterion for accepting or not the semivariographic adjustment. Thus, it was possible to compare all the information among all different sample arrangements.

The cross-validation process, which consists of the removal of each observation belonging to the dataset with subsequent value estimation by the interpolation method (ordinary kriging), was used to verify the reliability of the adjusted mathematical model. The model chosen was the one that best estimated the observed values, which is, the one that produced a linear regression equation between the observed values, as a function of the estimated values closer to the bisector (intercept equal to zero and angular coefficient = 1) (Isaaks & Srivastava, 1989).

The parameters models adjusted to the experimental semivariograms were used in the estimation of the attributes from non-sampled locations by means of the ordinary kriging technique. In this process the estimates were made based on the equation (Faria, 2013):

$z^{*}\left(x_0\right)={\sum }_{i=1}^N{}_i{z(x}_i)$ (7)

where Z^* is the value to be estimated at the non-sampled point x₀; N, the number of measured values Z(x_i) presented on the estimation and ( _i the weights associated with each measured value Z(x_i).

As a way of validating the maps results, the correlation coefficient of the numerical data generated in the kriging process of each attribute was analyzed. For this purpose, the Pearson correlation matrix was assembled using the Excel spreadsheet.

Considering that the spatial correlation between two variables can be verified and measured, and considering the difficulties involved in the evaluation of soil attributes in depth, the cross-semivariographic adjustment was used to test the feasibility of applying the co-kriging technique, according to the following expression (Guerra, 1988):

$\widehat{\mathrm{\gamma }}\left(h\right)=\mathrm{}\frac{1}{2N\left(h\right)}\mathrm{}\sum _{i=1}^{N\left(h\right)}\left\{\left[Z_1\left(X_{1i}+h\right)-Z_1\left(X_{1i}\right)\right]\mathrm{*}\left[Z_2\left(X_{2i}+h\right)-Z_2\left(X_{2i}\right)\right]\right\}$ (8)

where: Z₁ and Z₂ correspond to the values of two variables correlated; h is the distance between same variable samples; N (h) is the number of values of Z₁ and Z₂, separated by a vector of distance h.

Finally, to compare and validate the data generated by co-kriging in relation to the original mapping, was used once again the Pearson correlation matrix.

RESULTS AND DISCUSSION

Table 1 shows descriptive analysis of the attributes studied. According to the classification of Pimentel-Gomes & Garcia (2002), tree height (TH) and breast circumference (CBH) presented average data variability, in the same way as water content (WC) for all arrangements and depths evaluated.

Table 1:
Descriptive statistical analysis of tree height, circumference at breast height, water content and soil mechanical penetration resistance for different sample arrangements in a Brazilian Savannah Typic hapludox

Arrangements	Measures of Descriptive Statistics
	Mean*	Median	Value		Standard Deviation	Coefficients			Probability (a
			Min.	Max.		Variation(%)	Kurtosis	Asymmetry	Pr<w	FD
Tree Height (m)
GRID 122	23.2(A)	23.4	14.4	32.2	2.831	12.2	0.545	0.196	0.742	NO
GRID 84	23.4(A)	23.5	15.4	30.1	2.833	12.1	-0.071	-0.017	0.901	NO
GRID 48	22.9(A)	23.1	15.4	30.1	2.675	11.7	0.663	-0.134	0.764	NO
GRID 42	23.6(A)	23.6	18.5	28.6	2.387	10.1	-0.124	0.199	0.708	NO
Circumference at Breast Height (cm)
GRID 122	83.9(A)	83.2	65.4	113.4	9.952	11.9	0.329	0.751	<0.001	IN
GRID 84	84.7(A)	84.1	65.4	113.4	10.526	12.4	0.054	0.626	0.025	TN
GRID 48	85.4(A)	84.4	68.3	113.4	11.163	13.1	-0.043	0.596	0.086	NO
GRID 42	84.5(A)	84.1	65.4	113.4	12.412	14.7	-0.192	0.685	0.039	TN
Water Content (0 - 0.10 m)
GRID 122	11.8(A)	12.1	7.9	15.0	1.575	13.3	-0.275	-0.461	0.024	TN
GRID 84	11.6(A)	11.8	7.9	15.0	1.653	14.2	-0.462	-0.367	0.112	NO
GRID 48	11.6(A)	11.7	8.1	15.0	1.598	13.8	-0.402	-0.218	0.805	NO
GRID 42	11.6(A)	11.8	7.9	14.4	1.714	14.8	-0.208	-0.585	0.050	NO
Water Content (0.10 - 0.20 m)
GRID 122	12.1(A)	12.4	8.6	15.1	1.523	12.6	-0.701	-0.377	0.006	IN
GRID 84	11.9(A)	12.1	8.7	15.1	1.576	13.3	-0.865	-0.188	0.091	NO
GRID 48	11.9(A)	12.2	8.8	15.1	1.530	12.9	-0.768	-0.191	0.318	NO
GRID 42	11.8(A)	12.0	8.7	14.3	1.633	13.8	-0.947	-0.336	0.072	NO
Soil Mechanical Penetration Resistance (0 - 0.10 m)
GRID 122	6.74(A)	6.48	3.45	11.89	1.84	27.4	0.118	0.694	<0.001	IN
GRID 84	6.75(A)	6.37	3.45	11.89	1.98	29.4	0.110	0.786	<0.001	IN
GRID 48	6.89(A)	6.34	3.62	11.89	2.00	29.1	0.281	0.901	<0.001	IN
GRID 42	6.73(A)	6.40	3.62	11.35	1.81	26.9	0.235	0.767	<0.001	IN
Soil Mechanical Penetration Resistance (0.10 - 0.20 m)
GRID 122	9.99(A)	9.41	5.40	17.89	2.26	22.7	1.411	1.124	<0.001	IN
GRID 84	10.1(A)	9.42	5.40	17.89	2.45	24.3	1.236	1.106	<0.001	IN
GRID 48	10.0(A)	9.74	6.97	17.09	1.99	20.0	1.028	1.028	<0.001	IN
GRID 42	10.4(A)	9.86	6.97	17.09	2.36	22.2	1.071	1.053	<0.001	IN

• ^*
  ^a)FD = Frequency Distribution, being NO, TN, and IN respectively of the normal type, tending to normal and indefinite; *Averages followed by the same letter in the column do not differ from each other by F (p < 0.05).

The soil mechanical penetration resistance (PR) data indicated a high variation coefficient (VC) (Table 1). This high variability may be related mainly to external anthropic actions that influence the state of soil compaction. According to Silva et al. (2004), the PR is very influenced by conditions of soil management, whose action is not always homogeneous in the field, therefore, there is a high variability of PR directly in the soil. In general, soil compaction occurs regionally, with severe occurrence in sites with higher frequency of land management (Amado et al., 2007; Girardello et al., 2014), such as at the ends of the crop season (Silva et al., 2004), decreasing its intensity towards the center of the area, or even in areas with higher moisture content in the soil.

In addition, it was observed that the differences between the values of the VC of the studied attributes (Table 1) showed only a slight difference between them, indicating that even with the reduction of the number of observations within each sample Grid, the values were not very dispersed around the average value.

The frequency distribution for TH was normal in all Grids, although for the Circumference Breast Height, it varied between normal, tending to normal and indefinite (Table 1). For water content the normal type was present for most of the studied arrangements (WC1 Grid ₈₄ , Grid ₄₈ , Grid ₄₂ , WC2 Grid ₈₄ , Grid ₄₈ and Grid ₄₂ ) whereas only the data of the larger mesh showed distribution of tending to normal type (WC1 Grid ₁₂₂ ) and undefined (WC2 Grid ₁₂₂ ). On the other hand, the PR presented frequency distribution only of the indefinite class.

Despite these observations, when the other statistical parameters (means, and standard deviation) were evaluated, it was noted, even with different sample sizes and arrangements, the statistical values presented little oscillation, indicating, according to the test-F applied, that the mean data of all smaller grids would be similarly representing the population of the sampled area in relation with the largest grid (p < 0.05).

Thus, in Table 1, the mean values for tree height indicated values between 22.9-23.6 m standing within the average range for this species (Lorenzi et al., 2003). However, the circumference of the breast showed averages between 83.9 - 85.4 cm, slightly above the mean of 79.2 cm observed by Nakayama et al. (2017) for the same species.

The data of water content indicated average values for WC1 ranging between 11.6 - 11.8% and for WC2 between 11.8 - 12.1% while PR1 presented average values of 6.74 - 6.89 MPa and PR2 ratio between 9.99 - 10.38 MPa (Table 1). Values for PR between 2.0 and 4.0 MPa have been considered critical to the development of roots for several crops (Arshad et al., 1996; Suzuki et al., 2007), with their harmful effects amplified when the soil presents low WC (Tavares Filho & Tessier, 2009; Tavares et al., 2014).

The Table 2 presents the adjusted parameters to the experimental semivariograms, as well as the cross-validation parameters that attest the final model used for data interpolation by omnidirectional ordinary kriging. The TH showed special dependence for all sampled Grids, with effective limits (A) ranging from 65.3 - 93.2 m, which reflects on more homogeneous and continuous kriging maps for the area management. On the other hand, CBH showed special dependence only for the largest arrangement (Grid ₁₂₂ ), with a short effective range of 12.2 m, evidencing small spatial continuity, generally with mappings showing isolated "management islands". Therefore, due to the spatial characteristic of the CBH in the study area, the reduction of the number of points caused a loss in capture of the spatial dependence, that is, the other Grids denoted pure nugget effect (pne) and impossibility of mapping.

Table 2:
Models and parameters of experimental semivariograms for tree height, circumference at breast height, water content and soil mechanical penetration resistance for different sample arrangements in a Brazilian Savannah Typic hapludox

Arrangements	Semivariograms Adjustments Parameters
	Models (a)	C 0	C 0 +C	A(m)	R 2	SSR (B)	SDE (C)	Cross Validation
							%	a	B	r
Tree Height
Grid 122	exp	9.100X10-1	8.401	74.1	0.839	9.510	89.2	1.010	0.956	0.622
Grid 84	sph	3.600	8.088	82.3	0.953	0.627	55.5	-1.950	1.079	0.577
Grid 48	sph	2.130	7.078	65.3	0.898	0.735	69.9	-0.980	1.041	0.542
Grid 42	sph	1.710	5.936	93.2	0.700	2.630	71.2	1.660	0.927	0.499
Circumference at Breast Height
Grid 122	gau	1.000X10-1	90.800	12.2	0.895	676.0	99.9	16.060	0.815	0.381
Grid 84	pne	103.866	103.866	-	-	-	-	-	-	-
Grid 48	pne	115.357	115.357	-	-	-	-	-	-	-
Grid 42	pne	152.900	152.900	-	-	-	-	-	-	-
Water Content (0 - 0.10 m)
Grid 122	gau	6.400X10-1	4.517	180.8	0.975	7.170 X10-1	85.8	-0.460	1.041	0.844
Grid 84	gau	7.900 X10-1	4.631	181.7	0.978	6.220 X10-1	82.9	-0.630	1.053	0.848
Grid 48	gau	5.600 X10-1	3.739	148.8	0.985	1.920 X10-1	85.0	0.390	0.965	0.826
Grid 42	sph	1.300 X10-1	4.208	165.9	0.943	7.010 X10-1	96.9	-0.020	1.000	0.834
Water Content (0.10 - 0.20 m)
Grid 122	gau	4.500 X10-1	4.102	164.5	0.970	6.550 X10-1	89.0	-0.350	1.029	0.859
Grid 84	gau	6.400 X10-1	3.969	162.6	0.969	6.000 X10-1	83.9	-0.420	1.034	0.849
Grid 48	gau	4.300 X10-1	3.612	156.6	0.978	2.900 X10-1	88.1	0.530	0.954	0.832
Grid 42	sph	1.000 X10-2	4.107	180.9	0.940	8.010 X10-1	99.8	0.560	0.952	0.833
Soil Mechanical Penetration Resistance (0 - 0.10 m)
Grid 122	exp	1.132	3.483	57.6	0.793	8.132 X10-1	67.5	3.510	0.478	0.200
Grid 84	sph	1.798	3.597	99.3	0.811	5.540 X10-1	50.0	1.580	0.762	0.362
Grid 48	sph	1.205	3.613	108.2	0.706	1.270	66.6	1.100	0.834	0.456
Grid 42	pne	2.875	2.875	-	-	-	-	-	-	-
Soil Mechanical Penetration Resistance (0.10 - 0.20 m)
Grid 122	exp	8.040 X10-1	2.971	35.4	0.820	6.370 X10-1	72.9	0.350	0.966	0.406
Grid 84	sph	2.600 X10-1	6.409	31.2	0.857	6.790 X10-1	95.9	4.370	0.566	0.321
Grid 48	exp	7.560 X10-1	3.062	126.3	0.821	4.217 X10-1	75.3	1.950	0.802	0.397
Grid 42	pne	5.582	5.582	-	-	-	-	-	-	-

• ^*
  ^a)Adjusted models, where: gau = gaussian, exp = exponential, sph = spherical, pne = pure nugget effect; (b) SSR = sum of the squares of the residues; (d) SDE = spatial dependency evaluator.

It was observed that for all sample Grids WC presented spatial dependence with values ranging between 148.8 m and 181.7 m for the layer of 0 - 0.10 m and between 156.6 m and 180.9 m for the layer of 0.10 - 0.20 m (Table 2). These ranges (lower and higher) represented respectively 62.8% and 76.7% of the larger separation distance between points used for the construction of the experimental semivariograms, using TH as example, reflected in kriging maps with management belts more homogeneous and continuous, with little fragmentation or isolated management islands.

The spatial continuity values (A) found for TH (Table 2) corroborated with what was observed by Carvalho et al. (2012) which was working with the Eucalyptus camaldulensis in a sample grid of 120 points in an area of 75 x 285 m, under the same type of soil class, observed an extent of 75.7 m. On the other hand, these authors observed for the Circumference at Breast Height greater extents to what was observed in the present study. This fact may indicate, with respect to height, that this plant attribute has an important genetic influence, being less variable in function of the different physical conditions of the soil. Contrarywise, it is possible to speculate that the circumference at breast height is more influenced by soil conditions, and therefore, presented such difference.

The semivariographic parameters pointed to WC (Table 2) gives the possibility to diagnose that this attribute appears to have a fairly homogeneous spatial continuity for the studied area. Data observed by Gonçalves et al. (1999); Grego & Vieira (2005) also stated that the soil moisture was not distributed randomly in the area, having a well-defined spatial dependence on the superficial layers of the soil. This result may be related to an uniform spatial distribution of clay and organic matter content in the soil, which are the most important factors for water retention in the soil.

In turn, the PR (Table 2) presented pure nugget effect (pne) for the last arrangement with the lowest sample number, at both depths (PR1 and PR2), a result that may be related to a greater soil attribute variability and to its response in relation to the reduction of sampling points, especially in the sample model adopted in Grid ₄₂ (Figure 1d). Considering that the PR showed smaller extents when using a more densified grid (Grid ₁₂₂ ), the reduction of points caused loss in detail, affecting larger extents (Grid ₈₄ and Grid ₄₈ ) which, in a way, made the mapping more generalized. Thus, it is believed that the last arrangement (Grid ₄₂ ) due to drastically reducing the sampling and presenting a spatial model of collection differentiated from the others, resulted in the loss of the spatial dependence detection capacity.

This shows that the extent of spatial dependence depends on the variable studied, and soil water content can be sampled with larger spacing’s (smaller sample size), while point reductions, such as those performed in the sampling of resistance to penetration, may result in generalized spatial dependence, thus this attribute should be sampled at smaller distances (larger sample size). These results confirm what was observed by Molin et al. (2015), which recommends a higher sample density to detect spatial dependence of soil resistance to penetration mainly because such attribute presents high variability in the studied soil.

For the PR attributes (Table 2) that showed spatial dependence, we observed two distinct situations that repeat at both depths. In the first case, for the 0 - 0.10 m depth, in the Grid₁₂₂ which represents the largest arrangement, an effective extent of 57.6 m and when the refinement points were removed (Grid ₈₄ and Grid ₄₈ ) could be noted that the attribute continued to show spatial dependence, but with larger ranges (between 99.3 and 108.2 m). The same trend could be observed in the depth of 0.10 - 0.20 m, but on this layer this occurred only on the Grid₄₈.

The resistance to penetration is an attribute that presented the biggest variability (Table 1), because all localized compacting spots corroborate to its high values based on the fact that it is a well-known attribute widely affected by external interferences (Silva et al., 2004; Girardello et al., 2014). Reducing points in this case will result in a less detailed mapping. Resende & Coelho (2014) reinforces that poorly planned sampling may negatively affect soil management, because distortions in the soil attributes mapping can cause all related practices in the area to be disconnected from the real variability in this environment. Fact which was not observed in the present study for water content and for Tree Height.

For Tree Height and Water Content, independent of the sampling arrangement (Table 2), was observed that the adjustment parameters (R², SSR and SDE) presented excellent response. This fact is verified by the cross-validation method that showed high correlations between observed and estimated, attesting good mapping of the spatial patterns estimated by kriging in any of the arrangements (Grid ₁₂₂ , Grid ₈₄ , Grid ₄₈ , Grid ₄₂ ). Fact which is not observed for the soil mechanical penetration resistance and neither for CBH. Therefore, in order to evaluate the performance of the mappings, the kriging maps are presented in Figures 2, 3 and 4.

In Figure 2, it can be noted that for the Tree Height, both maps presented well-defined spatial patterns, evidencing the tallest trees regionalized in the central area of the plot, a fact that is justified mainly by competition for light, which is commonly observed in the field. Besides, it is reliable to say that TH monitoring would not require the use of 122 sample points (Grid ₁₂₂ ), since spatial detailing would not be impaired. Considering Table 2, due to best parameters adjustment, one could easily choose to use the Grid₈₄ arrangement without major problems.

Figure 2:
Kriging maps of tree heigth (m) and circumference at breast height (cm) for different sample arrangements in a Brazilian Savannah Typic hapludox.

For CBH (Figure 2), as discussed previously (Table 2), was observed mapping only for Grid ₁₂₂ , which presented a small spatial continuity, displaying small isolated areas.

The Figure 3 shows that, except for the mapping of the WC array in Grid ₄₂ (Figure 3d), the other mappings were quite similar.

Figure 3:
Kriging maps of water content (%) for different sample arrangements in a Brazilian Savannah Typic hapludox.

It should be noted that such similarities were not only in the spatial performance of the many specific management extents but were also evident in the minimum and maximum values of each estimated management range, indicating the possibility of working with a minimum number of sampling points without loss of quality (Figure 3). Likewise, Souza et al. (2014) observed that with a considerable reduction in sampling points (approximately 50%) for the evaluation of clay, sand, among other parameters for the chemical management of soil, continued to present similarity using the kriging estimation.

Thus, considering these information (Table 2), it is steady to state that the WC arrangement in Grid ₄₈ (Figure 1c) emerged as the best cost-benefit arrangement for the study area in the assessment of water content.

While in Figure 4, it is observed that PR was an attribute whose spatial variability tended to appear in short distances, when compared to those observed for WC.

Figure 4:
Kriging maps of soil mechanical penetration resistance (MPa) for different sample arrangements in a Brazilian Savannah Typic hapludox.

For the soil compaction, Molin et al. (2015) emphasizes the need for a large number of samples (low spatial dependence) and sub-samples (low sampling area for the cone index), since it is about an anthropogenic variability. Thus, to detect the actual variability in the area and to elaborate a more accurate mapping, a larger sample number is needed (Oliveira et al., 2011). It was verified that the reduction of points, despite detecting spatial dependence causing a reduction in scale detail.

As a way of validating the comparison between the maps, we analyzed the correlations of the data estimated by the kriging from the Pearson correlation matrix. For Tree Height, the values generated by ordinary kriging technique for Grids with 84, 48 and 42 sample points presented high correlation coefficients (r) when paired with the initial kriging generated from Grid ₁₂₂: r = 0.90 (Grid ₁₂₂ vs Grid ₈₄ ); r = 0.75 (Grid ₁₂₂ vs Grid ₄₈ ) and r = 0.78 (Grid ₁₂₂ vs Grid ₄₂ ). For water content, better performances were observed with the following correlation coefficients (r): in the 0-0.10 m layer: r = 1.00 (Grid ₁₂₂ vs Grid ₈₄ ); r = 0.98 (Grid ₁₂₂ vs Grid ₄₈ ); r = 0.94 (Grid ₁₂₂ vs Grid ₄₂ ); and in the 0.10-0.20 m layer: r = 1.00 (Grid ₁₂₂ vs Grid ₈₄ ); r = 0.98 (Grid ₁₂₂ vs Grid ₄₈ ), and r = 0.96 (Grid ₁₂₂ vs Grid ₄₂ ). This fact was not observed with the same magnitude within the PR data, probably due to the higher variation coefficient (VC) observed (Table 1) and the greater range of scopes observed among the different Grids for this attribute (Table 2).

Considering the practical application of the cokriging function, was analyzed the possibilities of determining the WC (0.10-0.20 m) from co-variable data of WC (0-0.10 m) and determining the WC (0.10-0.20 m) from co-variable data of PR (0-0.10m). Thus, in Table 3, are presented the parameters of the cross-linked estimates semivariograms of WC2 and PR2 and their respective correlations (Pearson correlation matrix) with their original kriging mappings from the 122 points (Grid ₁₂₂ ).

Table 3:
Crossed semivariograms parameters adjusted to water content (WC2) and to soil mechanical penetration resistance (PR2) on different sample arrangements in a Brazilian Savannah Typic hapludox

Arrangements (a)		Cross Semivariograms Adjustments Parameters								(e) Correlation Coefficient (n = 7315)
		Models (b)	C 0	C 0 +C	A(m)	R 2	SQR (c)	SDE (d)	Cross Validation
								%	r	r
Water Content	WC2 42= f.WC1 122	sph	1.000 X10-2	3.912	175.1	0.946	6.990 X10-1	99.7	0.676	0.842
	WC2 48= f.WC1 122	gau	2.600 X10-1	3.448	155.7	0.985	2.210 X10-1	92.5	0.794	0.942
	WC2 48= f.WC1 84	gau	2.600 X10-1	3.448	155.7	0.985	2.210 X10-1	92.5	0.825	0.944
	WC2 48= f.WC1 48	gau	2.300 X10-1	3.410	152.4	0.984	2.080 X10-1	93.3	0.846	0.944
	WC2 84= f.WC1 122	gau	4.700 X10-1	4.011	165.0	0.970	7.230 X10-1	88.3	0.846	0.957
	WC2 84= f.WC1 84	gau	4.700 X10-1	4.011	165.0	0.970	7.230 X10-1	88.3	0.811	0.962
Penetration Resistance	PR2 48 =f.PR1 122	sph	1.000 X10-1	1.398	37.3	0.273	0.379	99.9	0.163	0.297
	PR2 48 =f.PR1 84	exp	1.000 X10-1	1.225	75.3	0.376	0.388	99.9	0.467	0.792
	PR2 48 =f.PR1 48	sph	1.500 X10-1	1.302	107.8	0.589	0.491	88.5	0.453	0.588
	PR2 84 =f.PR1 122	exp	1.000 x10-3	2.596	55.5	0.683	0.604	100.0	0.122	0.435
	PR2 84 =f.PR1 84	sph	1.190	2.388	89.8	0.500	0.933	50.2	0.409	0.688

• ^(a)WC = water content, 1 and 2 refer respectively to the profile depth: 0 - 0.10 m and 0.10 - 0.20 m, attributes followed by ₁₂₂, ₈₄, ₄₈ and ₄₂ refer respectively to sample arrangements; ^(b)adjusted models, where: gau = gaussian, exp = exponential, sph = spherical; ^(c)SSR = sum of the squares of the residues; ^(d)SDE = spatial dependency evaluator; ^(e)correlation coefficient between the cokriging maps with the initial kriging maps (at depth 0.10 - 0.20 m in the Grid₁₂₂ arrangement), n = number of data analyzed, ^**significant at 1% probability.

In Table 3, was observed that all combinations provided great semivariograph parameters. However, the combination WC2₄₂ = f.WC1₁₂₂ deserves emphasis since it was the only adjustment where the choice of the number of interpolating neighbors denoted cross validation values (r) abruptly changing (from 0.05 to 0.68), corroborating to the understanding that the form factor, in this case, different from the other grids, presented as a negative factor for the use of the kriging/cokriging technique. The other combinations were satisfactory, based on the practical viewpoint, it was feasible to determine WC for the 0.10-0.20 m layer from only 48 WC2 samples and 48 WC1 samples, obtaining a result whose linear coefficient correlation (r) between the maps (WC2₄₈ = f.WC1₄₈ vs WC2₁₂₂) was 0.94.

Regarding the soil mechanical penetration resistance (Table 3), the performance of the cokrigings was not considerable in relation to the WC, as evidenced by the lower values of the spatial determination coefficients R² (from 0.27 to 0.68) and the cross validation (from 0.16 to 0.46). This fact may be a reflection of the sensitive differences in spatial continuity observed between 0 - 0.10 m (short distances) and 0.10-0.20 m (with greater continuity). Thus, the best cokriging was observed for PR2₄₈ = f.PR1₈₄, a fact that is justified because it is the data set that represented the best cross-validation and the best correlation of the estimated map compared with the original map (Grid ₁₂₂ ).

Therefore, Figure 5 shows the semivariogram of adjustment and the mapping of cokrigings among the best interactions for estimation of WC2.

Figure 5:
Crossed semivariogram and cokriging map to water content in the 0.10 - 0.20 m depth as a function of the WC1 in a Brazilian Savannah Typic hapludox.

CONCLUSIONS

Except for the circumference at breast height, the reduction of collection points did not expressively affect the parameters of the classical statistics of all other parameters evaluated (means, data distribution).

For tree height and water content, the reduction of collection points did not significantly affect the geostatistical parameters, allowing estimating mapping similar to the initial mapping (WC Grid ₁₂₂ ) with high correlation coefficients.

For the soil mechanical penetration resistance mapping, the reduction of points provides a loss of accuracy and detail, mostly because this attribute showed spatial dependence with short distance ranges, mainly in the superficial layer.

For the study area, the WC at the 0.10-0.20 m layer could be estimated from 48 samples of the water content in the first layer (0-0.10 m) and 48 samples in the layer of 0.10-0.20 m with a result whose correlation coefficient (r) was equal to 0.94 when compared to initial kriging from an arrangement with 122 sample points.

The soil mechanical penetration resistance in the layer of 0.10-0.20 m can be estimated from 84 samples in the 0-0.10m layer and only 48 samples in the 0.10-0.20 layer presenting an equal correlation coefficient equal to 0.79 when compared to the initial kriging from the array with 122 sample points.

ACKNOWLEDGEMENTS, FINANCIAL SUPPORT AND FULL DISCLOSURE

The authors declare that there is no conflict of interest in the conduct and publication of this work.

REFERENCES

• Almeida KSSA, Souza LS, Paz VPS, Silva FTS, Santos DN & Pereira JSL (2017) Variabilidade espacial da condutividade hidráulica do solo saturado em Latossolo Amarelo distrocoeso, no município de Cruz das Almas. Irriga, 22:259-274.
• Alves BGX & Reis EF (2017) Micro variabilidade espacial de atributos físicos do solo. Revista Espacios, 38:1-26.
• Amado TJC, Pontelli CB, Santi AL, Viana JHM & Sulzbach L (2007) Variabilidade espacial e temporal da produtividade de culturas sob sistema plantio direto. Pesquisa Agropecuária Brasileira, 42:1101-1110.
• Arshad MA, Lowery B & Grossman B (1996) Physical tests for monitoring soil quality. In: Doran JW & Jones AJ (Ed.) Methods for assessing soil quality. Madison, Soil Science Society of America. p.123-141.
• Carvalho MP, Mendonça VZ, Pereira FCBL, Arf MV, Kappes C & Dalchiavon FC (2012) Produtividade de madeira do eucalipto correlacionada com atributos do solo visando ao mapeamento de zonas específicas de manejo. Ciência Rural, 42:1797-1803.
• Cherubin MR, Santi AL, Eitelwein MT, Amado TJC, Simon DH & Damian JM (2015) Dimensão da malha amostral para caracterização da variabilidade espacial de fósforo e potássio em Latossolo Vermelho. Pesquisa Agropecuária Brasileira , 50:168-177.
• Cherubin MR, Santi AL, Eitelwein MT, Menegol DR, Da Ros CO, Pias OHC & Berghetti J (2014) Eficiência de malhas amostrais utilizadas na caracterização da variabilidade espacial de fósforo e potássio. Ciência Rural, 44:425-432.
• Donagema GK, Campos DVB, Calderano SB, Teixeira WG & Viana JHM (2011) Manual de métodos de análise de solo. 2^nd ed. Rio de Janeiro, Embrapa/CNPS. 212p.
• Faria GA (2013) Uso de um sistema de informação geográfica para análise geoestatística em solos. Master Dissertation. Universidade Federal de Lavras, Lavras, 113p.
• Gelain E (2016) Arranjo amostral para mapeamento de atributos do solo. Doctoral Thesis. Universidade Federal da Grande Dourados, Dourados. 87p.
• Gimenez LM & Zancanaro L (2012) Monitoramento da fertilidade de solo com a técnica de amostragem em grade. Informações Agronômicas, 138:19-25.
• Girardello VC, Amado TJC, Santi AL, Cherubin MR, Kunz J & Teixeira TG (2014) Resistência à penetração, eficiência de escarificadores mecânicos e produtividade da soja em Latossolo argiloso manejado sob plantio direto de longa duração. Revista Brasileira de Ciência do Solo, 38:1234-1244.
• Gonçalves ACA, Folegatti MV & Silva AP (1999) Estabilidade temporal da distribuição espacial da umidade do solo em área irrigada por pivô central. Revista Brasileira de Ciência do Solo , 23:155-164.
• Grego CR & Vieira SR (2005) Variabilidade espacial de propriedades físicas do solo em uma parcela experimental. Revista Brasileira de Ciência do Solo , 29:169-177.
• Grego CR, Oliveira RP & Vieira SR (2014) Geoestatística aplicada a agricultura de precisão. In: Bernardi ACC, Naime JM, Resende AV, Bassoi LH & Inamasu RY (Ed.) Agricultura de precisão: resultados de um novo olhar. Brasília, Embrapa, p.74-83.
• Guerra PAG (1988) Geoestatística operacional. Brasília, Ministério de Minas e Energia. 145p.
• Isaaks EH & Srivastava RM (1989) Applied geostatistics. Nova York, Oxford University Press. 413p.
• Landim PMB (2015) Sobre geoestatística e mapas. Terrae Didatica, 2:19-33.
• Lemos Filho LCA, Ferreira LLN & Lyra DL (2017) Variabilidade espacial de atributos do solo indicadores de degradação ambiental em microbacia hidrográfica. Revista Agroambiente On-line, 11:11-20.
• Lima CGR, Carvalho MP, Souza A, Costa NR & Montanari R (2016) Correlação entre componentes da produtividade da cana-de-açúcar com pH de um Argissolo Vermelho distrófico do noroeste paulista. Engenharia na Agricultura, 24:120-130.
• Lorenzi H, Souza HM, Torres MAV & Bacher LB (2003) Árvores exóticas no Brasil: madeiras, ornamentais e aromáticas. 1^st ed. Nova Odessa, Platarum. 368p.
• Molin JP, Amaral LR & Colaço AF (2015) Agricultura de precisão. São Paulo, Oficina de textos. 238p.
• Montanari R, Marques Júnior J, Pereira GT & Souza ZM (2005) Forma da paisagem como critério para otimização amostral de Latossolos sob cultivo de cana-de-açúcar. Pesquisa Agropecuária Brasileira , 40:69-77.
• Nakayama FT, Silva BE & Curi AS (2017) Avaliação de 12 espécies de eucalipto para reflorestamento para as condições edafoclimaticas da região da Alta Paulista. Revista Científica Anap Brasil, 10:27-35.
• Oliveira RB, Silva AF, Quartezani WZ, Lima JSS & Zimback CR (2011) Levantamento do tipo de malha amostral, tamanho de área e número de pontos utilizados em análise geoestatística. In: 2º Simpósio de Geoestatística Aplicada em Ciências Agrárias, Botucatu. Proceedings, FCA - Universidade Estadual Paulista. p.19-20.
• Pimentel-Gomes FP & Garcia CH (2002) Estatística aplicada a experimentos agronômicos e florestais. Piracicaba, FEALQ. 309p.
• Resende AV & Coelho AM (2014) Muestreo para mapeo y manejo de la fertilidad del suelo. In: Mantovani EC & Magdalena C (Ed.) Manual de agricultura de precisión. Montevideo, PROCISUR. p.49-59.
• Robertson GP (2004) GS⁺: geostatistics for environmental sciences. Michigan, Gamma Desing Software. 152p.
• Rolim GS, Camargo MBP, Lania DG & Moraes JFL (2007) Classificação climática de Köppen e de Thornthwaite e sua aplicabilidade na determinação de zonas agroclimáticas para o Estado de São Paulo. Bragantia, 66:711-720.
• Santos HG, Almeida JA, Oliveira JB, Lumbreras JF, Dos Anjos LHC, Coelho MR, Jacomine PKT, Cunha TJF & Oliveira VA (2018) Sistema brasileiro de classificação de solos. 5^th ed. Brasília, Embrapa Solos . 590p.
• Schlotzhaver SD & Littell RC (1997) SAS system for elementary statical analysis. 2^nd ed. Cary, SAS. 456p.
• Silva MCC, Andreotti M, Costa NR, Lima CGR & Pariz CM (2017) Soil physical attributes and yield of winter common bean crop under a no-till system in the brazilian Cerrado. Revista Caatinga, 30:155-163.
• Silva VR, Reichert JM & Reinert DJ (2004) Variabilidade espacial da resistência do solo à penetração em plantio direto. Ciência Rural , 34:399-406.
• Souza ZM, Souza GS, Marques Júnior JM & Pereira GT (2014) Número de amostras na análise geoestatística e na krigagem de mapas de atributos do solo. Ciência Rural , 44:261-268.
• Stolf R (2014) Penetrômetro de impacto Stolf- programa de manipulação de dados em Excel-VBA. Revista Brasileira de Ciência do Solo , 38:774-782.
• Stolf R (1991) Teoria e teste experimental de fórmulas de transformação dos dados de penetrômetro de impacto em resistência do solo. Revista Brasileira de Ciência do Solo , 15:229-235.
• Suzuki LEAS, Reichert JM, Reinert DJ & Lima CLR (2007) Grau de compactação, propriedades físicas e rendimento de culturas em Latossolo e Argissolo. Pesquisa Agropecuária Brasileira , 42:1159-1167.
• Tavares Filho J & Tessier D (2009) Compressibility of Oxisol aggregates under no-till in response to soil water potential. Revista Brasileira de Ciência do Solo , 33:1525-1533.
• Tavares UE, Montenegro AAA, Rolim MM, Silva JS, Vicente TFS & Andrade CWL (2014) Variabilidade espacial da resistência à penetração e da umidade do solo em Neossolo Flúvico. Water Resources and Irrigation Management, 03:79-89.
• United States Department of Agriculture - USDA. Soil survey staff, Keys to soil taxonomy. 12th ed. Washington, DC, USA, 2014.
• Van Groenigen JW, Siderius W & Stein A (1999) Constrained optimization of soil sampling for minimization of the kriging variance. Geoderma, 87:239-259.
• Vieira SR (2000) Geoestatística em estudo de variabilidade espacial do solo. In: Novais RF, Alvarez V & Schaefer CE (Ed.) Tópicos em ciência do solo. Viçosa, Sociedade Brasileira de Ciência do Solo, p.1-55.
• Yamamoto JK & Landim PMB (2013) Geoestatística: conceitos e aplicações. São Paulo, Editora Oficina de Letras. 216p.
• Zimback CRL (2001) Análise espacial de atributos químicos de solos para fins de mapeamento da fertilidade do solo. Professorship Thesis. Universidade Estadual Paulista, Botucatu. 114p

Publication Dates

• Publication in this collection
  03 Dec 2021
• Date of issue
  Nov-Dec 2021

History

• Received
  07 Aug 2019
• Accepted
  11 Sept 2021