The RothC 22 model was applied to predict climate change-driven SOC changes on global agricultural soils using the RothC version implementation in the R package soilassessment 42. Soil carbon stocks were estimated by running the model backwards, starting in 2018 and finishing in January 1919. The GSOCmap (FAO-ITPS, 2019, 1 km resolution) was used to specify the initial value of the soil carbon stocks for the 932k modelled points. In addition to carbon stocks, model driving data were used that are similar to the data compilation suggested by the FAO initiative “Global soil organic carbon sequestration potential” 21:
- clay contents in the 0–30 cm topsoil layer taken from soil texture maps of the OpenLandMap with a 1 x 1 km resolution (https://doi.org/10.5281/zenodo.1476854)
- monthly precipitation, mean temperature and potential evapotranspiration from the CRU TS v 4.03 raster dataset with a resolution of 0.5° 40
- land use of 2018 from the European Space Agency (ESA) Climate Change Initiative (CCI) – Copernicus Climate Change Service in a resolution of 300 m
- monthly land cover derived by NVDI from MODIS - MOD13A2 datasets (https://lpdaac.usgs.gov/products/mod13a2v006/), after the method suggested in 21
Prior to model runs, the model was initialised by spin-up runs to derive carbon input at equilibrium (Cinequi) and related pool distributions in 2018. These spin-up runs were done with an analytical solution of RothC 43 to minimise computational time.
Two model scenarios were run, and both explicitly ignored any changes in agricultural practices on both SOC decomposition and C inputs. Instead, only potential climate-driven changes in SOC were modelled. In scenario 1, a constant annual carbon input similar to the input at equilibrium was assumed. In scenario 2, the annual carbon input from 2018 to 1919 was derived by scaling the Cinequi using the ratio of the recent NPP (NPP(t)) and NPP in the reference period 1919–2018 (NPPref):
$$Cin \left(t\right)=Ci{n}_{equi}\times \frac{NPP \left(t\right)}{NP{P}_{ref}}$$
This is not exactly the same approach applied in 21, however both approaches produce identical C-input estimates.
NPP was estimated by the MIAMI model 44 based on annual precipitation (P) and annual mean temperature (T):
$$NP{P}_{T}=\frac{3000}{1+exp\left(1.315-0.0119T\right)}$$
$${NPP}_{P}=3000\left(1-exp\left(-0.0000664P\right)\right)$$
$$NPP=min\left(NP{P}_{T}, NP{P}_{p}\right)$$
Monthly soil water deficit, required to derive the rate modifying factor b 22 in the RothC model, was quantified in forward mode starting in January 1919. No land-use change and no change in the temporal pattern of soil coverage were assumed. Owing to computational limitations, scenario runs were done for a sample of 1% of all raster grids selected by random sampling, which resulted in about 932,389 runs given the 1 km x 1 km resolution of the underlying SOC map. Modelled SOC, carbon input, NPP, weather data and soil moisture deficit were stored for 10-year time intervals. These variables were aggregated according to the Köppen-Geiger climate classification map 45.
Raster maps of modelled soil carbon changes were spatially aggregated to a 0.1° resolution using the terra package in R 46.
Linear regression models were fitted to explain the variability of SOC stock changes (absolute and relative) across climatic zones. Initial SOC stocks, mean annual temperature (MAT), mean annual precipitation (MAP), water balance as well as changes in temperature, precipitation and water balance were used as explanatory variables. The best model was chosen based on the Akaike Information Criterion (AIC) and model residuals were visually checked for normal distribution using quantile-quantile plots.