Acessibilidade / Reportar erro

Mapping soil carbon, particle-size fractions, and water retention in tropical dry forest in Brazil

Mapeamento de carbono, frações granulométricas e água do solo em Floresta Tropical Seca no Brasil

Abstract

The objective of this work was to compare ordinary kriging with regression kriging to map soil properties at different depths in a tropical dry forest area in Brazil. The 11 soil properties evaluated were: organic carbon content and stock; bulk density; clay, sand, and silt contents; cation exchange capacity; pH; water retention at field capacity and at permanent wilting point; and available water. Samples were taken from 327 sites at 0.0-0.10, 0.10-0.20, and 0.20-0.40-m depths, in a tropical dry forest area of 102 km2. Stepwise linear regression models for particle-size fractions and water retention properties had the best fit. Relief and parent material covariates were selected in 31 of the 33 models (11 properties at three depths) and vegetation covariates in 29 models. Based on external validation, ordinary kriging obtained higher accuracy for 21 out of 33 property x depth combinations, indicating that the inclusion of a linear trend model before kriging does not necessarily improve predictions. Therefore, for similar studies, the geostatistical methods employed should be compared on a case-by-case basis.

Index terms:
caatinga; digital soil mapping; gamma radiometric survey; geostatistics; pedometrics

Resumo

O objetivo deste trabalho foi comparar krigagem ordinária com regressão-krigagem para mapear atributos do solo, em diferentes profundidades, em área de Floresta Tropical Seca no Brasil. Os 11 atributos do solo avaliados foram: conteúdo e estoque de carbono orgânico; densidade do solo; conteúdos de argila, areia e silte; capacidade de troca catiônica; pH; retenção de água na capacidade de campo e no ponto de murcha permanente; e água disponível. As amostras foram retiradas de 327 locais a 0,0-0,10, 0,10-0,20 e 0,20-0,40 m de profundidade, em área de Floresta Tropical Seca de 102 km2. Os modelos de regressão linear "stepwise" tiveram o melhor ajuste para as frações granulométricas e as propriedades de retenção de água. Foram selecionadas covariáveis de relevo e material de origem em 31 dos 33 modelos (11 atributos em três profundidades) e de vegetação em 29 modelos. Com base na validação externa, a krigagem ordinária obteve maior acurácia para 21 das 33 combinações atributo vs. profundidade, o que é indicativo de que a inclusão de um modelo linear de tendência antes da krigagem não necessariamente melhora as predições. Portanto, para estudos semelhantes, os métodos geoestatísticos empregados devem ser comparados caso a caso.

Termos para indexação:
caatinga; mapeamento digital de solos; levantamento gamarradiométrico; geoestatística; pedometria

Introduction

Tropical dry forests correspond to 42% of the world's tropical and subtropical forests (Murphy & Lugo, 1986MURPHY, P.G.; LUGO, A.E. Ecology of tropical dry forest. Annual Review of Ecology and Systematics, v.17, p.67-88, 1986. DOI: 10.1146/annurev.es.17.110186.000435.
https://doi.org/10.1146/annurev.es.17.11...
), of which about 54% occur in South America (Miles et al., 2006MILES, L.; NEWTON, A.C.; DEFRIES, R.S.; RAVILIOUS, C.; MAY, I.; BLYTH, S.; KAPOS, V.; GORDON, J.E. A global overview of the conservation status of tropical dry forests. Journal of Biogeography, v.33, p.491-505, 2006. DOI: 10.1111/j.1365-2699.2005.01424.x.
https://doi.org/10.1111/j.1365-2699.2005...
). In Brazil, and possibly other tropical countries, laws concerning tropical dry forests are ill-defined, grouping them with moister or drier forest types. This creates a conundrum, hindering the proper conservation of this important ecosystem that serves as a habitat for many endemic species, as well as a food and wood source for traditional peoples, but is still little studied. Soil maps can be used to provide information to support government and private sector decisions for promoting the use and conservation of this poorly-known, threatened ecosystem and for guiding future research (Miles et al., 2006MILES, L.; NEWTON, A.C.; DEFRIES, R.S.; RAVILIOUS, C.; MAY, I.; BLYTH, S.; KAPOS, V.; GORDON, J.E. A global overview of the conservation status of tropical dry forests. Journal of Biogeography, v.33, p.491-505, 2006. DOI: 10.1111/j.1365-2699.2005.01424.x.
https://doi.org/10.1111/j.1365-2699.2005...
).

Soils vary both horizontally and vertically in the landscape, and are continuously modified by internal and external factors, generating complex spatial patterns (Grunwald, 2009GRUNWALD, S. Multi-criteria characterization of recent digital soil mapping and modeling approaches. Geoderma, v.152, p.195-207, 2009. DOI: 10.1016/j.geoderma.2009.06.003.
https://doi.org/10.1016/j.geoderma.2009....
). According to Boucneau et al. (1998)BOUCNEAU, G.; VAN MEIRVENNE, M.; THAS, O.; HOFMAN, G. Integrating properties of soil map delineations into ordinary kriging. European Journal of Soil Science, v.49, p.213-229, 1998. DOI: 10.1046/j.1365-2389.1998.00157.x.
https://doi.org/10.1046/j.1365-2389.1998...
, soil spatial predictions can be obtained either by: soil property estimation, using mean and range values from a choropleth soil map, for example; or spatial interpolation, using kriging, for instance. While the former assumes uniformity within mapping units and abrupt changes at their boundaries, the latter assumes that all variation is gradual; however, in reality, soil properties may present both gradual and abrupt changes spatially.

Among spatial interpolation methods, kriging constitutes a best linear unbiased estimator and has been extensively used to map soil properties. One of its variants - regression kriging (RK) - has received special attention as a means to incorporate the variation of soil-forming environmental factors (so-called environmental covariates) into the variation of the target soil property (Knotters et al., 1995 KNOTTERS, M.; BRUS, D.J.; OUDE VOSHAAR, J.H. A comparison of kriging, co-kriging and kriging combined with regression for spatial interpolation of horizon depth with censored observations. Geoderma, v.67, p.227-246, 1995. DOI: 10.1016/0016-7061(95)00011-C.
https://doi.org/10.1016/0016-7061(95)000...
; Kravchenko & Robertson, 2007 KRAVCHENKO, A.N.; ROBERTSON, G.P. Can topographical and yield data substantially improve total soil carbon mapping by regression kriging? Agronomy Journal, v.99, p.12-17, 2007. DOI: 10.2134/agronj2005.0251.
https://doi.org/10.2134/agronj2005.0251...
; Vasques et al., 2010 VASQUES, G.M.; GRUNWALD, S.; COMERFORD, N.B.; SICKMAN, J.O. Regional modelling of soil carbon at multiple depths within a subtropical watershed. Geoderma, v.156, p.326-336, 2010. DOI: 10.1016/j.geoderma.2010.03.002.
https://doi.org/10.1016/j.geoderma.2010....
). Regression kriging integrates remote, field, and laboratory data, and statistical methods in a quantitative estimation framework for soil property or class mapping - digital soil mapping -, with applications ranging from precision agriculture (Grunwald et al., 2015 GRUNWALD, S.; VASQUES, G.M.; RIVERO, R.G. Fusion of soil and remote sensing data to model soil properties. Advances in Agronomy, v.131, p.1-109, 2015. DOI: 10.1016/bs.agron.2014.12.004.
https://doi.org/10.1016/bs.agron.2014.12...
) to global mapping (Hengl et al., 2014 HENGL, T.; JESUS, J.M. de; MACMILLAN, R.A.; BATJES, N.H.; HEUVELINK, G.B.M.; RIBEIRO, E.; SAMUEL-ROSA, A.; KEMPEN, B.; LEENAARS, J.G.B.; WALSH, M.G.; GONZALEZ, M.R. SoilGrids1km - global soil information based on automated mapping. PLoS ONE, v.9, p.e105992, 2014. DOI: 10.1371/journal.pone.0105992.
https://doi.org/10.1371/journal.pone.010...
).

Ordinary kriging (OK) and RK are commonly used interpolation methods in digital soil mapping, and have been compared in some studies for mapping various soil properties. Knotters et al. (1995) KNOTTERS, M.; BRUS, D.J.; OUDE VOSHAAR, J.H. A comparison of kriging, co-kriging and kriging combined with regression for spatial interpolation of horizon depth with censored observations. Geoderma, v.67, p.227-246, 1995. DOI: 10.1016/0016-7061(95)00011-C.
https://doi.org/10.1016/0016-7061(95)000...
obtained better results with RK using apparent electrical conductivity as a covariate than with OK or cokriging, when mapping the depth to soft layers (peat or unripe clay) in a 97-ha area in the Netherlands. Kravchenko & Robertson (2007) KRAVCHENKO, A.N.; ROBERTSON, G.P. Can topographical and yield data substantially improve total soil carbon mapping by regression kriging? Agronomy Journal, v.99, p.12-17, 2007. DOI: 10.2134/agronj2005.0251.
https://doi.org/10.2134/agronj2005.0251...
, mapping total carbon using topographic and crop yield as covariates, observed little or no improvement with RK in comparison to OK in 12 sites of 0.36 ha in Michigan, in the United States. In the same country, in Florida, RK outperformed OK for total carbon in three out of five depth intervals (Vasques et al., 2010 VASQUES, G.M.; GRUNWALD, S.; COMERFORD, N.B.; SICKMAN, J.O. Regional modelling of soil carbon at multiple depths within a subtropical watershed. Geoderma, v.156, p.326-336, 2010. DOI: 10.1016/j.geoderma.2010.03.002.
https://doi.org/10.1016/j.geoderma.2010....
). As reported in previous studies, preference for RK over OK is not granted nor expected, and there is no clear indication of when to use OK, RK, or any other geostatistical method for more accurate results.

The objective of this work was to compare ordinary kriging with regression kriging to map soil properties at different depths in a tropical dry forest area in Brazil.

Materials and Methods

The study was conducted at Parque Estadual da Mata Seca, a dry forest state park, in the north of the state of Minas Gerais, Brazil (43º58'48"W, 14º52'12"S) (Figure 1). The park is located in the ecotone of the Brazilian Atlantic Forest (tropical coastal moist forest), Cerrado (savanna), and Caatinga (xeric shrubland). It was created in 2000 and has been preserved since then. Previously, the land was mainly used for agriculture, cattle ranching, and logging for charcoal production. The study area has 102 km2, encompassing the part of the park prior to its expansion in 2009. Mean annual temperature and precipitation are around 24ºC and 827 mm, respectively.

Figure 1
Location of study area, that is, of Parque Estadual da Mata Seca, a dry forest state park in Brazil (A), as well as digital elevation model (B), and sampling design, soil suborders, and geologic units (C). CX, Cambissolos Háplicos (Haplustepts); CY, Cambissolos Flúvicos (Fluventic Haplustepts); GM, Gleissolos Melânicos (Argiaquolls); GX, Gleissolos Háplicos (Fluvaquents, Endoaqualfs); LA, Latossolos Amarelos (Haplustox); LV, Latossolos Vermelhos (Haplustox, Eutrustox); LVA, Latossolos Vermelho-Amarelos (Haplustox, Eutrustox); MX, Chernossolos Háplicos (Haplustolls, Argiustolls); RY, Neossolos Flúvicos (Ustifluvents, Ustorthents); and VX, Vertissolos Háplicos (Haplusterts) (Soil Survey Staff, 2014SOIL SURVEY STAFF. Keys to soil taxonomy. 12nd ed. Washington: United States Department of Agriculture, 2014.). Soil suborder map compiled from Coelho et al. (2013) COELHO, M.R.; DART, R. de O.; VASQUES, G. de M.; TEIXEIRA, W.G.; OLIVEIRA, R.P. de; BREFIN, M. de L.M.; BERBARA, R.L.L. Levantamento pedológico semi-detalhado (1:30.000) do Parque Estadual da Mata Seca, município de Manga-MG. Rio de Janeiro: Embrapa Solos, 2013. 264 p. (Embrapa Solos. Boletim de pesquisa e desenvolvimento, 217)..

Vegetation in the study area can be roughly divided into three main types, considering plant species, size, and seasonality: typical dry forest, comprised mostly of deciduous and thorny species that dominate the study area, covering about 69% of the area in the central part; riparian forest, consisting of evergreen species that occur in 25% of the area on richer and moister soils in the eastern part; and "carrasco" forest, characterized by a few thorny and contorted species adapted to poor, sandy, and excessively drained soils, occurring in less than 6% of the area.

Elevations range from 434 to 523 m (Figure 1 B), and slope gradients from 0 to 75%. Three geologic units are present in the area, along with their associated soil classes according to Coelho et al. (2013) COELHO, M.R.; DART, R. de O.; VASQUES, G. de M.; TEIXEIRA, W.G.; OLIVEIRA, R.P. de; BREFIN, M. de L.M.; BERBARA, R.L.L. Levantamento pedológico semi-detalhado (1:30.000) do Parque Estadual da Mata Seca, município de Manga-MG. Rio de Janeiro: Embrapa Solos, 2013. 264 p. (Embrapa Solos. Boletim de pesquisa e desenvolvimento, 217). and to Soil Survey Staff (2014)SOIL SURVEY STAFF. Keys to soil taxonomy. 12nd ed. Washington: United States Department of Agriculture, 2014. (Figure 1 C). The "Urucuia" Group takes up about 9% of the area and is mainly composed of sandstone from the Cretaceous Period. It has the poorest soils of the park, which are classified as Latossolos Amarelos (Embrapa, 2006), that is, Haplustox (Soil Survey Staff, 2014SOIL SURVEY STAFF. Keys to soil taxonomy. 12nd ed. Washington: United States Department of Agriculture, 2014.), and also the poorest vegetation, the "carrasco". The "Bambuí" Group was formed during the Neoproterozoic Era and is composed of siltstone, shale, marl, and limestone from the Serra de Santa Helena and Lagoa do Jacaré geologic formations. It covers about 66% of the area, where typical dry forest occurs, and encompasses the highest diversity of soils, including Latossolos (Oxisols), Cambissolos (Inceptisols), Chernossolos (Mollisols), and Vertissolos (Vertisols). Finally, the "Quaternary sediments" unit (25% of the area) comprises fluvial deposits from the São Francisco River, under riparian forest, where the predominant soils are Gleissolos Háplicos (Fluvaquents, Endoaqualfs), Cambissolos Flúvicos (Fluventic Haplustepts), and Neossolos Flúvicos (Ustifluvents, Ustorthents).

Soils were described and sampled at 327 sites, allocated according to three sampling designs, which were carried out in separate field campaigns. In the first design ("survey"), 222 locations were chosen by tacit knowledge for the purpose of soil survey, comprising 44 deep pits (>1 m), 83 shallow pits (<1 m), and 95 boreholes (<1 m). The second design ("cLHS") consisted of 41 shallow pits allocated by conditioned Latin hypercube sampling (Minasny & McBratney, 2006 MINASNY, B.; MCBRATNEY, A.B. A conditioned Latin hypercube method for sampling in the presence of ancillary information. Computers and Geosciences, v.32, p.1378-1388, 2006. DOI: 10.1016/j.cageo.2005.12.009.
https://doi.org/10.1016/j.cageo.2005.12....
), using elevation, normalized difference vegetation index, and land use as covariates. The third design ("random transect") was composed of 64 shallow pits randomly distributed along transects in different directions, cutting across relatively undersampled areas (Figure 1).

At all sites, soils were described and classified to the fourth hierarchical level (subgroup) of the Brazilian Soil Classification System (Embrapa, 2013). Sampling was done by horizon in deep pits and by layer in shallow pits and boreholes, at 0.0-0.10, 0.10-0.20, 0.20-0.40, 0.40-0.60, 0.60-0.80, and 0.80-100-m depths. Because soil sampling served different purposes in the experiment, in each field campaign, the shallow pits were sampled differently, with most samples collected at the first three layers to 0.40 m. Undisturbed samples for bulk density (BD) and water retention were obtained using 100-cm3 steel rings. Soil samples were analyzed in the laboratory according to the methods described by Donagema et al. (2011)DONAGEMA, G.K.; CAMPOS, D.V.B. de; CALDERANO, S.B.; TEIXEIRA, W.G.; VIANA, J.H.M. (Org.). Manual de métodos de análise de solos. 2.ed. rev. Rio de Janeiro: Embrapa Solos, 2011. 230p. (Embrapa Solos. Documentos, 132). for organic carbon (OC), clay, silt, and sand contents (g kg-1); BD (kg dm-3); pH; and water retention (volumetric percentage) at field capacity (WFC, -0.01 MPa) and at permanent wilting point (PWP, -1.5 MPa). Soil organic carbon stock (OCS, kg m-2) was derived by multiplying OC, BD, and layer thickness, corrected for the percent of coarse fragments (>2 mm), using the following equation: OCS = 0.01OC × BD × D × (1 - CF/1000), in which OCS is the organic carbon stock (kg m-2), OC is the organic carbon content (g kg-1), BD is bulk density (kg dm-3), D is layer thickness (cm), and CF is the coarse fragment (>2 mm) content (g kg-1). Available water (AW, volumetric percentage) was calculated by subtracting WFC from PWP.

Environmental data layers representing soil-forming factors were compiled in a geographic information system and extracted from the soil sampling sites to be used as model covariates (Table 1). A digital elevation model (DEM) was derived from an Ikonos (DigitalGlobe Inc., Westminster, CO, USA) stereo pair from December 2012 with 1-m spatial resolution, and then resampled to 10 m using bilinear interpolation (Figure 1 B). Relief covariates were calculated from the resampled DEM. Vegetation covariates were represented by two sets of RapidEye imagery (Planet Labs Inc., San Francisco, CA, USA) with 5-m spatial resolution, one from the dry and the other from the wet season, that is, from August 2012 and May 2013, respectively. The original five bands were radiometrically and atmospherically corrected to surface reflectance, resampled to 10 m by bilinear interpolation, and then used to derive vegetation indices. Parent material covariates included gamma-radiometric and magnetometric layers obtained from an aerial survey (Levantamento..., 2009 LEVANTAMENTO aerogeofísico do Estado de Minas Gerais: área 11A (Jaíba - Montes Claros - Bocaiúva). Rio de Janeiro: Serviço Geológico do Brasil, 2009.) with 125-m spatial resolution, which were resampled to 10 m by nearest neighbor resampling. The covariates were processed in the Saga geographic information system (Saga Development Team, University of Hamburg, Hamburg, Germany).

Table 1
Environmental data layers used as model covariates(1) (1) nir, near infrared; RE, RapidEye sensor; K, potassium; Th, thorium; U, uranium; kWh, kilowatt-hour; μR, microroentgen; ppm, parts per million; and nT, nanotesla. .

The 11 soil properties evaluated were mapped at the 0.0-0.10, 0.10-0.20, and 0.20-0.40-m depths, respectively. Horizon-based data from the 44 deep pits were recalculated for the three layers by depth-weighted averaging. Model training data was composed of all samples from the survey (222) plus cLHS (41) sets. External model validation data were composed of all 64 random transect samples from a separate field campaign. Descriptive statistics of the soil properties were made, and Pearson's linear correlations among them were derived.

Two geostatistical methods were compared to map soil properties, namely OK and RK. In RK, the soil properties were first predicted across the study area as a function of the 62 environmental covariates, using stepwise multiple linear regression (SMLR) as a global trend model. Second, the SMLR residuals were interpolated across the area by OK. Then, the SMLR predictions were added to the interpolated residuals to derive the final output maps. Semivariogram models - spherical, exponential, or Gaussian - were chosen visually and were fitted by ordinary least squares and/or manually.

Stepwise multiple linear regression was implemented using Akaike's information criterion (AIC) with a p-value of 0.05 for including or excluding variables. All linear regression assumptions were checked. Influential outliers were identified and removed by an outlier test consisting of comparing the Bonferroni corrected p-value from a t-test of the Studentized residuals against a threshold of 0.05. Multicollinearity was minimized by removing variables with variance inflation factors >10.

The goodness-of-fit of the models was assessed by the coefficient of determination (R2), whereas prediction uncertainty was determined by mean error (ME) for biasedness and by root mean square error (RMSE) for accuracy. The RMSE of the validation was used to select the best geostatistical method for each soil property at each depth. Data analyses were conducted in the R software (R Core Team, 2015R CORE TEAM. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing, 2015.).

Results and Discussion

Most soil properties had a close to normal distribution, as indicated by their mean, median, skewness, and kurtosis values (Table 2). On average, soils had more sand, followed by clay, then silt. However, particle-size distribution varied according to soil class. Latossolos Amarelos (Haplustox), Latossolos Vermelho-Amarelos (Haplustox, Eutrustox), and Neossolos Flúvicos (Ustifluvents, Usthorthents) had the largest sand content, which decreased with depth. Latossolos Amarelos and Latossolos Vermelho-Amarelos are formed on sand-rich parent material, mainly sandstone, whereas Neossolos Flúvicos are made up of sandy materials deposited on the floodplain of the São Francisco River. Clay content, however, increased with depth, with the largest values found in Vertissolos Háplicos (Haplusterts), Gleissolos Háplicos (Fluvaquents, Endoaqualfs), Cambissolos Háplicos (Haplustepts), and Latossolos Vermelhos (Haplustox, Eutrustox). These soils are formed from clay-richer rocks of the "Bambuí" Group.

Table 2
Descriptive statistics of soil properties(1) (1) N, number of observations; SD, standard deviation; OC, organic carbon content (g kg-1); BD, bulk density (kg dm-3); OCS, organic carbon stock (kg m-2); CEC, cation exchange capacity (cmolc kg-1); clay, silt, and sand contents (g kg-1); WFC, water retention at field capacity (volumetric percentage); PWP, water retention at permanent wilting point (volumetric percentage); and AW, available water (volumetric percentage). .

Soil OC, OCS, and CEC were strongly correlated (Table 3), and decreased with depth; for OCS calculation, the 0.20-0.40-m layer has double the thickness of the first two. The highest values were found in Chernossolos Háplicos (Haplustolls, Argiustolls) and Vertissolos Háplicos, which also had high pH. Both soils are formed from nutrient-richer parent material; however, Chernossolos Háplicos occur on footslopes and Vertissolos Háplicos on depressions. It should be highlighted that water retention properties were stable across depths and were highly correlated with particle-size fractions and CEC. Water retention was greater in Cambissolos Flúvicos (Fluventic Haplustepts), Vertissolos Háplicos, and Cambissolos Háplicos; the two latter with high clay contents. Latossolos Amarelos and Latossolos Vermelho-Amarelos had the lowest OC, OCS, pH, CEC, water retention, clay, and silt content values, since they were formed from poorer and sandier parent materials. These soil property values are consistent with those of a previous study in the same region (Oliveira et al., 1998OLIVEIRA, C.V.; KER, J.C.; FONTES, L.E.F.; CURI, N.; PINHEIRO, J.C. Química e mineralogia de solos derivados de rochas do Grupo Bambuí no norte de Minas Gerais. Revista Brasileira de Ciência do Solo, v.22, p.583-593, 1998. DOI: 10.1590/S0100-06831998000400003.
https://doi.org/10.1590/S0100-0683199800...
), in which relatively high OC, CEC, and pH were also found in Vertissolos (Vertisols).

Table 3
Pearson’s linear correlations among soil properties(1) (1) OC, organic carbon content (g kg-1); BD, bulk density (kg dm-3); OCS, organic carbon stock (kg m-2); CEC, cation exchange capacity (cmolc kg-1); clay, silt, and sand contents (g kg-1); WFC, water retention at field capacity (volumetric percentage); PWP, water retention at permanent wilting point (volumetric percentage); and AW, available water (volumetric percentage). .

The SMLR models achieved R2 varying from 0.23, for OC at 0-0.10 m, to 0.84, for sand content at 0.20-0.40 m (Table 4), with all linear regression assumptions met. Models included from 2 to 14 covariates at 5% probability. Relief (elevation and terrain derivatives) and parent material (gamma radiometric and magnetometric) covariates were selected in 31 of the 33 models (11 properties at three depths). Vegetation covariates were selected in 29 models, with general preference for data from the dry (26 models) over the wet (20 models) season. Among the available covariates, the intensity of the magnetic field (27 models) was by far the most selected variable, followed by diffuse insolation (15 models), ratio of thorium to potassium (13 models), and surface reflectance from RapidEye in the red-edge region (band 4, 690-730 nm) from the wet season (13 models). Red edge is the region where the reflectance of live plants increases from the red (chlorophyll absorption) to the near-infrared region, which has been strongly linked to chlorophyll content, leaf area index, and plant moisture (Filella & Peñuelas, 1994FILELLA, I.; PEÑUELAS, J. The red edge position and shape as indicators of plant chlorophyll content, biomass and hydric status. International Journal of Remote Sensing, v.15, p.1459-1470, 1994. DOI: 10.1080/01431169408954177.
https://doi.org/10.1080/0143116940895417...
). The global trend SMLR models of all soil properties are presented in Table 5.

Table 4
Stepwise multiple linear regression model fit and external validation results, with the best geostatistical method in bold for each property(1) (1) nir, near infrared; RE, RapidEye sensor; K, potassium; Th, thorium; U, uranium; kWh, kilowatt-hour; μR, microroentgen; ppm, parts per million; and nT, nanotesla. .
Table 5
Stepwise multiple linear regression models for soil properties(1) (1) R2, training coefficient of determination; OC, organic carbon content (g kg-1); BD, bulk density (kg dm-3); OCS, organic carbon stock (kg m-2); CEC, cation exchange capacity (cmolc kg-1); clay, silt, and sand contents (g kg-1); WFC, water retention at field capacity (volumetric percentage); PWP, water retention at permanent wilting point (volumetric percentage); and AW, available water (volumetric percentage). Covariates are described in Table 1. .

These findings show the strong parent material control of the soil property patterns in the area, but also the effect of sun heat and vegetation. This makes sense in this region where plant patterns are controlled both by climatic conditions, which regulate evapotranspiration, leaf fall, and phenology (Pezzini et al., 2014PEZZINI, F.F.; RANIERI, B.D.; BRANDÃO, D.O.; FERNANDES, G.W.; QUESADA, M.; ESPÍRITO-SANTO, M.M.; JACOBI, C.M. Changes in tree phenology along natural regeneration in a seasonally dry tropical forest. Plant Biosystems, v.148, p.965-974, 2014. DOI: 10.1080/11263504.2013.877530.
https://doi.org/10.1080/11263504.2013.87...
), and by soil conditions, which regulate water and nutrient availability (Souza et al., 2007SOUZA, J.P. de; ARAÚJO, G.M.; HARIDASAN, M. Influence of soil fertility on the distribution of tree species in a deciduous forest in the Triângulo Mineiro region of Brazil. Plant Ecology, v.191, p.253-263, 2007. DOI: 10.1007/s11258-006-9240-2.
https://doi.org/10.1007/s11258-006-9240-...
). Oliveira-Filho et al. (1998)OLIVEIRA-FILHO, A.T.; CURI, N.; VILELA, E.A.; CARVALHO, D.A. Effects of canopy gaps, topography, and soils on the distribution of woody species in a central Brazilian deciduous dry forest. Biotropica, v.30, p.362-375, 1998. DOI: 10.1111/j.1744-7429.1998.tb00071.x.
https://doi.org/10.1111/j.1744-7429.1998...
identified available light, that is, canopy gaps, as a critical factor controlling the plant species' abundance distribution in a tropical dry forest in Central Brazil.

The spatial patterns of the soil properties were related, either positively or negatively (Table 3 and Figures 2 and 3). In the long range, they responded to the variation of geology, which, to some extent, also controls the variation of vegetation. In the short range, they reflected the spatial distribution of soil classes and their properties. The sampling density of about 3.2 samples per square kilometer was enough to derive better OK than RK maps using many (and some costly) environmental covariates for most cases; however, the effect of sampling density on prediction quality was not tested. Since soil properties were correlated, all OK maps look similar. Despite this, the RK maps show differences related to the covariates used in the global trend models, even though they globally resemble the other maps; in addition, the short-range variation of these maps is more detailed according to the variation of the environmental covariates.

Figure 2
Final maps of soil properties based on their best geostatistical methods: A, organic carbon content (OC) at 0.0-0.10-m depth by regression kriging (RK), and at 0.10-0.20 and 0.20-0.40 m by ordinary kriging (OK); B, bulk density (BD) at 0.0-0.10, 0.10-0.20, and 0.20-0.40 m by OK; C, organic carbon stock (OCS) at 0.0-0.10, 0.10-0.20, and 0.20-0.40 m by OK; D, pH at 0.0-0.10, 0.10-0.20, and 0.20-0.40 m by RK; and E, cation exchange capacity (CEC) at 0.0-0.10 and 0.10-0.20 m by RK, and at 0.10-0.20 m by OK.

Figure 3
Final maps of soil properties based on their best geostatistical methods: A, clay content at 0.0-0.10, 0.10-0.20, and 0.20-0.40 m by ordinary kriging (OK); B, silt content at 0.0-0.10 m by OK, at 0.10-0.20 m by regression kriging (RK), and at 0.10-0.20 m by OK; C, sand content at 0.0-0.10, 0.10-0.20, and 0.20-0.40 m by OK; D, water retention at field capacity (WFC) at 0.0-0.10, 0.10-0.20, and 0.20-0.40 m by RK; E, water retention at permanent wilting point (PWP) at 0.0-0.10, 0.10-0.20, and 0.20-0.40 m by OK; and F, available water (AW) at 0.0-0.10 and 0.10-0.20 m by RK, and at 0.10-0.20 m by OK.

The western and northwestern portions of the study area are on sand-rich parent material. These portions have: the sandiest and poorest soils in terms of fertility, that is, low carbon and CEC; low water holding capacity, that is, low WFC, PWP, and AW; and the poorest vegetation of the park, the "carrasco". The eastern portion constitutes the São Francisco River floodplain, with soils with fluvic characteristics, good water holding capacity, and medium fertility. The central portion corresponds to the area of the "Bambuí" Group, with variable lithology and the presence of calcareous materials, occasionally as rock outcrops. The most common soils in the area are Latossolos Vermelhos and Cambissolos Háplicos, each occupying 31% of the area. Soils with largest carbon stocks, that is, Vertissolos Háplicos and Chernossolos Háplicos, appear in this portion and together occupy less than 3% of the study area.

Total soil OCS at 0.40 m was determined by adding the values obtained at the three depths derived using OK. The study area stores approximately 5.65x105 tons of soil OC in the first 40 cm, of which 34% are found in Cambissolos Háplicos and 30% in Latossolos Vermelhos, due to their large coverage. The "Urucuia" Group, "Bambuí" Group, and "Quaternary sediments" geologic units hold 6, 69, and 24% of the OCS at 0.40 m, respectively, with the largest mean (5.79 kg m-2) in the "Bambuí" Group.

Considering external validation, OK outperformed RK for 21 out of 33 property x depth combinations, though, for most cases, OK and RK had similar prediction quality (Table 4). The preference for OK over RK did not follow a clear pattern. In fact, in some cases, the best method varied among depths for the same soil property. Nonetheless, this preference still indicates that the local variation of soil properties was, in general, a stronger predictor than environmental covariates. This is regardless of the quality of the SMLR models with RK, which, in some cases, had R2 >0.60 but still performed worse than OK, as observed for clay and sand contents, and PWP. As a rule, the quality of soil predictions depends on soil property, number and distribution of samples, environmental covariates, and soil-landscape correlations, among other factors. It should be noted that neither the number of samples nor the pool of environmental covariates (soil-forming factor groups) was tested in the present study.

Semivariogram parameters varied among soil properties (Table 6). Sand content and PWP had the highest ranges, whereas pH and BD had the lowest ones. Across depths, the same semivariogram model was chosen and similar ranges were found for the same property, indicating similar spatial dependence structures. The strength of spatial continuity, as evidenced by the nugget:sill ratio, was as follows: higher for the particle-size fractions and water retention properties; intermediate for the chemical properties, such as pH and CEC; and lower for OC, OCS, and BD. After global trend modeling by SMLR, partial sill and total sill variances decreased in all semivariograms, when compared to the original properties, whereas nugget variance and range changed inconsistently. These results are expected since the SMLR models explained part of the total variance of soil properties, but not necessarily their short- (nugget) or long- (range) distance variations. In any case, changing the semivariogram parameters affected kriging weights and derived spatial patterns, as well as kriging variances.

Table 6
Semivariogram parameters for original soil properties and residuals from stepwise multiple linear regression (SMLR)(1) (1) OC, organic carbon content (g kg-1); BD, bulk density (kg dm-3); OCS, organic carbon stock (kg m-2); CEC, cation exchange capacity (cmolc kg-1); clay, silt, and sand contents (g kg-1); WFC, water retention at field capacity (volumetric percentage); PWP, water retention at permanent wilting point (volumetric percentage); and AW, available water (volumetric percentage). .

Previous studies have also shown little or no improvement of RK over OK, at plot (Kravchenko & Robertson, 2007 KRAVCHENKO, A.N.; ROBERTSON, G.P. Can topographical and yield data substantially improve total soil carbon mapping by regression kriging? Agronomy Journal, v.99, p.12-17, 2007. DOI: 10.2134/agronj2005.0251.
https://doi.org/10.2134/agronj2005.0251...
), farm/catchment (Zhu & Lin, 2010ZHU, Q.; LIN, H.S. Comparing ordinary kriging and regression kriging for soil properties in contrasting landscapes. Pedosphere, v.20, p.594-606, 2010. DOI: 10.1016/S1002-0160(10)60049-5.
https://doi.org/10.1016/S1002-0160(10)60...
), or watershed scale (Vasques et al., 2010 VASQUES, G.M.; GRUNWALD, S.; COMERFORD, N.B.; SICKMAN, J.O. Regional modelling of soil carbon at multiple depths within a subtropical watershed. Geoderma, v.156, p.326-336, 2010. DOI: 10.1016/j.geoderma.2010.03.002.
https://doi.org/10.1016/j.geoderma.2010....
). Kravchenko & Robertson (2007) KRAVCHENKO, A.N.; ROBERTSON, G.P. Can topographical and yield data substantially improve total soil carbon mapping by regression kriging? Agronomy Journal, v.99, p.12-17, 2007. DOI: 10.2134/agronj2005.0251.
https://doi.org/10.2134/agronj2005.0251...
found that RK produced only a modest improvement in accuracy compared to OK and performed poorly in data sets with strong spatial correlation in the target variable, even when the regression model was relatively strong. The same behavior was observed in the present study. In two contrasting landscapes (agricultural versus forested) in the United States, RK was superior to OK only when the spatial structure could not be well captured by point-based observations or when a strong relationship existed between the target soil property and the covariates (Zhu & Lin, 2010ZHU, Q.; LIN, H.S. Comparing ordinary kriging and regression kriging for soil properties in contrasting landscapes. Pedosphere, v.20, p.594-606, 2010. DOI: 10.1016/S1002-0160(10)60049-5.
https://doi.org/10.1016/S1002-0160(10)60...
). In the same country, in Florida, Vasques et al. (2010) VASQUES, G.M.; GRUNWALD, S.; COMERFORD, N.B.; SICKMAN, J.O. Regional modelling of soil carbon at multiple depths within a subtropical watershed. Geoderma, v.156, p.326-336, 2010. DOI: 10.1016/j.geoderma.2010.03.002.
https://doi.org/10.1016/j.geoderma.2010....
observed that the preference for RK over OK depended on the depth of the measurement and on the regression method used to estimate soil total carbon in a 3,585-km2 watershed. Studies in tropical dry forest with similar extents or sampling densities were not found to compare with the results obtained in the present study.

Conclusions

  1. Soil property patterns are linked to each other, as well as to topographic, geologic, and vegetation patterns in a tropical dry forest area in Brazil.

  2. Ordinary kriging is a relatively simple and efficient method to create soil property maps, but it is conditioned to a good sampling design, since, contrary to regression kriging, it only relies on soil samples and not on environmental covariates.

  3. There is no straightforward decision on whether to use ordinary kriging or regression kriging for mapping soil properties, because the quality of predictions depends on soil property, number and distribution of samples, environmental covariates, and soil-landscape correlations, among other factors, which indicates that, for similar studies, the employed geostatistical methods should be compared on a case-by-case basis.

Acknowledgements

To Dr. Mário Marcos do Espírito Santo, Dr. Ricardo Luís Louro Berbara, Dr. Maria de Lourdes Mendonça Santos Brefin, to Instituto Estadual de Florestas (IEF) of the state of Minas Gerais, Brazil, to the staff of Parque Estadual da Mata Seca, and to Serviço Geológico do Brasil (CPRM), for support; and to Embrapa (grant No. 03.10.06.013.00.00) and Inter-American Institute for Global Change Research-United States National Science Foundation (grant No. GEO-04523250), for financial support.

References

  • BOUCNEAU, G.; VAN MEIRVENNE, M.; THAS, O.; HOFMAN, G. Integrating properties of soil map delineations into ordinary kriging. European Journal of Soil Science, v.49, p.213-229, 1998. DOI: 10.1046/j.1365-2389.1998.00157.x.
    » https://doi.org/10.1046/j.1365-2389.1998.00157.x
  • COELHO, M.R.; DART, R. de O.; VASQUES, G. de M.; TEIXEIRA, W.G.; OLIVEIRA, R.P. de; BREFIN, M. de L.M.; BERBARA, R.L.L. Levantamento pedológico semi-detalhado (1:30.000) do Parque Estadual da Mata Seca, município de Manga-MG Rio de Janeiro: Embrapa Solos, 2013. 264 p. (Embrapa Solos. Boletim de pesquisa e desenvolvimento, 217).
  • DONAGEMA, G.K.; CAMPOS, D.V.B. de; CALDERANO, S.B.; TEIXEIRA, W.G.; VIANA, J.H.M. (Org.). Manual de métodos de análise de solos 2.ed. rev. Rio de Janeiro: Embrapa Solos, 2011. 230p. (Embrapa Solos. Documentos, 132).
  • FILELLA, I.; PEÑUELAS, J. The red edge position and shape as indicators of plant chlorophyll content, biomass and hydric status. International Journal of Remote Sensing, v.15, p.1459-1470, 1994. DOI: 10.1080/01431169408954177.
    » https://doi.org/10.1080/01431169408954177
  • GRUNWALD, S. Multi-criteria characterization of recent digital soil mapping and modeling approaches. Geoderma, v.152, p.195-207, 2009. DOI: 10.1016/j.geoderma.2009.06.003.
    » https://doi.org/10.1016/j.geoderma.2009.06.003
  • GRUNWALD, S.; VASQUES, G.M.; RIVERO, R.G. Fusion of soil and remote sensing data to model soil properties. Advances in Agronomy, v.131, p.1-109, 2015. DOI: 10.1016/bs.agron.2014.12.004.
    » https://doi.org/10.1016/bs.agron.2014.12.004
  • HENGL, T.; JESUS, J.M. de; MACMILLAN, R.A.; BATJES, N.H.; HEUVELINK, G.B.M.; RIBEIRO, E.; SAMUEL-ROSA, A.; KEMPEN, B.; LEENAARS, J.G.B.; WALSH, M.G.; GONZALEZ, M.R. SoilGrids1km - global soil information based on automated mapping. PLoS ONE, v.9, p.e105992, 2014. DOI: 10.1371/journal.pone.0105992.
    » https://doi.org/10.1371/journal.pone.0105992
  • KNOTTERS, M.; BRUS, D.J.; OUDE VOSHAAR, J.H. A comparison of kriging, co-kriging and kriging combined with regression for spatial interpolation of horizon depth with censored observations. Geoderma, v.67, p.227-246, 1995. DOI: 10.1016/0016-7061(95)00011-C.
    » https://doi.org/10.1016/0016-7061(95)00011-C
  • KRAVCHENKO, A.N.; ROBERTSON, G.P. Can topographical and yield data substantially improve total soil carbon mapping by regression kriging? Agronomy Journal, v.99, p.12-17, 2007. DOI: 10.2134/agronj2005.0251.
    » https://doi.org/10.2134/agronj2005.0251
  • LEVANTAMENTO aerogeofísico do Estado de Minas Gerais: área 11A (Jaíba - Montes Claros - Bocaiúva). Rio de Janeiro: Serviço Geológico do Brasil, 2009.
  • MILES, L.; NEWTON, A.C.; DEFRIES, R.S.; RAVILIOUS, C.; MAY, I.; BLYTH, S.; KAPOS, V.; GORDON, J.E. A global overview of the conservation status of tropical dry forests. Journal of Biogeography, v.33, p.491-505, 2006. DOI: 10.1111/j.1365-2699.2005.01424.x.
    » https://doi.org/10.1111/j.1365-2699.2005.01424.x
  • MINASNY, B.; MCBRATNEY, A.B. A conditioned Latin hypercube method for sampling in the presence of ancillary information. Computers and Geosciences, v.32, p.1378-1388, 2006. DOI: 10.1016/j.cageo.2005.12.009.
    » https://doi.org/10.1016/j.cageo.2005.12.009
  • MURPHY, P.G.; LUGO, A.E. Ecology of tropical dry forest. Annual Review of Ecology and Systematics, v.17, p.67-88, 1986. DOI: 10.1146/annurev.es.17.110186.000435.
    » https://doi.org/10.1146/annurev.es.17.110186.000435
  • OLIVEIRA, C.V.; KER, J.C.; FONTES, L.E.F.; CURI, N.; PINHEIRO, J.C. Química e mineralogia de solos derivados de rochas do Grupo Bambuí no norte de Minas Gerais. Revista Brasileira de Ciência do Solo, v.22, p.583-593, 1998. DOI: 10.1590/S0100-06831998000400003.
    » https://doi.org/10.1590/S0100-06831998000400003
  • OLIVEIRA-FILHO, A.T.; CURI, N.; VILELA, E.A.; CARVALHO, D.A. Effects of canopy gaps, topography, and soils on the distribution of woody species in a central Brazilian deciduous dry forest. Biotropica, v.30, p.362-375, 1998. DOI: 10.1111/j.1744-7429.1998.tb00071.x.
    » https://doi.org/10.1111/j.1744-7429.1998.tb00071.x
  • PEZZINI, F.F.; RANIERI, B.D.; BRANDÃO, D.O.; FERNANDES, G.W.; QUESADA, M.; ESPÍRITO-SANTO, M.M.; JACOBI, C.M. Changes in tree phenology along natural regeneration in a seasonally dry tropical forest. Plant Biosystems, v.148, p.965-974, 2014. DOI: 10.1080/11263504.2013.877530.
    » https://doi.org/10.1080/11263504.2013.877530
  • R CORE TEAM. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing, 2015.
  • SOIL SURVEY STAFF. Keys to soil taxonomy 12nd ed. Washington: United States Department of Agriculture, 2014.
  • SOUZA, J.P. de; ARAÚJO, G.M.; HARIDASAN, M. Influence of soil fertility on the distribution of tree species in a deciduous forest in the Triângulo Mineiro region of Brazil. Plant Ecology, v.191, p.253-263, 2007. DOI: 10.1007/s11258-006-9240-2.
    » https://doi.org/10.1007/s11258-006-9240-2
  • VASQUES, G.M.; GRUNWALD, S.; COMERFORD, N.B.; SICKMAN, J.O. Regional modelling of soil carbon at multiple depths within a subtropical watershed. Geoderma, v.156, p.326-336, 2010. DOI: 10.1016/j.geoderma.2010.03.002.
    » https://doi.org/10.1016/j.geoderma.2010.03.002
  • ZHU, Q.; LIN, H.S. Comparing ordinary kriging and regression kriging for soil properties in contrasting landscapes. Pedosphere, v.20, p.594-606, 2010. DOI: 10.1016/S1002-0160(10)60049-5.
    » https://doi.org/10.1016/S1002-0160(10)60049-5

Publication Dates

  • Publication in this collection
    Sept 2016

History

  • Received
    27 Aug 2015
  • Accepted
    27 Jan 2016
Embrapa Secretaria de Pesquisa e Desenvolvimento; Pesquisa Agropecuária Brasileira Caixa Postal 040315, 70770-901 Brasília DF Brazil, Tel. +55 61 3448-1813, Fax +55 61 3340-5483 - Brasília - DF - Brazil
E-mail: pab@embrapa.br