DIAGNOSIS AND SPATIAL VARIABILITY OF SOIL FERTILITY AND CROP PRODUCTION IN A TEAK AREA IN EASTERN PARÁ STATE

This study aims to generate the diagnosis and spatial variability maps of soil fertility attributes, as well as teak (Tectona grandis L.f.) production through geostatistics. The study was conducted on a commercial farm in the municipality of Capitão Poço, state of Pará, Brazil. Soil and yield data were collected from 155 sampled georeferenced points and 143 were used in the study (after outlier removal). The collected soil samples were submitted to laboratory analysis to obtain the values of the following variables: pH, organic matter (OM), calcium (Ca), magnesium (Mg) and potassium (K), available phosphorus (P), aluminum (Al), Base Sum (BS), cation exchange capacity at pH 7.0 (T), effective cation exchange capacity (t), base saturation (V%) and aluminum saturation (m %). Thereafter a diagnosis, linear correlation univariate and geostatistical analysis were applied in the resulting data. The univariate statistics showed that normal distribution is not required when evaluating spatial variability of chemical and production variables. Soil fertility diagnosis showed K, OM and P as the most limiting parameters in the commercial teak plot and the importance to fertilize forest areas. Positive correlation was found between P, OM, K with volume per tree. All soil fertility and teak volume variables showed spatial dependence, which enabled the production of spatial variability maps. The variability maps showed to be used complementary with univariate statistics to enable more precise interventions in a teak production area, as it showed the shortage of K in the area and the relation of P and volume per tree. v.26 n.1 2020 DIAGNOSIS AND SPATIAL VARIABILITY OF SOIL FERTILITY AND CROP PRODUCTION IN EASTERN PARÁ STATE


INTRODUCTION
Teak (Tectona grandis L.f.) is a forest species in great demand worldwide. Its timber is attractive due to its quality and durability (Keogh, 2013). In Brazil, it is planted mainly in the states of Mato Grosso, Pará, Tocantins, Rondônia and São Paulo. In 2017 the state of Pará, second higher producer state, was responsible for 12% of teak production (IBGE, 2017). In order to increase the representation of the state on teak production, researches on fertility and silvicultural management practices are needed.
Studies that aimed to find better management practices in teak areas have already shown that the use of techniques such as thinning and greater spacing affect positively on its production (Silva et al., 2016). The management of soil fertility, in temporal or spatial character needs special attention in order to ensure a better understanding of the nutritional requirements of the species, such as the sustainability of the use of correctives and fertilizers (Fernández-Moya et al., 2013;Fernández-Moya et al., 2017). With the aim of further understanding the response of teak on fertilization Fernández-Moya et al. (2017) have shown positive results of the tree after N-P-K application, the authors also pointed the need to take in consideration environmental and silvicultural factors in future researches, with the objective to avoid loss of investments and pollution.
All the studies mentioned above show only the average behavior of the variables of interest, which means that the management of these areas are carried out in a uniform way, that is, the same amount of input throughout the area, which can lead to waste of inputs, resulting in financial and environmental damage, since there are variability in the soil attributes and in growth behavior in teak tree production areas (Alvarado et al., 2014).
In order to ensure a better use of fertilizers and a more suitable management of timber products generated from cultivation in planted forests, it is necessary to use modern management techniques, that take into account the variability of the area, such as precision forestry (Santos et al., 2017a;Gil et al., 2018).
Geostatistics is a precision tool, through which is possible to make spatial inferences about unobserved values through observed values of important attributes in forest production and, thus, to produce maps of spatial variability that facilitate the management of planted forests (Yamamoto and Landim, 2013;Santos et al., 2017a). The application of precision forestry in teak areas has proved to be a feasible alternative for management of such areas in a localized manner in order to use inputs more efficiently and rationally (Gil et al., 2018). Besides that, the use of precision forestry enables the determination of a relation between soil and yield attributes in order to achieve better yield rates (Pelissari et al., 2012;Pita, 2012).
The objective of this study was to generate the diagnosis and spatial variability maps of soil fertility attributes, as well as teak production through geostatistics in order to show the applicability of this tool in their management practices.

MATERIAL AND METHODS
The study was carried out on São Luiz Farm, a property owned by Tietê Agrícola Ltda, located in the municipality of Capitão Poço, northeastern Pará State. The climate in the municipality is Ami type, according to the Köppen's classification. The teak farm has a reforested area of 833.28 ha. The samples for the study were collected from June to July 2017, the seventh year of cultivation of the area, in field 07 of the farm, with an area of 50 ha (Figure 1).
Soil in most of the area is an Argisol (Santos et al., 2018). The area also has two spots of plinthosol. Prior to the implantation of teak trees, the area was cultivated with black pepper (12 ha) and mostly pasture (40 ha).
In its implementation year (2010), the teak plantation had a 4 x 4 spacing in single lines which represents a population of 31,250 trees. In its fourth year, 45% the total population was thinned, resulting in a population of 17,187 teak trees. The seedlings were originated by clones A1, A2 and A3 from Proteca.
Regarding the management of the area, glyphosate herbicide was first applied for the desiccation of spontaneous vegetation and pasture, followed by two heavy harrows, subsoiling and light harrowing. The soil was corrected by applying dolomitic lime at a rate of 4,000 kg ha -1 to the total area and subsequent application of 3000 kg ha -1 of Arad rock phosphate (32% total P 2 O 5 ) also to the whole plot. For planting fertilization, 300 kg ha -1 N-P-K of formulation 08-28-16 was applied once and 100 kg ha -1 of KCl in the first two years. The planting depth was on the first 15 cm of the center of subsoiling furrow. After thinning, in the fourth year, 500 g per tree of phosphate fertilizer (monoammonium phosphate) was applied.
For the execution of this study, soil samples and production data were collected by means of a Garmin navigation GPS. All samples were geographically spatialized, and data for a total of 155 points, distributed in the study area, were collected. MACEDO NETO et. al

Soil sampling
At each soil collection point, five single samples were taken, in the surface soil layer (0-20 cm), considering the GPS error radius, to form a composite sample (Molin et al., 2015). After collection, the samples were sent to a soil laboratory for determination of pH, organic matter (OM), calcium (Ca), magnesium (Mg) and potassium (K), available phosphorus (P) and aluminum (Al). Base sum (BS), cation exchange capacity at pH 7.0 (T), effective cation exchange capacity (t), base saturation (V%) and aluminum saturation (m %) were all calculated using the data results from the variables of the laboratory analysis (EMBRAPA, 2017).

Forest inventory
The inventory was compiled as follows: at each of the 155 central soil collection points, a 6-m radius was plotted and all trees within that area were measured for DBH (diameter at breast height) and total height (Pita, 2012). After collecting DBH and total height, it was possible to determine the volume per tree (which was the production variable used in this study) according to the equation described by Figueiredo et al. (2005).

Univariate statistics and geostatistics
Following sample collection, all data were submitted to descriptive statistical analysis, where the Shapiro-Wilk normality test was first applied to the variables in the study. Those variables that did not present normal distribution were submitted to boxcox transformation. Variables that did not show normal distribution even after this transformation but presented reduced asymmetry were used, as reducing asymmetry helps to approximate the distribution of a variable towards a normal one (Hair et al., 2009).
The coefficient of variation (CV) was also one of the statistical variables taken as reference for data distribution, taking into account the criteria of Pimentel- Gomes and Garcia (2002).
Also, data that were very discrepant in the normal distribution analysis were removed. The removal of these points was performed by plotting them in a GIS software and, since it had the spatial coordinates of each sample, it was possible to perform a spatial analysis in the study area. The data considered anomalous, in relation to its neighbors were removed from the database (Molin et al., 2015). As the Pearson correlation requires the same amount of observations for its realization, the entire line of that discrepancy was eliminated from the database of the present study. A total of 12 samples showed discrepant values in comparison to the rest of the data, therefore, they were eliminated. In total, 143 samples were used in this study.
After performing descriptive statistics, the fertility status of each of the 143 soil samples was evaluated, where, according to the Alvarez et al. (1999) and Raij  MACEDO NETO et. al (2017), they were categorized as very low, low, medium, good (or high), and very good (or very high). The results of these diagnostics are presented as a percentage for each of the categories.

CERNE
A correlation matrix between soil chemical attributes and volume was produced. As some of the variables did not show normal distribution, Spearman's correlation, which is for non-parametric data, was used (Pontes et al., 2018).
The development of spatial variability maps started with the calculation of a semivariogram, which is the mathematical model that enables the characterization of the spatial dependence among the variables under study. Semivariogram modeling was performed using the equation described by Vieira et al. (1983). The models tested were linear, spherical, exponential and Gaussian (Vieira et al., 1983;Liebhold et al., 1993).
In order to select the model that best represents the spatial distribution of each evaluated variable, the following parameters were used: coefficient of determination (R²) which, when presenting values close to one (1) exhibit good adjustments to the model according to the criterion of Downing (1986) and (2) cross-validation represented according to the value of the linear coefficient (a) and angular coefficient (b), in which the best adjustments are represented by a values near or equal to zero and b values close to one (Vieira et al., 2010).
By using the appropriate semivariogram model, it was possible to perform kriging, which is the production of the spatial variability map which, through spatial dependence modeling, allows the determination of the value of any variable in the study area (Dionisio et al., 2015).
Software GS + was used for semivariogram modeling and for the production of spatial variability maps (Kriging). The variability maps produced were exported to Surfer 11.0 for editions. ArcGIS 10 was also used for the development of the location map.

Production potential
The volume of teak ranged from 0.21 m³ to 0.59 m³ per tree and 50% of the evaluated sample points had an average value of 0.37 m³ per tree (Figure 2). The difference between the minimum and maximum values of the production potential shows that growth varied greatly in the evaluated plot. Farias et al. (2003) and Oliveira (2007) showed the variability of potential yield in a citrus orchard and, through the use of geostatistics they were able to produce spatially located maps, which can be used as tools in the treatment of inconsistency in a particular production area. This practice avoids the waste of inputs, of either financial or environmental character. Pelissari et al. (2012) by correlating chemical attributes with development variables in a teak area they were able to estimate, through geostatistics, the spatial variability of DBH and total height over a period of 8 years, which allowed the possibility of targeted interventions, aiming at maximizing wood production. Silva et al. (2016) e Rossi et al. (2011) when assessing the growth of teak areas in the municipalities of Cáceres and Monte Dourado, respectively, obtained average volume per teak tree values lower than the present study. Silva et al. (2016) who evaluated areas with 11 and 16 years, in different spacing, obtained a maximum individual volume value of 0.082 m³ at 11 years old and 0.147 m³ at 16 years old, which were located in their areas of greatest spacing. Rossi et al. (2011) who evaluated areas between 1 and 5 and another 26 years old, managed to estimate an average value, through models of hypsometric relations, of individual volume of 0.116 m³ for teak at 7 years of age.
In the area of the municipality of Cáceres there was no thinning, and this may have influenced the low dendrometry values found (Silva et al., 2016). As can be seen in the history of the study area of the são luiz farm, thinning was carried out in its fourth year of planting, which may have contributed to the satisfactory values found in the study. This fact reiterates the importance of applying management practices in forest areas, with the objective of greater gains.

Descriptive statistics
Soil chemical attributes, pH, K, T, t, OM, V% and volume of the trees presented a coefficient of

DIAGNOSIS AND SPATIAL VARIABILITY OF SOIL FERTILITY AND CROP PRODUCTION IN EASTERN PARÁ STATE
variation (CV) of less than 30%. The remaining variables presented CV above this value, which constitutes high dispersion of data (Table 1). The smallest value of CV, for soil attributes, was found for OM, showing that this is the most homogeneous variable in the soil. Al had the biggest CV found and, consequently, that is the most heterogeneous variable in soil. Pelissari et al. (2012) and Gil et al. (2018) also found the Al as the most heterogeneous variable in their studies, with CV values for this variable close or higher than 100%. Studying the variability and spatial correlation of micronutrients and OM with macadamia nut production Paris et al., 2020 also found the OM as the most homogeneous variable in their studies, with a CV value of 16%.
The variables pH, V% and volume presented normal distribution, while Ca, BS, T and t had to be transformed by box-cox transformation to be normalized. K, Al, OM and m% were not normalized even after the box-cox transformation, but were maintained as this showed a reduction in asymmetry. Mg did not present normal distribution, but it did not need to be submitted to box-cox transformation, as its asymmetry values were reduced after the removal of outliers.
According to Paris et al. (2020) data normality is not a requisite to the application of geostatistics, and it is recommended the use of this tool when the soil shows heterogeneity, so it is possible to see how each variable is distributed in the study area.
A series of papers did not find normal distribution on their soil attributes or production variables in teak areas (Pelisari et al.,2012;Gil et al., 2018) or areas with other cultures (Paris et al., 2020;Zonta et al., 2014) applied geostatistics in their studies and had no limitation with the fact that part or all their variables presented a non-normal distribution.

Soil fertility diagnosis
Following descriptive statistics, the 143 soil samples were classified according to the interpretation criteria determined by Alvarez et al. (1999) and Raij (2017). Soil pH ranged from 4.2 to 5.9. Moreover, 13% of the total samples were considered very low, 81% low and 6% good (Figure 3). Gil et al. (2018) found mean values of pH in a teak area around 5.22, which were a slightly acidic value and limiting for the growth of the species.
Calcium ranged from 0.8 to 7.8 cmol c dm -³, with 8% of all samples being low, 45% medium, 38% good and 8% very good (Figure 3) (Alvarez et al., 1999). These values are less than those found by Fernández-Moya et al. (2013), who analyzed the fertility of 23 different teak areas in Central America and obtained an average content of Ca of 13.22 cmol c dm -³. In comparison to the work of Rosa et al. (2015), who found average Ca values of 0.5 cmol c dm -³ in a teak area in Mato Grosso, the exchangeable Ca concentrations in the soil found in this study were higher.
Mg and K content ranged from 0.2 to 1.8 cmol c dm -³ and 8 to 60 mg . dm -³, respectively. For Mg, 12% of samples were in the low class, 45% in the medium,  For those results of better Mg fertility (above 0.9 mg dm -³, which refers to a status categorized as good), the values found in this paper are similar to those found by Behling (2009), who observed Mg values around 1.3 cmol c dm -³ in tree rows in the top layer (first 20 cm) of a teak growing area, located in Mato Grosso. In relation to K, which had values that ranged from 8.00 mg dm -³ to 60.00 mg dm -³, almost 100% of the observations (97%) are of very low or low status (<41.00 mg dm -³) which shows the deficiency of this nutrient in almost the entire production area (Figure 3). These values were below those found by Behling (2009) who observed K values around 100 mg dm -3 in tree rows in the topsoil.
The base sum (BS) ranged from 1.1 to 9.4 cmol c dm -³, with 7% of sample observations being categorized as low, 52% medium, 35% good, and 6% very good ( Figure 3). Regarding V%, the variation ranged from 27% to 81%, with 15% of observations of low status, 52% average, 30% good and 2% very good, which shows that most of the charges in the soil are comprise exchangeable bases (Alvarez et al., 1999). Values of P ranged from 1.0 to 14.4 mg dm -³ where 22% of the observations showed very low status, 54% low, 20% medium, 3% high ( Figure  3) (Raij, 2017). Results that showed better fertility for P (greater than 8.0 mg dm -³, which refers to good status) had greater values than those found in the literature, such as those by Fernández-Moya et al. (2013) and by Rosa et al. (2015), who observed P concentrations in teak areas around 3 and 2.4 mg dm -³, respectively.
OM values ranged from 10 to 31 g kg -1 and 85% of the observations were within the low status and 15% in the medium (Fig. 3). The samples with values greater than 20 g kg -1 had the best classification of OM (Alvarez et al., 1999). Pelissari et al. (2014) found OM values in a nine-year-old teak plantation located in Mato Grosso that ranged from 15 to 35 g kg -1 , which are values similar to those found in this paper.
Regarding Al and m%, they showed satisfactory values (low or very low) in almost 100% of the study area, which are excellent results for teak production, since low values of these variables are related with better development of this species. (Alvarado and Fallas, 2004).
Soils suitable for teak production must be deep, of medium texture, well drained and with good water holding capacity (Alvarado 2012;Alvarado and Mata, 2013;Jerez-Rico and Coutinho, 2017). Regarding soil physical-chemical characteristics pH is one of the most important in the production of teak, as it is closely linked to the base saturation and aluminum saturation of the soil, that is, its acric and aluminum character (Alvarado and Fallas, 2004;Santos et al., 2018). The application of lime can be a good alternative in reducing acidity in areas of teak production. According to Alvarado and Fallas (2004), the reduced growth, in height, of teak trees less than 25 years old is closely linked to aluminum saturation values above 3%.
The soil of the study area is an argisol, in its vast majority, this type of soil has an aluminum character, which are normal soil characteristics in the soils of the tropic humid Santos et al., 2018) which makes it unsuitable for teak production. The application of lime at the time of planting may have been the reason for finding BS results characterized as medium and high in its vast majority and low or very low values of Al in almost 100% of the study area, which shows the importance of soil management in teak areas. According to Alvarado (2012), the requirement for nutrients for teak increases in plantations over the age of 9 years, and this shows the possibility of this crop's response to fertilization in adult ages.

Correlation between soil and yield variables
Spearman's correlation showed that the volume had significant and positive correlation with OM, P and K (Table 2).
Organic matter has a great influence on the physical attributes of soil, such as better particle aggregation, increased microporosity, better water infiltration rate and increased oxygen availability. Better soil physical characteristics of soil promote teak development, as aeration favors root development, therefore tree development (Wehr et al., 2017).
The better soil P availability favors root development and photosynthetic efficiency of trees. Phosphorous makes a major contribution to the absorption of other nutrients and energy production (Hawkesford et al., 2012). When evaluating the efficiency of nutrient utilization in teak root development, Behling et al. (2014) observed that P, together with sulfur, are the most relevant nutrients to the development of fine and medium roots. In an evaluation on the availability of phosphate sources in Eucalyptus plantations, Dias et al. (2017) observed the availability of this nutrient contributed to root development.
In relation to plants, potassium is linked with stomatal opening and closing, osmotic regulation and enzyme activation (Hawkesford et al., 2012). In teak plantations, it is among the most required nutrients. In a 7.5-year-old teak crop, K was the nutrient with the second highest concentration in the plant parts, in leaves it was the third largest concentrated (Behling, 2009).
Fertilization of teak stands results in benefits for its production. Fernández -Moya et al. (2017), when evaluating NPK application in six-and seven-year-old teak stands observed an increase in production variables (DBH, Height and volume) after fertilization, specially of phosphorus, where about 90 g of triple superphosphate per tree was applied. The results obtained were 13 to 15 meter in height, 15 to 17 cm of DBH and 0.11 to 0.15 m³ volume per tree, values lower than those found in this study for volume (Table 2).

Geostatistics
The variables evaluated in this work showed better adjustments to the exponential, spherical and Gaussian models (Table 3). For soil variables, the exponential model was the most frequent one. Similar results were found by Zanão Junior et al. (2010), who evaluated the spatial variability of soil macronutrients in two no-tillage areas in the state of Minas Gerais and obtained a better semivariogram adjustment for the exponential model for most of the evaluated nutrients.
The evaluated parameters of the models adjusted to semivariograms presented R² values ranging from 0.80 to 0.97, thus indicating that they present reliability in the adjustments (Table 3). Regarding cross-validation, most variables had a value of around 0.0 which is their ideal value. The pH and V% presented values far distant from 0.0, 1.46 and 2.34 respectively, while the Mg, SB, Al and T presented negative values (Table 3).
According to Vieira et al. (2010) when values of a that detach from zero in a positive way, it is because of an overestimation of small values and an underestimation of high values, in spatial variability; values of a are negative because the opposite occurs. Regarding b, all variables showed satisfactory values, close to 1.0. Positive a values were also found by Pelissari et al. (2012) and Gil et al. (2018). Freddi et al. (2017) obtained negative values of a when applying geostatistics on main data components in a Latosol in the Cerrado Amazonian.
The spatial dependence indexes were classified as strong and moderate (table 3). Spatial dependence shows how much of spatial variability is related to data sampling and chance, that is, the degree of randomness in the survey (Dionisio et al., 2015;Santos et al., 2017b). The pH, Ca, K, P, OM and m% showed strong SDI with values ranging from 0.00 (no sampling error) to 0.25 (or 25% sampling error). On the other hand, Mg, SB, Al, T, t and V% and volume showed moderate SDI, showing SDI values ranging from 0.30 to 0.49 which shows that they presented up to 49% of sample error. Similar to that observed by Zonta et al. (2014), the spatial dependence of this study showed that the pattern of regionalized variables was not random and that the sample mesh used was sufficient to study the variability in the area.

Variability of soil chemical attributes and volume
After adjusting the semivariograms of the chemical attributes of the soil and the production, it was possible to perform the kriging of the variables in order to observe their spatial variability in the study area (Figure 4).
Besides showing the local pattern of soil chemical attributes in the studied plot, spatial variability also enabled spatial correlation of each of the soil variables with volume ( Figure 4).
It can be observed in the spatial variability maps that there was no tendency for the highest concentrations of most nutrients and other indicators of soil fertility (Ca, Mg, pH, V%, T and t) to be in the places where the highest volumes were found. In relation to P, it can be observed that the places where their highest concentrations were found showed the largest volumes (Figure 4). It was not possible to observe spatial  MACEDO NETO et. al tendency between K and OM values, variables related to production through Spearman's correlation (Figure 4). The fertilization practices that were performed in the plot were carried out homogeneously, where the same amount of fertilizer was applied to the total area. Such management may result in uneven soil fertility, generating areas with better fertility than others since soil heterogeneity is not taken into account (Alvarado et al., 2014). This uneven management may be responsible for the different results of production values that were found in the study area. It is possible to observe that the best volume results were found in the center of the plot, where the best phosphorus values were obtained.
Correlation with phosphorus may be responsible for the largest volume values in the center of the field, since this area was previously cultivated with black pepper (Piper nigrum L.). Phosphate fertilizers for previous pepper cultivation may have had a residual effect on the production of the plot. The residual effect of phosphorus in plantations is related to the slow availability of this nutrient provided by low solubility of phosphate fertilizers (Caione et al., 2013).
It is important to highlight the great variability of teak volume found in the experiment (Figure 4), therefore highlighting geostatistics as a relevant tool to analyze the development behavior in plantations. This information may assist in monitoring forest development, as done by Pelissari et al. (2012) for teak and by Santos et al. (2017a) for eucalyptus, as well as carrying out forest inventories and thinning operations.

CONCLUSION
The soil fertility diagnosis showed K, OM and P as the most deficient parameters in a commercial teak plot.
All soil fertility and teak volume variables showed spatial variability and SDI were considered strong and moderate.
The increase in P, K and OM values is correlated with the increase in teak volume and the variability maps clearly show the distribution of the largest volumes in the areas with the best available P concentrations, highlighting the importance of this nutrient to the productivity of these trees.
The variability indicates the wide variation in the fertility variables and the volume of teak trees, which allows precise interventions in the cultivation, an important factor for the planning of localized fertilization operations and thinning of individuals in this species that has a cultivation cycle of over 20 years. ¹Nugget effect; ²Spatial variance; ³Range (meters); 4Spatial Dependence Index calculated by C 0 /(C 0 +C 1 ); 5Coefficient of determination; 6Linear coefficient; 7Angular coefficient. a¹ = range; SDI= Spatial dependence index; R²= coefficient of determination; a= linear coefficient; b= angular coefficient.