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

Submitted on August 7, 2019 and accepted on September 1 1, 2021. 1 Universidade Estadual Paulista, Campus Ilha Solteira, Faculdade de Engenharia, Ilha Solteira, São Paulo, Brazil. cesar.lima@unesp.br; maa.carvalho@hotmail.com; nelsongiovanini@outlook.com; pauloteodoro@agronomo.eng.br 2 Universidade Estadual Paulista, Campus Jaboticabal, Faculdade de Ciências Agrárias e Veterinárias, Jaboticabal, São Paulo, Brazil. alan.panosso@unesp.br 3 Universidade Estadual Paulista, Campus Botucatu, Faculdade de Ciências Agrárias, Botucatu, São Paulo, Brazil. nidiarcosta@gmail.com *Corresponding author: cesar.lima@unesp.br Sample arrangements and spatial variability characterization of dendometrics parameters of Eucalyptus camaldulensis and physical soil attributes


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 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): (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 122 (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 84 (Figure 1b 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): (2) where: N (h) is the number of observed experimental pairs of Z (xi) and Z (xi + h) separated by a distance h.
being the maximum distance in which the variogram was defined, C 0 corresponds to the nugget effect; C 0 + 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 2 ), and c) the highest spatial dependence (SDE), which was evaluated as described in Zimback (2001): where: SDE is the indication for space dependency; C is the structural variance; C + C 0 is the threshold.
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). Rev. Ceres, Viçosa, v. 68, n.6, p. 597-608, nov/dec, 2021 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): where Z * is the value to be estimated at the non-sampled point x 0 ; 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): (8) where: Z 1 and Z 2 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 1 and Z 2 , 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. 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.

RESULTS AND DISCUSSION
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 84 , Grid 48 , Grid 42 , WC2 Grid 84 , Grid 48 and Grid 42 ) whereas only the data of the larger mesh showed distribution of tending to normal type (WC1 Grid 122 ) and undefined (WC2 Grid 122 ). 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 crossvalidation 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 arrange-ment (Grid 122 ), 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.
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)   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 superfi-cial 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 42 ( Figure  1d). Considering that the PR showed smaller extents when using a more densified grid (Grid 122 ), the reduction of points caused loss in detail, affecting larger extents (Grid 84 and Grid 48 ) which, in a way, made the mapping more 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.

Models (a) C 0 C 0 +C A(m) R 2 SSR (b)
Arrangements generalized. Thus, it is believed that the last arrangement (Grid 42 ) 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 122 which represents the largest arrangement, an effective extent of 57.6 m and when the refinement points were removed (Grid 84 and Grid 48 ) 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 48 .
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 2 , SSR and SDE) presented excellent response. This fact is verified by the crossvalidation 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 122 , Grid 84 , Grid 48 , Grid 42 ). 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. Rev. Ceres, Viçosa, v. 68, n.6, p. 597-608, nov/dec, 2021  Rev. Ceres, Viçosa, v. 68, n.6, p. 597-608, nov/dec, 2021 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 122 ), since spatial detailing would not be impaired. Considering Table 2, due to best parameters adjustment, one could easily choose to use the Grid 84 arrangement without major problems.
For CBH (Figure 2), as discussed previously (Table 2), was observed mapping only for Grid 122 , 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 42 (Figure 3d), the other mappings were quite similar.
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 48 ( 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.
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 . 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 122 ).
In Table 3, was observed that all combinations provided great semivariograph parameters. However, the combination WC2 42 = f.WC1 122 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, Rev. Ceres, Viçosa, v. 68, n.6, p. 597-608, nov/dec, 2021 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 48 = f.WC1 48 vs WC2 122 ) 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 2 (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 48 = f.PR1 84 , 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 122 ). Therefore, Figure 5 shows the semivariogram of adjustment and the mapping of cokrigings among the best interactions for estimation of WC2.

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