Mapping of soil organic carbon stocks for spatially explicit assessments of climate change mitigation potential

Data from a total of four sites, each covering an area of 100 km² (10 × 10 km) from three east African countries were included in the current study (figure 1). Dambidolo (35.0° East; 8.6° North) and Mbinga (35.2° East; 11.1° South) represented higher rainfall agricultural areas with a mean annual precipitation (MAP) (Hijmans et al 2005) of about 1300 mm, while Kipsing (37.1° East; 0.5° North) and Mega (38.3° East; 4.2° North) represented pastoral rangeland landscapes with drier climates (600–850 MAP). All sites were located in highland areas, with altitudes ranging from 1000 to 1500 m above mean sea level. These sites were selected to test the methods presented in the current study as they represent rangeland and smallholder agricultural systems with contrasting environments, climates and different soil types.

The field methods used in the current study are referred to collectively as the Land Degradation Surveillance Framework (LDSF) (Vågen et al 2012), which is a standardized set of procedures for the establishment of land health baselines and for monitoring of change. The aboveground components of the framework are based on the FAO Land Cover Classification System (LCCS) (Di Gregorio and Jansen 1998).

In short, within each 10 × 10 km site, 16 cluster centroids were stratified into 2.5 × 2.5 km tiles and their locations within the tiles were randomized, but buffered to avoid overlapping with neighboring tiles. Around each cluster centroid 10 sampling plots, each a 1000 m² circle, were randomly located to fall within a circular area of 1 km² using a 564 m radius from the cluster center-point. Each plot consisted of four 100 m² subplots.

At the plot level, basic site characteristics were recorded, including name, ID, georeferencing (coordinates) of the center-point, altitude, date and a photograph was taken for later reference. The data collected included information on slope and landform, vegetation cover types and strata, land use, land ownership and primary current use. Other information collected included presence/absence of soil and water conservation structures, and descriptions if present. At the subplot level, signs of visible erosion were recorded, including erosion types and rock/stone/gravel cover on the soil surface was scored. High erosion prevalence was determined if three of the four subplots had observations of visible erosion. Woody- and herbaceous cover ratings were made using a Braun-Blanquet (Braun-Blanquet 1928) vegetation rating scale from 0 (bare) to 5 (>65% cover). Woody plants (shrubs (1.5–3 m height) and trees (>3 m height)) were counted, and distance-based measurements were made using the T-square method (Krebs 1989) to determine vegetation distribution.

Soil samples were collected at each subplot, then composited into one top- and one subsoil sample for each plot. Standard top- and subsoil samples were collected at 0–20 cm and 20–50 cm depth increments, respectively. In the absence of depth restrictions, 160 topsoil and 160 subsoil samples were hence collected from each site. Any auger depth restrictions were noted (in cm) if present during soil sampling.

2.2.1. Soil cumulative mass sampling.

The soils sampled were analyzed at the World Agroforestry Centre's Spectral Diagnostics Laboratory in Nairobi, Kenya using dry combustion for SOC measurements on 10% of the samples and mid-infrared (MIR) scanning on all samples collected (Vågen et al 2010). Infrared spectroscopy (IR) is an established technology for rapid, non-destructive characterization of the composition of materials based on the interaction of electromagnetic energy with matter (Ben-Dor and Banin 1995, Shepherd and Walsh 2002, Viscarra Rossel et al 2006, Vågen et al 2006, Brown 2007). The MIR analysis was done according to methods and procedures described in Terhoeven-Urselmans et al (2010). Particle size analysis was done using a Horiba Model LA 950A2 particle size analyzer, which can analyze particles in the size range 0.04–3000 μm. Samples used in this study were subjected to dispersion with Calgon for 4 min before particle size analysis.

All analysis and modeling was conducted in R-statistics version 2.15.1 (R Development Core Team 2011).

Models for estimating the stocks of SOC to 30 cm were developed by first calculating the mass of soil carbon (mSOC) at each depth increment in the soil profile and taking the cumulative sum by depth as follows;

$$\begin{eqnarray} &&m S O{C}_{i}=\sum _{l=1}^{i}S O{C}_{l}\cdot {\mathrm{soilmass}}_{l}. \end{eqnarray} \tag{ 1 }$$

The mSOC estimates were scaled to represent a standard area of 1 m², and the change in the resulting SOC stock (kg  m⁻²) with depth was estimated by using a linear mixed effects model from the 'nlme' package in R (Bates and Pinheiro 1998, Pinheiro and Bates 2000) with three nested levels of grouping;

$$\begin{eqnarray} &&\begin{array}{@{}l@{}} \displaystyle \fl {\mathrm{SOCstock}}_{i j k}={\beta }_{0}+{X}_{i j k}{\beta }_{1}+Z i,j k{b}_{i}\\ +~{Z}_{i j,k}{b}_{i j}+{Z}_{i j k}{b}_{i j k}+{\varepsilon }_{i j k},\\ i=1,\cdot 4 j=1,\cdot 6,k=1,\cdot 0,\\ \displaystyle \fl {b}_{i}\sim \mathcal{N}(0,{\phi }_{1}),{b}_{i j}\sim \mathcal{N}(0,{\phi }_{2}),{b}_{i j k}\sim \mathcal{N}(0,{\sigma }^{2}I), \end{array} \end{eqnarray} \tag{ 2 }$$

where the intercept is represented by β₀, the depth slope by β₁, and b₁, b_ij and b_ijk are the random effects at site (i), cluster within site (j) and plot within cluster within site (k) random effects, while denotes the within group errors. The above model was used to predict SOC stocks to 30 cm (SOC₃₀) for each soil profile (n = 640).

Finally, an ensemble modeling approach using the 'randomForest' library in R (Breiman 2001) was used to predict SOC₃₀ as a function of Landsat ETM+ reflectance for each sampling plot. A random forest model represents a form of bagging (Breiman 1996), where large sets of de-correlated trees are built and the resulting predictions are averaged (Prasad et al 2006). An 'out-of-bag' test sample is held out and used to estimate model error and for the calculation of variable importance. The resulting model was applied to a set of Landsat imagery from the GLS 2005 archive for the generation of maps of predicted SOC₃₀. A data set consisting of 67% of the plots in the data set was used to develop the prediction model, while the remaining 33% was used for validation.

Differences in rainfall between sites explain about 20% of the variation in SOC₃₀ stocks, therefore there must be other factors that strongly determine SOC stocks in these areas. A closer examination of the texture of the soils included in the study shows that average sand content is higher in Mbinga (21%) than in Mega (14%), and highest in Kipsing with about 47%. The interaction between sand content and rainfall explains 47% of the variations in SOC₃₀ stocks between sites. This is most likely part of the explanation for the similar SOC₃₀ stocks on average observed in Mega and Mbinga (figure 2), despite the much higher rainfall in Mbinga. One implication of this is that SOC stocks are strongly determined by inherent soil properties such as soil texture, which are not influenced much by management. Similar relationships between sand content, MAP and soil carbon were reported by Jobbágy and Jackson (2000), which are due to the much lower protection of SOC against microbial degradation in sandy soils than in soils that are clay rich.

The data used in this study had a wide range of land cover types, ranging from open grasslands, bushlands and woodlands to croplands and forest, although forested ecosystems have a small sample size in this data set. The only site with some forest was Mbinga, where estimated SOC₃₀ in forested areas was 7 kg  m⁻² on average. Croplands were found mainly in the humid Dambidolo and Mbinga sites, and estimated SOC stocks varied strongly within these areas. The cultivated areas also had highly variable histories of conversion from native vegetation to cultivation. In Dambidolo, for example, most of the currently cultivated areas were converted from grasslands in the last 25 yr. In the drier ecosystems, grasslands had the highest estimated SOC₃₀, while shrublands and woodlands had lowest SOC₃₀ stocks overall (figure 3).

Zoom In Zoom Out Reset image size

While the effects of agricultural management practices on SOC stocks have been widely studied, the effects of soil erosion on SOC storage within land management practices are generally not well understood (Ogle et al 2005). The findings of the current study show that soil erosion affects SOC stocks strongly in all sites and across different land cover types. Hence, land management practices that lead to the control and reduction of soil erosion will be important mitigation strategies in these landscapes. Model estimates suggest that SOC stocks are lower by 0.9 kg m⁻² (p < 0.01) on average in eroded versus non-eroded areas. This effect is stronger in woodland and grassland systems in Mbinga, where SOC₃₀ stocks are lower by between 2 and 3 kg m⁻² in eroded areas (summarized in figure 3).

In the drier ecosystems, the majority of the shrublands and woodlands had a high prevalence of soil erosion, possibly due to overgrazing leading to low herbaceous cover in these sites. Also, semi-arid ecosystems generally have aggressive rainfall events that can be highly erosive. The results suggest that the impact of grazing was high in both Mega and Kipsing. Impact scores were high for grazing in more than 90% of the plots surveyed in these sites. Although preliminary, these findings suggest that maintaining woody cover in these systems is not enough to reduce soil erosion and effective climate change mitigation strategies hence need to focus on the maintenance of a continuous herbaceous cover in order to reduce soil erosion and increase the potential for carbon storage in these ecosystems.

Soil organic carbon stocks to 30 cm were predicted using Landsat ETM+ ground calibrated reflectance (middle-right panel in figure 2), with results indicating that the random forest models used have a reasonable level of prediction performance, with a multiple R-squared of 0.67 and 0.65 for calibration and validation data sets, respectively. These models are expected to become more robust in future assessments as a wider range of sites are included, including humid tropical forests.

The maps generated from these models at 30 m resolution (figure 4), clearly show the between-site differences in SOC₃₀ stocks and also the high level of variability within each individual site (figure 4). These maps provide estimates of SOC stocks that are spatially resolved enough for assessing the complexity of SOC storage in these landscapes and for the identification of hotspots where management interventions may be targeted. When combined with assessments of the spatial distribution of erosion prevalence, spatially explicit recommendations can be made about areas where mitigation potential is high. Further efforts to develop the models used for generating the maps presented here are focusing on including sites also from the humid tropics, as well as near-real time predictions of SOC stocks for monitoring purposes.

Zoom In Zoom Out Reset image size