Data sources
We collected soil C stock outputs from 24 Earth system models in the Coupled Model Intercomparison Project phase 5 (CMIP5) 29 and phase 6 (CMIP6) 30. Results from historical experiments, RCP 8.5 in CMIP5 models, and SSP585 in CMIP6 models were analyzed. The historical simulations run from 1850 to 2005 for CMIP5 models and from 1850 to 2014 for CMIP6 models. Eight models reported results in both CMIP5 and CMIP6. Note that GFDL-ESM2G and HadGEM2-ES simulations in CMIP5 start at 1861 and 1860, respectively. RCP8.5 and SSP585 are high radiative forcing scenarios with a prescribed atmospheric CO2 mole fraction ranging from 378 to 935 ppm. We used model runs with prescribed CO2 concentrations for consistent comparison of the strength of the CO2 fertilization effects among models.
Biome definition
A biome map of 1º × 1º resolution was used for consistent biome-level analysis21. The original 0.05° × 0.05° land cover map was obtained from the MODIS (Moderate Resolution Imaging Spectroradiometer) Land cover MCD12Q1 product 43. One of 16 land cover types was assigned to each 1º × 1º grid cell by taking the most common land cover from the original underlying 0.05° × 0.05° land cover map. Each 1º × 1º grid cell was assigned to one of 8 biomes: tundra, boreal forest, tropical rainforest, temperate forest, desert and shrubland, grasslands and savannas, cropland and urban, or permanent wetland21. Note that snow and ice are excluded for biome assignment. Briefly, tundra includes grassland and open shrubland north of 55°N. Boreal forest includes evergreen needleleaf forest, deciduous needleleaf forest, woody savanna north of 50°N, mixed forest north of 50°N, and closed shrublands north of 50°N. Tropical rainforest includes evergreen broadleaf forest between 25°N and 25°S. Temperate forest consists of deciduous broadleaf, evergreen broadleaf outside of 25°N–25°S, and mixed forest south of 50°N. Desert and shrubland includes barren or sparsely vegetated, open shrubland south of 55°N, and closed shrubland south of 50°N. Grassland and savanna includes woody savanna south of 50°N, savanna, and grasslands south of 55°N. Cropland and urban includes cropland, urban, and mosaic land cover. Permanent wetlands are unmodified.
All model outputs were regridded to 1° × 1° for biome- and grid-scale comparisons. Our regridding approach assumed conservation of mass. Note that we used native grids for global total SOC.
Random forest surrogate models for ESMs
We used a machine learning random forest algorithm to construct a surrogate model for each ESM. The feature variables include decadal mean of contemporary (1996–2005) net primary productivity (NPP), surface air temperature (Ta), precipitation (P), soil C content (cSoil), soil C turnover time (τ) and their changes over the 21st century except for cSoil which is the target variable. In the random forest, we used 300 trees, maximal depth of 20, and minimal sample split of 4. In total, we built 24 RF surrogate models using the regridded 1° × 1° model outputs.
We computed relative importance (RI) for the feature variables in each RF model and across land biomes. Relative importance of a given feature variable was defined as the normalized increase in mean square error (MSE) when this variable was excluded in the RF algorithm \(RI =\varDelta {MSE}_{i}/\left({\sum }_{j=1}^{n}\varDelta {MSE}_{j}\right)\) where i is a given feature and j are the other features.
Major uncertainty sources for model predictions of SOC change are model structure and external inputs. Given the high computational cost of running an Earth system model, it is difficult to assess the relative contributions of the two factors to the modeled SOC change by running factorial combinations of the two factors. Here, we built a RF surrogate model for each ESM. By forcing each of the 24 RF-models with all 24 sets of features, we obtained 576 (24×24) sets of SOC change. The predicted SOC change by the RF-models accounted for uncertainties from the model structure and external inputs.
We statistically assessed the uncertainties introduced by variations in feature variables and model structure.
$${SS}_{model}=n\times \left(\sum _{i=1}^{n}{({\overline{X}}_{model}-\overline{X})}^{2}\right)$$
1
$${SS}_{feature}=n\times \left(\sum _{i=1}^{n}{({\overline{X}}_{feature}-\overline{X})}^{2}\right)$$
2
\({SS}_{interaction}=\sum _{i=1}^{n\times n}{(X-{\overline{X}}_{model}-{\overline{X}}_{feature}+\overline{X})}^{2}\) (3),
where SS is the sum of square error and n is the number of surrogate models. \({SS}_{model}\), \({SS}_{feature}\) and \({SS}_{interaction}\) are the sum of square errors across models, features and their interactions, respectively. \(\overline{X}\), \({\overline{X}}_{model}\) and \({\overline{X}}_{feature}\) are the overall mean, the mean value across models and features, respectively. The variation explained by RF models is \({SS}_{model}/({SS}_{model}+{SS}_{feature}+{SS}_{interaction})\) and by model features is \({SS}_{feature}/({SS}_{model}+{SS}_{feature}+{SS}_{interaction})\). The remaining variation is attributed to the interaction between RF models and feature variables.