Science of the Total Environment 526 (2015) 58–69

Contents lists available at ScienceDirect

Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv

Spatial variability of soil carbon stock in the Urucu river basin, Central Amazon-Brazil Marcos Bacis Ceddia a,⁎, André Luis Oliveira Villela b, Érika Flávia Machado Pinheiro a, Ole Wendroth c a b c

Department of Soil, Institute of Agronomy, Universidade Federal Rural do Rio de Janeiro (UFRRJ), Seropédica, RJ 23890-000, Brazil Colégio Técnico da UFRRJ, RJ, Seropédica 23890-000, Brazil Department of Plant & Soil Sciences, University of Kentucky, College of Agriculture, Lexington, KY, USA

H I G H L I G H T S • • • •

The SOC stocks across 30 and 100 cm depth were 3.28 and 7.32 kg C m− 2, respectively. SOC stocks were 34 and 16%, respectively, lower than those observed in other studies. Regions with waterlogging soils presented the lowest SOC stock. Heterotopic cokriging, using CTI, improved the accuracy of SOC stock maps.

a r t i c l e

i n f o

Article history: Received 16 December 2014 Received in revised form 27 March 2015 Accepted 28 March 2015 Available online xxxx Editor: Charlotte Poschenrieder Keywords: Içá Formation Isotopic cokriging Heterotopic cokriging CTI

a b s t r a c t The Amazon Forest plays a major role in C sequestration and release. However, few regional estimates of soil organic carbon (SOC) stock in this ecoregion exist. One of the barriers to improve SOC estimates is the lack of recent soil data at high spatial resolution, which hampers the application of new methods for mapping SOC stock. The aims of this work were: (i) to quantify SOC stock under undisturbed vegetation for the 0–30 and the 0–100 cm under Amazon Forest; (ii) to correlate the SOC stock with soil mapping units and relief attributes and (iii) to evaluate three geostatistical techniques to generate maps of SOC stock (ordinary, isotopic and heterotopic cokriging). The study site is located in the Central region of Amazon State, Brazil. The soil survey covered the study site that has an area of 80 km2 and resulted in a 1:10,000 soil map. It consisted of 315 field observations (96 complete soil profiles and 219 boreholes). SOC stock was calculated by summing C stocks by horizon, determined as a product of BD, SOC and the horizon thickness. For each one of the 315 soil observations, relief attributes were derived from a topographic map to understand SOC dynamics. The SOC stocks across 30 and 100 cm soil depth were 3.28 and 7.32 kg C m−2, respectively, which is, 34 and 16%, lower than other studies. The SOC stock is higher in soils developed in relief forms exhibiting well-drained soils, which are covered by Upland Dense Tropical Rainforest. Only SOC stock in the upper 100 cm exhibited spatial dependence allowing the generation of spatial variability maps based on spatial (co)-regionalization. The CTI was inversely correlated with SOC stock and was the only auxiliary variable feasible to be used in cokriging interpolation. The heterotopic cokriging presented the best performance for mapping SOC stock. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Models of carbon cycle require accurate estimations of the carbon storage in different reservoirs. Regarding the soil compartment, global C pools are difficult to estimate because of the limited knowledge about specific properties of soil types (Sombroek et al., 1993; Batjes, 1996), the considerable spatial variability of soil C even within the same soil mapping unit (Cerri et al., 2000), and the different effects of factors controlling soil organic C dynamics (Pastor and Post, 1986; ⁎ Corresponding author. E-mail address: [email protected] (M.B. Ceddia).

http://dx.doi.org/10.1016/j.scitotenv.2015.03.121 0048-9697/© 2015 Elsevier B.V. All rights reserved.

Parton et al., 1987). Thus, regional studies are necessary to refine large-scale estimations obtained by aggregation of regional estimates (Bernoux et al., 2002). Estimates of SOC stocks for the Brazilian Amazon are presented in Table 1. Cerri et al. (2000), using pedotransfer functions, estimated 41 Pg C stocked in the top 1 m in the Brazilian Amazon, which represents about half of all soil carbon stocked in Brazil, emphasizing the importance of the Amazon Forest in carbon storage. The potential total soil C stock of the Brazilian Amazon under native vegetation, estimated by Bernoux et al. (2003) and Batjes (2005) was 22.7 and 23.9–24.2 Pg C in the 0–30 cm layer, respectively. For the Brazilian Amazon (about 500 Mha), Batjes and Dijkshoorn (1999) observed that organic C

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

59

Table 1 An overview of the papers published about SOC stock in the Amazon region (5,217,423 km2). Source

De Moraes et al. (1995)

SOC (Pg)

N

0–30 cm

0–100 cm

21

47

Batjes and Dijkshoorn (1999) 25

46.5

Cerri et al. (2000)

23.4

41

Bernoux et al. (2002)

22.7

–

Batjes (2005) Fidalgo et al. (2007)

23.9–24.2 42.3–43.8 23,9–24,2 –

This study

–

–

Methodology

1162 1:5 M soil map linked to national profile data — mean SOC content per soil unit. Missing bulk density (BD) data replaced with the mean measured BD of the corresponding soil group. Source of dates: EMBRAPA (1981); RADAMBRASIL (1973–1983). 1828 1:5 M SOTER for Latin America — mean SOC content per FAO soil unit. Estimates of missing BD data. Source of dates: SOTER-LAC 1:5 M soil map linked to national soil profile dataset — medians for SOC content per soil type. Missing BD data from linear regression. 3969 Soil-Vegetation Associations (SVA), using national soil data (1:5 M); Calculation of representative C stock values for each SVA category producing a map of their distribution, and extraction of part of the Brazilian Amazon portion from the stock produced for the whole of Brazil. Source of dates: (EMBRAPA, 1981; RADAMBRASIL, 1973–1983; IBGE, 1988). 584 Updated 1:5 M SOTER database for Brazil (Batjes et al., 2004) and simulation of phenoforms (Batjes, 2005). 363 IPCC/UNEP/OECD/OEA (IPCC, 1997). Soil-vegetation association (1:5 M). Source of dates: IBGE (2004). 96 Data were obtained from a detailed soil survey (1:10,000). SOC stock for each soil profile was calculated by summing C stocks by horizon, determined as a product of BD, SOC and the horizon thickness. In each one of the 315 soil observations, relief attributes were derived from a topographic map. Maps of SOC stock were derivate using geostatistical techniques.

stock was 25 Pg C (0–30 cm) and 46.5 Pg C (0–100 cm), similar to the values reported by de Moraes et al. (1995). Despite the quality and the importance of these SOC stock estimations for the amazon region, it is important to highlight some common aspects of these results presented in Table 1, such as: a) directly or indirectly, the data source used by different authors has the same origin (RADAMBRASIL), which was generated a long time ago (from 1973–1984); b) the maps used by different authors were published at a coarse scale and the data available for SOC calculations were very sparse along the Amazon region (~ 0.7–3.5 soil profiles for each 10,000 km2 — Table 1); c) the Legal Amazon region is a huge territory covering an area of 5,217,423 km2 (59% of the total Brazilian territory), encompassing different ecosystems. The latter implies that a model based on low-density data information, probably fails to predict SOC stock for the entire Amazon Region; d) the estimates of SOC stocks have been made through soil mapping or soil-vegetation units linked to representative soil profile data; e) in general, authors do not clearly reveal the mechanistic relation between the SOC stock, relief and vegetation. Some decisions during the mapping process are not completely explained, and consequently, maps become hardly reproducible. The aspects presented above demonstrate the necessity to generate new data at a higher spatial resolution, as well as to apply appropriate methods for mapping SOC stock. A broad range of statistical methods has been applied towards digital soil mapping of SOC stock. Most commonly, multiple and linear regressions have been used for spatial quantifications of carbon stock (Moore et al., 1993; Powers and Schlesinger, 2002; Thompson et al., 2006). Based on the variability structure, geostatistics can be used to estimate the spatial distribution of carbon stock through kriging with or without external drift or cokriging (Simbahan et al., 2006). Some authors have applied regression kriging to predict SOC stock or soil organic matter (Bhatti et al., 1991; Hengl et al., 2004). Considering the demand for reliable and reproducible maps of SOC stock based on data that were collected at a high spatial resolution in the Amazon region, we conducted a study in an area located in the Central part of the Amazon State. The aims of this work were: (i) to quantify SOC stock under undisturbed vegetation for the 0–30 cm and 0–100 cm reference layers in the Amazon Forest; (ii) to correlate the SOC stock with soil mapping units and relief attributes; (iii) to evaluate three geostatistical techniques to generate spatial variability maps of SOC stock (ordinary kriging, isotopic cokriging and heterotopic cokriging). We hypothesized that relief is the dominating process responsible for soil and vegetation variability along the study site. As a soil attribute, SOC stock is a consequence of these relationship and consequently

some relief attributes can be used to adequately estimate spatial variability of SOC stock using cokriging. 2. Methods 2.1. Study area The study site is located in the central region of the Amazon State, nearby the Urucu River, in the municipality of Coari, Brazil (Fig. 1). The study site belongs to Içá Formation, and encompasses an area of 80 km2 in the Amazon Forest. The region, about 640 km from Manaus (State capital), can be accessed only by boat or airplane. The climate is equatorial (Af — Köppen climate classification) where the temperature of the coldest month is higher than 20 °C, with no pronounced dry period and a mean annual precipitation of 2500 mm. According to RADAMBRASIL (1978), the soils of the study site were formed from sediments of Içá Formation. The sediments of Içá Formation cover an area of 563.264 km2 (36% of the Amazon State) and were deposited in the Tertiary–Quaternary Period. The Içá Formation consists of fine to medium sandstone and siltstone, locally with clay conglomerates and yellow-red colors (Galvão et al., 2012). The Holocene alluvium of the Quaternary Period deposits is related to the current Amazonian drainage networks. 2.2. Soil database Between the years of 2008 and 2009, a soil survey was conducted in the Oil Province of the Urucu River, called Geologo Pedro de Moura. The work resulted in the generation of a soil map, which covers an area of 80 km2 (Fig. 1), along with its respective report (Villela, 2013). Throughout the soil survey, 315 field observations were performed, consisting of 96 soil profiles and 219 soil boreholes (Fig. 2). Due to the limitations imposed by native vegetation, the 315 field observations were restricted to the vicinity of access roads. The soils were classified as suggested by the Brazilian Soil Classification System (Embrapa, 1999). The soil-mapping units of the study site, as well as the number of profiles and area of occurrence are shown in Table 2. All the data used for SOC stock calculation were obtained from the report of the soil survey. In each horizon of the 96 soil profiles, the following soil attributes were determined: soil organic carbon (SOC), soil bulk density (BD), sand, silt and clay content. SOC was measured by wet combustion, according to the methodology proposed by Walkley and Black (1932). BD was measured using undisturbed soil core samples (Kopeck rings: 4.2 cm height and 4.0 cm diameter). Soil particle

60

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

Fig. 1. The location of the study area in the Central Amazon, Brazil (upper). The soil map of the study site is shown in the bottom of the figure. The dashed line shows where the soil map is located in the Central Amazon Region.

size distribution (sand, silt and clay content) was measured according to the Pipette method (Embrapa, 1997). 2.3. Soil organic carbon stock calculation For each of the 96 soil profiles the calculation of the SOC stock (also called soil organic carbon storage, soil organic carbon pool

or soil organic carbon density) was performed in the layers of 0–30 and 0–100 cm depth. Commonly, these two depths are the most used in soil SOC studies (Batjes, 2000), although the IPCC (1997) methodology recommends SOC stock only at 0–100 cm depth. The classical way of calculating C densities (C mass per area) for a given depth consists of summing C stocks by horizon, determined as a

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

61

Fig. 2. Spatial distribution of soil profiles and boreholes across the study site (315 observations).

product of BD, SOC concentration, and horizon thickness, according to Bernoux et al. (2002): SOC stock ¼ ðSOC  BD  TÞ where: SOC stock soil organic carbon stock (kg C m−2); SOC soil organic carbon (g kg−1); BD soil bulk density (Mg m−3); T horizon thickness (m). In the soil survey, the soil profiles were divided into A-, B-, and Chorizons. In most cases, the calculations concerned two horizons where the first horizon typically was entirely above 30 cm, and the second one crossed this 30 cm or 100 cm depth. When a horizon crossed the 30 or 100 cm boundary, simply the portion of the horizon that was above that depth was used to calculate its SOC stock. 2.4. Topography attributes Topography has the potential to explain large parts of the variation of the SOC stock. Thus, models accounting for terrain attributes can provide more realistic estimates of SOC stock. Topographic maps were based on Radar image (X Band, 0.5 m resolution, HH polarization) throughout the area and its surroundings at the 1:5000 scale resulting

in 2-m contour line intervals and a hydrographic network. A total of 4 terrain attributes (slope, aspect, curvature and compound topographic index) were calculated in order to assess the relationship between SOC stock and relief. The topographic map in “dxf” format was exported to ArcGIS software version 9.3 (ESRI, 2006) where the map was converted to a raster at a spatial resolution of 5 m. Then, the digital elevation model (DEM) was generated. Using the surface module of the spatial analyst tool, slope, aspect and curvature maps were generated from the DEM. The compound topographic index (CTI), often referred to as the steady-state wetness index, is a quantification of catenary landscape position (Gessler et al., 1995). It is defined as:
CTI ¼ ln

 As tanβ

where: As is the specific catchment area (m2) per unit width orthogonal to the flow direction [(flow accumulation + 1) ∗ (pixel area in m2)], and β is the slope angle (radians). The relief attributes generated in ArcGis (slope, curvature, aspect and CTI) were extracted for each soil profile (96 points — soil attributes and relief attributes with the same locations) and Borehole (219 points — only soil morphology information and relief attributes). Cross-semivariograms were computed for the extracted values of relief attributes and SOC stock in order to evaluate their feasibility as auxiliary variables to improve the accuracy of SOC stock estimation.

Table 2 Mean soil organic carbon stock (SOC) in each soil mapping unit at 0–30 cm and 0–100 cm (kg C m−2) and percentage of soil organic carbon stocks (% SOC) in the upper layer (0–30 cm). MU

MU1 MU2 MU3 MU4 MU5 MU6 MU7 Total

SOC (30 cm)

SOC (100 cm)

SOC

N

Area (ha)

%

Mean

Min.

Max.

CV

Mean

Min.

Max.

CV

%

10 21 33 5 15 7 5 96

1631.7 1703.9 1491.7 297.9 1073.1 1394.9 373.3 7966.5

20.5 21.4 18.7 3.7 13.5 17.5 4.7 100

3.08 3.63 3.54 3.97 3.32 2.82 2.53 3.28⁎

1.90 2.07 1.54 2.97 1.84 1.74 1.67 –

6.42 5.74 6.44 5.41 5.15 3.72 5.48 –

34 31 30 29 31 27 51 –

6.03 7.82 7.40 8.73 7.90 8.02 5.01 7.32⁎

4.26 5.16 3.74 7.81 3.64 5.05 3.26 –

8.02 10.40 11.93 11.26 11.87 13.05 8.45 –

16 21 25 17 31 34 42 –

51 46 47 45 42 35 52 45⁎

Min. — minimum value; Max. — maximum value; CV — coefficient of variation. N — number of observations; % SOC — mean percentage of soil organic carbon stocks at 0–30 cm depth. MU — mapping unit; MU1 — consociation CXbd; MU2 — complex PVAal–CXal–CXbd; MU3 — Consociation CXal; MU4 — complex CXal–PVAa; MU5 — complex CXal–Pad; MU6 — complex PAal–Pad; MU7 — PACd. PVAal/PVAa: typic hapludults; CXbd: fluvaquentic distrudepts; CXal: typic dystrudepts; PAal — Pad: oxyaquic hapludults; MU7 — PACd: typic endoaquults. ⁎ is the average value.

62

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

Fig. 3. Spatial distribution of the training and validation dataset used for geostatistical analysis.

2.5. The Soil-Relief-Vegetation (SRV) landscape unities along the study site The interpretation of the SOC stock distribution along the study site demanded the creation of SRV (landscape unity). These unities are important, once the SOC stock is a consequence of the combined effect of soil formation factors. Furthermore, organizing the carbon stock data in the form of landscape unit, it is possible to obtain a more integrated vision and evaluate the existence of regional patterns to assist in the construction of prediction models. The SRV unities were developed based on the co-occurrence of soil map unities (MU), relief forms, vegetation classes and geology. The

Soils-MU

Vegetation class* Fac Fda Fda Fdb

Geology* I.F.- H.S. I.F.- P.S. I.F.- P.S. I.F.- P.S.

Relief Forms* APf C11 T21 EP2

information of relief forms, vegetation classes and geology were obtained from RADAMBRASIL project (1978). The SRVs created in this work followed the criteria below: SRV1 — The co-occurrence of soil map unity MU1-(CXbd) + relief form APf — River plains + vegetation class Fac — Flooded Lowland Open Tropical Rainforest. This landscape unity has 10 soil profiles with SOC stock measured. SRV2 — The co-occurrence of soil map unity MU1 + MU2 + MU3 + MU4-(PVAal–CXal–CXbd, CXal, CXal–PVAa, respectively) + relief form C11 — Dried out areas on flat-topped + vegetation class

Carbon Stock Lower CS and concentrated in the surface layers. Carbon source: leaves and tree trunks, and surface roots. Higher CS throughout the soil profile. Carbon source: leaves and tree trunks, and roots distributed from the surface until deeper layers. Higher CS throughout the soil profile. Carbon source: leaves and tree trunks, and roots distributed from the surface until deeper layers. Lower CS and concentrated in the surface layers. Carbon source: leaves and trunks of treetops and surface roots. N SOC 100 SOC30 %SOC Clay A Clay B CTI Slope

Soil-Relief-Vegetation Unity- SRV bc a a c c 10 SRV1 (MU1+Fac+HS+Apf) 51 179 6.03 3.08 163 a a a a 59 7.93ª 3.60 46 355 372 SRV2 (MU2+MU3+MU4+Fda+PS+C11) ab a a b b 22 299 7.60 3.39 47 278 SRV3 (MU5+MU6+Fda+PS+T21) bc bc c a a 5 5.01 2.53 52 196 184 SRV4 (MU7+Fdb+PS+EP2) Average values with the same letter, in each column, are statistically equals (Tukey Test- P< 0.05). N- Number of soil profiles per SRV unity.

a

9.9 c 7.6 c 8.1b ab 9.0

b.

2.5 a 11.1 ab 4.6 b 2.6

Fig. 4. A schematic representation of the relationship between soil carbon stocks, soil mapping unity, vegetation classes, geology and relief forms. I.F. — Içá Formation; H.C. — Holocene Sediments; P.S. — Pleistocene Sediments; Fac — Flooded Lowland Open Tropical Rainforest; Fda — Upland Dense Tropical Rainforest; Fdb — Upland Open Tropical Rainforest; APf — River plains; C11 — Dried out areas on flat-topped; T21 — Tabular Interfluves; EP2 — Bi-plained superficies-flatlands. CS — Carbon Stocks. * Source: RADAMBRASIL (1978).

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

63

Table 3 Descriptive statistics of SOC stocks and relief attributes of the 96 soil profiles and 315 data, respectively. Attributes

Na

Mean

Median

Min.b

Max.c

Kurtosis

Skwewness

CV (%)

SOC 30 cm (kg C m−2)⁎ SOC 100 cm (kg C m−2) Sand A (g kg−1)d Silt A (g kg−1)e Clay A (g kg−1)f Sand + silt A (g kg−1)g Clay B (g kg−1)h Slope (%)⁎

96 96 96 96 96 96 96 96 96 96 96 315 315 315 315

3.46 7.49 501.7 347.8 150.5 849.5 314.4 7.51 192.7 0.07 9.00 7.46 174.8 0.07 9.11

3.23 7.55 490 333.5 142.5 857.5 312 7.4 209.4 0.09 8.53 6.8 168.6 0.09 8.70

1.54 3.26 172.0 26.0 21.0 678.0 13.0 0.40 0.77 −1.48 6.37 0.2 0.77 −1.41 6.32

6.44 13.05 918.0 687.0 322.0 979.0 563.0 22.10 359.60 0.86 15.02 28.9 359.7 1.0 18.72

−0.10 −0.19 −0.73 −0.70 −0.57 −0.57 0.13 −0.06 −1.38 1.12 1.66 0.83 −1.30 1.47 5.49

0.62 0.14 0.34 0.11 0.38 −0.38 −0.33 0.75 −0.19 1.21 1.37 0.95 0.10 −0.52 1.99

32.5 27.9 35.3 44.2 48.8 8.6 34.4 72.3 57.9 469.4 18.9 71.6 61.2 514 19.7

Aspect Curvature⁎ CTIi,⁎ Slope (%)⁎ Aspect Curvature⁎ CTIi,⁎ a

N — number of data. Min. — minimum value. c Max. — maximum value. d Sand A — sand content in A horizon. e Silt A — silt content in A horizon. f Clay A — clay content in A horizon. g Sand + Silt A — sand silt content in A horizon. h Clay B — clay content in B horizon, except BA. i CTI — compound topographic index. ⁎ Dates were not normally distributed (p b 0.05, Kolmogorov–Smirnov test). b

Fda — Upland Dense Tropical Rainforest. This landscape unity has 59 soil profiles with SOC stock measured. SRV3 — The co-occurrence of soil map unity MU5 + MU6-(CXal + PAd, PAal + PAd) + relief form T21 — Tabular Interfluves + vegetation class Fda — Upland Dense Tropical Rainforest. This landscape unity has 22 soil profiles with SOC stock measured. SRV4 — The co-occurrence of soil map unity MU7-(PACd) + relief form EP2 — Bi-plained superficies-flatlands + and vegetation class Fdb — Upland Open Tropical Rainforest. This landscape unity has 5 soil profiles with SOC stock measured. For each SRV unity, the average weight value of SOC stock up to 100 and 30 cm, texture fraction (Clay, Silt and Sand), and relief attributes (slope, aspect, curvature and CTI) were calculated. The average values of these attributes were weighted according to the number of soil profile belonging to each SRV. The statistical analysis comparing the average weight values of these attributes are shown in Fig. 4 (bottom).

2.6. Geostatistical estimation methods In order to generate spatial variability maps of SOC stock (target variable or primary variable) and take the advantage of the available topographic attributes, three geostatistical techniques were compared: Table 4 SOC stock estimation for the study site in relation to the literature (kg C m−2). Authors

SOC 30 cm

SOC 100 cm

Ceddia (this paper) De Moraes et al. (1995)a Batjes and Dijkshoorn (1999)a Cerri et al. (2000)a Bernoux et al. (2002)a Batjes (2005)a Fidalgo et al. (2007)a Average deviation (%)

3.28 4.02 (+23%) 4.79 (+46%) 4.48 (37%) 4.35 (+33%) 4.61 (+41%) 4.01 (+22%) 34

7.32 9.01 (+23%) 8.91 (+22%) 7.86 (+7%) – 8.25 (+13%) – 16

a SOC was converted from Pg C to kg C m−2, considering the Legal Amazonia area of 5,217,423 km2. Values in parentheses mean the percentage of deviation in relation to this study.

ordinary kriging (OK), isotopic cokriging (ICOK) and partially heterotopic cokriging (HCOK). OK is the most widely used kriging method. It serves to estimate a value at a point of a region for which a variogram is known, using a data in the neighborhood of the estimation location (Wackernagel, 2003). For the study site, SOC stock is the primary variable Z1, measured at sampled locations u to estimate SOC stocks at unsampled locations (Z⁎ok(u)). 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 is written as a linear combination of the n(u) data Z(uα) with a single unbiasedness constraint: 

ZOK ðuÞ ¼

nðuÞ X

λa ðuÞ½Z1 ðuÞ; with

a¼1

nðuÞ X

OK

λa ¼ 1:

a¼1

The unknown local mean m(u) is filtered from the linear estimator by forcing the kriging weights to sum to 1. Cokriging (COK) is a multivariate extension of kriging in which the auxiliary information (in this case, relief attributes) is incorporated in

Table 5 Soil texture and relief attributes related to soil mapping units. MU

Sanda

Silta

Claya

Sand + silta (g kg−1)

Clayb

Slope (%)

Asp

Curv

CTI

MU1 MU2 MU3 MU4 MU5 MU6 MU7

593 552 504 619 438 349 403

320 251 352 216 395 522 491

88 196 145 165 168 129 106

913 804 855 835 832 871 894

156 351 326 416 358 313 191

1.6 10.2 9.5 11.3 5.6 3.5 2.3

230.74 166.84 186.92 281.51 168.06 198.96 241.74

0.12 0.31 0.25 0.36 0.43 0.12 0.05

11.3 8.4 8.6 8.0 8.7 9.1 11.1

MU — mapping unit; MU1 — Consociation CXbd; MU2 — complex PVAal–CXal–CXbd; MU3 — Consociation CXal; MU4 — complex CXal–PVAa; MU5 — complex CXal–Pad; MU6 — complex PAal–Pad; MU7 — PACd. PVAal/PVAa: typic hapludults; CXbd: fluvaquentic distrudepts; CXal: typic dystrudepts; PAal — Pad: oxyaquic hapludults; MU7 — PACd: typic endoaquults. Asp. — aspect; Curv. — curvature; CTI — compound topographic index. a Sand, silt, clay and sand + silt in the A-horizon. b Clay in the B-horizon.

64

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

B

A

C

Fig. 5. Vegetation type in the study area. A) Fda — Upland Dense Tropical Rainforest, B) Fac — Flooded Lowland Open Tropical Rainforest and C) Fdb — Upland Open Tropical Rainforest.

the estimation at unsampled locations by accounting for spatial correlations of the primary variable Z1 with secondary attributes Zv (Goovaerts, 1998). Experimental auto- and cross-semivariograms of all target and auxiliary variables were computed and modeled in the form of a linear model of coregionalization (Goulard and Voltz, 1992). The COK estima⁎ (u)), for an unsampled location u is then: tor (Zcok 

ZCOK ðuÞ ¼

nðuÞ X

λa ðuÞ½Z1 ðuÞ þ

Nv X

λv ðuÞ½Zv ðuÞ;

v¼1

a¼1

with: nðuÞ X

2.7. Evaluation of geostatistical performance

λa ¼ 1

a¼1

and, nv X

When all the variables (target and secondary variables) have been measured at the same sampling locations it is called isotopic cokriging (from the Greek: iso = same and topos = location). Otherwise, when some sample locations are not shared by all variables (the target variable is under-sampled), which is called partially heterotopic cokriging. With partially heterotopic data it is advisable, whenever possible, to infer the cross-semivariogram or the covariance function model on the basis of the isotopic subset of the data (Wackernagel, 2003).

λv ¼ 0

v¼1

Cokriging does not require that the secondary information is available at all locations u to be estimated. The influence of the secondary information on estimating Z depends on (i) the correlation between primary and ancillary variables, (ii) the spatial continuity of the attributes, and (iii) the sampling density and spatial configuration of primary and ancillary variables (Simbahan et al., 2006).

Both kriging and cokriging, like any other method of estimation, involve an error. This error is due to the fact that the variable at the location to be estimated is generally somewhat different from the estimated value (Journel and Huijbregts, 1978). Cross-validation is a statistical method for validating a predictive model, which involves the creation of two subsets of the data (training and validation subsets). The training subset is used to fit the prediction model, while the validation subset is held out for model evaluation. Averaging the quality of the predictions across the validation sets yields an overall measure of prediction accuracy. In order to run and compare the accuracy of the prediction models, the whole dataset containing the 96 data of SOC and relief attributes was randomly split 75/25 into training and validation sets. Thus, the dataset used for model development (training) consisted of 72 points, while the dataset used for validation consisted of 24 points (validation). The distribution of

Table 6 Matrix of linear correlation coefficients. Attributes

SOC (100 cm)

SOC (30 cm)

Sanda

Silt + claya

Claya

Sand + silta

Clayb

Slope

Asp

Curv

CTI

SOC 100 cm SOC 30 cm Sanda Silta Claya Sand + silta Clayb Slope Asp Curv CTI

1.00 0.53 −0.05 −0.01 0.14 −0.14 0.30 0.11 −0.04 0.25 −0.27

1.00 −0.07 0.02 0.13 −0.13 0.25 0.04 −0.03 −0.01 −0.02

1.00 −0.91 −0.51 0.51 −0.39 0.14 0.05 −0.05 −0.07

1.00 0.10 −0.10 0.14 −0.30 0.02 0.02 0.23

1.00 −1.00 0.64 0.28 −0.16 0.09 −0.31

1.00 −0.64 −0.28 0.16 −0.09 0.31

1.00 0.28 −0.08 0.13 −0.40

1.00 0.01 0.04 0.60

1.00 0.05 0.18

1.00 −0.38

1.00

Asp. — Aspect; Curv. — Curvature; CTI — compound topographic index. Pearson Correlation values are significant in bold (alpha = 0.050). a Sand, silt + clay, clay and sand + silt in the A-horizon. b Clay in the B-horizon.

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

65

Table 7 Descriptive statistics of the training and validation dataset for soil organic carbon stock across the 0–100 cm soil depth. Dataset

n

Mean

SD.

C.V.

Min.

Max.

Skewness

Kurtosis

KS (p)

Whole set Training Validation

96 72 24

7.49 7.47 7.01

2007 2.06 1.83

27.29 27.58 26.17

3.26 3.26 3.84

13.05 11.93 11.87

0.14 −0.08 0.65

−0.19 −0.44 0.77

– 0.9971a 0.2106b

KS (p) — Kolmogorov–Smirnov normality test probability. (p) values lower than 0.05 are not normal distributed. a Both training and whole dataset share a common normal distribution. b Both validation and training dataset share a common normal distribution.

both the training and the validation dataset along the study area is presented in Fig. 3. Both ordinary kriging and isotopic cokriging were performed using a data file composed by 72 collocated data of SOC stock up to 0–30 and 0– 100 cm depth, and the relief attributes (slope, aspect, curvature and CTI). A second data file containing 315 relief attributes data (72 collocated in soil profile with SOC stock and 243 with only relief attributes) were used to predict SOC stock (partially heterotopic cokriging). Important to highlight that, as recommended by Wackernagel (2003), the same cross semivariogram model fitted on the basis of the isotopic subset of the data (72 soil profiles) was used to predict SOC stock using partially heterotopic cokriging (dataset with 315 data). The predicted values (Pi) of SOC stock using OK, ICOK and HCOK were compared with the 24 observed values (Oi) of the validation dataset (Fig. 4). The differences between the 24 predicted and observed values (Pi − Oi) were used to calculate the following error measurements: Mean Error (ME), Mean Absolute Error (MAE) and Root Mean Square Error (RMSE), as follows: ME ¼

n 1X ðP −Oi Þ n i¼1 i

MAE ¼

n 1X jP −Oi j n i¼1 i

" RMSE ¼

n 1X 2 ðP −Oi Þ n i¼1 i

#0:5 :

The ME is used for determining the degree of bias in the estimates, often referred to as the bias (Isaaks and Srivastava, 1989). According to Willmott (1982), MAE and RMSE are among the best overall measures of model performance, as they summarize the mean difference in the units of observed and predicted values. Compared to MAE, the RMSE is more sensitive to extreme values. All of descriptive statistical and geostatistical analyses were conducted using the software Geoestat (Vieira et al., 1983). The comparison between the average values of SOC, clay content and relief attributes according to different SRV unities was performed by the R Development Core Team (2014) (R software (version 3.1.1). Kriged and cokriged SOC stocks maps were generated using the program ARCGIS 9.3 (ESRI, 2006). 3. Results and discussion 3.1. The SOC stock and its relationship to landscape attributes Descriptive statistical parameters at depths of 0–30 and 0–100 cm of SOC stock, and all the results of particle size fractions and relief attributes for the 96 soil profiles and the 315 survey locations, respectively, are presented in Table 3. In general, the soil textural compositions show a predominance of sand and silt fractions, reflecting the characteristics of the parent material of the soil. With the exception of the data of sand + silt contents of the A horizon, the coefficients of variation are high, especially those related to relief (slope, aspect and curvature). Regarding the frequency distribution, all granulometric attributes and SOC

stocks up to 0–100 cm deep followed a normal distribution, while the relief attributes (slope, curvature and CTI) and the SOC stock up to 0– 30 cm deep presented a log-normal distribution. The SOC stock values in the layer at 0–30 cm, ranged from 1.54 to 6,44 kg m−2, with an average value of 3.46 kg m− 2. Increasing the layer thickness from 30 to 100 cm soil doubles the carbon stock, once the minimum and maximum values changed to 3.26 and 13.05 kg m−2, respectively, with average values of 7.49 kg m−2. Based on these average values, around 46% of the SOC stocks are concentrated in the upper 0–30 cm soil layer, which is more sensitive to changes in land use and deforestation than the deep soil layers (Six et al., 2002). Commonly the SOC stock for the Amazon region is presented in association with the main soil classes and soil mapping units. In order to compare our results with those observed in other studies we also present the SOC stock according to soil mapping units (Table 2). The SOC stock values in the layer at 0–30 cm, ranged from 2.53 (Consociation– PACd) to 3.97 kg m−2 (complex CXal–PVAa), with an average areaweighted value of 3.28 kg m−2. The same mapping unities (Consociation–PACd and complex CXal–PVAa) also presented the lowest and the highest SOC stock at 0 to 100 cm soil depth (5.01 kg C m− 2 and 8.73 kg C m−2 respectively), with an average weighted value of 7.32 kg m−2. If these average weighted values were used to estimate the SOC stock across the upper 30 and 100 cm soil depth for all Içá Formation, the total SOC stock would be of 1.85 Pg and 4.12 Pg, respectively. However, the Içá Formation covers an area representing about 36% of the Amazon State. Therefore, any extrapolation considering data from site-specific research must be executed carefully because the validity of the assumption that the conditions of the study area would hold over this huge territory would need to be proven based on more field work along the Içá Formation. In order to compare our results with the results from the literature, we converted the average values presented in Table 1 to kg C m−2 (Table 4). Irrespective of the authors, the results presented from other studies were systematically higher than our results (for 30 and 100 cm soil depth). The average overestimation for the 0–30 cm soil depth was 34%, and decreased to 16% for the 0–100 cm soil depth. An overestimation of 16% is quite acceptable, considering that the SOC stocks estimated in the cited studies were based on different sources of soil and vegetation data and at exploratory scale. Besides, the dry bulk density was not measured but estimated using different models. Considering that, our results were based on measurements of carbon content and bulk density, the models used in other studies seem to perform better for the 0–100 cm compared to the 0–30 cm depth. On the other hand, the systematic overestimation of 34% in the 0–30 cm soil depth calls the attention towards the necessity of achieving a better understanding of the relationship between SOC stock, soil types, relief and vegetation along the amazon region as well as the impact of estimation protocols on the result. Our results indicate that across both depths considered here, well drained soils (complex CXal–PVAa, for example) exhibit higher SOC stocks than poorly drained soils (Consociation– PACd). This behavior can explain the causes of the differences between our results and of the other authors. According to Stevenson (1982), SOC stock is a function of five factors of formation: parent material, climate, relief, organisms and time. Hence, an analysis of the factors of soil formation in the area may help explain the data shown in Tables 2 and 5 and the differences with

66

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

35

a

30

0.2 0.15

25 20

0.1

15 10

0.05

5

Slope 72 Data

Log Slope 315 Data

0

0 0

2000

4000

6000

8000

10000 12000 14000

1.4

b

1.2

0 0.0003

2000

4000

6000

8000

10000 12000 14000

g

0.00025

1.0

0.0002

0.8

0.00015

0.6

0.0001

0.4 0.2

0.00005

Log Curvature 72 Data

0.0

Log Curvature 315 Data

0

0 16000

2000

4000

6000

8000

10000 12000 14000

c

14000

Semivariance

f

0

2000

4000

6000

8000 10000 12000 14000

0.25

h

0.20

12000 10000

0.15

8000 0.10

6000 4000

0.05

2000

Aspect 72 data

Log Aspect 315 Data

0

0.00 0

2000

4000

6000

0

8000 10000 12000 14000

1.4

d

1.2

4000

6000

8000 10000 12000 14000

1.4

i

1.2

1.0

1.0

0.8

0.8

0.6

0.6

0.4

2000

0.4

0.2

0.2

Log CTI 72 Data

0.0

Log CTI 315 Data

0.0 0

2000

4000

6000

8000

10000 12000 14000

7.0

e

CS100 72 data CS30 72 Data

6.0

0 0 0.00

5.0

-0.02

4.0

-0.04

3.0

-0.06

2.0

2000 4000 6000 8000 10000 12000 14000 2000 4000 6000 8000 10000 12000 14000 CS100 x Log CTI 72 Data

-0.08

1.0

-0.10

0.0 0

2000

4000

6000

8000

10000 12000 14000

-0.12

j

Distance (m) Fig. 6. Experimental semivariograms of relief attributes (slope, curvature, aspect, Log CTI 72 data, 6a, b, c, d, and Log slope, curvature, Log aspect, Log CTI 315 data, 6f, g, h, i, respectively), SOC stock at 0–30 cm and 0–100 cm (6e) and cross-semivariograms of SOC stock 0–100 cm × Log CTI-72/315 data (6j), with their linear model of coregionalization fitted.

other studies presented in Table 4. According to RADAMBRASIL (1978), the relief is the main factor influencing the variability of soil types and its attributes, as well as the vegetation in the Içá Formation. All soils of the Içá Formation were developed from Tertiary–Quaternary sediments and are composed of clay, silt and very fine sand. During the periods of

Pleistocene and Holocene climate fluctuations exhibited wet to dry well marked seasons, culminating in the current climate. The current climate (temperature higher than 20 °C, with no pronounced dry periods and a mean annual precipitation of 2500 mm) excludes any water deficit for plants. In fact, depending on the characteristics of the relief, the soils

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

67

Table 8 Models and their parameters fitted for semivariogram of SOC stock, relief attributes and cross-semivariograms between SOC stock and CTI. Atributes

Models

C0

C1

Sill

Range (m)

C0/(C0 + C1) %

SOC 30 cm 72 data SOC 100 cm 72 data Slope-72 data Log slope-315 data Log curvature-72 data Log curvature-315 data Aspect-72 data Log aspect-315 data Log CTI-72 data Log CTI-315 data SOC 100 cm × Log CTI 72 data

Pure nugget Spherical Pure nugget Spherical Spherical Pure nugget Pure nugget Pure nugget Spherical Spherical Spherical

1.45 3.00 28.4 0.14 0.70 0.00023 12,909 0.20 0.85 0.80 −0.025

– 1.70 – 0.038 0.38 – – – 0.21 0.25 −0.045

– 4.70 – 0.178 1.08 – – – 1.06 1.05 −0.070

– 6800 – 6800 6800 – – – 6800 6800 6800

– 64 – 79 65 – – – 80 76 36

C0 — nugget effect; C1 — contribution; Sill — C0 + C1; C0/(C0 + C1) — proportion of nugget effect in relation to the total semivariance.

of some regions become temporarily or frequently saturated with water, influencing the type and diversity of vegetation, and consequently the supply of carbon to the soil. According to Eglin et al. (2008) in waterlogging zones, the SOC stocks are lower in humus horizons and higher in the understorey vegetation. This behavior occurs because in these zones, the tree density is low, the litterfall is small and the understorey vegetation receives more sun radiation. Fig. 4 shows a scheme of the combined effects of relief, soil and vegetation on soil carbon stocks. The soil developed in relief forms classified as C11 (Dried out areas on flat-topped) and T21 (Tabular Interfluves), are well drained, with no water deficits for plants and with higher clay content (MU 2, 3, 4, 5 and 6). The photos of two soil profiles in the left side of Fig. 4 represent these environments. The soil colors are yellowish and reddish, denoting a well-drained soil. The Upland Dense Tropical Rainforest (Fda) is the main vegetation type covering these soils and the input of carbon is higher not only for the leaves and tree trunks but also due to roots growing into deeper layers of the soil. In Fig. 5A, we present the picture of this type of vegetation, which has a higher plant density, and also larger height and trunk diameter. On the other hand, the soils observed in relief forms classified as APf (River plains) and EP2 (Bi-plained superficies-flatlands), are poorly and imperfectly drained, respectively (MU1 and MU7). The photos of two soil profiles in the right side of Fig. 4, shows grayish soil colors, typical for a poorly drained soil. The Flooded Lowland Open Tropical Rainforest (Fac) and Upland Open Tropical Rainforest (Fdb) are the main vegetation covering these soils, once its species are more adapted to soil aeration restriction, mainly at subsurface. Consequently, the carbon stocks are lower, since both the input of carbon from leaves, tree trunks and roots are lower and concentrated at surface layers. In the lowlands, near rivers and water courses, the highest density of palm trees is common (Fig. 5B.). In the Upland Open Tropical Rainforest (Fig. 5C), both the density and trunk diameter of the plants are lower. Yet in Fig. 4 (bottom), we present the statistical comparison of the weighted mean values of the attributes SOC, clay content, slope and CTI, considering different landscape units, which were defined by the combination of the Solo-Relief-Vegetation units (SRV). The landscape units SRV1 and SRV4, presented the lowest SOC stock at 0–100 cm, as well was the lowest clay content (both A and B horizons) and slope. These two SRV units also present the highest CTI value. On the other hand, the SRV2 and SRV3 units present the highest SOC stock down to 100 cm, the highest clay content and slope, and the lowest CTI values. The results compiled in Fig. 4, in the form of landscape units (SRV) provide subsidies for a better understanding of how the carbon stocks Table 9 Cross-validation parameters of prediction models for SOC (100 cm). Prediction methods

ME

MAE

RMSE

Ordinary kriging Isotopic cokriging Heterotopic cokriging

0.299 0.237 0.209

1.595 1.598 1.538

1.986 1.995 1.913

in the upper 100 cm soil depth relate to the topography and the vegetation of the study area. In the present study it was not possible to introduce quantitative information of the vegetation (such as NDVI and EVI), which could further elucidate the understanding of the carbon stock distribution in the region. Considering the available auxiliary data on relief, we highlight the relationship between SOC stock and CTI, since this index already reflected the slope in the calculation and may be largely determined throughout the study area. The soil profiles with lower CTI exhibit higher values of SOC stocks and higher clay content. This pattern is confirmed by the significance of the correlation coefficients shown in Table 6. Even though the correlation coefficient is relatively low, the potential contribution of CTI index for carbon stock prediction along the area was explored in the next section of this study. 3.2. The spatial variability of SOC stock and the performance of geostatistical interpolation methods As the dataset with 96 SOC stock values was split into 75/25% for training and validation, respectively, it is important to assure that both subsets share the same distribution function. According to the results presented in Table 7, both subsets showed the same normal distribution function. Fig. 6 and Table 8 display the experimental semivariograms and cross-semivariograms, as well as, list their fitted parameters for the SOC stock (0–30 and 0–100 cm), relief attributes, and crosssemivariograms between the variables CTI and SOC stock in the 0– 100 cm depth. The SOC stock at 0–100 cm exhibited spatial dependence (Fig. 6e), while that in 0–30 cm varied randomly in space (pure nugget effect). The random behavior of the SOC stock in the upper 30 cm may be a consequence of the carbon accumulation in the upper soil layer being rather strongly associated with factors external to the soil, such as specific local variations of vegetation cover. However, when calculating the SOC stock across greater depths, for example SOC stock between 0 and 100 cm, the influence of pedogenetic soil characteristics dominated the spatial carbon dynamics manifested in a structured semivariance behavior. The structured variability allows generating a map of SOC stocks across the upper 100 cm soil depth by ordinary kriging. It also shows the importance of monitoring SOC over depths greater than 30 cm. Ordinary cokriging (isotopic and heterotopic) is more demanding than other algorithms inasmuch as two semivariograms (primary and secondary variable) and their cross-semivariograms need to be known and a common model of co-regionalization needs to be defined. In other words, to take the advantage of cokriging, besides the semivariogram of SOC stock up to 100 cm, the definition of the semivariogram of one of the relief attributes and the crosssemivariograms between both variables becomes necessary. Considering this demand, only CTI was used as a secondary variable in cokriging (isotopic — 72 data and heterotopic — 315 data), as it exhibited spatial structure (Fig. 6d and i, respectively) and covariance structure with SOC up to 100 cm (Fig. 6 j).

68

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

Fig. 7. Maps of SOC stock (kg C m−2) spatial variability at 0–100 cm by ordinary kriging (a), isotopic cokriging (b), heterotopic cokriging (c), kriging variance (d), isotopic cokriging variance (e), and heterotopic cokriging variance (f).

In all experimental semivariograms and cross-semivariograms, the model that best fitted was spherical with a range of 6800 m. Comparing the parameters of the semivariogram models for SOC stocks up to 100 cm with the cross-semivariogram (SOC stock data × Log CTI), it became clear that the use of CTI, as a secondary variable, reduced the proportion of the nugget effect in relation to the sill. The lower the nugget effect proportion, the higher is the effect of the distance of the data locations on the cokriging estimation. The lower nugget effect proportion in the cross-semivariogram means that the cokriging variance is lower than the kriging variance. The consideration of performance of ICOK and HCOK in comparison to OK, as shown by ME, MAE and RMSE in Table 9, indicates that only HCOK shows an improvement in prediction of SOC stock. The HCOK prediction presented lower bias (lower ME) and higher accuracy (MAE and RMSE), simultaneously. The increase of the CTI information (from 72 to 315 data) improves the quality of the estimation. These results are in accordance with Goovaerts (1998) and Odeh et al. (1994, 1995), who also found that heterotopic cokriging was superior over isotopic cokriging and ordinary kriging. Fig. 7 shows the maps of spatial variability of SOC stock by kriging, isotopic and heterotopic cokriging (Fig. 7a, b and c, respectively), as well as kriging and cokriging variances (Fig. 7d, e and f), respectively. Overall, the maps of SOC stock spatial variability have similarities in estimated values; however, the major difference is apparent when comparing the maps of kriging variance with those of cokriging variance. According to Nielsen and Wendroth (2003), both kriging and cokriging variances have to be considered cautiously with respect to the magnitude of the variance. Comparing the relative distribution of estimation

variance for different sampling strategies is the main benefit of (co-) kriging variance maps. 4. Conclusion - The SOC stocks across the upper 30 cm and 100 cm soil depth are 3.28 and 7.32 kg C m−2, respectively, and at about 45% of total SOC stock was stored in the upper 30 cm, which is the most susceptible to land use change caused by clear cutting and agricultural or grazing use. - The SOC stocks down to 30 and 100 cm soil depth, were in average, 34 and 16%, respectively, lower than those observed in other studies. This difference is probably a consequence the different mapping procedures applied in other studies, which were based on low-density data, coarse soil and vegetation maps, and lack of measured dry bulk density. Consequently, it did not capture the locally specific relations between SOC stock, relief and vegetation occurring in this study site. - Our hypothesis that the relief is the dominating processes responsible for soil and vegetation variability along the study site was confirmed. The SOC stock is higher in soils developed in relief forms exhibiting well-drained soils, which are covered by Upland Dense Tropical Rainforest. In this Soil Relief Vegetation (SRV) landscape unit, the input of carbon is higher not only through leaves and tree trunks, but also through root growth into deeper layers of the soil. On the other hand, SOC stock is lower in soils developed in relief forms with high waterlogging, which are covered by Flooded Lowland Open Tropical Rainforest or Upland Open Tropical Rainforest.

M.B. Ceddia et al. / Science of the Total Environment 526 (2015) 58–69

In this SRV landscape unity, the carbon stocks are lower, since both the input of carbon from leaves, tree trunks and the roots are lower and concentrated at surface layers. - SOC stock across the upper 100 cm exhibited spatial dependence allowing the generation of spatial variability maps based on spatial (co)-regionalization. However, reliable maps could not be produced for SOC in the surface layer (0 to 30 cm). The lack of spatial structure in this layer is probably caused by random determinants, such as specific variations of carbon sources and biological activity. - The CTI was inversely correlated with SOC stock and was the only auxiliary variable feasible to be used in cokriging interpolation. The heterotopic cokriging presented the best performance for mapping SOC stock. - Finally, the results show that the generation of new field data in a detailed scale encompassing various landscapes of the amazon region can significantly improve both the understanding of the relations between SOC stock and relief and vegetation attributes, and also the modeling and mapping of SOC stock. The quality of the digital SOC stock mapping has great potential to improve since the availability of hyperspectral satellite imagery and digital elevation models with higher spatial resolution are increasing.

Acknowledgments The authors acknowledge Petrobrás (Petróleo Brasileiro-SA) for the technical and financial support to perform this study (ContractPetrobras/UFRRJ/FAPUR, No. 0050.0036944.07.2). The authors thank the Conselho Nacional de Desenvolvimento Cientifico e Tecnológico (CNPq: 249594/2013-7) and Coordenação de Aperfeiçoamento de Pessoal de Nivel Superior (CAPES: 1574-14-0) for the Pos Doctoral Scholarship granted to Drs. Ceddia and Pinheiro. Finally, the authors would also like to thank the reviewers for their valuable comments and suggestions. References Batjes, N.H., 1996. Total carbon and nitrogen in the soils of the world. Eur. J. Soil Sci. 47, 151–163. Batjes, N.H., 2000. Effects of mapped variation in soil conditions on estimates of soil carbon and nitrogen stocks for South America. Geoderma 97, 135–144. Batjes, N.H., Bernoux, M., Cerri, C.E.P., 2004. Soil data derived from SOTER for studies of carbon stocks and change in Brazil (ver. 1.0). Report 2004/03. ISRIC – World Soil Information Wageningen, Netherlands. Batjes, N.H., 2005. Organic carbon stocks in the soils of Brazil. Soil Use Manag. 21, 22–24. Batjes, N.H., Dijkshoorn, J.A., 1999. Carbon and nitrogen stocks in the soils of the Amazon region. Geoderma 89, 273–286. Bernoux, M., Carvalho, M.D.S., Volkoff, B., Cerri, C.C., 2002. Brazil's soil carbon stocks. Soil Sci. Soc. Am. J. 66, 888–896. Bernoux, M., Volkoff, B., Carvalho, M.D.S., Cerri, C.C., 2003. CO2 emissions from liming of agricultural soils in Brazil. Glob. Biogeochem. Cycles 17, 1049. Bhatti, A.U., Mulla, D.J., Frazier, B.E., 1991. Estimation of soil properties and wheat yields on complex eroded hills using geostatistics and thematic mapper images. Remote Sens. Environ. 37, 181–191. Cerri, C.C., Bernoux, M., Arrouays, D., Feigl, B., Piccolo, M.C., 2000. Carbon pools in soils of the Brazilian Amazon. In: Lal, R. (Ed.), Global Climate Change and Tropical Ecosystems. Advances in Soil Science. CRC Press, Boca Raton, pp. 33–50. de Moraes, J.L., Cerri, C.C., Melillo, J.M., Kicklighter, D., Neil, C., Skole, D.L., Steudler, P.A., 1995. Soil carbon stocks of the Brazilian Amazon Basin. Soil Sci. Soc. Am. J. 59, 244–247. Eglin, T., Walter, C., Nys, C., Follain, S., Forgeard, F., Legout, A., Squividant, H., 2008. Influence of waterlogging on carbon stock variability at hillslope scale in a beech forest (Fougèrest forest-West France). Annals of Forest Science 65. Springer Verlag, Germany, p. 20. Embrapa 1981. Mapa de Solos do Brasil, escala 1:5,000,000. Servico Nacional de Levantamento e Conservação de Solos, Rio de Janeiro. Embrapa. Centro Nacional de Pesquisa em Solos, 1997. Manual de métodos de análises de solo. 2nd ed. Documentos 1. Embrapa CNPS, Rio de Janeiro (212 pp.).

69

Embrapa. Centro Nacional de Pesquisa em Solos, 1999. Sistema Brasileiro de Classificação de Solos. 2nd ed. Embrapa, Brasília (412 pp.). ESRI, 2006. Environmental Systems Research Institute, ArcInfo 9.3. Redlands, CA. Fidalgo, E.C.C., et al., 2007. Estoque de Carbono nos Solos do Brasil. Rio de Janeiro: Boletim de Pesquisa e Desenvolvimento. Embrapa Solos, p. 121. Galvão, P.H.F., Demétrio, J.G.A., Souza, E.L., Pinheiro, C.S.S., Baessa, M.P. Martins, 2012. Hidrogeologia e geometria dos aquíferos das formações cretáceas Içá e Solimões, Bacia Paleozoica do Solimões, na região de Urucu, Amazonas. Rev. Bras. Geogr. 42 (Suppl. 1), 142–153 (Dezembro. www.sbgeo.org.br). Gessler, P.E., Moore, I.D., Mckenzie, N.J., Ryan, P.J., 1995. Soil-landscape modelling and spatial prediction of soil attributes. Int. J. Geogr. Inf. Syst. 9, 421–432. Goovaerts, P., 1998. Ordinary cokriging revisited. Math. Geol. 30, 21–42. Goulard, M., Voltz, M., 1992. Linear coregionalization model: tools for estimation and choice of cross-variogram matrix. Math. Geol. 24 (3), 269–286. Hengl, T., Heuvelink, G.B.M., Stein, A., 2004. A generic framework for spatial prediction of soil variables based on regression-kriging. Geoderma 120, 75–93. IBGE, 1988. Mapa de vegetac¸a˜o do Brasil, escala 1:5,000,000. Fundacao Instituto Brasileiro de Geografia e Estatistica, Rio de Janeiro. IBGE, 2004. Instituto Brasileiro de Geografia e Estatistica. Mapa de Vegetacao do Brasil, Rio de Janeiro Escala 1: 5 000 000 Avaiable: ftp://geoftp.ibge.gov.br/mapas/tematicos/ mapas_murais/. IPCC/UNEP/OECD/IEA, 1997. Revised 1996 IPCC Guidelines for National Greenhouse Gas Inventories: Reporting Instructions (Vol. Intermediary Calculation to Derive These Annual Fluxes 1); Workbook (Vol. 2); Reference Manual (Vol. 3). Intergovernmental Panel on Climate Change, United Nations Environment Programme, Organization for Economic Co-Operation and Development, International Energy Agency, Paris. Isaaks, E., Srivastava, R., 1989. An Introduction to Applied Geostatistics. Oxford University Press, New York (561 pp.). Journel, A.G., Huijbregts, Ch.J., 1978. Mining Geostatistics. Academic Press, London. Moore, I.D., Gessler, P.E., Nielsen, G.A., 1993. Soil attributes prediction using terrain analysis. Soil Sci. Soc. Am. J. 57, 443–452. Nielsen, D.R., Wendroth, O., 2003. Spatial and Temporal Statistics: Sampling Field Soils and Their Vegetation. Catena Verlag, Reiskirchen (398 pp.). Odeh, I.O., McBrateney, A.B., Chittleborough, D.J., 1994. Spatial prediction of soil properties from landform attributes derived from a digital elevation model. Geoderma 63, 197–215. Odeh, I.O., McBrateney, A.B., Chittleborough, D.J., 1995. Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging. Geoderma 67, 215–226. Parton, W.J., Schimel, D.S., Cole, C.V., Ojima, D.S., 1987. Analysis of factors controlling soil organic matter levels in Great Plains. Soil Sci. Soc. Am. J. 51, 1173–1179. Pastor, J., Post, W.M., 1986. Influence of climate, soil moisture and succession on forest carbon and nitrogen cycles. Biogeochemistry 2, 3–27. Powers, J.S., Schlesinger, W.H., 2002. Relationships among soil carbon distributions and biophysical factors at nested spatial scales in rain forests of northeastern Costa Rica. Geoderma 109, 165–190. R Development Core Team, 2014. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Version 3.1.1. 3-90005107-0 (URL http://www.R-project.org). RADAMBRASIL, 1973. Levantamento de Recursos Naturais (Geologia, Geomorfologia, Pedologia, Vegetacao, Usa Potencial da Terra). Ministerio das Minas e Energia (DNPM), Rio de Janeiro Brazil. Simbahan, G.C., Dobermann, A., Goovaerts, P., Ping, J.L., Haddix, M.L., 2006. Fine resolution mapping of soil organic carbon based on multivariate secondary data. Geoderma 132, 471–489. Six, J., Conant, R.T., Paul, E.A., Paustian, K., 2002. Stabilization mechanisms of soil organic matter: implications for C-saturation of soils. Plant Soil 241, 155–176. Sombroek, W.G., Nachtergaele, F., Hebel, A., 1993. Amounts, dynamics and sequestering of carbon in tropical and subtropical soils. Ambio 22 (7), 417–426. Stevenson, F.J., 1982. Humus Chemistry. Wiley, New York. Thompson, J.A., Pena-Yewtukhiw, E.M., Grove, J.H., 2006. Soil landscape modeling across a physiographic region: topographic patterns and model transportability. Geoderma 133, 57–70. Vieira, S.R., Nielsen, D.R., Biggar, J.W., 1983. Geostatistical theory and application to variability of some agronomical properties. Hilgardia 51 (3), 1–75. Villela, A.L.O., 2013. Mapeamento Digital de Solos da Formação Solimões, Urucu-AM. (PhD Thesis), CPGA-CS. UFRRJ, Seropedica, RJ-Brazil (101 pp.). Wackernagel, H., 2003. Multivariate Geostatistics: An Introduction With Applications. 3rd edition. Springer-Verlag, Berlin (387 pp.). Walkley, A., Black, I.A., 1932. An examination of the Degtsjareff method for determining soil organic matter and proposed modification of the cromic acid titulation method. Soil Sci. 29–38. Willmott, C.J., 1982. Some comments on the evaluation of model performance. Bull. Am. Meteorol. Soc. 63, 1309–1313.