The use of Pedotransfer functions and the estimation of carbon stock in the Central

Computer models have been used to assess soil organic carbon (SOC) stock change. Commonly, models require to determine soil bulk density (Db), a variable that is often lacking in soil data bases. To partly overcome this problem, pedotransfer functions (PTFs) are developed to estimate Db from other easily available soil properties. However, only a few studies have determined the accuracy of these functions and quantified their effects on the final quality of the spatial variability maps. In this context, the objectives of this study were: i) to develop one PTF to estimate Db in soils of the Brazilian Central Amazon region; ii) to compare the performance of PTFs generated with three other models generally used to estimate Db in soils of the Amazon region; and iii) to quantify the effect of applying these PTFs on the spatial variability maps of SOC stock. Using data from 96 soil profiles in the Urucu river basin in Brazil, a multiple linear regression model was generated to estimate Db using SOC, pH, sum of basic cations, aluminum (Al+3), and clay content. This model outperformed the three other PTFs published in the literature. The average estimation error of SOC stock using our model was 0.03 Mg C ha−1, which is markedly lower than the other PTFs (1.06 and 1.23 Mg C ha−1, or 15 % and 17 %, respectively). Thus, the application of a non-validated PTF to estimate Db can introduce an error that is large enough to skew the significant difference in soil carbon stock change.


Introduction
Pedotransfer functions (PTFs) are predictive models of certain soil properties using data from soil surveys (Bouma, 1989).These functions fill the gap between the available soil data and the properties that are more useful or required for a particular model or quality assessment.In this work, PTFs are used as physical-mathematical models that allow the estimation of D b from soil data, which are needed to convert carbon contents from percentage of dry weight to carbon mass per unit of area (Howard et al., 1995;Benites et al., 2007).
The lack of D b data is a limiting factor in regions such as the Brazilian Amazon.With a few exceptions, the lack of detailed maps along the region is partly due to: a) accessibility (often only possible by boat or airplane), b) a vast territory extension is covered by the Tropical Amazon Forest, and c) lack of more detailed ancillary maps, including topographic, geologic, and geomorphologic maps.Bernoux et al. (1998) and Tomasella and Hodnett (1998) provided the first baseline to predict D b from databases.These authors used data from the RADAMBRASIL project (RADAMBRASIL, 1978).More recently, Benites et al. (2007) used data from the Soil Archives of Embrapa (Brazilian Corporation of Agriculture Research) to develop a PTF to estimate D b for most Brazilian biomes.These three PTFs are frequently applied to predict D b of soils in Brazil (Bernoux et al., 2002).
Typically, PTF-estimated attributes are used directly in numerical modeling, although the accuracy of the attribute estimates and their effect on modeling results are not often investigated or they are simply ignored (Deng et al., 2009).Despite the importance of these SOC stock estimations for Brazilian soils in the Amazon region, it is important to highlight that most studies that have applied the PTFs above to estimate D b do not present quantitative information about the error caused by the application of these PTFs.Bernoux et al. (1998), reported that D b estimation using their PTF could lead to an error of ≤ 10 % in the final calculation of carbon stock (CS) in 323 horizons used in the study.
In this context, the purpose of this study was to predict D b from readily available soil properties of Içá Formation in the Brazilian Amazon region, considering the minor effort and uncertainty principles.We have also compared the performance of the models generated in this study with those in the literature (Bernoux et al., 1998;Tomasella and Hodnett, 1998;Benites et al., 2007).Finally, we quantified the error due the application of different PTFs to estimate D b in the spatial variability of soil carbon stock.

Study site, soils, and the database
The study site is located in the central region of the Amazon State near the Urucu River in the municipality of Coari, Brazil, more specifically between the geographic coordinates 4º 45'S and 65º 25'W and the average elevation is 60 m above the sea level (Figure 1).The region is located about 640 km from the state capital Manaus, and it can only be accessed by boat Amazon region or airplane.The climate is equatorial (Af-Köppen Climate Classification) where the temperature in the coldest month is higher than 20 °C, with no pronounced dry period and a mean annual precipitation of 2500 mm.
The sampled area covers about 80 km 2 , spanning longitudinally across the Içá Geologic Formation, which is part of the Solimões Geologic Domain.The soils of the study site were formed from sediments of Içá Geologic Formation (RADAMBRASIL, 1978).The Içá Formation consists of fine to medium sandstone and siltstone with clay conglomerates and yellow-red.The Holocene alluvium of the Quaternary Period deposits is related to the current Amazonian drainage networks.The sediments of Içá Formation cover an area of 563,264 km 2 (36 % of the Amazon State) and were deposited in the Tertiary-Quaternary Period.
A soil survey was conducted in the Oil Province of Urucu River (named Geólogo Pedro de Moura) between the years of 2008 and 2009.This work resulted in the generation of a soil map along with its respective report (Ceddia et al., 2015), which covers an area of 79.665 km 2 (Figure 1).Throughout the soil survey, 96 soil profiles were described and sampled by horizon, totaling 483 horizons/samples.Due to the limitations imposed by the native vegetation, the 315 field observations were restricted to the vicinity of access roads and only data from 96 soil profiles were used for this study.
The soils were classified based on the Brazilian Soil Classification System (Embrapa, 1999).The soilmapping units of the study site as well as the number of profiles and area of occurrence are shown in Table 1.
Undisturbed soil sampling for D b was done by the core method using standard steel cylinders of 53 cm 3 volume (h = 42 mm, d = 40 mm).In each of the 96 soil profiles, the steel cylinders were inserted into the center of each soil horizon (perpendicular to the surface).The soil-filled cylinder was carefully removed from the ring holder and the soil extending beyond both cylinder ends was trimmed flush using a sharp knife.Protective plastic covers were used to prevent samples from drying out.The samples were transported to the laboratory and were oven dried (105 °C) until constant weight (24-48 h).

Development of the pedotransfer function
The first step in the PTFs generation was the selection of the model development (training) and validation of the dataset.The 96 soil profiles were randomly split (80/20) into training and validation sets.Thus, the dataset used for training consisted of 77 soil profiles (adding up to 378 horizons), while the dataset used for validation consisted of 19 soil profiles (adding up to 105 horizons).The spatial distribution of the training and validation soil profiles along the study area are presented in Figure 1.
The stepwise multiple regression routine was used for explanatory analysis relating D b to soils attributes.The Akaike Information Criterion (AIC) with a p-value of 0.05 was used to include or exclude variables.All linear model assumptions were checked (multicollinear-  ity, homoscedasticity and the normality of regression residues).Multicollinearity was minimized by removing variables with variance inflation factors > 4. To verify the premise of homoscedasticity was performed the Breusch-Pagan test and regression normality residues was performed by the Kolmogorov-Smirnov (KS).The stepwise regression was performed by the R software (R Development Core Team, version 3.1.1).
The explanatory analysis was conducted considering two possibilities: a) construction of a unique regression model to estimate soil bulk density for all soil depths, and b) construction of two regression models, one for surface horizons (A and AB) and another for subsurface horizons (BA, B, BC, C).The final choice of the model to be used was based on the evaluation of the indices AIC, R 2 and standard error (SE) of the stepwise regression.Therefore, the best model presented the highest R 2 and the lowest SE and AIC.

The evaluation of the pedotransfer function performance
This study compares the applicability of the proposed multiple regression model with three existing models in the literature (Bernoux et al., 1998;Tomasella and Hodnett, 1998;Benites et al., 2007), which are presented in Table 2.
We simulated a situation where the database of the study area lacked D b data.Thus, to estimate the carbon stock, it would be necessary to use a pedotransfer function available in the literature.Commonly, when D b data are not available in soils of the Amazon region, researchers choose one of these three PTFs, although the criteria to choose a specific PTF is subjective and not clearly explained.Bernoux et al. (1998) and Tomasella and Hodnett (1998) generated the first two PTFs to predict D b data from properties for soil across the Amazon basin and both used the soil data set generated by the RADAMBRASIL project (RADAMBRASIL, 1978).More recently, Benites et al. (2007) generated a model to predict D b from readily available soil properties of Brazilian soils found in most biomes.The latter study constructed a database from the Soil Archives of Embrapa Solos in Rio de Janeiro, Brazil.
The predicted values (yi) of D b, using different PTFs were compared with the 104 observed values (ŷi) in the validation dataset.The differences between the 104 predicted and observed values (ŷ i -y i ) were used to calculate the following error measurements: the mean prediction error (MPE), Eq. (1); the root mean squared prediction error (RMSPE), Eq. (2); the mean absolute error (MAE), Eq. (3); and the prediction coefficient of determination (R 2 ), Eq. ( 4).The results were also evaluated graphically by the ratio 1:1 of the observed vs predicted values.
COV y y Var y Var y = (4) where: ŷ i is observed D b of ith soil sample; y i is predicted D b of ith soil sample, and n is total number of observations.
The MPE (accuracy error) enables the evaluation of an average tendency for overestimation (positive values) or underestimation (negative values).The closer to zero MPE is, the greater the model accuracy.MAE is an accuracy indicator, but it does not reveal the trend to over-or underestimation.This is because it uses an absolute value because of the difference between the observed and predicted data.The RMSPE should be zero, when a perfect fit between the observed and predicted data is achieved.

Spatial variability map of soil organic carbon stock using measured and estimated values of D b by PTFs
For each of the 96 soil profiles, the calculation of the SOC stock was performed at 0-100 cm depth.The classical way of calculating SOC stock (C mass per area) for a given depth consists of summing C stocks by horizon, determined as a product of D b , SOC content, and horizon thickness (Eq.5), according to Bernoux et al. (2002): where: SOC stock is soil organic carbon stock (kg C m −2 ); SOC is soil organic carbon (g kg −1 ); D b is the soil bulk density (Mg m −3 ) and, T is the horizon thickness (m).
Table 2 − Soil bulk density estimate models from the literature (Benites et al., 2007).In the soil survey, the soil profiles were divided into horizons (A, B, and C).In most cases, the calculations concerned two horizons where the first horizon was typically entirely above 100 cm, and the second one crossed this 100 cm depth.When a horizon crossed the 100 cm boundary, only the horizon portion above that depth was used to calculate its SOC stock.
Considering the measured data , as well as the four estimates of D b (model generated in this work and three other published in the literature (Bernoux et al., 1998;Tomasella and Hodnett, 1998;Benites et al., 2007), a total of five dataset of SOC stocks at the layers 0-100 cm depth were generated.For each of the five dataset of SOC stocks, experimental semivariograms were calculated for spatial dependence evaluation and a theoretical model that best represented the data variability was set.The experimental semivariogram, g(h), of n spatial observations Z(xi), i = 1, ... n, was calculated using equation 6: where: N(h) = the number of observations separated by distance h; Z xi = the soil attribute value measured at a specific point (x 1 ) of the grid; Z xi+h = the soil attribute value measured at a specific neighbor point apart by distance h (xi+h).
The theoretical model was validated by the Jack-knife tool (self-validation or cross-validation).Thereafter, five spatial variability maps of SOC stocks were generated by the ordinary kriging (OK) method.OK only uses primary data such as SOC stocks measured at sampled locations u to estimate SOC stocks at unsampled locations (Wackernagel, 2003).For the study site, SOC stock is the primary variable Z i (u), measured at sampled locations u to estimate SOC stocks at unsampled locations (Z u OK * ( )).The stationarity of the mean is assumed only within a local neighborhood W(u), centered at the location u being estimated.Here, the mean is deemed to be a constant but unknown value, i.e., m(u´)=constant but unknown, ∀u´∈ W(u).The OK estimator (Eq.7) is written as a linear combination of the n(u) data Z i (u) with a single unbiasedness constraint (Eq.8), as below: The unknown local mean m(u) is filtered from the linear estimator by forcing the kriging weights (l) to sum to 1 (Eq.8).The weights l are chosen so that the estimate Z u OK * ( ) is unbiased, and that the estimation variance σ OK u ( ) 0 2 (Eq.9) is less than any other linear combination of the observed values.The minimum variance of Z u OK * ( ) is given by: and is obtained when where: g(u i , u j ) = the semivariance of z between the sampling points u i and u j ; g(u 1 , u 0 ) = the semivariance of z between the sampling point u i and the unvisited point u 0 ; Both quantities g(u i , u j ) and g(u 1 , u 0 ) are obtained from the theoretical model fitted to the experimental semivariogram; μ = the Langrange multiplier required for the minimization.
The semivariogram calculation, Jack-knife and ordinary kriging procedures were conducted using the software Geoestat (Vieira at al., 1983).The kriged files were exported to software ArcGIS (ESRI, version 9.3), where the spatial variability maps of SOC were made.

The evaluation of the spatial variability error of SOC stock from measured and estimated values of D b by PTFs
In order to compare the maps generated from different data of D b , spatial analysis through map algebra (ArcGIS Raster Calculator Function) was performed.The spatial variability map of carbon stocks generated from the measured data (96 SOC stock values using measured values of carbon content and D b ) was considered the most appropriate map (reference map).Then, we calculated how much of each PTF used to estimate the D b would over-or underestimate the soil carbon stock.This procedure allowed determining the carbon stock residual (up to 1 m) from the use of each PTF.The residuals were calculated subtracting (pixel by pixel) the carbon stock values in the reference map, in relation to the others generated from the application of PTFs.For this study, each pixel has an area of 1 ha (resolution of 100 by 100 m).

Descriptive statistics and the correlation between D b and other soil attributes
Descriptive statistics of soil properties are shown in Table 3.The 483 soil horizons covered a wide range of soil textural classes, although sandy and franco silty were the most predominant classes.We highlight the relatively high silt content (mean value of 339, reaching up to 721 g kg −1 ), since this granulometric fraction is not commonly higher than 200 g kg −1 in the main classes of Brazilian soils.However, along the Central Amazon Region, the soils are formed over sediments of Içá and Solimões formations (consisted of fine to medium sandstone and siltstone).High values of silt were already reported by Tomasella and Hodnett (1998).The authors found silt values in their dataset that reached up to 800 g kg −1 .Considering the chemical attributes, the soils Sci.Agric.v.74, n.6, p.450-460, November/December 2017 Içá Formation, as well as the type and diversity of vegetation and consequently the carbon supply to the soil (RADAMBRASIL, 1978;Ceddia et al., 2015).The coefficient of variation for D b was 16 %, slightly higher than that previously reported for the Amazon Basin (Moraes et al., 1995;Bernoux et al., 1998).Moraes et al. (1995) observed CV values of 7 % for Alfisols and Ultisols and 13 % for Oxisols, whereas others (Bernoux et al., 1998) have reported a CV of 10 % and 11 % for Oxisols and Alfisols and Ultisols, respectively.

The regression model generated and the comparison with previously published models
The results of the stepwise analysis to select the regression models and respective predictors for D b for all soil depths, as well as for the surface and sub surface horizons, are presented in Table 4.All models generated showed a normal distribution of the regression residuals (Kolmogorov-Smirnov test).With the exception of the regression model for soil surface layer, the models (all horizons and subsurface) fulfilled the homoscedasticity premise as confirmed by the Breusch-Pagan test (Table 4).
The D b varies according to soil depth (Harrison and Bocock, 1981;Leonaviciute, 2000) however, our results show no improvement in the prediction capacity of the models after doing a separate calibration of superficial horizons from the sub-superficial horizons.Similarly, De Vos et al. (2005) did not find any significant enhancement in the predictive capacity when datasets from top and subsoil layers in forest soils of the Flanders region (Belgium) were separated.Other authors working with different regions in the world, such as Han et al. (2012) in China, Sequeira et al. (2014) and Heuscher et al. (2005) in the United States, also did not find any improvement in the prediction of D b generating separated regression model for surface and subsurface soil layers.The PTFs published in the literature, which are recommended to be applied to the Amazon region (Bernoux et al., 1998;Tomasella and Hodnett, 1998) and to most Brazilian biomes (Benites et al., 2007), are used for all soil depths.However, more recently the performance of models was evaluated in soils of most Brazilian biomes presented high acidity and aluminum toxicity and low sum of basis (SB).As characteristics of the study site, the mean Al +3 content was 4.03 cmol c dm −3 (ranging from 0.25 to 12.00 cmol c dm −3 ).All chemical properties, except for pH measurements, had a CV > 45 %.The SOC contents ranged from 0.10 to 60.30 g kg −1 and a CV of 74 %.
The skewness and kurtosis coefficient could be used to infer about the normal data distribution (symmetric histogram).A zero value for both coefficients means that the attribute presents a normal distribution.Although there is no clear-cut guidelines, most studies consider data to be approximately normal in shape if the skewness and the kurtosis values range from -1.0 to +1.0 (Huck, 2012).Considering this range for skewness and kurtosis, with the exception of H + and SOC content, all soil properties follow within a normal distribution.
The mean D b measured was 1.25 Mg m −3 , with minimum and maximum values of 0.49 and 1.67 Mg m −3 , respectively.The lowest D b value is due to the highest amount of soil organic matter (SOM) in well drained soils (complex CXa1-typic dystrudepets and PVAa-typichapludults) as opposed to poorly drained soils (consociation PACd-typicendoaquults) in the Amazon Forest (Ceddia et al., 2015).The relief is considered the main factor influencing the variability of soil types and their attributes in   1.00 1.00 Surface horizons refer to A and AB horizons; Sub-surface horizons refer to BA, B, BC and C horizons; Int.= Intercept; SOC = soil organic carbon; AIC = Akaike Information Criterion; SE = Standard error of the estimation; KS (p value) -values higher than 0.05 mean that the regression residues are normally distributed (KS = Kolmogorov-Snirnov test); BP (p value) -values higher than 0.05 means that the regression residues are homocedastic (BP = Breusch-Pagan test); SB = sum of basic cations; VIF = Variance Inflation Factor; As a reference, when the VIF value of the soil attribute is lower than 4, it means that the respective attribute does not present colinearity with any other predictor variable.
Sci. Agric.v.74, n.6, p.450-460, November/December 2017 after separation for surface (to 30 cm depth) and subsurface (below 30 cm) (Benites et al., 2007).The authors did not found any advantage of partitioning the dataset into groups of soil depth and soil orders (Benites et al., 2007).Here, we highlight the model number 1 in Table 4 (Eq.11) that could be used for all soil depths.This model explained 67 % of the variance and presented the lowest value of AIC and the highest value of R 2 .The SOC was the main predicted variable followed by pH H2O , SB, Al +3 and Clay.D b = 1.171 -0.0237(SOC) + 0.0622(pH) -0.0230(SB) -0.0124(Al +3 ) + 0.0002(Clay) (11) Consistent with the findings of many authors, SOM clearly plays a dominant role (Adams, 1973;De Vos et al., 2005;Heuschek et al., 2005;Han et al., 2012) because of its much lower density than mineral soil particles and its aggregation effect on soil structure.The D b strongly correlates with SOM content and soil texture (Adams, 1973;Manrique and Jones, 1991).Barros and Fearnside (2015), working with oxisoils, also in the Central Amazon region, found that the clay content accounts for about 70 % of the variation in soil bulk density.However, the importance of each variable depends on the study site.For example, De Vos et al. (2005) found that the addition of texture as a predictor has a minor effect on D b estimation of forest soils.In fact, the authors observed that SOM accounted for 55 to 57 % of the total variation in D b , whereas soil texture explained only 20 to 26 %.
In the literature, the relationship between D b and chemical attributes (such as pH, SB and Al +3 ) is scarcely reported and a direct physical link is not clearly presented.In most cases, these attributes are inserted into the model due to their availability in most datasets, because chemical attributes are determined at low costs (Bernoux et al., 1998).For instance, Bernoux et al. (1998) used pH in water and Benites et al. (2007) used SB to estimate D b.However, these authors did not provide a physical explanation regarding the relationship between either pH or SB with D b .In tropical soils, the main clay minerals are kaolinite, goethite, hematite and gibbsite, which are colloids with variable charges depending on the pH.The main clay minerals in the study site are kaolinite and goethite, consequently, as the pH is very low (4.58 on average) and organic colloids are mainly responsible for the CEC of the soil, cation bridges are formed between soil particles.Cation flocculants (Al +3 , Fe +3 and H + ) promote the approximation of the colloidal particles, which is the first step for the aggregation formation.The next step of aggregate formation is cementation, where SOC has an important role, acting as an agent that cements flocculated particles (Tisdall and Oades, 1982).Finally, soil macro-aggregates are formed, which improves the soil porosity and consequently reduces D b .In uncultivated soils, like most forest soils, the OM has a dominant effect over D b (Adams, 1973) and, naturally, becomes the main predictor variable (De Vos et al., 2005).Considering that the main colloids of the study site in charges that depend on the pH, we hypothesize that Al +3 and H + availability increases as the pH is lowered.Therefore, colloids flocculate and increase soil porosity and decrease soil bulk density.
Compared to other PTFs, the proposed model (Eq.11) presented the best agreement (Figures 2A, B, C  and D).Considering the MPE index, our proposed model presented the lowest bias, once it reached the ideal value (MPE = 0).The model proposed by Tomasella and Hodnett (1998) tended to underestimate D b (MPE = -0.03kg dm −3 ), whereas the models by Benites et al. (2007) and by Bernoux et al. (1998) 5).Barros and Fearnside (2015) also developed PTFs for Oxisols in the Central Amazon region and compared the performance of their model with those presented by authors reported in this work (Benites et al., 2007;Bernoux et al., 1998;Tomasella and Hodnett, 1998).The authors found that the application of these three PTFs overestimate soil bulk density, which is in agreement with our findings.However, we highlight that the soil bulk density determined by Barros and Fearnside (2015) is significantly lower than what we found in the region of the Urucu River (average and median values of 0.66 and 0.62, respectively).This can explain why even the model developed by Tomasella and Hodnett (1998) also overestimated soil bulk density in the Oxisols evaluated by Barros and Fearnside (2015).
Therefore, the model developed using local data of the study site to predict D b outperformed the three others, both in terms of accuracy and precision.These results confirm that it is very difficult for one particular PTF to be precise and accurate to predict a soil attribute in a vast territory.This is particularly true in the Amazon region that encompasses different ecosystems within an area of 5,020,000 km 2 .The low performance of the other three models is probably a consequence of the low-density data used to generate the regression models (Table 2).These models did not capture adequately the high variability of the Amazon soils.The better performance of the model generated in this study was expected, since local models can capture more effectively the specificities of the relationship among soil attributes.This is also why PTFs should not be generalized to wider regions that do not follow the assumptions underlying the relationship between the target and predicting variables.Considering these aspects, it is important to highlight that the PTF generated in this work should only be used in regions belonging to the Içá Formation and that are covered with forest.

Applying the measured and predicted values of D b to model the spatial variability of SOC stock
As the measurement of D b is essential to predict SOC stock, the impact of using PTFs in the spatial variability map of SOC stock was evaluated.Thus, at this stage, we compared the semivariogram parameters and maps of SOC stock performed from the measured data with those generated when applying the four PTFs.The application of geostatistical techniques requires that the spatial dependence between the observations of SOC stock be proved.For this, the experimental semivariogram should be calculated in order to set a theoretical model that best represents data variability (Figures 3A, B, C and D).
The SOC stock at 0-100 cm exhibited spatial dependence.In fact, for all experimental semivariograms, both the lag distance (h) and semivariance g (h) increased until they reached an approximately constant value called sill variance (known as the priori variance of the random variable (Table 6)).The same theoretical model (spherical) was fitted for all the experimental semivariograms, which differ from each other in the nugget parameters (C 0 ), structural variance (C 1 ) and range (a).The semivariogram of reference (Figure 3A) is the one calculated from the measured data of D b and soil organic carbon, as these are the best data we have to evaluate the SOC stock variability.In general, the parameters of the semivariogram using the pedotranfer generated in this study were close to those of the reference semivariogram (Figure 3B).The experimental semivariogram calculated using D b estimated from the PTF proposed by Tomasella and Hodnett (1998) presented the highest range (Figure 3D).The range is the lag distance at which the semivariogram reach its sill.This is the spatial dependence and beyond it, the variance bears no relation to the separation distance (Webster and Oliver, 1990).On the other hand, the experimental semivariogram calculated using D b estimated using the PTF proposed by Benites et al. (2007) and Bernoux et al. (1998), presented   the lowest ranges (spatial dependence of 2180 and 2400 m, respectively) (Figures 3C and E).
The semivariogram with the highest nugget effect value (3.30) was the one generated from the application of the PTF proposed by Tomasella and Hodnett (1998), while the lowest nugget effect value (2.38) was the one developed by Benites et al. (2007).This implies that the SOC stock has a higher initial variability when using the model proposed by Tomasella and Hodnett (1998).
The structured variability using ordinary kriging allowed the generation of the spatial variability map of SOC stocks up to 100 cm of soil depth (Figures 4A, B, C,  D and E).The SOC stock map generated with measured values of D b and SOC (reference map) ranged from 4.89 up to 9.93 kg C m −2 (Figure 4A).These values are similar to those found in the literature for the Amazon region.Ceddia et al. (2015), using different results published in the literature, compared the average values of SOC stock up to the 100 cm soil depth for the same region and found that the estimative of SOC stock ranged from 7.32 to up to 9.01 kg cm −2 .
The SOC stock variability maps performed with D b estimated by our PTF and the PTFs presented by Tomasella and Hodnett (1998) (Figures 4B and D, had amplitude values closer to those in the reference map.The SOC stock variability maps using the D b estimated by PTFs proposed by Benites et al. (2007) and Bernoux et al. (1998) tended to overestimate these values (Figures 4C and E).These maps showed the highest upper limit values of SOC stock (13.96kg cm −2 and 12.81 kg cm −2 , respectively).Bernoux et al. (2002) argue about the uncertainty sources in estimating SOC stock.The first source related to different database information used in the estimates.The second source referred to the error associated with the estimation of D b by PTFs.However, according to the authors, the most important source of uncertainty originates from the SOC analytical methods.Therefore, they reaffirm the need for a more complete documentation (metadata) about the database used to generate PTFs.According to the authors, in addition to the statistical data information, the database must contain which methods were used to obtain the attributes used as predictors.The authors refer to (Garten and Wullschleger, 1999) as traditional method in soil science and claim that this method is not completely accurate and is considerate the main source of uncertainty in estimating SOC.We also observed in this study that variations in D b through ap-

The residuals caused by applying different PTFs on the spatial variability maps of SOC stock
The spatial variability maps of the residuals as well as the statistics of the residuals are shown in Figures 5A,  B, C and D and Table 7, respectively.The residual classes shown in Figures 5A and C confirmed that the PTFs generated in this study and those developed by Tomasella and Hodnett (1998) presented the lowest residual and tend to underestimate even more the SOC stock, mainly when compared to the maps generated by Benites et al. (2007) and Bernoux et al. (1998).These results are in accordance with the values of the MPE index previously discussed in item 3.2.We highlight the overestimation of SOC stock when applying the PTFs generated by Benites et al. (2007) and Bernoux et al. (1998).Using these PTFs, the overestimation of SOC stock could reach up to 4.20 and 2.95 kg cm −2 , respectively (Figures 5B and D).The mean value of the SOC stock applying the PTF function generated in this study is practically the same of that in the reference map (7.49Mg C ha −1 ).On the other hand, applying the PTFs developed by Benites et al. (2007) and Bernoux et al. (1998), resulted in mean values of SOC stock of 8.56 and 8.73 Mg C ha −1 , respectively.When applying these PTFs, the mean residual value of SOC stock in relation to the reference map is 1.06 and 1.23 Mg C ha −1 , a mean overestimation of 15 and 17 %, re-  spectively (Table 7).Considering the residuals found in this work, the error in the estimation of SOC stocks, using different PTFs to estimate D b could be higher than 10 %, as previously reported by Bernoux et al. (1998).
The results here presented can also give support to answer a significant and generic problem with the estimation of changes in terrestrial biospheric carbon, which is the smallest detectable change.The smallest difference in SOC stock that could be detected after 5 years under an herbaceous bioenergy crop was about 1 Mg C ha −1 in the southeastern United States (Garten and Wullschleger, 1999).Considering the average residual of SOC stock caused by applying the PTFs developed by Benites et al. (2007) and Bernoux et al. (1998) in the study site, it is possible to conclude that the simple application of a PTF to estimate D b could introduce an error enough to offset the smallest detectable change of the SOC stock estimate.

Conclusions
A linear regression model was generated to estimate D b for soils of the Central Amazon region, Brazil.The developed PTF used soil attributes easily found in soil survey reports such as SOC, pH in water, sum of basic cations, Al +3 , and clay content as predictor of D b .The PTF generated should only be used in regions belonging to the Içá Geological Formation and covered with forest.
Compared to the three most used PTFs to estimate D b in Amazon soils (Bernoux et al., 1998;Tomasella and Hodnett, 1998;Benites et al., 2007), our PTF outperformed the other models, showing not only better accuracy, but also the lowest bias.
We also quantified the effect of using PTFs to estimate the D b on spatial variability of SOC stock.The results showed that the uncertainty caused by the estimative of D b using PTFs could be much higher than the values reported in the literature.The PTFs of Benites et al. (2007) and Bernoux et al. (1998) caused an overestimation of 1.06 and 1.23 Mg C ha −1 in the SOC stock, which represented 15 % and 17 %, respectively.This means that the simple application of a PTF to estimate D b could introduce an error large enough to skew the significant difference in soil carbon stock change.

Figure 1 −
Figure 1 − The location of the study site in the Central Amazon (Brazil) and the spatial distribution of the soil profiles used for training and validation of the pedotransfer function.
tended to overestimate D b (0.11 and 0.15kg dm −3 , respectively).The MAE value ranged between 0.09 and 0.19 kg dm −3 , and the highest value was observed in the Bernoux et al. (1998) model, followed by the Tomasella and Hodnett (1998) and the Benites et al. (2007) models.The Bernoux et al. (1998) model presented the highest RMSPE value and our model the lowest (Table

Figure 3 −
Figure 3 − Carbon stock semivariograms using measured and estimated data of D b .A) Semivariogram generated with the measured data of D b ; B) Semivariogram generated applying the D b estimated using the PTF proposed in this work; C) Semivariogram generated applying the D b estimated using the PTF proposed by Benites et al. (2007); D) Semivariogram generated applying the D b estimated using the PTF proposed by Tomasella and Hodnett (1998); E) Semivariogram generated applying the D b estimated using the PTF proposed by Bernoux et al. (1998).Shp = Spherical model (nugget effect, contribution and range values); PTFs = Pedotransfer functions; SOC = soil organic carbon.

Figure 4 −
Figure 4 − Spatial variability maps of SOC stock using ordinary kriging.A) Map generated with the measured data of D b ; B) Map based on the D b estimated using the PTF generated in this work; C) Map based on the D b estimated using the PTF generated by Benites et al. (2007); D) Map based on the D b estimated using the PTF generated by Tomasella and Hodnett (1998); E) Map based on the D b estimated using the PTF generated by Bernoux et al. (1998).SOC = Soil Organic Carbon; PTFs = Pedotransfer functions.

Figure 5 −
Figure 5 − Spatial variability maps of the SOC stock residual.A) Residual based on the D b estimated by the PTF generated in this work; B) Residual based on the D b estimated using the PTF generated by Benites et al. (2007); C) Residual based on the D b estimated of the PTF generated using Tomasella and Hodnet (1998); D) Residual based on the D b estimated using the PTF generated by Bernoux et al. (1998).SOC = Soil Organic Carbon; PTFs = Pedotransfer functions.

Table 1 −
Distribution of soil class in the data set.

Table 3 −
Descriptive statistics of the soil properties.

Table 4 −
Candidate models to predict soil bulk density for all soil depth and dividing the dataset into surface and subsurface horizons.

Table 5 −
Evaluation indices of proposed and existing regression models for D b estimation.

Table 6 −
Models and their parameters fitted for SOC stock's semivariograms.

Table 7 −
Statistics of the SOC stock maps and the residual caused by the application of different PTFs.