Global‐Scale Convergence Obscures Inconsistencies in Soil Carbon Change Predicted by Earth System Models

Soil carbon (C) responses to environmental change represent a major source of uncertainty in the global C cycle. Feedbacks between soil C stocks and climate drivers could impact atmospheric CO2 levels, further altering the climate. Here, we assessed the reliability of Earth system model (ESM) predictions of soil C change using the Coupled Model Intercomparison Project phases 5 and 6 (CMIP5 and CMIP6). ESMs predicted global soil C gains under the high emission scenario, with soils taking up 43.9 Pg (95% CI: 9.2–78.5 Pg) C on average during the 21st century. The variation in global soil C change declined significantly from CMIP5 (with average of 48.4 Pg [95% CI: 2.0–94.9 Pg] C) to CMIP6 models (with average of 39.3 Pg [95% CI: 23.9–54.7 Pg] C). For some models, a small C increase in all biomes contributed to this convergence. For other models, offsetting responses between cold and warm biomes contributed to convergence. Although soil C predictions appeared to converge in CMIP6, the dominant processes driving soil C change at global or biome scales differed among models and in many cases between earlier and later versions of the same model. Random Forest models, for soil carbon dynamics, accounted for more than 63% variation of the global soil C change predicted by CMIP5 ESMs, but only 36% for CMIP6 models. Although most CMIP6 models apparently agree on increased soil C storage during the 21st century, this consensus obscures substantial model disagreement on the mechanisms underlying soil C response, calling into question the reliability of model predictions.


Introduction
Future carbon (C) storage in soil and its potential for climate mitigation remain uncertain in all global change scenarios, especially at a global scale.On one hand, soils could gain soil organic C (SOC) due to CO 2 fertilization via plant photosynthesis (Terrer et al., 2021) and land management (Bossio et al., 2020).On the other hand, empirical evidence points to large C losses from warming (Nottingham et al., 2020;Soong et al., 2021) and associated disturbance (Pellegrini et al., 2018;Walker et al., 2019).Earth system models (ESMs) are the main tool for predicting global SOC change under future climate scenarios.Although society relies on ESMs for designing climate policy and mitigation strategies, ESM predictions remain uncertain.Global SOC changes predicted by ESMs in the fifth Coupled Model Intercomparison Project (CMIP5) range from a 270 Pg C gain to a 50 Pg C loss (Todd-Brown et al., 2014).Consequently, researchers have sought to improve the reliability of SOC models in ESMs.Reliable models are both accurate-their encoded mechanisms and outputs match the real world-and certain, meaning that their predictions are consistent.Accurate models are not very useful in decision-making if their predictions are highly uncertain.Models that are certain but inaccurate can be misleading because they make biased, incorrect predictions with apparently high confidence.Current ESM projections in soil carbon have low confidence due to outdated knowledge about soil C turnover (Todd-Brown et al., 2014).Addressing structural uncertainties, rather than refining parameters, is crucial for improving confidence in these projections.Initially, this approach may widen the range of soil carbon projections, but it will ultimately enhance confidence in managing climate change.
Given the complexity of mechanisms represented by ESMs, uncertainty in SOC predictions should be expected (Bradford et al., 2016).For example, whether tropical forests will gain or lose C is an open question.Productivity in these forests can respond positively to elevated CO 2 (Schimel et al., 2015), while SOC loss may be less sensitive to warming temperatures such as Q 10 (Koven et al., 2017), potentially allowing tropical forests to store more SOC in the future.Alternatively, tropical forests may experience regional warming-induced drought (Tao et al., 2022), which may cause large areas of forest mortality and consequently SOC loss.In permafrost regions, soils may lose C with climate change due to more intense warming (Rantanen et al., 2022), vulnerable C stocks (Mishra et al., 2021;Natali et al., 2021), and high temperature sensitivity (Davidson & Janssens, 2006;Koven et al., 2017).Permafrost thaw can also create thermokarst landscapes and expose deep SOC to microbial decomposition (Schuur et al., 2015).On the other hand, release of mineral nitrogen due to warming may enhance plant growth (Melillo et al., 2002) and lead to SOC gain in permafrost soils (Koven et al., 2015).Moreover, the dominance of these mechanisms may vary with time.A recently model analysis of the arctic-boreal carbon sink revealed a tipping point in 2050-2080 when the carbon sink strength starts to decrease (Braghiere et al., 2023).
To assess model reliability within the current generation of ESMs, we evaluated the SOC change under highemission scenarios in the CMIP Representative Concentration Pathway RCP 8.5 or Shared Socioeconomic Pathways SSP 585.We analyzed 24 available ESMs used in CMIP phase 5 (CMIP5) and phase 6 (CMIP6) to quantify the change in SOC stocks from 1996 to 2100.Based on prior analysis of CMIP5 models, we expected SOC stocks in CMIP6 models to increase over this time (Todd-Brown et al., 2013).We evaluated key mechanisms underlying the SOC changes, including the strength of the CO 2 fertilization effect on ecosystem productivity, the warming effect on soil C turnover, and local conditions such as climate and soil C stocks.We determined the relative importance (RI) of these mechanisms.In addition, we built random forest (RF) surrogate models for SOC dynamics in these ESMs to quantify the uncertainty in SOC change contributed by model structure and model inputs.

Data Sources
We collected soil C stock outputs from 24 ESMs (Table S1 in Supporting Information S1) based on the availability of the model outputs in CMIP5 (Arora et al., 2013) and CMIP6 (Arora et al., 2020).Results were analyzed from historical experiments with RCP 8.5 in CMIP5 models and SSP585 in CMIP6 models.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 CO 2 mole fraction ranging from 378 to 935 ppm.We used model runs with prescribed CO 2 concentrations for consistent comparison of the strength of the CO 2 fertilization effects among models.

Biome Definition
A biome map of 1°× 1°resolution (Figure S1 in Supporting Information S1) was used for consistent biome-level analysis (Todd-Brown et al., 2013).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, grassland and savanna, cropland and urban, or permanent wetland (Todd-Brown et al., 2013).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 no regridding was done for computations with global total SOC.

Random Forest Surrogate Models for Soil Carbon Dynamics in ESMs
We used a machine learning RF algorithm to construct a surrogate model for the soil carbon dynamics in each ESM (Figure S2 in Supporting Information S1).The predictors included decadal mean of contemporary (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(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 RFs, 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.By forcing each of the 24 RF-models with all 24 sets of predictors, we obtained 576 (24 × 24) sets of SOC change.The predicted SOC change by the RF-models accounted for uncertainties from the RF model structure and RF predictors.
We computed RI for the predictors 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 when this variable was excluded in the RF algorithm RI = ΔMSE i /( ∑ n j=1 ΔMSE j ) where i is a given feature and j are the other features.
We statistically assessed the uncertainties introduced by variations in predictors and RF model structure.The uncertainty analysis is for all land pixels without dividing them into biomes.
2 ) (1) 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.X, X model , and 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.

SOC Change at a Global Scale
We analyzed 24 fully coupled ESMs to investigate the predicted SOC change under the high emission scenarios RCP 8.5 in CMIP5 and SSP585 in CMIP6.Consistent with our expectation that global SOC would increase under global change, the majority (17 out of 24) of these models predicted a net gain in SOC, with a mean of 43.9 Pg [95% CI: 9.2-78.5Pg] C gain (Figures 1a and 1b) by 2100.Among the 17 models, 11 predicted more than 50 Pg C

SOC Change at a Biome Scale
The models predicted diverse biome-level SOC changes (Figure 2, Figure S3 in Supporting Information S1).Due to temperature constraints on decomposition, we expected large SOC losses with global change in high latitude biomes but SOC gains in warm regions.In contrast to this expectation, the biomes in cold regions (i.e., tundra and boreal forest) showed divergent patterns across the models with both C gains and losses simulated.At the biome level, models did not predict consistent changes in SOC between CMIP5 and CMIP6 (Figure 2, Figure S3 in Supporting Information S1).Except for MPI, all models present in both CMIPs predicted opposing SOC changes in several biomes.CanESM predicted opposing SOC changes in all biomes between the two model versions.The two MPI models predicted the same direction of SOC change, but with a large difference in magnitude.

Predictors of SOC Change
At the grid scale, the spatial distribution of SOC change was highly predictable using a RF approach.The predictors included contemporary mean annual temperature (Ta) and precipitation (P), NPP, turnover rate (τ) and their changes at the end of the 21st century, as well as contemporary mean soil C stock (cSoil) (Figure 3, Figures S4 and S5 in Supporting Information S1).Among the predictors, Ta and NPP showed consistent increases during the 21st century (Figures S6-S10 in Supporting Information S1).On average, the nine predictors accounted for more than 80% of the variation in SOC change on a grid level and reached as high as 93% in several models (Figure 3).
The most important predictors varied among models.Change in NPP (ΔNPP) was the most important feature in nine out of the 24 models and was the second most important feature in 11 models (Figures S4 and S5 in Supporting Information S1).Contemporary mean soil C stock and air temperature were the two most dominant variables in five models, and mean annual NPP was the dominant variable in three models.When cSoil was the dominant predictor, the models usually simulated SOC loss in the high latitudes.However, when ΔNPP was the most important variable, there was not a clear trend toward either C loss or gain for a given region.Precipitation and its change were the least important predictors for nearly all models.Interestingly, ΔNPP and Ta were the dominant predictors in the three models with the largest C gains (i.e., HadGEM2-ES, MPI-ESM-LR, and bcc-csm1-1).The large C loss in GISS-E2-R was associated with cSoil and Δτ as dominant predictors.In addition, the dominant predictors changed in five out of the eight models that were compared across CMIP5 and CMIP6.For example, the dominant predictor was ΔNPP in bcc-csm1-1 but cSoil in BCC-CSM2-MR.Similarly, the dominant predictor changed from mean annual NPP to cSoil across the two versions of CESM.In the MIROC and UK ESMs, the dominant predictors remained the same with ΔNPP and cSoil for MIROC, and Ta and ΔNPP for UK ESMs.
At a biome scale, the dominant predictors for SOC change varied greatly among biomes and across models (Figure 4, Figures S11 and S12 in Supporting Information S1).The dominant factors always differed between the two versions of the same model.For example, cSoil determined soil C change in the tundra biome in MIROC-ESM, but not in MIROC-ESM2.Across the two cold biomes (i.e., tundra and boreal forest), the dominant predictors were not always consistent.For example, in both versions of CESM, cSoil determined soil C change in tundra regions, whereas mean annual temperature was more important than cSoil in the boreal forest.On the other hand, in several models, ΔNPP was the most frequent dominant factor across biomes.For example, ΔNPP was overwhelmingly important for almost all biomes in models such as bcc-csm1-1, HadGEM2-ES and MPI-ESM1-2-LR.

Uncertainty in RF Surrogate Models for Soil Carbon Dynamics in ESMs
We used RF surrogate models for soil carbon dynamics to quantify the uncertainties in SOC change from RF model structure and predictors.Including the uncertainties from surrogate model structures and predictors, soils gained 43.7 Pg C (95% CI: 37.0-50.4Pg C) during the 21st century (Figure S13 in Supporting Information S1).
The variation was mainly determined by model structure (57.7%).RF predictors only accounted for 20.2% of the variation.The remaining variation was explained by the model-predictor interaction.There was a clear difference in explained variation between CMIP5 and CMIP6 RF models (Figure 5).In CMIP5-RF results, the model structure accounted for more than 60% of the variation whereas it only accounted for 36% in CMIP6-RF results, even lower than the RF predictors (39%).Furthermore, there was a significant statistical difference in SOC change between CMIP5-RF models and CMIP6-RF models (df = 286, F = 2.06, p = 0.04) (Figure 5).

Discussion
Earth system models in CMIP5 and CMIP6 predicted an average global SOC gain under high emission scenarios by 2100.The variation in SOC change was significantly lower in CMIP6 compared with CMIP5 results.Similarly, previous studies showed increasing convergence of predicted contemporary SOC in CMIP6 relative to CMIP5 (Braghiere et al., 2023;Varney et al., 2022).However, we conclude that the apparent convergence was a result of inherent divergence in dominant mechanisms and processes in the models.In other words, the low variation (convergence) in SOC change in CMIP6 was achieved by revising different processes (divergence) among models.
The tighter clustering of modeled SOC change in CMIP6 likely reflects a "regression toward the mean" phenomenon rather than a real increase in model performance.When we examined the key drivers of variation in SOC change, it was clear that different models arrived at similar answers for very different reasons.Model developers may have adjusted model parameters and structures to reduce large deviations from data-driven estimates of global soil C stocks following CMIP5 analyses (Todd-Brown et al., 2013, 2014).Still, a high risk of model bias remains due to missing process representations and limited diversity in model structure.Additionally, more SOC observations are also beneficial for model development.If the mean SOC itself contains large uncertainty, "regression toward the mean" would be less meaningful.
Increasing confidence in ESM predictions requires reducing spread among models, but only after fully representing key Earth system processes and their associated uncertainties in the model structure (Bradford et al., 2016).Notably, none of the ESMs except GFDL represents explicit microbial processes (e.g., direct uptake of soil carbon by microbes such as through Michaelis-Menten kinetics), yet models including microbial processes predict rather different C dynamics (Allison et al., 2010;Wieder et al., 2013) and sink capacity (Hararuk et al., 2015;Shi et al., 2018).Many CMIP6 models omit or underrepresent other fundamental mechanisms critical for long-term C dynamics, such as nutrient constraints on vegetation production under global change; large-scale vegetation mortality risk under warming-associated long-term drought and fire or insect outbreak in forest ecosystems; permafrost thaw leading to large and fast carbon release from tundra regions; land-use changes and management practices such as deforestation, biomass burning and forest degradation.Soil carbon saturation is another mechanism which may prevent continuing soil carbon gain.

Negative Feedback of Soil Carbon to Global Change
Our analysis calls into question current ESM predictions of SOC balance.In most CMIP6 system models, the net balance of SOC change was positive, resulting in a negative soil carbon-climate feedback during the 21st century.
In these models, higher C inputs due to elevated CO 2 and warming more than compensated for C losses triggered by climate change.Moreover, the rate of SOC increase was not predicted to decline, meaning that these models are projecting soils to remain C sinks beyond the 21st century (Koven et al., 2021).
At the same time, many empirical studies conflict with these projections.A recent meta-analysis of manipulative field experiments shows that soil did not gain C under the combined treatment of warming and CO 2 fertilization (Bai et al., 2023).Further evidence shows that elevated CO 2 stimulates the accumulation of soil C in the short term, but the SOC gains do not persist in the longer term (van Groenigen et al., 2017), potentially due to warminginduced SOC loss (Bai et al., 2023).Synthesizing >100 CO 2 enrichment experiments, Terrer et al. (2021) reported a trade-off as CO 2 concentrations rise, with larger increases in plant biomass associated with lower soil carbon gains due to priming effects.Additionally, at larger scales, empirical studies have shown declining ecosystem productivity in tropical forest due to land use and climate change (Gatti et al., 2023), indicating that soil carbon loss may follow.
The large differences in soil carbon-climate feedback between the two CMIP phases are striking.None of the eight models included in both CMIPs showed similar predictions of SOC change across the versions.Several models even predicted opposite responses of global total soil C. The C sinks significantly shrank in CMIP6 relative to CMIP5 for the three models that predicted the largest SOC gains.Particularly, the two largest SOC gains with more than 200 Pg in CMIP5 models (HadGEM2-ES and MPI-ESM-LR) shrank to less than 100 Pg due likely to the newly added nitrogen cycle (Jones et al., 2016).Indeed, we found that the increase of NPP under global change is much less in the two models of CMIP6 than those of CMIP5.However, the shrunken SOC gain can be caused by other mechanisms.For example, neither bcc-csm1-1 nor BCC-CSM2-MR had nitrogen cycles in the land components.However, it managed to reduce SOC gain from more than 150 Pg to nearly 0. It was realized by greater sensitivity of high-latitude SOC to temperature change.These findings indicate significant changes in ESM structure and parameterization, which were likely intended to improve model realism.Still, disparate predictions of SOC change at the biome scale suggest that attempts at improving the realism of CMIP6 models may not yet have improved accuracy.

Biome-Level SOC Change
In many CMIP6 models, divergent biome-level responses led to global convergence in predicted SOC change.In half of the models, most or all of the biomes responded with small increases in soil C with global change, contributing to convergence.However, these small, positive responses were associated with different processes across biomes and models.Some models (e.g., CanESM5) simulated greater SOC due to increased productivity to global change, whereas other models predicted SOC increase in the areas with high local productivity.Furthermore, the opposite responses between biomes in cold regions and those in other regions in some models also contributed to the global convergence-SOC losses in cold regions were compensated by gains in other regions (e.g., BCC-CSM2-MR).
We expected large SOC losses with global change in high latitude biomes due to a high magnitude of warming and the destabilization of permafrost (Natali et al., 2021;Schuur et al., 2022).Contrary to this expectation, most of the CMIP models predicted SOC gains in tundra and boreal regions.We saw large C losses in only three models (GFDL-ESM2G, GISS-E2-R, and BCC-CSM2-MR) resulting from different mechanisms.Faster turnover in a warmer climate led to the large loss in GISS-E2-R whereas large amounts of soil C in the high latitudes were vulnerable to the increased temperature in BCC-CSM2-MR.A C loss up to 50 Pg was also modeled in GFDL-ESM2G related to high SOC storage.Alternatively, high latitude soils might store more C if climate warming alleviates temperature limitation of plant growth.In the Hadley and MPI models, temperature was indeed the most dominant factor for SOC change, and both models showed greater productivity.
Besides the cold biomes, the biomes in hot or dry regions also showed divergent responses.In the tropical forest, we saw either a large increase (BCC-CSM2-MR) or a decrease (MIROC-ESM) in soil C, driven by changes in NPP.Empirical studies have shown inconsistent responses of productivity to elevated CO 2 and warming in tropical forest ecosystems.These drivers could increase productivity (Schimel et al., 2015), but tropical trees may also be susceptible to mortality caused by climate warming and drought.In dry biomes, elevated CO 2 can increase plant water use efficiency, potentially boosting plant growth and SOC storage (Morgan et al., 2011).Indeed, SOC increased by 80 Pg in the shrubland and desert ecosystems simulated by MPI-ESM-LR.However, we did not see consistently large increases in the dry regions of other models.

Relative Importance in Model Predictors
Surprisingly, nine predictors could account for most of the spatial variation in SOC change (∼93%) using our RF technique.We expected to see that changes in NPP and turnover rate would be two dominant factors determining SOC change (Todd-Brown et al., 2014).Although we found dominance of ΔNPP in many models, increased NPP did not always translate into an increase in SOC.In some models, an greater NPP meant less soil C loss (see CanESM2 and MIROC-ESM in Figures 3 and 5).Furthermore, the same two factors, air temperature together with ΔNPP, predicted more than 100 Pg SOC gain in three models.These results are consistent with the idea that global change removes temperature constraints on plant growth in cold regions.In contrast to ΔNPP, changes in turnover rate rarely explained spatial variation in SOC change (but see GISS-E2-R in Figure 3).In addition, changes in climate (Ta and P) had minimal direct impact in most models, though climate may have affected NPP and turnover time indirectly.The RF approach represents nonlinear relationships among predictors that likely affected SOC change.
Contemporary SOC was also a good predictor in some of the CMIP6 models.We chose this predictor under the assumption that large SOC pools would release more C under warming, potentially dominating over the CO 2 fertilization effect on productivity.Indeed, several models showed net C losses in cold regions with high C stocks.
We caution that incomplete model spin-up could affect current SOC and its future trajectory.For example, in cold regions where models might not reach equilibrium, soils may gain C due to the internal dynamics of the models.However, it is difficult to separate and quantify the spin-up effect with the data currently available in CMIPs.

Conclusion
Our analysis refutes the idea that convergence in model-predicted SOC change is due to promising improvements in model performance.Clearly it is possible to reduce the spread in model predictions of global soil C change.However, the underlying discrepancies in biome-scale changes and driving mechanisms both within and across models are a cause for concern.Even though ESM predictions now cluster more tightly around the mean soil C gain, they all changed something different.These underlying issues suggest that the models have not in fact reached a consensus and that predictive confidence is low for soil C feedbacks to climate change.Moreover, uncertainty has not yet peaked because the models are still missing key processes.Some models may be getting the right answers for the wrong reasons, or all the models may be getting the same wrong answer (a bias).There is a risk that the current modeled potential for SOC storage is overestimated.If future models include additional processes such as microbial constraints, thermal and water dynamics in permafrost, and updated fire regimes in tropical and tundra regions, soils could turn from projected sinks into sources of CO 2 with serious consequences for climate change mitigation.

•
Earth system models (ESMs) predicted global soil C gains under the high emission scenario during the 21st century • Dominant processes underlying soil C change at global or biome scales differed among models • Model structure accounted for more than 60% variation of the global soil C change predicted by Coupled Model Intercomparison Project phase 5 ESMs, but only 36% for Coupled Model Intercomparison Project phase 6 models Supporting Information: Supporting Information may be found in the online version of this article.
gain.A gain over 150 Pg C was found in three CMIP5 models (HadGEM2-ES, MPI-ESM-LR, and bcc-csm1-1).In contrast, GISS-E2-R predicted a large SOC loss of about 150 Pg C and CanESM2 predicted a loss over 50 Pg C. SOC change varied much less in CMIP6 models compared to CMIP5 models (Figure1c).Most of the CMIP6 models simulated C gain, except for a slight loss in three of the 12 models.The largest gain of nearly 100 Pg C was predicted by CNRM-ESM2-1.Models such as MPI-ESM, MIROC, NorESM, and CNRM in CMIP6 also predicted high C gains between 50 Pg and 100 Pg C. The three models that predicted moderate SOC losses (<30 Pg C) were ACCESS, BCC and TaiESM1.Moreover, the eight models available in both CMIP5 and CMIP6 simulated rather different SOC changes between the two model versions.For example, the UK model (HadGEM2-ES) predicted a nearly 300 Pg C gain in CMIP5, but little change in CMIP6 (UKESM1-0-LL).Similarly, the MPI model in CMIP5 (MPI-ESM-LR) predicted a gain of more than 200 Pg C, which shrank to less than 100 Pg C in CMIP6 (MPI-ESM1-2-LR).The BCC model predicted a gain of more than 150 Pg C in CMIP5 (bcc-csm1-1), but little net change in the CMIP6 model (BCC-CSM2-MR).
For example, C gains were over 100 Pg C in HadGEM2-ES and over 50 Pg C in bcc-csm1-1 and MPI-ESM-LR.On the contrary, BCC-CSM2-MR lost 100 Pg C and GISS-E2-R lost 70 Pg C from cold regions.In most cases, SOC change in boreal and tundra biomes showed the same direction of response to global change within a given model.We also found model differences in biome-level SOC change in warm biomes.BCC-CSM2-MR predicted a sink of 60 Pg C in tropical forests, whereas MIROC-ESM predicted a loss of 30 Pg C there.Shrubland and desert gained 60 Pg C in MPI-ESM-LR and grasslands gained 40 Pg C in HadGEM2-ES, but other models showed little change or significant losses in these biomes.SOC change in temperate forest was less variable than for other biomes (Figure 2, Figure S3 in Supporting Information S1).Large SOC gains predicted in HadGEM2-ES, MPI-ESM-LR, and bcc-csm1-1 were contributed by different biomes.Dry regions were much larger C sinks in MPI-ESM-LR, cold regions sequestered more C in the

Figure 1 .
Figure 1.Changes in soil organic carbon under global change during the 21st century.The 12 dashed lines are the Earth system models in Coupled Model Intercomparison Project phase 5 (CMIP5) and the 12 solid lines are those in Coupled Model Intercomparison Project phase 6 (CMIP6).The first eight models reported results in both CMIP5 and CMIP6.The cumulative change of soil organic C (SOC) was tracked from 1996 to the end of the 21st century (a).The probability distribution of the cumulative change showed an average C gain of 43.9 Pg (95% CI: 9.2-78.5Pg) (b).The variation in SOC change declined significantly from CMIP5 models to CMIP6 models (c).

Figure 2 .
Figure 2. Biome-level change in soil organic carbon.Coupled Model Intercomparison Project phase 5 models are shown in blue text and Coupled Model Intercomparison Project phase 6 models are shown in red text.SOC change varied between biomes in the same model and across models.See Figure S1 in Supporting Information S1 for the biome delineations.

Figure 3 .
Figure 3. Global distribution of soil organic carbon change for eight models present in both Coupled Model Intercomparison Project and CMIP6.The SOC change was derived as the decadal difference between mean soil C in 2091-2100 and 1996-2005.The inset figures show the relative importance of each predictor for the soil C change generated by the random forest approach.R 2 and mean square error are provided for each model.Note that the color bars share the same scale, with units of kg C m 2 .Refer to Supplementary Figures 2 and 3 for the rest of the models.

Figure 4 .
Figure 4. Relative importance (RI) in biome-level predictors for soil organic carbon change.The figure shows the eight models present in both Coupled Model Intercomparison Project and CMIP6.A random forest approach was applied in each of the eight biomes to determine the RI of the nine predictors.Refer to Figures S9 and S10 in Supporting Information S1 for the rest of the models.

Figure 5 .
Figure 5. Uncertainty of soil carbon change in random forest surrogate models.CMIP5-RF predicted a soil C gain of 54.3 Pg (95% CI: [37.6-71.1 Pg C]), and CMIP6-RF predicted a soil C gain of 34.2 Pg (95% CI: [24.6-43.8Pg C]).The importance of model structure declined from Coupled Model Intercomparison Project phase 5 to Coupled Model Intercomparison Project phase 6.