Spatial variability of physical and chemical attributes of some forest soils in southeastern of Brazil

Capão Bonito forest soils, São Paulo state, Brazil, have been used for forestry purposes for almost one century. Detailed knowledge about the distribution of soil attributes over the landscape is of fundamental importance for proper management of natural resources. The purpose of this study was to identify the variability and spatial dependence of chemical and physical attributes of Capão Bonito forest soils. A large soil database of regional land was raised and organized. Most of the selected variables were close to the lognormal frequency range. Soil texture presented a higher range in the A horizon, and the nugget effect and sill were greater in the B horizon. These differences are attributed to the parent material of the region (Itararé Geologic Formation), which presents uneven distribution of sediments. Chemical attributes related to soil fertility presented a higher spatial dependence range in the B horizon, probably as a result of more intensive management and erosion history of the superficial soil layer. Maps for some attributes were interpolated. These had specific areas of occurrence and a wide distribution along the perimeter of the Capão Bonito District Forest, allowing a future site-specific soil management.


Introduction
Forest soil refers to a soil that has been developed by and/ or sustains some forest vegetation (Gonçalves, 2002).Capão Bonito forest soils, São Paulo state, Brazil, have been used for forestry purposes for almost one century.Pedologic survey history from this region presents three phases: (i) soil map from the São Paulo state (SNPA, 1960), (ii) pedologic map from the São Paulo state (Oliveira, et al., 1999), (iii) semi-detailed pedologic surveys (Gonçalves, 1997;Rizzo, 2001).This evolution aimed to achieve a greater level of detail, but actually its interpretation is what allows a precise identification of forest management units (Gonçalves, 1988).The method used nowadays for soil cartographic representation is characterized by abrupt boundaries of mapping units.It is crucial to consider the spatial range of soil properties verifying other ways to represent them and providing a better final quality on generated maps (Burrough, 1991;Webster, 2008).Soil attributes must be considered as regional variables (Matheron, 1963) as numeric spatial functions that vary from one region to another with apparent extension whose variability cannot be represented by a simple mathematic function (Guerra, 1988).
Geostatistics has been widely applied in studies on forest soil variability.Several areas of forestry and ecology used the semivariogram as a tool to detect and understand the physical and chemical soil variability related to forest productivity (Bognola et al., 2008), the soil survey (Novaes Filho et al., 2007), and the forest management effects on soil (Lima et al., 2008).
These studies fulfilled a great demand for mapped information of edapho-pedologic attributes, mainly applying to agricultural productivity modeling, fertilization management and delimitation of management units.Therefore, geostatistic results should be used as a routine method to allow more precise recommendations (Vieira, 2000) for in-site soil management (Molin, 2001).This study was carried out with the aim of characterizing the variability and spatial dependence of chemical and physical attributes of Capão Bonito forest soils.

Materials and Methods
The study was carried out on Eucalyptus stands from Capão Bonito (23°55' S and 48°22' W; 604-732 m altitude range), south of the São Paulo state (southeast Brazil).This region is located predominantly on sedimentary rocks of the Itararé Geologic Formation (IPT, 1981), and a smaller area in the far south is located on crystalline rocks represented by the Canastra Geologic Group and Pilar Geologic Complex (IPT, 1981) originating poorly developed soils.There are even records of fluvial sediments occurring during the Pleistocene period along Paranapanema riverbanks (Rizzo, 2001) and Almas riverbanks (Setzer, 1949), influencing pedogenesis of local soils.Local relief indicates that the area belongs to the Paraná Sedimentary drainage basin and to a small region, at the extreme south, to the Atlantic Orogenic Belt (Ross and Moroz, 1997).The Koeppen climate classification of the region is Cfa (maritime temperate climate), i.e., humid and hot summer, according to Setzer (1946).
One hundred and eighty-three soil points sampled by Gonçalves (1997) and Rizzo (2001) were used, including 139 pedologic profiles and 44 samples collected with a Dutch auger.These data were obtained from printed reports, georeferenced and compiled in a soil database for the studied area.From 2006 to 2008, 244 more georeferenced samples collected with a Dutch auger were performed for all the planted areas aiming to review the pedologic mapping, in addition to increasing soil sampling, leading to a total of 427 soil point samples (Figure 1).The soil database from the Capão Bonito forest district includes localization information of each sample, soil and pedologic horizon type, horizon depth, sampling method and physical and chemical soil attribute values.All points are located in forest stand areas (21,883 ha).Among the sampling years 1997 and 2008 only a small percentage of this area suffered intervention of tillage.The techniques of tillage followed the precepts of the minimum soil preparation (Gonçalves et al., 2004).Soil samples have the same distribution of classes of pedological surveys.The soil types Rhodic Haplodux, Typic Hapodlux, Rhodic Paleudult, Rhodic Kandiudalf and Dystrochrept present, respectively, 82.5 %, 5.3 %, 1.4 %, 0.4 % and 10.4 % occurrence of pedological surveys (Gonçalves, 1997;Rizzo, 2001).
Exploratory data analysis was performed by the statistical program SYSTAT v.11 (Wilkinson, 2004).The exploratory step is crucial because it provides identification of possible diverse values and shapes of variable distribution.Position (mean, median and mode), dispersion (minimum, maximum and standard deviation) and distribution measurements (coefficients of variation, skewness and kurtosis) were performed.Normality hypothesis was verified according to the W test at 5 % (Shapiro and Wilk, 1965) in which normal or logarithmic distribution tendencies were verified.
Experimental omnidirectional semivariograms were adjusted by the geostatistic program GS+ v.9.In this program, theoretical model selection (semivariogram) is performed based on the smallest reduced sum of squares (RSS) and on the greatest de- termination coefficient (R 2 ) (Robertson, 2008).In some cases, semivariograms were manually adjusted or "at feeling", originating more realistic and satisfactory results in relation to automatic adjustment by the GS+ v.9 program.Experimental semivariograms were determined until approximately 50 % of the geometric camp (Figure 1), since after this value the semivariogram did not seem to be correct (Guerra, 1988), i.e., its accuracy was reduced due to the small number of possible pairs to calculate the semivariance at this distance.A geometric camp of 12,000 m with partition groups (lags) of 1,200 m was considered, as these lags are the estimators of the experimental semivariograms (Deutsch and Journel, 1998).Although the study area is divided into two subareas, it was considered as one area, because it presents a small separation distance, and the subareas have the same topography and soil groups.Theoretical models considered, such as spherical, exponential, Gaussian and linear, were described by Guerra (1988), Vieira (2000) and Andriotti (2003).Only this theoretical semivariogram group was considered because it usually covers the general dispersion situation of soil science spatial events (Burrough and McDonnell, 1998;Soares, 2006).Through GS + v.9 cross validation, the correlation coefficient of selected models was obtained.
Scaled experimental semivariograms were calculated according to Vieira et al. (1997) by dividing their values by data variances, obtaining semivariograms with a sill close to 1.0.It was possible to plot many soil attributes on one graph along with their proximity when elected, indicating a similarity of range along landscapes (Vieira et al., 1997).The spatial dependence index (SDI) was used according to Zimback (2001), which measures a sample's structural variance effect on the total variance (sill).SDI comprises the following interpretation break: weak SDI ≤ 25 %, moderate SDI between 25 % and 75 % and strong SDI ≥ 75 %.This index is a complement of the traditional method recommended by Cambardella et al. (1994) in which the nugget weight effect (randomness) on total variance is evaluated.Through structural parameters obtained from experimental semivariograms, maps of some properties were created with the geographic information system ArcMap v.9.3 (ESRI, 2008).A punctual ordinary kriging estimator was used for geostatistic interpolation.To assist map interpretation, subtitles were composed with similar amounts by dividing the distribution into four groups that enclose the same occurrence area.Thus, each subtitle group represents 25 % of the mapped area.

Results and Discussion
Capão Bonito forest soils presented distinct textures (Table 1).Samples of the A horizon had 70 to 870 g kg -1 of total sand.Medium, coarse and very coarse sand presence was small (< 20 g kg -1 ); few areas reached values of 441 g kg -1 (coarse sand in the A horizon).This feature became evident by high coefficients of skewness, low coefficient of kurtosis (platykurtic) and by closer tendency to a lognormal distribution of medium, coarse and very coarse sand.Total, very fine and fine sand presented distributions close to normal, positive skewness and low kurtosis coefficient (leptokurtic).Oxidic soils from Capão Bonito presented an average of 115 g kg -1 of silt in the A horizon with positive skewness and platykurtic distribution tending to be lognormal.The mean clay content in the B horizon (310 g kg -1 ) was slightly more than that of the A horizon (284 g kg -1 ) due to the presence of the Rhodic Kandiudalf in the area.Rhodic Hapludox with clayey texture was predominant, thus demonstrating a slight positive skewness for clay content; leptokurtic distribution was more similar to the lognormal.Similar results were obtained with the silt/ clay ratio with a mean value equal to 0.4.This ratio presented positive skewness and platykurtic distribution and was similar to the lognormal.The mean content of natural clay was 113 and 106 g kg -1 in the A and B horizons, respectively, presenting positive skewness and a tendency to a lognormal distribution.Thus, these data clearly represent areas with a predominance of a Hapludox in which low clay dispersion and mobility are usual (EMBRAPA, 2006).
SOM presented a mean content of 23.5 g kg -1 in the A horizon and 13 g kg -1 in the B horizon (Table 1).SOM tended to the normal distribution and reduced positive skewness.Available phosphorus presented high skewness and kurtosis coefficients with mean contents of 4.4 and 3 mg kg -1 in the A and B horizons, respectively.Potassium, calcium and magnesium presented strong positive skewness, high kurtosis coefficients (platykurtic) and distributions tending to be lognormal.Mean values of potassium (0.7 mmol c kg -1 ), calcium (2.3 mmol c kg -1 ) and magnesium (2 mmol c kg -1 ) in the A horizon were lower than the soil critical point (Gonçalves, 1995).The mean sodium content was 0.2 mmol c kg -1 with strong skewness distribution to the left because most of the soils are poor in this element.Aluminum + hydrogen followed distributions tending to normal.Data of pH H 2 O and pH KCl in the A horizon presented a mean value of 3.8 with distributions tending to normal.In the B horizon, the mean was 3.9 and the distribution tended to be lognormal.The maximum pH value recorded was 5.0 (pH H 2 O) and the minimum value was 2.6 (pH KCl).Delta pH values presented negative skewness and distribution close to normal.The most weathered soils were found at areas with delta pH values equal to 0.3 (positive balance charge).
The sum of bases presented high positive skewness and platykurtic distribution (high C.V.) tending to be lognormal.The mean of this attribute was very low, only 5.1 and 3 mmol c kg -1 in the A and B horizons, respectively.The mean values of CEC were 24.7 mmol c kg -1 (effective) and 85.7 mmol c kg -1 (pH7) for the A horizon.Their distributions presented negative positive skewness behavior only for CEC pH7 at the superficial layer.Most of the attributes, in both soil layers, presented distributions tending to be lognormal.In other studies, this was also the best adjustment for most of the edaphic properties (Cambardella et al., 1994;Zimback, 2001).The variables that required transformation, and with p-value close to zero by the normality test (Shapiro and Wilk, 1965), are presented in Table 1.
Only coarse sand (Figure 2b), natural clay (Figure 2d), phosphorus (Figure 2d), exchangeable aluminum and sodium (Figure 2f), pH KCl (Figure 2g) and effective CEC (Figure 2h), in the B horizon, did not present spatial dependence, with the SDI = 0 (Table 2).These traits were represented by a linear model (or potential, with α = 1), which usually represents non-steady events (Andriotti, 2003).For these variables, spatial correlation did not exist and the event could be understood as random  at this sampling intensity.This indicated that the sampling intensity was insufficient to detect the spatial dependence of these variables, considering the scale of a usual semi-detailed pedologic survey.
Very fine sand, medium sand, coarse sand, total sand, silt, clay and natural clay content in the A horizon were adjusted to exponential functions, whereas very fine sand, total sand, silt and clay content were adjusted to spherical models.Other texture variables, such as fine sand, very coarse sand and silt/ clay ratio in both soil layers were also modeled with exponential functions.Both theoretical semivariogram models presented a rapid rise at the origin in which a linear behavior could be seen   ( Vieira, 2000), although exponential models faster increase at the origin than spherical model (Soares, 2006).This finding can be seen in Figures 2a, 2b and 2c for the traits of very fine sand, total sand and clay in which curves in the A horizon showed slightly higher growth (greater slope) at the origin than that in the B horizon.This means that there was greater spatial continuity as the distance from B to A horizon increased.
In general, semivariogram structural traits of very fine sand, very coarse sand, total sand, silt and clay were superior in the A horizon (Table 2), indicating that in this layer these samples were correlated to greater distances in relation to the sub-superficial layer.At areas with Rhodic Hapludox, Corá et al. (2004) and Souza et al. (2004) also found greater spatial dependence for soil texture variables at the 0-20 cm layer.The nugget effect and sill were superior in the B horizon for fine sand, medium sand, very coarse sand, silt and silt/clay ratio (Table 2).A greater standard deviation and/or coefficient of variation for medium sand, total sand, silt, clay, silt/clay ratio and natural clay occurred in the B horizon (Table 1) These ratios detected by semivariograms demonstrated a greater range in the sub-superficial layer.This may be a result of the raw parent material heterogeneity, because approximately 98.5 % of the area presented sedimentary rocks from the Itararé Geologic Formation (IPT, 1981), which have shown an uneven distribution of grain sizes.
SOM had greater spatial continuity close to the origin in the A horizon, adjusting better by the Gaussian model (Figure 2d).It presented weak spatial dependence.In the B horizon, an exponential model with moderate SDI was adjusted.Spatial dependence range was inferior in the A horizon, demonstrating that the natural dispersion of SOM in the B horizon was more homogeneous.More intense soil management in the past and erosion history might have influenced SOM dispersion in the superficial layer.The spatial distribution of phosphorus in the 0-30 cm layer was similar to SOM (Figure 2d and Table 2).Likewise, Ortiz et al. (2006) did not find spatial dependence for SOM in forest soils under Eucalyptus stands in the Paraibuna region (SP state).Boruvka et al. (2005) and Regalado and Ritter (2006) found low spatial dependence for SOM in the superficial layer.
The exchangeable potassium, calcium and magnesium cations presented moderate spatial dependence with R 2 higher than 0.7 (Table 2).The spatial dependence of these traits indicates that in the sub-superficial layer the sample dispersion was superior in relation to the A horizon.This fact can be verified through experimental semivariograms (Figure 2e) and descriptive statistics (Table 1).These cations also presented high dependence continuity in the B horizon related to the superficial layer.Figure 2e shows a range decrease of them in the A horizon since the faster the semivariogram rises, the more discontinued the regionalized variable is (Guerra, 1988).Similar to SOM, events that happen on soil surface, such as management operations and erosion, influence the natural dispersion of these chemical elements and lead to a greater apparent homogeneity (compared to the surface layer) in the soil sub-superficial layer.On the other hand, Corá et al. (2004) found that the constant fertilizer and limestone applications and intensive soil preparation are the main reasons for greater continuity and spatial dependence in the topsoil layer.
In the topsoil layer, exchangeable aluminum presented a Gaussian semivariogram with moderate SDI and range higher than the exchangeable aluminum + hydrogen (Figure 2f), adjusted to the exponential model (Table 2).With higher sampling intensity compared to the present study, Rufino et al. (2006) did not find spatial dependence for exchangeable aluminum, calcium, magnesium, pH and sum of bases in soils under Eucalyptus stands in Luís Antônio (SP state, Brazil).Exchangeable sodium did not present spatial dependence.Variables of pH presented greater nugget effect and sill in the A horizon (Table 2) which at greater distance, presented a greater variance dispersion in relation to the B horizon (Figure 2g).The spatial dependence range was greater in the A horizon for pH H 2 O, which was adjusted to an exponential model with R 2 equal to 0.96.In the A horizon, the delta pH was fitted to an exponential model, presenting a R 2 equal to 0.73 and a strong SDI.Rahman et al. (1996) found low spatial dependence for pH in forest soils, adjusting semivariograms to a linear model with undetermined range.In soils covered by Pinus nigra stands, Basaran et al. (2006) found a pure nugget effect for pH and low SDI for SOM.
The sum of bases and CEC presented moderate spatial dependence (Table 2).In the top soil layer, the effective CEC and CEC pH7 were represented by an exponential model (R 2 > 0.8), and the sum of bases presented spatial continuity according to a spherical theoretical variogram with a high adjustment index (R 2 = 0.97).In the B horizon, the sum of bases was adjusted to an exponential model and CEC pH7 was represented by a Gaussian function.
Soil attribute maps were interpolated by ordinary kriging, which produces a mediator effect (Figure 3).This means that the kriging estimator tends to overestimate low values and underestimate high values.The higher the dispersion around the calculated mean is, the higher is the relief.
The legends of the maps were developed by the criterion of equal areas, i.e., each class of the legend has the same area of occurrence.Thus, some soil properties from each region can be easily seen.More than 50 % of the area presented soils with clayey texture (Figure 3c), about 25 % of the area presented values greater than 0.7 for silt/clay ratio, suggesting limitations of the current Rhodic Hapludox distribution (Figure 3d).SOM amounts above 30 g kg -1 occurred in about 25 % of the area (Figure 3f).Phosphorus content was considered low (van Raij et al., 1996) in approximately 25 % of the area (Figure 3g).More than 50 % of the soils presented exchangeable potassium content under the critical point (Gonçalves, 1995) (Figure 3h).Exchangeable calcium and magnesium contents show that most regions need limestone input (Figures 3i and j).This kind of maps provides an overview of the area and confirms the existence of distinct management zones, which deserve special attention.Moreover, the maps obtained will adjust ecophysiological spatial models of wood productivity (Landsberg and Gower, 1997;Almeida et al., 2010), for classifying the quality of forest sites (Fisher and Binkley, 2000), for the estimation and mapping of erosion (Brady and Weil, 2008) and for management of nutrients budgets in areas of forest plantations.

Conclusions
The most frequent distribution was lognormal comprehending 63 % of the variables.Sampling density was insufficient to identify spatial variability of the following soil attributes: coarse sand (B horizon), natural clay (B horizon), phosphorus (A horizon), exchangeable sodium (A horizon), aluminum (B horizon), pH KCl (B horizon) and effective CEC (B horizon).Soil texture had higher spatial dependence in the A horizon whereas nugget and sill were superior in the B horizon.The soil nutrients presented wider range in the B horizon likely due to intensive soil management and topsoil erosion.

Figure 1 -
Figure 1 -Map identifying the range of pedologic collection; points recorded and soil type along the Capão Bonito forest district.The legend of soils is made according to EMBRAPA (2006), FAO (1998) and Soil Survey Staff (1999).

Figure 2 -
Figure 2 -Scaled omnidirectional experimental semivariograms of the Capão Bonito forest soil chemical and physical variables.

Figure 3 -
Figure 3 -Maps of some physical and chemical variables of Capão Bonito forest soils.Except for the silt/clay ratio map, all of the maps refer to A horizon.

Table 1 -
Results of descriptive exploratory analyses for soil attributes at 0 -30 cm depth (A horizon) and at 30 -80 cm depth (B horizon).

Table 2 -
Models, parameters and quality of experimental semivariograms adjusted to Capão Bonito forest soil data base in the 0 -30 cm layer (A horizon) and 30 -80 cm layer (B horizon).