Observational benchmarks inform representation of soil organic carbon dynamics in land surface models

Representing soil organic carbon (SOC) dynamics in Earth system models (ESMs) is a key source of uncertainty in predicting carbon–climate feedbacks. Machine learning models can help identify dominant environmental controllers and establish their functional relationships with SOC stocks. The resulting knowledge can be integrated into ESMs to reduce uncertainty and improve predictions of SOC dynamics over space and time. In this study, we used a large number of SOC field observations (n=54 000), geospatial datasets of environmental factors (n=46), and two machine learning approaches (namely random forest, RF, and generalized additive modeling, GAM) to (1) identify dominant environmental controllers of global and biome-specific SOC stocks, (2) derive functional relationships between environmental controllers and SOC stocks, and (3) compare the identified environmental controllers and predictive relationships with those in models used in Phase 6 of the Coupled Model Intercomparison Project (CMIP6). Our results showed that the diurnal temperature, drought index, cation exchange capacity, and precipitation were important observed environmental predictors of global SOC stocks. While the RF model identified 14 environmental factors that describe climatic, vegetation, and edaphic conditions as important predictors of global SOC stocks (R²=0.61, RMSE = 0.46 kg m⁻²), current ESMs oversimplify the relationships between environmental factors and SOC, with precipitation, temperature, and net primary productivity explaining > 96 % of the variability in ESM-modeled SOC stocks. Further, our study revealed notable disparities among the functional relationships between environmental factors and SOC stocks simulated by ESMs compared with observed relationships. To improve SOC representations in ESMs, it is imperative to incorporate additional environmental controls, such as the cation exchange capacity, and refine the functional relationships to align more closely with observations.

1 Introduction

Soil is the largest actively cycling carbon pool in terrestrial ecosystems and stores almost twice the amount of carbon currently present in the atmosphere (Lal, 2016). Even a small change in soil carbon stocks can lead to large changes in the atmospheric CO₂ concentration, influencing future climate change trajectories. Additionally, soils play a crucial role in capturing atmospheric CO₂ through the storage of soil organic carbon (SOC) (Hinge et al., 2018). Therefore, the sequestration, protection, and sustainable management of SOC stocks can be a promising climate mitigation strategy (Lal, 2020). Accurately representing global SOC storage and its environmental controllers is essential for predicting realistic SOC changes under different land use and climate change scenarios. However, there is currently no consensus among Earth system models (ESMs) on how to represent the spatial distributions of global SOC storage and its fate under future climate change scenarios (Friedlingstein et al., 2014; Arora et al., 2020).

Multiple environmental variables, including climatic and topographic factors, land use history, and edaphic properties, have been identified as possible controllers of SOC storage (Luo et al., 2021; Mishra et al., 2022). However, current ESMs use only a limited number of environmental factors to represent SOC storage and dynamics. A recent study that compared SOC stocks from multiple ESMs with observations indicated a large knowledge gap in both ESMs and observations (Georgiou et al., 2021). Therefore, it is crucial to compare ESM simulations with global SOC observations to evaluate model performance and identify key environmental controllers of global SOC storage.

Benchmarking ESM simulations against observed data is a common approach for model evaluation (Luo et al., 2012; Collier et al., 2018). By comparing model simulations with observations, model strengths, deficiencies, and needed improvements can be identified. The insights gained from SOC benchmarking can lead to new ESM land model structures (by identifying key processes) and new parameterizations (by quantifying key relationships between SOC and environmental variables). Therefore, benchmarking analysis of ESMs is an effective tool for reducing uncertainties in predicting SOC dynamics and can provide more realistic information for managing SOC under changing climate conditions (Lauer et al., 2017).

Currently, ESMs predict SOC stocks primarily with model representations that depend on soil temperature, moisture, and belowground net primary production (Todd-Brown et al., 2013). ESMs capture the positive correlation between net primary productivity (NPP) and precipitation, leading to high SOC stocks in regions with high NPP and sufficient moisture (Sun et al., 2016). Higher temperatures increase soil respiration, which, in the short term, reduces SOC storage. Over the longer term, increased soil respiration can release nutrients, promoting plant growth and increasing belowground carbon inputs, thereby contributing to SOC stocks; however, the balance of these effects can take centuries to manifest (Mekonnen et al., 2022). In ESMs, the temperature sensitivity of soil respiration is often defined using the Q₁₀ or Arrhenius equations (Wynn et al., 2006), although modifications for extreme temperatures are likely needed (Jiang et al., 2013; Azizi-Rad et al., 2022).

Table 1 Descriptive statistics of global soil organic carbon stocks at the 0–100 cm depth interval.

TS represents tropical and subtropical.

Among environmental factors, soil moisture plays a crucial role in plant growth, microbial activity, carbon inputs, and the decomposition of litter and SOC. Global soil carbon stocks correlate with mean annual precipitation, emphasizing the significance of water availability in SOC dynamics. The relationship between soil moisture and microbial activity follows a curve, reaching a maximum at the optimal moisture content. Variations in soil moisture can either inhibit or enhance microbial activity, impacting SOC decomposition rates and carbon cycling (Moyano et al., 2013, 2018; Wieder et al., 2018; Davidson et al., 2012). This nonlinear soil moisture function is crucial for predicting SOC turnover, although its exact form varies among models (Sierra et al., 2015). Different measures of soil moisture are vital for understanding water availability across different scales, serving as indicators for soil–water relationships and ecosystem functioning. Previous studies suggest various functional forms, such as linear, quadratic, or asymptotic, to capture the impact of moisture on soil microbial activity, with relative water saturation proving to be a reliable predictor across diverse soil types (Moyano et al., 2013). The temperature function in soil carbon models represents the sensitivity of SOC decomposition to temperature, as well as the availability of soluble substrates that drive carbon substrate decomposition (Davidson et al., 2012). Based on the Q₁₀ equation, a 10 °C temperature increase roughly doubles the rate of soil respiration, reflecting increased microbial activity that accelerates organic matter decomposition and increases CO₂ emissions from the soil. Recent research emphasizes that the variability in the temperature sensitivity of SOC decomposition is influenced by the microbial community composition. A comprehensive understanding of the temperature function requires accounting for microbial community dynamics in soil carbon models, as these dynamics are intricately linked to temperature. This includes both accelerated microbial activity and faster O₂ depletion, which together affect soil oxygen dynamics. Despite this complexity, the empirical relationship between soil respiration and temperature, as represented by the Q₁₀ relationship, remains essential for predicting the impacts of temperature change on soil carbon dynamics and understanding its global implications for carbon cycling (Lloyd and Taylor, 1994).

In a previous US continental-scale study, we derived empirical nonlinear relationships between SOC and environmental factors that achieved comparable prediction accuracy to a random forest (RF) machine learning approach (Mishra et al., 2022). In this study, we apply a similar approach using both global field observations and ESMs to (1) identify key observed environmental controllers of global SOC stocks and their functional relationships and (2) evaluate ESMs using these observational benchmarks. Simulated SOC stocks from three ESMs used in Phase 6 of the Coupled Model Intercomparison Project (CMIP6) – the Community Earth System Model (CESM; Hurrell et al., 2013), the UK Earth System Model (UKESM; Sellar et al., 2019), and the Beijing Climate Center model (BCC; Xiao-Ge et al., 2019) – were benchmarked with over 50 000 SOC profile observations across the globe. We used a machine learning (i.e., random forest) approach with 46 environmental factors to identify the key environmental controllers of SOC stocks at the global scale. We then applied a generalized additive model (GAM) to derive the predictive relationships between these key environmental factors and SOC stocks in both observations and ESM-simulated SOC stocks. The specific objectives of this study were to (1) identify dominant environmental controllers of SOC stocks in field observations and CMIP6 ESM simulations, (2) derive the observed and ESM-modeled functional relationships between environmental factors and SOC stocks, and (3) analyze these functional relationships to inform needed improvements in ESM representations of SOC dynamics.

2 Materials and methods

2.1 Soil organic carbon stock observations

We used two datasets of SOC stocks: one for the upper 30 cm of soil (i.e., 0–30 cm) and one for the upper meter of soil (i.e., 0–100 cm). We note that, by limiting our analysis to these depth intervals, we may not be accounting for the total SOC stocks, as large SOC stocks can be found at much deeper depths in some soils (e.g., peatlands). The World Soil Information Service (WoSIS) compiled SOC profiles across the globe after quality assessment. The 2019 snapshot of the WoSIS dataset contained 111 380 soil profiles with SOC content information (given in grams of carbon per gram of soil) at different soil depths (Batjes et al., 2020). We estimated the SOC stock (g C m⁻²) at different soil layers using the following expression:

$$\begin{array}{}\text{(1)}& \mathrm{SOC}\mathrm{stock}=\mathrm{SOC}\mathrm{content}\times \left(\mathrm{1}-{\displaystyle \frac{G}{\mathrm{100}}}\right)\times \mathrm{BD}\times D,\end{array}$$

where G is the coarse fragment fraction (%), BD is the bulk density of soil (g m⁻³), and D is the soil layer depth (m).

When the measured bulk density value was absent from the dataset, we used a pedotransfer function (Yigini et al., 2018) to estimate the soil bulk density:

$$\begin{array}{}\text{(2)}& \mathrm{BD}=\mathit{\alpha }+\mathit{\beta }\times \exp (-\mathit{\gamma }\times \mathrm{OM}).\end{array}$$

Here, OM is the organic matter, equivalent to SOC × 1.724, with SOC content in percent (%); α, β, and γ are fitting parameters. We found α=0.32, β=1.30, and γ=0.0089 after fitting WoSIS data to this equation.

Another dataset that we used in this study was compiled from Mishra et al. (2021). This dataset contained 2546 soil profiles with SOC stock (g C m⁻³) information from permafrost regions in North America, northern Eurasia, and the Qinghai–Tibet Plateau. In total, we used 113,926 soil profile observations from these two data sources. SOC stocks of different soil layers were then summed to SOC stocks at the 0–30 cm and 0–100 cm depth intervals. Because not all of these soil profiles covered the whole 0–30 cm or 0–100 cm intervals, we used a total of 54 000 soil profiles that included SOC stock information for both depth intervals. The geographical distributions of soil profiles used in this study are shown in Fig. 1. Because SOC stock values across the globe were highly skewed, we used a natural logarithm transformation in this study.

Table 2 Prediction accuracies of random forest models across biomes and at the global scale with respect to predicting SOC stocks.

TS represents tropical and subtropical. RMSE is the root-mean-square error.

2.2 Environmental predictors of SOC stocks

The storage and cycling of SOC are controlled by multiple environmental factors. In this study, we used observations of 46 environmental variables, which represented major soil-forming factors (McBratney et al., 2003). A total of 21 of the 46 environmental variables were climatic variables, including the annual average temperature, precipitation, evapotranspiration, drought severity index, and statistics for different temporal scales (e.g., during the wettest and driest quarter in a year). A total of 13 of the 46 variables described soil properties (e.g., clay content, sand content, silt content, soil texture, pH, and cation exchange capacity). Six variables represented topographic factors (e.g., elevation and soil depth). Six variables represented land use and land cover types. All of the categorical variables were converted to integer variables, and the environmental variables were resampled to a common 1 km resolution. The environmental factors, their original spatial resolution, and data sources are provided Table S1 in the Supplement.

2.3 Selection of dominant environmental controllers of SOC stocks

We used RF to select dominant environmental predictors of SOC stocks within biomes and at the global scale in both observations and ESMs. RF is an ensemble learning method that is an extension of the classical classification and regression trees (CART) approach. Building a collection of uncorrelated CARTs through bootstrapping the samples and applying the random subspace method at each branch of the trees, RF improves prediction performance (Breiman, 2001; Wiesmeier et al., 2011; Mishra et al., 2020). RF is well known for its strength in modeling highly nonlinear relationships between the predictors and is robust to overfitting (Chagas et al., 2016). Moreover, RF is not very sensitive to the choice of the hyperparameters, making it one of the most popular off-the-shelf models for many classification and regression problems.

In this study, we trained the RF model using SOC content as a response variable and environmental factors as predictors. The model performance was evaluated using the coefficient of determination (R²) and root-mean-square error (RMSE). A 10-fold cross-validation was used to compute the R² and RMSE values. Biome-specific analyses were conducted on a subset of the global dataset. For biome classification, we used the International Geosphere–Biosphere Programme (IGBP) land classes (Loveland and Belward, 1997). The “randomForest” package in R was used to train an RF model using all of the observed environmental factors in the dataset and to identify dominant environmental controllers of SOC stocks. Prior to fitting into the final model, we performed a potential collinearity test among the environmental variables by calculating pairwise correlations and the variance influence factor (VIF) values. Predictors showing a VIF value greater than 10 were omitted, leaving 14 uncorrelated environmental predictors of SOC stocks in the observations.

Figure 1 Spatial and statistical distributions of 54 000 soil organic carbon profiles used in this study.

Figure 2 Box plot of the soil organic carbon (SOC) stock for each biome analyzed in this study. The horizontal line in the middle of the boxes is the median, while their lower and upper limits correspond to the respective first and third quartiles. TS represents tropical and subtropical.

2.4 Generalized additive model

The generalized additive model (GAM) approach is an extension of generalized linear models that employs spline functions to model nonlinear relationships between predictor and response variables (Arnold et al., 2013). In GAMs, the relationship between the predictor and response variable can be modeled as follows (Hastie and Tibshirani, 1987):

$$\begin{array}{}\text{(3)}& Y=C+{\sum }_{i=\mathrm{1}}^p{f}_i\left(X_i\right).\end{array}$$

Here, Y is the response variable (SOC), C is a constant, X_i represents the environmental controller variables, f_i is a spline function for X_i, and p is the total number of environmental controllers. We used the “mgcv” package in R to build nonlinear relationships between environmental factors and SOC stocks using GAMs for both the SOC field observations and CMIP6 ESM-simulated SOC data (Arnold et al., 2013). The performance of GAMs was evaluated using the R² and RMSE.

2.5 Earth system model outputs

We downloaded and aggregated the ESM-simulated SOC and environmental controller data from three ESMs that participated in CMIP6: the Community Earth System Model (Hurrell et al., 2013), the UK Earth System Model (Sellar et al., 2019), and the Beijing Climate Center model (Xiao-Ge et al., 2019). These ESMs included most of the environmental factors used by CMIP6 ESMs. ESMs did not report depth-dependent soil carbon projections, making direct comparison with depth-dependent SOC observations difficult. The majority of land models used in the ESMs were designed to simulate topsoil carbon for topsoil depth; thus, we assumed that the simulated soil carbon is contained within the top 1 m of the soil profile to simplify comparison with observations.

3 Results

3.1 Descriptive statistics of SOC observations

The average global SOC stock within the 0–1 m depth interval was 13.5 kg C m⁻², ranging from 0.14 to 435.3 kg C m⁻². Our results indicate substantial variability in global-scale SOC observations, as the standard deviation (18.2 kg C m⁻²) is greater than the average SOC stock. Summary statistics of SOC stocks at the global scale and within different biomes are presented in Table 1. Boreal forests and temperate forests exhibited higher SOC stocks compared with other biomes, whereas tundra and the tropical and subtropical broadleaf forests showed lower and relatively similar average SOC stocks. Tundra and the tropical and subtropical grasslands and savannas exhibited similar and lower standard deviations with respect to SOC stock values. Conversely, boreal forests and shrublands showed higher standard deviations, indicating a broader range in SOC stock values. The distributions of total SOC stocks across different biomes are presented in Fig. 2.

3.2 Dominant environmental controllers of SOC stocks in observations and ESMs

At the global scale, we found that the diurnal temperature, drought severity index, annual temperature, and cation exchange capacity are the dominant environmental controllers of SOC stocks in the observations (Fig. 3). By including all environmental controllers, the RF model explained 61 % of observed global spatial variation in SOC stocks. The R² values ranged from 43 % in boreal forests to 62 % in croplands (Table 2), with the importance of key environmental controllers varying across biomes (Fig. 4). In croplands, precipitation, drought, diurnal temperature, and cation exchange capacity were identified as the dominant controllers of SOC stocks. In grasslands, annual temperature, cation exchange capacity, and sand content were the dominant controllers. In forests, the cation exchange capacity, precipitation, and temperature were the dominant controllers, while in shrublands, annual temperature, soil pH, and cation exchange capacity were the most important controllers. In savannas, soil-related variables, temperature, and precipitation were the most important controllers. Across all land cover types, we found that the cation exchange capacity and seasonal climatic variables were the dominant environmental controllers of SOC stocks.

In contrast, the RF model with eight environmental variable predictors made near-perfect predictions of ESM-simulated SOC stocks (average R²=0.95; the R² values for the UKESM, CESM, and BCC models were 0.99, 0.89, and 0.98, respectively). In contrast to the results obtained from the observed SOC stocks, the dominant controllers of ESM-simulated SOC stocks were annual temperature, NPP, and annual precipitation (Fig. 5), with NPP being by far the most dominant predictor of SOC stocks in the UKESM.

Figure 3 Importance of different environmental factors for predicting the global SOC stocks in observations.

Figure 4 Strength and importance of environmental factors for predicting observed SOC stocks within different biomes.

Figure 5 Importance of different environmental factors for predicting global soil organic carbon stocks in three CMIP6 Earth system models.

3.3 Predictive relationships between environmental factors and SOC stocks

The dominant environmental controllers of observed SOC stocks, as identified by the RF model, were used in a GAM to derive predictive relationships. We retrieved explicit analytical expressions by fitting splines derived from the GAM to the observation dataset. Despite its role as the sole carbon source for soil, our results did not show NPP as a strong controller of observed SOC stocks (Fig. 6a). In contrast to field observations, all ESMs showed a significant dependence (exponential increase) of SOC stocks on NPP. Our results also showed that observed SOC stocks increased almost linearly with observed annual precipitation (Fig. 6b). In contrast, ESMs showed different relationships between SOC and precipitation. We found a nonlinearly increasing SOC with precipitation in the CESM, an initial sharply increasing and then decreasing relationship in the UKESM, and a decreasing relationship in the BCC ESM. Regarding the relationship of SOC storage with soil texture and elevation, ESMs did not capture the observed relationships. Our results indicated that observed SOC stocks decreased with clay content in the interval between 0 % and 20 %, whereas they increased with clay content above 20 % (Fig. 6c). Similarly, observed SOC stocks increased with silt content up to 55 % and decreased thereafter (Fig. 6d).

Figure 6 Predictive relationships between environmental factors and SOC stocks in observations (black line) and CMIP6 Earth system models (different colors; see legend). Red circles are computed from fitted curves. The shading around the solid line indicates the 95 % confidence interval.

SOC stock functional relationships differed between the three ESMs, and, in many cases, they also differed from the relationships that we derived from observations. In terms of the effects of annual temperature on modeled SOC storage, we found that SOC stocks decreased with increasing annual temperature and were most sensitive to temperature within the range of 0 and 10 °C (Fig. 6e). However, while the three ESMs captured the general negative relationship between SOC storage and temperature, none of them correctly described the varying sensitivity of SOC stocks across different temperature ranges, especially in extreme temperature ranges (< 0 and > 20 °C). With respect to representing the control of elevation on SOC storage, only the UKESM showed patterns consistent with observations: SOC storage remained stable at elevations below 2000 m and decreased at elevations above 2000 m (Fig. 6f).

4 Discussion

Previous studies have suggested that the spatial variation in SOC depends on multiple environmental factors, such as climatic and edaphic variables, geography, and vegetation. In this study, we found that climatic variables (i.e., temperature and precipitation) are the most important controllers of global SOC stocks, followed by edaphic variables (i.e., cation exchange capacity), topography (i.e., elevation), and vegetation (i.e., NPP). Using boosted regression trees, Luo et al. (2021) studied edaphic and climatic controls on SOC dynamics at different soil depths and found that soil type and climatic variables are the most important variables for explaining SOC stocks (Luo et al., 2021). In this study, we found that seasonal climatic variables, such as diurnal temperature range and precipitation seasonality, are among the most important environmental controllers of the spatial variation in SOC stocks. This result indicates the critical role of seasonal and interannual climatic variables in understanding SOC dynamics.

The importance of climatic variables on global SOC storage emerges from their close links with processes that affect ecosystem productivity and soil microbial activity. Consistent with our findings, Wiesmeier et al. (2014) reported climatic variables (temperature and precipitation) as significant controllers of SOC stocks up to 1 m depth in German soils under oceanic climate conditions. Similarly, Sreenivas et al. (2014) used RF to predict SOC variability across semiarid and humid areas of India in the top 30 cm of soil and found that the top three environmental controllers were land cover, the mean temperature of the hottest months, and mean annual precipitation. In our analysis, the overall relative importance of climatic variables was significantly higher than other variables at both the global and biome scales.

Soil properties were identified as the second most important controllers of global SOC stocks. Soil properties influence various processes that govern soil carbon dynamics. For example, soil properties impact microbial activity, porosity, and oxygen availability in the soil profile, which directly or indirectly control soil water dynamics, plant growth, and SOC stocks. Consistent with our findings, Luo et al. (2021) reported that sand content, silt content, and soil pH were significant controllers of SOC stocks in all soil depths globally.

The Palmer drought index used in this study explicitly includes factors such as evapotranspiration, precipitation, and temperature. This index accounts for controls on the soil water balance, making it a better predictor than precipitation alone. Drought severity and duration play crucial roles in influencing the extent of soil carbon losses through microbial respiration (Borken and Matzner, 2009). A decline in soil CO₂ efflux is observed as both the quantity and frequency of precipitation events decrease (Harper et al., 2005). In the initial phases of drought, elevated soil CO₂ emissions occur due to the rapid response of plants and microorganisms to environmental stress (Ru et al., 2018). However, as drought intensifies, overall CO₂ emissions diminish due to reduced root growth and microbial CO₂ efflux caused by increasing soil dryness (Hasibeder et al., 2015). Similar to the impacts of drought duration, drought intensification results in a decrease in total CO₂ emission by suppressing soil microbial activity and associated soil CO₂ fluxes (Harper et al., 2005; Hu et al., 2020).

Cation exchange capacity is a soil property that indicates the active soil surface to which SOC can be adsorbed, with polyvalent metal cations playing a significant role in SOC stabilization by binding organic compounds to mineral surfaces (O'Brien et al., 2015; Solly et al., 2020). O'Brien et al. (2015) found that exchangeable soil Ca²⁺ is a significant predictor of SOC stocks. This relationship is supported by the mechanism in which Ca²⁺ and Mg²⁺ promote clay flocculation and bind organic matter to clay surfaces. Solly et al. (2020) reported that SOC and cation exchange capacity are significantly related in both topsoil and subsoil with a strong positive relationship.

After climatic factors and cation exchange capacity, topography and vegetation (NPP) were important controllers of observed global SOC stocks. The effect of NPP on observed SOC stocks was found to be relatively small (∼ 6 % in soils at 0–100 cm depth). Similar to our findings, Luo et al. (2021) reported that NPP explained about 10 % of the variation in SOC stocks. NPP provides the primary inputs of carbon to the soil, and NPP generally increases with moisture, temperature, and CO₂ levels, up to a certain limit. NPP also depends on the availability of soil nutrients. Most ESMs tend to overestimate the increase in SOC pools in response to NPP increases (Todd-Brown et al., 2013). The effect of NPP on SOC also varies depending on biome type and soil depth (Luo et al., 2021; Georgiou et al., 2021). The contribution of NPP to SOC stocks mostly depends on how much NPP ends up in the soil and how it is translocated to different soil depths. Georgiou et al. (2021) reported a saturating relationship of SOC stocks with increasing NPP in a global observational dataset. However, Chen et al. (2018) reported high SOC stocks associated with increasing productivity and soil water-holding capacity (Chen et al., 2018).

The three CMIP6 ESMs that we analyzed predicted SOC stocks mostly as a function of temperature, precipitation, and NPP. These ESMs simulated positive correlations between SOC stocks and NPP (Fig. 6a), resulting in higher SOC stocks in areas with a high NPP (Shi et al., 2013; Sun et al., 2016). In these ESMs, the effects of temperature and precipitation on SOC stocks are driven by soil respiration. Most current ESMs simulate the response of soil respiration to temperature using either the Q₁₀ or Arrhenius equations (Wynn et al., 2006), where higher temperatures lead to increased soil respiration, which, all else being equal, eventually reduces SOC stocks (Fig. 6e). Our results showed diverse controls of precipitation on SOC stocks in different ESMs. Todd-Brown et al. (2013) showed that ESMs simulate soil respiration as either (1) increasing monotonically with precipitation or (2) rising to a plateau under optimal precipitation and then experiencing a subsequent decline with further increases in precipitation. Consistent with those results, the ESMs that we analyzed in this study showed differences in the functional relationship between SOC storage and precipitation (Fig. 6b).

In this study, we found that while current ESMs consider key environmental controllers, such as soil temperature and moisture, with respect to regulating SOC storage, they show large inter-model variations in representing the functional relationships between these factors and SOC at the global scale. Meanwhile, none of the three ESMs investigated in this study align with the patterns that emerged from observations. These results signify potential unrealistic parameterization or missing critical processes in model representation. Moreover, our results highlight the importance of including other environmental factors in simulating global SOC storage. The observed global SOC stocks are influenced by processes beyond those currently considered in ESMs, such as temperature, precipitation, and NPP. Effects of other environmental factors, such as the drought severity index and cation exchange capacity should also be considered in future representations of SOC dynamics in ESMs. Our results highlight the critical role of observational data in benchmarking ESM simulations and informing model structures and parameterization. Although our findings cannot be directly used to develop model parameterizations, they can (1) point to categories of functional forms for controllers, (2) guide where efforts may be best applied to improve model functional forms (e.g., for the dominant controllers), and (3) provide information for modelers regarding where their models may differ significantly from observed functional forms for emergent relationships.

We note that although ecosystem-specific environmental factors (e.g., the presence or absence of certain species, such as earthworms or termites, or indicators of anthropogenic land management, such as the use of fertilizers or conservation agriculture practices) are not represented in the current generation of ESMs, incorporating these factors – specific to ecosystems like croplands or forests – may improve the accuracy of SOC stock predictions in comparison to observations.

5 Conclusions

Our results document a disagreement between the environmental controllers of SOC stocks in observations and CMIP6 ESM simulations. Specifically, while global SOC observations indicate that the diurnal temperature, drought index, annual temperature, cation exchange capacity, and other soil-related variables are critical for controlling SOC stocks at the global scale, ESMs overemphasize the roles of NPP, annual temperature, and annual precipitation in simulating SOC stocks. Moreover, the functional relationships between key environmental factors and SOC stocks show significant uncertainty among ESMs, with no agreement with those derived from observations. Optimizing current model parameterizations and developing new model structures that incorporate more processes in the soil carbon cycle are crucial for improving the simulation of global SOC storage in future ESMs. Our results highlight the importance of benchmarking ESMs against observations to improve the mechanistic understanding of the global soil carbon cycle.