The Eastern Plains of Colombia is located in the northeastern departments of Arauca, Casanare, Meta, and Vichada, covering an area of approximately 26 million hectares (Figure 1). The tropical and isothermal climate is characterized by distinct wet and dry seasons, and mean annual temperatures range between 24°C and 28°C in the plains and between 13°C and 21°C in the mountains. The region has a unimodal precipitation regime with an average annual precipitation of 2,991 mm. Boarded by the Andean Eastern Cordillera in the west and the Amazon in the south, the Llanos includes major tributaries of the Orinoco River such as the Arauca, Meta, and Vichada rivers. The Eastern Plains is divided into six landscapes -piedmont, alluvial terraces, alluvial overflow plain, aeolian plain, alluvium, and the high plains-with an elevation that rarely exceeds 300 m. Predominant soils in the Eastern Plains are Oxisols and Ultisols that originate from alluvial deposits transported from erosion of the Eastern Cordillera during the late Pleistocene (Goosen, 1971). Native vegetation is typical of savannah, although forest vegetation can be found along the major rivers and their tributaries.

This study examines legacy data collected by the International Center for Tropical Agriculture (CIAT) in collaboration with AGROSAVIA (before Corpoica), the Agustin Codazzi Geographical Institute (IGAC), and the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM) for Colombia’s National Forest Inventory (IFN). Data obtained from CIAT and AGROSAVIA consisted of 158 point observations with three replicates from the topsoil (0–30 cm) sampled in 2013, and 64 soil profiles from IGAC, the Colombian institute responsible for producing and hosting soil information (http://ssiglwps.igac.gov.co:8888/siga_sig/Agrologia.seam). Following the IFN sampling strategy (Peña et al., 2014) IDEAM has collected 107 soil samples from 30 sites (conglomerates of subsamples) distributed in the Arauca, Meta, and Vichada departments. Sites are characterized by five sub-parcels that were arranged in the form of a square wherein four sub-parcels were collected 80 m in each cardinal direction of a central sample. Data collection and analysis by the IDEAM and the IFN are still on going, therefore many of the sites contain a range of one to five sub-parcels. Soil organic carbon content of data provided by CIAT-AGROSAVIA and IGAC were determined by the Walkley-Black methodology (Walkey and Black, 1934), and dry combustion by the IDEAM.

We acknowledge that the models are created using SOC content determined using different laboratory methods, without the use of a correction factor to collaborate findings. However, as measurements were very similar in areas of proximity we proceeded with the full data set. In total, 653 topsoil point samples were used in this study.

To predict the amount and distribution of SOC quantities according to the scorpan (a mnemonic for factors used in the prediction of soil attributes) framework developed by McBratney et al. (2003), we used a range of primary and secondary geospatial datasets to build a stack of environmental covariates. A 90 m resolution digital elevation model (DEM) was used to derive terrain parameters using Arc-GIS and SAGA-GIS (Jarvis et al., 2008). The DEM is a product of the United States National Aeronautics and Space Administration’s (NASA) Shuttle Radar Topographic Mission (SRTM), and the data voids were filled by Jarvis et al. (2008). Prior to use for digital soil mapping, the DEM was corrected to remove depressions and sinks and to ensure a regular flow of the landscape surface. A basic terrain analysis was conducted to derive representative terrain attributes using SAGA-GIS (analytical hill shading, aspect, slope, wetness index, profile and plan curvature, convergence index, and catchment area).

Other environmental covariates included:

• A topographic wetness index (TWI) of the DEM was derived in SAGA GIS as a representative of topographic control on hydrological processes.

The amount and distribution of SOC content in the Eastern Plains of Colombia was determined by the random forest approach (Breiman, 2001) using the R environment (R Core team, 2016). Prior to the generation of the model, each environmental covariate was resampled to the resolution of the DEM and a principal component analysis (PCA) was used to reduce the number of environmental covariates and determine the primary covariates that were most representative of the landscape. The multivariate PCA was conducted in SAGA-GIS.

The exploratory analysis of the data showed that the dataset was not normally distributed (Anderson-Darling test <0.05), negatively skewed (coefficient of skew <0), and had a high degree of peakedness (Kurtosis >1). Therefore, a Box-Cox transformation was used to normalize data prior to input for spatial modeling.

Random forest modeling is a non-parametric technique that fits classification trees to the dataset, and aggregates similarities among the observations to fit decision trees. This method consists of the selection of a random sample from the dataset to build a tree. This sample was then randomly replaced by others in the dataset (bootstrapped) and added to the training dataset. The trees were then split at nodes that were defined by a random subset of the predictor variables. The size of the predictor subset was determined by the square root of the total number of predictor variables, and nodes were split based the predictor that minimizes the regression error. A “forest” was then created as this process was repeated several times, until no further improvements in error could be achieved. Data that were not used are called “out-of-bag” samples and used to determine the regression errors for the trees. In this study, a random forest model with 500 trees, and a node size of three was ran with nine predictors using the randomForest R-package (Liaw and Wiener 2002). Although there is still debate about the best digital modeling methodology, this model has several advantages when compared to other digital mapping strategies as the bootstrap samples and random selection across the many predictors helps to reduce variance among the predictors (Breiman, 2001).

A cross-validation analysis was conducted to assess the predicative performance of the random forest model. Due to the small amount of data points, 90% of the dataset was used in the training dataset to develop the model, and the remaining 10% was used for model testing. The model was cross-validated 100 times, and the performance matrix of each repetition was averaged. A diagnostic comparison of the predicted and observed SOC values was performed to illustrate the predictive performance of the random forest model. The coefficient of determination (R ²), mean error (ME), and root-mean-squared error (RMSE) were used to determine the model performance.

Soil organic carbon stock (ton/ha) was calculated using the SOC content, bulk density, and coarse rock fragment maps of the Eastern Plains according to Eq. (1). The SOC content map was the product of the random forest methodology, and the bulk density and coarse fragment maps were retrieved from the International Soil Reference and Information Center’s (ISRIC) SoilGrids automated soil mapping database. ISRIC SoilGrids maps were used for this project due to the lack of bulk density and coarse fragments data within this region. The SoilGrids maps are created using state-of-the-art machine learning methods that utilize soil profiles and environmental covariate data to globally map soil properties at a spatial resolution of 250 m. Prior to the determination of SOC stock, the bulk density and coarse fragment maps were clipped to the boundary of the Eastern Plains using GRASS-GIS. S O C s t o c k = 10 ∗ ( S O C c o n t e n t 1000 ∗ 30 100 ∗ B L D ∗ ( 100 − C R F ) 100 ) (1) Where SOC_stock is the SOC stock (t ha⁻¹), SOCcontent the soil organic carbon content (g kg⁻¹), BLD is the bulk density (units) and CRF the coarse rock fragments (units). The equation is multiplied by a factor of 10 to convert from kg m⁻² to t ha⁻¹.

The Taylor series method was used to assess SOC stock predictions throughout the region. This method was chosen in lieu of other uncertainty analyses because it allows for the identification of the main source of error in prediction estimates. Global accuracy measures for the three soil properties were used since bulk density and coarse rock fragment data either do not exist or could not be obtained for this specific study area. While global measures may not apply to the study area, they are the best estimates for this understudied region. The following standard deviations of the SoilGrids prediction error for SOC content, bulk density, and coarse rock fragment were used to calculate the standard deviation of the SOC stock map over 0–30 cm: sd (SOC) = 32.8 mg/kg, sd (bulk density) = 164.7 kg/m3, sd (coarse rock fragment) = 10.9%. The standard deviation of SOC stock predictions was determined as the square root of the prediction variance, and the relative error was calculated as the ratio of the standard deviation and the SOC stock prediction.

In total, the following nine environmental covariates were used to predict SOC in the Eastern Plains of Colombia: analytical hill shading, aspect, slope, catchment slope, TWI, PEI, geology, geomorphon, and land use. The continuous property map of SOC content (g kg⁻¹) in the Llanos region is presented in Figure 2. Soil organic carbon content was generated to the resolution of the DEM (90 m), and spatial predictions are shown to represent the amount of SOC contained in the topsoil (0–30 cm). The predicted SOC content was similarly distributed in each department, with values ranging between 3.6–35.6 (g kg⁻¹). The amount and distribution of SOC content reflects the differences in geology, and thus elevation and slope, across the landscape.

Model assessment was performed as the 100-fold cross-validation of 65 data points. Based on the cross-validation results, 50.28% of the variation within SOC content was explained by the model. RMSE is 0.461 and ME in the predictions of SOC was approximately 0.038. Since the ME value is close to zero, it indicates that there was a very low bias in the predictions and that there was no systematic over or under-prediction of the SOC content in the Eastern Plains. However, a plot of the predicted versus observed values indicated that the model under-predicted low values and over-predicted high values (Figure 3). Therefore, high and low SOC content on the predicted map should be treated with caution.

The total predicted SOC stock in the topsoil of the Eastern Plains is approximately 1.2 Gt ha⁻¹ (Figure 4). This indicates that the amount of carbon stored in the study area in tons of CO₂ equates to 3.7 Gt, which is equivalent to approximately 14 times of total Colombia emissions (258.35 Mt CO₂-eq) reported in the National Inventory of Greenhouse Gas Emissions (Ideam et al., 2016) and which represents 0.42% of the global emissions for the year 2012.

Based on the Taylor series method the standard deviation of the uncertainty propagation ranged between 30.8–171.2 (t ha⁻¹) (Figure 5). Of the three soil properties, uncertainty about the SOC content is the main source of uncertainty (99.69%), and this source of uncertainty is magnitudes higher than that of bulk density (0.16%) and coarse rock fragment (0.15%). Lower ranges of uncertainty are concentrated in areas closes to the Eastern Cordillera, while the largest amount of uncertainty were found in the southern region of the Plains and the northern regions of the Dissected High Plains. Furthermore, the mean relative error of the SOC stock predictions was 325.42%. Although not ideal, large estimates were expected since the SOC stock predictions were small in most areas of the Eastern Plains. As the ratio of the prediction error standard deviation and the prediction itself, since the standard deviation was large, the relative error of the estimates was large as well.