To calculate the spatial prediction of SOCs at 30 cm depth, we generated synthetic profiles of 0–30 cm depth from the soil local database for SOCc and SOCs data. The aggregation of the horizons was carried out using the equal-area spline technique through mass-preserving spline (‘mpspline’) function (Bishop et al., 1999). The estimation of SOCs in each soil profile were calculated as:
SOCs (Kg∙m^2 )=SOC (g⁄(Kg)∙BD(Kg∙m^3)∙[1-(CRFVOL/100)]∙HSIZE(cm))
were BD is bulk density, CRFVOL is percentage of coarse fragments (above 2 mm in diameter), and HSIZE is thickness of the horizons. Due to data gaps, BD was estimated by means of a pedotransfer function adapted from a regional study (Barahona and Santos, 1981). We used the R package GSTAT for the stock estimates and ‘mpspline’ function, where the propagated error was estimated by the Taylor Series Method (Hengl and Mendes de Jesus, 2016; Heuvelink, 1998; Malone et al., 2009).
Our DSM approach was based on the SCORPAN conceptual model using the soil forming environmental factors as soil spatial prediction function (McBratney et al., 2003). We generated a covariate stack based on 34 environmental factors to predict SOC
Prior to predictive model building, a regression matrix was performed including the best correlated environmental covariates with SOC local data. To select them, we considered a balance among higher Pearson coefficient of multiple linear regression, lower error (RMSE), and lower variance inflation factor (VIF) to identify statistical redundancy (Heiberger et al., 2005). We used the Akaike information criterion (AIC) to determine the best compromise between model accuracy and model parsimony (Rossel and Behrens, 2010).
The quantile regression forest (qrf) model was performed to predict SOCc and SOCs, since estimates an approximation of the full conditional distribution of the response variable, the inferred conditional quantiles to build prediction intervals were estimated as surrogate of the value of uncertainty associated with the response variable (Meinshausen, 2006). We used the ‘quantregForest’ package ‘implemented in R software for statistical computing (Meinshausen, 2006).
==================== Data Sources =========================
- Local soil database:
The local database was derived from the LUCDEME Project generated between years 1986-2004 by the “Ministerio de Medio Ambiente de España“ and the support of “Dirección General de Medio Ambiente de la Región de Murcia” (Alias and Ortiz, 1986). This legacy database consists of a sampling of 255 soil profiles representative of soil typologies over a topographic range of 0-1700 m in altitude. For each profile, there are morphological and analytical data for each horizon (903 horizons).
- Spataial Environmental Covariates:
We used dynamic and static variables as predictors for SOC. The static variables were 16 topographic parameters derived from a local digital elevation model (DEM) using the Terrain Analyst functions included in SAGA GIS software (Conrad et al., 2015). The DEM is available from Geographic Information National Centre (Spain), resulting from interpolation of LiDAR national images, with a 25 m pixel resolution, re-sampled into 100 m.
The dynamic variables included climatic variables (precipitation and temperature) (Ninyerola et al., 2005); land cover (IGN, 2012) reclassified into 13 classes; forest structural variables, aboveground biomass of forest trees cover from LiDAR data (Durante et al., 2019); and vegetation indexes (VIs). The calculated VIs were The Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), which are associated to ecosystem functional attributes related to seasonal dynamics of net primary productivity. These indices were derived from mean annual time series images (2001-2016) of MODIS-Terra images satellite using Google Earth Engine as describe in (Arenas-Castro et al., 2019).
===========================================================
