Modeling biochar effects on soil organic carbon on croplands in a microbial decomposition model (MIMICS-BC_v1.0)

Biochar (BC) application to croplands aims to sequester carbon and improve soil quality, but its impact on soil organic carbon (SOC) dynamics is not represented in most land models used for assessing land-based climate change mitigation; therefore, we are unable to quantify the effects of biochar application under different climate or land management conditions. Here, to fill this gap, we implement a submodel to represent biochar in a microbial decomposition model named MIMICS (MIcrobial-MIneral Carbon Stabilization). We first calibrate and validate MIMICS with new representations of the density-dependent microbial turnover rate, adsorption of available organic carbon on mineral soil particles, and soil moisture effects on decomposition using global field-measured cropland SOC at 285 sites. We further integrate biochar in MIMICS by accounting for its effect on microbial decomposition and SOC sorption/desorption and optimize two biochar-related parameters in these processes using 134 paired SOC measurements with and without biochar addition. The MIMICS-biochar version can generally reproduce the short-term (≤ 6 years) and long-term (8 years) SOC changes after adding (mean addition rate of 25.6 t ha⁻¹) biochar (R²= 0.79 and 0.97, respectively) with a low root-mean-square error (RMSE = 3.73 and 6.08 g kg⁻¹, respectively). Our study incorporates sorption and soil moisture processes into MIMICS and extends its capacity to simulate biochar decomposition, providing a useful tool to couple with dynamic land models to evaluate the effectiveness of biochar application with respect to removing CO₂ from the atmosphere.

1 Introduction

Soil organic carbon (SOC) is the largest terrestrial carbon pool, and increasing soil respiration in response to global warming can cause large carbon emissions to the atmosphere (Bond-Lamberty et al., 2018). On the other hand, SOC sequestration through improved land management practices has the potential to mitigate climate change by increasing soil carbon accumulation, as is the focus of initiatives such as the “4 per 1000” project (Minasny et al., 2017). Due to the limited temporal and spatial coverage of field SOC measurements, soil biogeochemical models have been widely applied to simulate SOC and its response to climate change and human activities (Eglin et al., 2010). Soil carbon models are evolving from first-order kinetics-based models with a simple representation of pool sizes and their turnover rates to microbial models with explicit representation of microbial roles in SOC decomposition and stabilization (Manzoni and Porporato, 2009; Sulman et al., 2018). For example, the MIcrobial-MIneral Carbon Stabilization (MIMICS) model is a process-based soil carbon model with explicit representations of nonlinear SOC decomposition dynamics related to microbial physiology, substrate quality, and the physical protection of SOC (Wieder et al., 2014, 2015). This model has been calibrated with global SOC data, can represent the current understanding of SOC decomposition and formation well (Wieder et al., 2015), and outperforms conventional first-order decomposition models with respect to simulating the spatial variation in SOC stocks in forest ecosystems on the continental scale (Zhang et al., 2020). However, the model has not been evaluated for agricultural sites or misses processes that theoretically should influence SOC dynamics.

Microbial interactions at the community level (e.g., competition) play a crucial role in controlling SOC dynamics, but they are usually omitted in microbial models (Georgiou et al., 2017), resulting in the unrestricted growth of the microbial community with increasing carbon input, which is unrealistic (Buchkowski et al., 2017; Wieder et al., 2013). In addition, field experiments show that physicochemical adsorption plays a more important role in controlling dissolved organic carbon (DOC) fluxes than the biodegradation process (Kalbitz et al., 2005). Although the adsorption mechanism is complex, depending on various factors such as pH, clay content, and the destruction and formation of soil aggregates (Mayes et al., 2012), some soil carbon models have implemented dynamic adsorption and desorption processes controlled by the DOC concentration and available mineral surface sites for binding (G. Wang et al., 2020, 2013). The availability of SOC is influenced by the adsorption process (Michalzik et al., 2003). Some adsorption kinetic equations, such as the Langmuir isotherm, have commonly been employed to depict the adsorption/desorption process. However, the MIMICS model lacks consideration of the adsorption process, thus not effectively elucidating its role in stabilizing SOC. Furthermore, the effect of soil moisture on SOC cannot be ignored, as soil moisture controls microbial activity, substrate availability, and further influences soil respiration and nitrogen mineralization (Manzoni et al., 2012; Schimel et al., 2007). A set of empirical functions for the soil moisture effects have been proposed for use in Earth system models (ESMs) (Moyano et al., 2013; Camino-Serrano et al., 2018), and a mechanistic moisture function that incorporates physicochemical and biological processes has also recently been developed (Yan et al., 2018). In previous MIMICS versions, an implicit or explicit density-dependent microbial turnover has been introduced (Wieder et al., 2015; Kyker-Snowman et al., 2020; Zhang et al., 2020; Georgiou et al., 2017), which causes an increase in microbial biomass turnover with increasing microbial community size, reflecting increasing pressure from competition for resources other than carbon (e.g., space) and from viral infections (Jansson and Wu, 2023), and a water scalar has been used to represent the soil moisture effects (Wieder et al., 2019). The inclusion of a density-dependent microbial turnover rate improved the accuracy of SOC prediction at the global scale, compared with MIMICS without this inclusion, and eliminated the correlation between simulated biases and the input of annual litterfall (Zhang et al., 2020). MIMICS with soil water modifications showed comparable predicted global soil carbon stocks compared to other models, but the extent to which soil water influences SOC turnover remains uncertain (Wieder et al., 2019). Therefore, based on these theories and model limitations, it is necessary to integrate the three aspects (density-dependent microbial turnover rate, adsorption processes, and soil moisture impacts) into one model version to improve the prediction accuracy of SOC dynamics. For agricultural lands, modeling the SOC decomposition processes is more challenging due to management practices, such as tillage and fertilization, that can significantly interrupt the carbon cycle and require specific parameterizations.

Biochar application to croplands as a soil amendment can improve the soil quality and increase crop production (Smith, 2016; Woolf et al., 2010). Moreover, as biochar is produced from biomass through pyrolysis processes and is recalcitrant to decomposition, it is also considered to be a promising negative-emission technology (NET) with respect to climate change mitigation (Fuss et al., 2018; Minx et al., 2018). The carbon dioxide removal (CDR) potential of biochar is estimated to be 0.5–2 Gt CO₂e yr⁻¹ (where CO₂e denotes CO₂ equivalent) (Fuss et al., 2018; Minx et al., 2018). However, biochar application affects SOC mineralization through various processes (Palansooriya et al., 2019; Luo et al., 2017), resulting in positive or negative priming effects (PEs, changes in native SOC mineralization) (Zimmerman et al., 2011). A recent meta-analysis showed that biochar induced negative PEs on average (−3.8 %), but the 95 % confidence interval (CI) of −8.1 % to 0.8 % also covered positive values (Wang et al., 2016). Biochar may induce positive PEs via the stimulation of microbial activity by providing additional nutrients for soil microbes (El-Naggar et al., 2019; Li et al., 2019). Positive PEs usually occurred over the shorter term (< 2 years) and then decreased or became negative over the longer term (Luo et al., 2011; Singh and Cowie, 2014; Ding et al., 2017). For example, biochar can reduce the SOC available for microbes by enhancing soil aggregate stability via associations between soil minerals and biochar (Zheng et al., 2018). Its porous structure and high surface area with a strong adsorption affinity for SOC can, thus, cause negative PEs (Zimmerman et al., 2011; Lehmann et al., 2021). PEs are also impacted by the properties of biochar (e.g., feedstock type and pyrolysis temperature) and soil climate (e.g., soil moisture) (Ding et al., 2017). Therefore, soil moisture could be closely related to the adsorption capacity of biochar and, thus, needs to be included in the model for predicting the PEs of biochar on SOC changes. The biochar decomposition and impacts on native SOC via PEs are important for the CDR potential of biochar, but these processes are not represented in most land carbon models (Lehmann et al., 2021), precluding the model capacity to fully assess the effectiveness of the large-scale application of biochar as a NET as well as its environmental impacts.

2 Materials and methods

2.1 Modifications of the MIMICS model

2.1.1 The default version of MIMICS (MIMICS-def)

There are seven carbon pools in MIMICS: two litter pools, two microbial biomass pools, and three SOC pools (Fig. 1). The litter inputs (LIT) are divided into metabolic (LIT_m) and structural (LIT_s) pools according to the litter quality (f_met, i.e., fraction of litter to LIT_m), which is linearly related to the lignin-to-nitrogen ratio (lignin:N; Table S1 in the Supplement). Microbial growth efficiency (MGE) determines the carbon fluxes from the two litter pools and the available SOC pool (SOC_a) for microbial biomass pools and heterotrophic respiration. The turnover of microbial biomass (τ) depends on the functional types of soil microbes (MIC_r and MIC_k for the r- and k-strategy, respectively). Three SOC pools represent the available (SOC_a), physically protected (SOC_p), and chemically recalcitrant SOC (SOC_c). SOC in the protected pools (i.e., SOC_p and SOC_c) is released to the available SOC pool (SOC_a) over time. A more detailed description of the model parameters and carbon fluxes can be found in Table S1 in the Supplement and in Wieder et al. (2015). The carbon decomposition rate (mg C cm⁻³ h⁻¹) of the litter and SOC pools is based on temperature-sensitive Michaelis–Menten kinetics (Allison et al., 2010; Schimel and Weintraub, 2003):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{s}}}{\mathrm{d}t}}=\mathrm{MIC}\times {\displaystyle \frac{V_{max}\times C_{\mathrm{s}}}{K_{\mathrm{m}}+C_{\mathrm{s}}}},\end{array}$$

where C_s (mg C cm⁻³) is the size of a substrate carbon pool (LIT or SOC) and MIC (mg C cm⁻³) is the size of the microbial carbon pool (MIC_r or MIC_k). V_max and K_m are the microbial maximum reaction velocity (mg C (mg MIC)⁻¹ h⁻¹) and the half-saturation constant (mg C cm⁻³), respectively, which depend on temperature, T (in °C).

$$\begin{array}{}\text{(2)}& {\displaystyle }{V}_{\max }=e^{V_{\mathrm{slope}}\times T+V_{\mathrm{int}}}\times a_v\times V_{\mathrm{mod}}\text{(3)}& {\displaystyle }{K}_{\mathrm{m}}=e^{K_{\mathrm{slope}}\times T+K_{\mathrm{int}}}\times a_k\times K_{\mathrm{mod}}\end{array}$$

Here, V_mod and K_mod represent the modifications of V_max and K_m, respectively, based on their dependence on litter quality, microbial functional type, and soil texture; a_v and a_k are the tuning coefficients of V_max and K_m, respectively; V_slope and K_slope are the regression coefficients; and V_int and K_int are the regression intercepts.

The turnover of MIC_r and MIC_k (MIC_τ, mg C cm⁻³ h⁻¹) at each time step depends on their specific turnover rate (k_mic, h⁻¹), annual total litter input (LIT_tot, g C m⁻² yr⁻¹), and f_met:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathrm{MIC}}_{\mathit{\tau }}& =a_{\mathit{\tau }}\times k_{\mathrm{mic}}\times e^{c\times f_{\mathrm{met}}}\\ & \times max(min\left(\sqrt{{\mathrm{LIT}}_{\mathrm{tot}}},\mathrm{1.2}\right),\mathrm{0.8})\times \mathrm{MIC},\end{array}\end{array}$$

where a_τ (=1.0, dimensionless) is the tuning coefficient of k_mic and c is the regression coefficient of MIC_r (0.3) and MIC_k (0.1). The carbon inputs from microbial biomass to SOC pools are determined by the microbial biomass turnover.

The carbon transfer from SOC_p to SOC_a (D, mg C cm⁻³ h⁻¹) represents the desorption of SOC_p from mineral surfaces or the breakdown of aggregates, calculated as a function of soil clay content (f_clay):

$$\begin{array}{}\text{(5)}& D=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}\times k_{\mathrm{d}}\times e^{-\mathrm{1.5}\times f_{\mathrm{clay}}},\end{array}$$

where k_d (=1.0, dimensionless) is a tuning coefficient of the desorption rate. The parameter values of the default MIMICS version can be found in Table S1.

2.1.2 MIMICS considering the density-dependent microbial turnover rate (MIMICS-T)

Similar to the logistic growth model in population ecology, various regulatory mechanisms (e.g., competition and viruses) can limit the microbial population size (Buchkowski et al., 2017; Jansson and Wu, 2023). The absence of restrictions other than carbon on population size results in predictions of microbial biomass that increase indefinitely with carbon input. Consequently, the response of predicted SOC to changes in carbon inputs is close to zero, which contradicts field observations (Georgiou et al., 2017). A density-dependent microbial turnover rate with β>1 was adopted to regulate the responses of soil microbial biomass to external environmental variations, such as carbon input, and thus SOC dynamics in previous microbial models (Georgiou et al., 2017; Zhang et al., 2020). We incorporated the density-dependent microbial turnover rate into MIMICS following Zhang et al. (2020). In the MIMICS-T version, we modified Eq. (4) to represent the increased microbial turnover rate with growing microbial biomass density (MIC, mg C cm⁻³):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathrm{MIC}}_{\mathit{\tau }}& =a_{\mathit{\tau }}\times k_{\mathrm{mic}}\times e^{c\times f_{\mathrm{met}}}\times \\ & max(min\left(\sqrt{{\mathrm{LIT}}_{\mathrm{tot}}},\mathrm{1.2}\right),\mathrm{0.8})\times {\mathrm{MIC}}^{\mathit{\beta }},\end{array}\end{array}$$

where β is the density-dependent exponent.

2.1.3 MIMICS-T with the additional representation of sorption (MIMICS-TS)

Although the MIMICS model can simulate the desorption process (yellow arrow from SOC_p to SOC_a in Fig. 1), the adsorption process is still missing. In the original version of MIMICS, fixed fractions of litter and microbial turnover are transferred to the physically protected SOC pool (SOC_p in Fig. 1); the SOC_p is then unprotected from mineral surfaces or the breakdown of aggregates using a desorption rate which is a function of the clay fraction. Therefore, we do not think that the original MIMICS model actually simulates sorption as a process, as sorption is dependent on the substrate concentration; therefore, the sorption rate should vary with the dissolved organic carbon concentration, rather than being proportional to the microbial carbon turnover rate as assumed in the original MIMICS. Thus, we further added the adsorption of available SOC to MIMICS, following Wang et al. (2013) and Mayes et al. (2012). The MIMICS-TS version includes a new sorption process (purple arrow from SOC_a to SOC_p in Fig. 1), but the original desorption process (i.e., yellow arrow from SOC_p to SOC_a in Fig. 1) remains unchanged. The sorption capacity of SOC_a (Q_max) increases with increasing clay content, and the carbon flux of the sorption process is calculated as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }{F}_{\mathrm{ads}}=K_{\mathrm{ads}}\times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{SOC}}_{\mathrm{p}}}{Q_{max}}}\right)\times {\mathrm{SOC}}_{\mathrm{a}},\text{(8)}& {\displaystyle }{K}_{\mathrm{ads}}=k_{\mathrm{d}}\times k_{\mathrm{ba}},\text{(9)}& {\displaystyle }{Q}_{max}={\mathrm{10}}^{(c_{\mathrm{1}}\times \log (\mathit{\%}\mathrm{clay})+c_{\mathrm{2}})},\end{array}$$

where F_ads is the carbon flux from SOC_a to SOC_p (mg C cm⁻³ h⁻¹); k_ba is the binding affinity; K_ads is the sorption rate of SOC_p, which is associated with the desorption rate (k_d); Q_max is the maximum sorption capacity of SOC_p (mg C cm⁻³ soil); and c₁ and c₂ are the coefficients for estimating Q_max from Mayes et al. (2012).

2.1.4 MIMICS-TS with soil moisture effects (MIMICS-TSM)

Finally, based on MIMICS-TS, we added soil moisture effects on decomposition into MIMICS. We tested two empirical functions for soil moisture used in the CENTURY model (Eq. 10; Parton et al., 2000) and the ORCHIDEE-SOM model (Eq. 11; Camino-Serrano et al., 2018), respectively. We also attempted to implement one mechanism-based function that captures the main physicochemical and biological processes of soil moisture in regulating soil respiration from Yan et al. (2018) (Eq. 12). The three functions of soil moisture are illustrated in Fig. S1.

$$\begin{array}{}\text{(10)}& {\displaystyle }{f}_{\mathrm{m}\mathrm{1}}\left(w\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+p_{\mathrm{1}}\times e^{(p_{\mathrm{2}}\times w)}}},\text{(11)}& {\displaystyle }{f}_{\mathrm{m}\mathrm{2}}\left(\mathit{\theta }\right)=max(\mathrm{0.25},min\left(\mathrm{1},k_{\mathrm{1}}\times {\mathit{\theta }}^{\mathrm{2}}+k_{\mathrm{2}}\times \mathit{\theta }+k_{\mathrm{3}}\right),\text{(12)}& {\displaystyle }{f}_{\mathrm{m}\mathrm{3}}\left({\displaystyle \frac{\mathit{\theta }}{\mathit{\phi }}}\right)=\left\{\begin{array}{l}\frac{K_{\mathit{\theta }}+{\mathit{\theta }}_{\mathrm{op}}}{K_{\mathit{\theta }}+\mathit{\theta }}\times (\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{op}}}{)}^{(\mathrm{1}+a\times n_s)},\mathit{\theta }<{\mathit{\theta }}_{\mathrm{op}}\\ (\frac{\mathit{\phi }-\mathit{\theta }}{\mathit{\phi }-{\mathit{\theta }}_{\mathrm{op}}}{)}^b,\mathit{\theta }\ge {\mathit{\theta }}_{\mathrm{op}},\end{array}\right.\end{array}$$

where f_mi ($i=\mathrm{1},\mathrm{2},\mathrm{3}$, unitless value in range from 0 to 1) is the response factor to soil moisture; w is the soil moisture indicator (AI, mm mm⁻¹); p₁ and p₂ are empirical parameters of the soil moisture scalar with values of p₁=30 and $p_{\mathrm{2}}=-\mathrm{8.5}$, respectively (Parton et al., 2000); θ is soil moisture (m³ m⁻³); k₁, k₂, and k₃ are soil moisture coefficients with values of 1.1, 2.4, and 0.29, respectively (Camino-Serrano et al., 2018); φ is the soil porosity related to soil bulk density; $\mathit{\theta }/\mathit{\phi }$ is the relative water content in soil pores; θ_op is an optimum soil moisture content parameter at which the heterotrophic respiration rate peaks; K_θ is moisture constant depending on organic–mineral associations; n_s is a saturation exponent depending on soil structure and texture; and a and b are respective SOC–microbial collocation and oxygen supply restriction factors (Yan et al., 2018).

We assumed that the kinetic parameters V_max and K_m respond to soil moisture, similarly to temperature in the Michaelis–Menten equation, by affecting enzyme activity and enzyme–substrate affinity, respectively. The soil enzyme–substrate affinity was found to increase with soil moisture due to the increased diffusion and movement of substrate, but the affinity may also decrease due to decreased substrate concentrations (Zhang et al., 2009). Thus, we translated the impacts of soil moisture on the enzyme–substrate affinity to changes in K_m. In MIMICS-TSM, the effects of soil moisture on the SOC decomposition rate are represented by multiplying the response factor by V_max and K_m as follows (Eqs. 13, 14):

$$\begin{array}{}\text{(13)}& {\displaystyle }{V}_{max}=e^{V_{\mathrm{slope}}\times T+V_{\mathrm{int}}}\times a_v\times V_{\mathrm{mod}}\times f_{\mathrm{m}i},\text{(14)}& {\displaystyle }{K}_{\mathrm{m}}=e^{K_{\mathrm{slope}}\times T+K_{\mathrm{int}}}\times a_k\times K_{\mathrm{mod}}\times f_{\mathrm{m}i}.\end{array}$$

The MIMICS models with three soil moisture functions – f_m1 (Eq. 10), f_m2 (Eq. 11), and f_m3 (Eq. 12) – are referred to as MIMICS-TSM_a, MIMICS-TSM_b, and MIMICS-TSM_c, respectively. The modifications in all MIMICS versions are summarized in Table 1.

2.1.5 Adjusted parameters for cropland SOC

Crop net primary productivity (NPP) at each site was used as the litter input to soil, but different crop types (e.g., maize, rice, and wheat) were not specified in the model. The leaf, root, and stem litter values were assumed as a fixed fraction of crop NPP. The ratio of carbon to nitrogen (C:N) and the ratio of lignin to carbon (lignin:C) of the leaf, root, and stem (Table S2) were used to calculate the metabolic fraction in the total crop litter (f_met); it was calculated as the mean metabolic fractions in the leaf, root, and stem weighted by the NPP in the three parts. In order to adapt MIMICS to simulate cropland SOC, we modified C:N and lignin:C in the three parts based on field measurements of the main crop types (Table S2; Abiven et al., 2005). A harvest index (HI) of 0.45 (Hicke and Lobell, 2004) was also applied to remove the harvested part of crop and obtain the litter input to soil (calculated as the crop aboveground NPP × (1-HI)).

2.2 Implementing biochar modeling in MIMICS

When applying biochar in croplands, a fraction of biochar was assumed to be lost during application (f_loss=2 %; Archontoulis et al., 2016). Although biochar is recalcitrant to decomposition with a long turnover time (556 ± 484 years) in general, it contains a labile fraction (108 ± 196 d), and its stability varies with different biochar feedstocks, pyrolysis temperatures, and soil properties (Wang et al., 2016). Because the sizes of the SOC_p and SOC_c pools in MIMICS are not measured directly in the field studies, the 98 % remaining fraction of added biochar is partitioned into three MIMICS SOC pools by assuming that 60 % goes to SOC_p, based on the measured proportions of added biochar within aggregates (Yoo et al., 2017); 20 % goes to SOC_a, according to the labile C portion in biochar (Roberts et al., 2010); and 20 % goes to SOC_c (Fig. 1). Note that biochar is not treated as a separate carbon pool but is rather assumed to mix with other carbon in the existing pools (Fig. 1). In addition to the increase in total SOC, some important processes controlling SOC accumulation and decomposition are affected by biochar addition. Thus, we modified the parameters related to decomposition and desorption of SOC (Fig. 1). The associated rationales, equations, and parameters are described in the following sections.

Figure 1 Framework of the MIMICS model with biochar addition (MIMICS-BC; adapted from Wieder et al., 2015). The turnover of microbial biomass (τ, blue arrows) is modified with the density-dependent microbial turnover rate (MIMICS-T; Eq. 6). The adsorption process of SOC_p to SOC_a (purple arrow) is newly added and is associated with the adsorption rate (K_ads) and the maximum sorption capacity (Q_max) (MIMICS-TS; Eqs. 7–9). The carbon decomposition processes (red arrows) are modified further with three soil moisture scalers that are applied to the microbial maximum reaction velocity (V_max) and the half-saturation constant (K_m) (MIMICS-TSM_a, MIMICS-TSM_b, and MIMICS-TSM_c; Eqs. 10–12). When biochar is added to soil, the biochar (BC) carbon with an assumed fraction loss (f_loss) is partitioned into SOC_p, SOC_a, and SOC_c based on f_bp, f_ba, and f_bc, respectively (purple arrows from the BC to SOC pools). The desorption process (orange arrow from SOC_p to SOC_a) is modified through changes in the desorption rate of SOC_p (D^′) with biochar addition. The carbon decomposition processes (red arrows) are modified by adjusting the microbial maximum reaction velocity ($V_{max}^{\prime }$) with biochar addition.

The negative PEs of biochar addition on SOC may be caused by (1) the inhibition of microbial activity via changes in the soil environment due to biochar or (2) by SOC protection against microbial utilization through mineral adsorption or aggregates (Zimmerman et al., 2011). We assumed that biochar addition decreased the mineralization of native SOC (negative PE) because of its porous structure and strong adsorption affinity to organic matter (Kasozi et al., 2010), which was reported to have significantly contributed to the negative PE mechanism due to biochar addition (Zheng et al., 2018; Zimmerman et al., 2011). A desorption coefficient (f_d, ha t⁻¹ C) was defined as a function of the biochar application rate (Rate_BC) based on Woolf and Lehmann (2012) and Archontoulis et al. (2016); thus, Eq. (5) was modified as follows:

$$\begin{array}{}\text{(15)}& D^{\prime }=D\times (\mathrm{1}+f_{\mathrm{d}}\times \mathrm{Rate}\mathrm{\_}\mathrm{BC}\times \mathrm{BC}\mathrm{\_}\mathrm{C}),\end{array}$$

where D^′ (mg C cm⁻³ h⁻¹) is the new desorption rate of SOC_p with biochar addition, and a negative value of f_d indicates a negative PE. The Rate_BC is the application rate of biochar (t BC ha⁻¹) and BC_C is the carbon content in biochar (tons of carbon per ton of BC). Because the adsorption and desorption of SOC are interrelated dynamic process, modification of the desorption process with biochar addition also impacts the adsorption process. Therefore, we only modified f_d in Eq. (15) to represent the negative PE of biochar.

We also assumed that biochar stimulated microbial growth and activity through its nutrient input, inducing a positive PE to SOC (El-Naggar et al., 2019). We defined a new decomposition rate coefficient (f_v, ha t⁻¹ C), which is a function of Rate_BC, and included it in MIMICS by modifying Eq. (2):

$$\begin{array}{}\text{(16)}& V_{\max }^{\prime }=V_{max}\times (\mathrm{1}+f_v\times \mathrm{Rate}\mathrm{\_}\mathrm{BC}\times \mathrm{BC}\mathrm{\_}\mathrm{C}),\end{array}$$

where $V_{max}^{\prime }$ is the new microbial maximum reaction velocity (mg C (mg MIC)⁻¹ h⁻¹) with biochar addition.

Biochar may also have a positive PE on SOC by increasing the degradation rate of available SOC by microbes (i.e., SOC_a in MIMICS). Therefore, we added a test through modifying the V_max as a function of the biochar addition rate only in the fluxes from SOC_a to MIC_r and MIC_k, instead of in all fluxes of decomposition (Eq. 16; red arrows in Fig. 1).

2.3 Model calibration and evaluation

2.3.1 Observational data collection

We collected 387 paired field measurements of SOC concentrations (g kg⁻¹) in croplands with or without biochar (BC) addition from 58 locations (see the site map in Fig. 2) from the published literature. Soil properties (clay content, Clay; bulk density, BD; and soil moisture, SM), climatic conditions (mean annual temperature, MAT, and mean annual precipitation, MAP), a biological variable (net primary productivity, NPP), and biochar-related characteristics (application rate, Rate_BC; the interval between biochar application and soil sampling, Age_BC; feedstock type, Feedstock_BC; and pyrolysis temperature, Temp_BC) were also collected when available. Auxiliary information (e.g., location, management, and crop type) and more detailed information can be found in Han et al. (2021).

Because some sites have multiple biochar addition experiments (e.g., pyrolysis temperature × aging time of biochar), the control SOC concentrations at the same site were averaged, and the SOC concentrations with biochar addition for a given rate (Rate_BC) were also averaged, omitting other characteristics of the biochar (like pyrolysis temperature). In total, 134 paired SOC data were used for model calibration and validation (Fig. 2). The depth of soil sampled varies among sites but is less than 30 cm in general. The biochar application rate varies from 0.9 to 120 t ha⁻¹ with a median value of 20 t ha⁻¹ (Fig. S2a). Most biochar addition experiments are short-term experiments, with a median Age_BC of 1.2 years (Fig. S2b). The main types of cultivated crop are maize, rice, and wheat.

There are SOC measurements on cropland sites from 58 control treatments (no BC application) and 134 measurements from biochar treatments at the 58 sites. One control treatment may correspond to multiple biochar treatments with different applied biochar rates at a single site. Considering that the 58 site observations may be inadequate to constrain all of the new features in the revised model, we also collected SOC data on croplands (no biochar addition) from three other published global datasets (227 sites in total; Sun et al., 2020; Geisseler et al., 2017; M. Zhou et al., 2017). Therefore, a total of 285 sites were used to calibrate and evaluate the model performance for simulating cropland SOC without biochar addition (Fig. 2). The MIMICS model can run for each site; however, to be consistent with the model input resolution of daily temperature in the transient simulation (Sect. 2.3.3), the resolution of 0.5° was used for site aggregation. In detail, all sites within a given 0.5° × 0.5° grid cell were aggregated on average, and the average value was used to compare the model result in this grid cell.

The stabilizing processes of SOC after biochar addition are usually slow, and the long-term impacts of biochar addition on SOC may be different from the short-term impacts (Ding et al., 2011). Therefore, long-term SOC observations with biochar addition are needed to validate possible mechanisms and evaluate the model performance of simulating SOC stability with biochar addition. However, the duration of most biochar addition experiments is short (74.2 % of data < 3 years); thus, we extended our collected control SOC data to 8 years according to the decomposition curve of biochar in soil fitted by a double first-order exponential decay model (Fig. S4; Wang et al., 2016). Note that the double exponential decay function is only applied to the observational records of measurement data, and this function is not used in the MIMICS model. Specifically, the 8 years of SOC data with biochar addition is the sum of field control SOC observations (short term) and the residual biochar carbon in soil after 8 years. These extended long-term data were also used for model calibration and model evaluation.

Soil properties that are not reported in the literature were extracted from gridded datasets using the coordinates of the sites: clay content was taken from the Global Soil Dataset for Earth System Models (GSDE; Shangguan et al., 2014), whereas SM was obtained from observations of the Soil Moisture Active Passive (SMAP) satellite (Entekhabi et al., 2010). Missing soil BD information in control treatments was filled according to the relationship between SOC and bulk density based on 4765 cultivated soil data from the second national soil survey (Song et al., 2005), and a decrease of 7.6 % (Omondi et al., 2016) from the control soil BD was assumed to fill the missing BD values in the biochar addition experiments. The climate variable MAT was extracted from WorldClim (Fick and Hijmans, 2017), whereas the mean annual aridity index (AI, i.e., the ratio of precipitation to potential evapotranspiration) used in the soil moisture equation (Eq. 10) was obtained from the Global Aridity Index and Potential Evapotranspiration Database (Zomer et al., 2022). The biological variable (i.e., NPP) was taken from the MODIS NPP dataset (Zhao and Running, 2010).

Figure 2 Locations of the field cropland SOC measurements with or without biochar addition collected in this study and SOC measurements without biochar addition from Sun et al. (2020), Geisseler et al. (2017), and M. Zhou et al. (2017). The number of sites is also shown in the legend. Note that one site may have multiple paired SOC data due to the various experimental conditions of biochar addition at our 58 collection sites. The cropland area percentage in each 10 km × 10 km grid cell is derived from EarthStat (http://www.earthstat.org, last access: 7 June 2024; Ramankutty et al., 2008).

2.3.2 Calibration and validation for MIMICS versions without biochar

All field SOC observations in the control treatments (without biochar) from the paired measurements and SOC from the other three global datasets (Fig. 2) were assumed to be in a steady state: under present climate and continuous input of crop NPP after 45 % removal of grain with a specific crop litter quality (Sect. 2.1.5, Table S2). SOC pools in MIMICS reached an equilibrium state after about 200 years of simulation (Fig. S3). To accelerate this process, we used the Newton–Raphson method (Press et al., 2007) to obtain the steady SOC state with the site-level inputs of annual mean crop NPP, MAT, Clay, SM, and BD in the parameter optimization. This approach is constructed based on the fundamental principles governing biogeochemical cycle processes in terrestrial ecosystems (e.g., respiration and carbon distribution). A set of ordinary differential equations were built to express the dynamics of carbon flows in soil over time, and these equations can be solved numerically to obtain steady carbon pool sizes. The solver of “mnewt” is used to solve equations by iteratively calculating the values of model function “modelx” and its Jacobian matrix “modeljacx” (see the codes for details in the “Code and data availability” section). The Shuffled Complex Evolution Algorithm (SCE-UA) has been proven to be a robust method for parameter optimization (Duan et al., 1994; Muttil and Jayawardena, 2008), and the SCE-UA method from the “spotpy” package in Python (Houska et al., 2015; https://pypi.org/project/spotpy/, last access: 7 June 2024) was applied here. Parameters are optimal when the root-mean-square error (RMSE; Eq. 17) between the simulated SOC and observed SOC concentrations is minimized. The Akaike information criterion (AIC; Eq. 18; Akaike, 1974), which considers both the model error and the number the model parameters, was also calculated to evaluate different MIMICS versions.

$$\begin{array}{}\text{(17)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n({\mathrm{SOC}}_{\mathrm{obs},i}-{\mathrm{SOC}}_{\mathrm{sim},i}{)}^{\mathrm{2}}}{n}}}\text{(18)}& {\displaystyle }\mathrm{AIC}=n\times \ln \left({\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n({\mathrm{SOC}}_{\mathrm{obs},i}-{\mathrm{SOC}}_{\mathrm{sim},i}{)}^{\mathrm{2}}}{n}}\right)+\mathrm{2}p\end{array}$$

Here, SOC_obs,i and SOC_sim,i are the observed and simulated SOC at each i site, respectively; n is the number of observations; and p is the number of model parameters to be optimized.

We randomly separated 80 % of the 285 sites for the MIMICS versions' (MIMICS-def, MIMICS-T, MIMICS-TS, and MIMICS-TSM) calibration, and we used the remaining 20 % for model validation. The R², RMSE, and AIC were calculated by comparing simulated SOC with the observed SOC in training and test datasets. The parameters optimized in different MIMICS versions are shown in Table S3. Soil depth was not explicitly considered in this study, and we assumed that the soil carbon concentrations (g kg⁻¹) are similar within the top 30 cm. Note that the parameters of soil moisture functions (Eqs. 10–12) are directly derived from the original literature (Parton et al., 2000; Camino-Serrano et al., 2018; Yan et al., 2018) and not optimized in MIMICS-TSM. We calibrated the models against our datasets including SOC and auxiliary information (Fig. 2) for the main crop types (maize, rice, and wheat), and the relationships between SOC in these crop types and model input variables (i.e., NPP, MAT, and Clay) were analyzed.

To explore possible soil moisture effects on SOC, we also tried a test assuming that soil moisture affects the microbial growth rate by regulating microbial growth (V_max) and turnover (τ) of MIC_r and MIC_k (Wieder et al., 2019); thus, we added a soil moisture factor (i.e., f(θ) in Eq. 11) on V_max and τ. We also conducted a sensitivity test of MIMICS input variables (i.e., MAT, Clay, NPP, SM, and BD) with four perturbation levels of −50 %, −25 %, 25 %, and 50 % to explore the effects of possible underrepresented processes on the cropland steady SOC.

2.3.3 Calibration and validation for MIMICS versions with biochar (MIMICS-BC)

For the version of MIMICS with biochar addition, for each site, we ran simulations with a control (without biochar addition) and experimental (with biochar addition) simulation for Age_BC year at hourly time steps, restarted from the previous SOC equilibrium. Note that these simulations for biochar addition are transient runs; thus, SOC is not in a steady state. In order to meet the daily time step of transient runs required by MIMICS, the two model runs are forced by 6 h surface temperature from Climatic Research Unit and Japanese reanalysis data (CRU-JRA; Kobayashi et al., 2015; Harris et al., 2014). Therefore, the climate forcing data are different from those from WorldClim used for the steady runs (Sect. 2.3.2), as hourly climate data are not available in WorldClim. The soil-related inputs of Clay, SM, and BD were assumed to be invariant in time and consistent with input data for the steady SOC runs. The absolute SOC changes (ΔSOC, g kg⁻¹; Eq. 19) in the simulated and observed SOC concentrations were compared after BC addition. The RMSE between simulated and observed ΔSOC was minimized using SCE-UA for parameter optimization. AIC and the slopes of regression lines between the simulated and observed SOC changes were analyzed.

$$\begin{array}{}\text{(19)}& \mathrm{\Delta }\mathrm{SOC}=X_t-X_{\mathrm{c}},\end{array}$$

where X_t and X_c are the observed (or simulated) SOC concentrations with and without biochar addition, respectively.

The 134 paired observations were randomly split into training samples for parameter optimization (80 % of the data) and test samples for model validation (20 % of the data). Four tests were conducted to evaluate the performance of MIMICS_TSMb-BC with respect to simulating SOC changes after biochar addition using the optimized parameter values in MIMICS-TSM_b (i.e., a_v, a_k, k_d, β, k_ba, c₁, and c₂; Table S3):

1. without biochar-related parameters (MIMICS_TSMb-BC_def);

2. with only one new biochar-related parameter (i.e., the desorption coefficient, f_d; Eq. 15) optimized (MIMICS_TSMb-BC_D);

3. with a further one new biochar-related parameter (i.e., the decomposition rate coefficient, f_v; Eq. 16) optimized and f_v included in all decomposition processes (MIMICS_TSMb-BC_DV);

4. same as (3) but f_v only included in the fluxes from SOC_a to the MIC pools (MIMICS_TSMb-BC_DV-SOCa).

Although MIMICS-TSM_b is not the model with the highest R² and lowest RMSE and AIC, the differences in the R², RMSE, and AIC among various versions are relatively small (Fig. S5). The new processes (density-dependent processes, sorption, and soil moisture scalars) are based on a theoretical understanding and have been shown to improve predictions of soil carbon in previous studies (Zhang et al., 2020; Liang et al., 2019; Abramoff et al., 2022). Thus, this version was used for further development of biochar processes in MIMICS. As an alternative model version, we also tested the implementation of biochar processes in MIMICS-T, which has the highest R² and lowest RMSE and AIC in model validation (Fig. S5b). The model versions and simulation settings are shown in Table 1 and Fig. 3, and the optimized parameters values in these tests are shown in Table S3.

The relationships between observed ΔSOC and model input variables and the partial correlations between biases (simulated minus observed ΔSOC) from the four tests and model input variables (soil-, climate-, biology-, and biochar-related variables) were also analyzed to detect the possible missing processes. The availability of carbon in biochar may affect the magnitude of PEs. Thus, we tested the MIMICS_TSMb-BC versions using a different initial value of the partitioning coefficient from biochar carbon to SOC_a (f_ba) (2 % according to Lychuk et al., 2014, compared with the original value of 20 % in Sect. 2.2). The partitioning coefficients of f_bp and f_ba were also optimized to evaluate the model performance. Considering the uncertainties in the MIMICS-BC parameters, we conducted a sensitivity test of biochar-related parameters (i.e., f_d, f_v, f_bp, and f_ba), microbe-related parameters (MGE and τ), and input variables (i.e., Rate_BC, Age_BC, NPP, Clay, and SM) by perturbing one parameter at a time while keeping all others unchanged for each site. Four perturbation levels of −50 %, −25 %, 25 %, and 50 % were set.

Table 1 Modifications in various MIMICS versions.

Figure 3 Diagram of field measurement SOC data and the model simulation settings. The simulated or observed ΔSOC is equal to SOC with the biochar addition treatment minus that in the control treatment (without biochar addition). Note that one control treatment may correspond to multiple BC treatments with different applied BC rates at one single site.

3 Results of model calibration and validation

3.1 Performance of different MIMICS versions without biochar for simulating cropland SOC

3.1.1 Calibration and validation for MIMICS versions without biochar

We split the data for calibration (80 %) and validation (20 %) of MIMICS versions without biochar-related processes. When using 20 % of the data for the independent model validation, MIMICS-T also performs best, with the highest accuracy (R²=0.56), the lowest RMSE (4.82 g kg⁻¹), and the lowest AIC (187.2) among all model versions (Figs. 4, S5b). Compared with MIMICS-def (R²=0.51, RMSE = 4.97 g kg⁻¹, AIC = 188.8) (Figs. 4e, S5b), MIMICS-TS and MIMICS-TSM_b have a better correlation (R²=0.52 and 0.52) but a higher RMSE (RMSE = 5.01 and 5.05 g kg⁻¹) and AIC (AIC = 197.7 and 198.6) between the observed and simulated cropland SOC concentration. The R² of the MIMICS-TSM versions ranges from 0.46 to 0.52, and R² of MIMICS-TSM_b is the highest among the model versions. We also evaluated the performance of the MIMICS-TSM_b version calibrated with cropland SOC data under different crop types. The model performance varies among different crop types (i.e., maize, rice, and wheat). The R² between the simulated SOC concentrations by MIMICS-TSM_b and observations is higher for maize and wheat (0.84 and 0.74, respectively; Fig. S6a, c) but lowest for rice (0.38; Fig. S6b). This is probably because the flooded condition in the paddy field limits SOC decomposition, which is partly supported by the weaker correlation between SOC and NPP for rice (R²=0.06; Fig. S7d) than that for maize and wheat (R²=0.77 and 0.54; Fig. S7a, g). Although MIMICS-TSM_b has the soil moisture processes (Sect. 2.1.4), the flooding conditions in the paddy field are not available and are not explicitly forced in these rice sites.

To align with the resolution of climate input variables used in the transient simulations (Sect. 2.3.1), we tested MIMICS after aggregating cropland SOC sites within each 0.5° × 0.5° grid cell instead of using each site directly, and the model can reproduce about 45 %–55 % of the SOC spatial variation (Fig. S8), slightly lower than that using site-specific data (R²=0.51–0.56; Fig. 4e). This shows that site-specific data are more accurate for model simulation, but the small difference also suggests that using grid climate data (0.5° × 0.5°) to drive the model has a relatively small impact on predicting SOC in the transient simulations.

Figure 4 Comparison between the observed and simulated SOC concentrations by (a) MIMICS-def, (b) MIMICS-T, (c) MIMICS-TS, and (d) MIMICS-TSM_b. Blue and red dots in panels (a)–(d) represent observation sites for model calibration (80 % sites) and validation (20 % sites), respectively. Panel (e) presents the R², root-mean-square error (RMSE), and Akaike information criterion (AIC) from the model calibration (left) and validation (right) for the four MIMICS versions. The values for the bar graphs in Fig. 4e can be found in Fig. S5. Relationships for the other MIMICS versions can be found in Fig. S9.

3.1.2 Sensitivity analyses for MIMICS versions without biochar

To explore the possible effects of soil moisture on SOC dynamics by affecting different processes (Sect. 2.3.2), we assumed that the microbial turnover (τ) was also affected by soil moisture in addition to the microbial reaction velocity (V_max). The model with the soil moisture effects does not predict SOC concentrations more accurately (R²=0.46, RMSE = 5.06 g kg⁻¹, AIC = 198.9; Fig. S10b) than the MIMICS-TSM_b version in which V_max and K_m are affected (R²=0.52, RMSE = 5.05 g kg⁻¹, AIC = 198.6; Figs. 4d, S5b). This may be because the inclusion of soil moisture effects on V_max and τ complicates the model processes, and other microbe-related observations, such as soil microbial carbon and soil heterotrophic respiration, are needed to further constrain these processes in MIMICS.

We analyzed the responses of the MIMICS model to changes in input variables, in order to identify the important variables and explore possible missing processes related to these variables (Sect. 2.3.2). The perturbation of MIMICS input variables shows that the size of the steady SOC pool is positively correlated with NPP and Clay but negatively correlated with MAT and BD. The responses of steady SOC to the perturbation of BD, MAT, and NPP are relatively large (Fig. 5). Improving processes associated with these variables in MIMICS will enhance the prediction accuracy of the model.

Figure 5 Sensitivity analysis of responses of the steady SOC simulated by MIMICS to the input variables (a) MAT, (b) Clay, (c) NPP, (d) SM, and (e) BD with different perturbation levels. The yellow lines and the dotted green lines in the box plots are the respective median and mean values of the output steady SOC changes at calibrated sites. The average SOC changes at calibrated sites for the four perturbation levels are shown in panel (f).

3.2 Calibration and evaluation of MIMICS-BC

3.2.1 Calibration and validation for MIMICS-BC versions

For the model validation using short-term observation data that were not used for calibration (Figs. 6, 7, S11, S12), the performance of MIMICS_T-BC_DV-SOCa (R²=0.80, RMSE = 3.38 g kg⁻¹, AIC = 69.8; Fig. 6e–g) is slightly better than MIMICS_T-BC_D (R²=0.79, RMSE=3.43 g kg⁻¹, AIC = 68.5) and MIMICS_T-BC_DV (R²=0.76, RMSE = 3.66 g kg⁻¹, AIC=74.1), except that the AIC (69.8) is higher than that of MIMICS_T-BC_D (68.5) (Fig. 6). By comparison, the performance of MIMICS_T-BC_def is poorer than these three versions. Among the MIMICS_TSMb-BC versions, MIMICS_TSMb-BC_DV performs best with respect to reproducing SOC changes with biochar addition, with the highest R² (0.79), the lowest RMSE (3.73 g kg⁻¹), and AIC (75.0) (Fig. 7e–g).

Figure 6 Relationships of short-term (≤ 6 years; black) and long-term (i.e., extended to 8 years; red) SOC changes after biochar addition (ΔSOC) between observations and models in the validation dataset. The MIMICS versions used include MIMICS_T-BC_def (a), MIMICS_T-BC_D (b), MIMICS_T-BC_DV (c), and MIMICS_T-BC_DV-SOCa (d). Comparisons of R² (e), the root-mean-square error (RMSE, f), and the Akaike information criterion (AIC, g) among the four MIMICS_T-BC versions are shown separately. For information on the model versions, the reader is referred to Table 1.

For the model validation using the long-term (extended to 8 years based on biochar decomposition curve; Wang et al., 2016) SOC changes after biochar addition, MIMICS_T-BC_def and MIMICS_TSMb-BC_def underestimate the extrapolated observations of SOC change (Figs. 6a, 7a). MIMICS_T-BC_D shows the best performance, with the lowest RMSE (3.84 g kg⁻¹) and AIC (74.7), among all of the MIMICS_T-BC versions (Fig. 6). Compared with MIMICS_TSMb-BC_def (R²=0.88, RMSE = 9.35 g kg⁻¹, slope = 0.08, AIC =120.7; Fig. 7a, e, f, g), predictions of MIMICS_TSMb-BC_D, MIMICS_TSMb-BC_DV, and MIMICS_TSMb-BC_DV-SOCa are more accurate, with a smaller RMSE (8.12, 6.08 and 6.78 g kg⁻¹; Fig. 7f), a smaller AIC (115.1, 101.5 and 107.4; Fig. 7g), a linear slope closer to 1 (0.29, 1.68 and 1.74; Fig. 7a–d), and a reasonable R² accuracy (0.45, 0.97 and 0.94; Fig. 7e). Among the different MIMICS_TSMb-BC versions, MIMICS_TSMb-BC_DV shows the best performance (Fig. 7). When assuming that biochar produces a PE only through affecting the utilization rate of SOC_a by microbes (MIMICS_TSMb-BC_DV-SOCa), the model accuracy is slightly poorer, with a lower R² (= 0.94), a higher RMSE (= 6.78 g kg⁻¹), and a higher AIC (= 107.4) than MIMICS_TSMb-BC_DV which assumes that all decomposition processes are affected (Fig. 7).

Figure 7 Same as Fig. 6 but for models based on MIMICS-TSM_b instead of MIMICS-T.

3.2.2 Error analysis for MIMICS-BC versions

The biases between the simulated and observed short-term SOC changes with biochar addition are significantly correlated with Rate_BC and Clay (p<0.05), but they only vary marginally with SM, MAT, and NPP when additional parameters are optimized (Fig. S13). The biases between long-term observations and simulations by MIMICS_TSMb-BC_def are significantly correlated with Rate_BC ($r=-\mathrm{0.81}$) (Fig. 8), suggesting that the model may underrepresent processes related to Rate_BC. By considering biochar effects on the SOC desorption (MIMICS_TSMb-BC_D), the correlations of model biases with Rate_BC, BD, SM, and NPP become weaker (Fig. 8). MIMICS_TSMb-BC_DV, which incorporates the biochar impacts on the microbial decomposition rate, further reduces the correlations between model biases and the variables of Rate_BC, Age_BC, and BD, but it increases biases related to SM. MIMICS_TSMb-BC_DV-SOCa, which includes the impacts on microbial decomposition rate only in the fluxes from the SOC_a to MIC pools, can also reduce the correlations between model biases and the variables of Rate_BC and BD, but the correlations change little with Clay and Age_BC (Fig. 8).

Figure 8 Correlations between the MIMICS_TSMb-BC biases (i.e., simulated long-term ΔSOC minus observed ΔSOC) and input soil-related (Clay, BD, and SM), climate-related (MAT), biology-related (NPP), and biochar-related (Rate_BC and Age_BC) variables for MIMICS_TSMb-BC_def, MIMICS_TSMb-BC_D, MIMICS_TSMb-BC_DV, and MIMICS_TSMb-BC_DV-SOCa. Asterisks indicate significant correlations (p<0.05).

3.2.3 Sensitivity analyses for MIMICS-BC versions

The test of the partitioning coefficient from biochar to SOC_a (f_ba) with 2 % for the MIMICS_TSMb-BC model (Sect. 2.3.3) shows a similar R² (0.35–0.79; Fig. S14) to that from f_ba=20 % in the short term (0.25–0.79; Fig. 7). The optimization of the partitioning coefficient from biochar carbon to SOC_p (f_bp) and f_ba shows a better performance (R²=0.80, RMSE = 3.44 g kg⁻¹, AIC = 66.7; Fig. S15) than MIMICS_TSMb-BC_DV without these two parameters optimized (Fig. 7).

The responses of ΔSOC to parameter perturbations show that f_v and f_d affect ΔSOC with biochar addition in the opposite directions and that ΔSOC is more sensitive to the partitioning coefficient f_bp than to f_d, f_v, and f_ba (Fig. 9). ΔSOC is more sensitive to Rate_BC than Age_BC. The sensitivity of the ΔSOC to τ is greater than that to MGE (Fig. S16), and the two parameters may be influenced by biochar but are not considered in the current MIMICS-BC versions. The sensitivity tests for the model input variables, i.e., crop net primary productivity (NPP), soil clay content (Clay), and soil moisture (SM), show that Clay is very important to the model outputs, while the impacts of NPP and SM are relatively smaller (Fig. S16).

Figure 9 Sensitivity analysis of MIMICS-BC model parameters of (a) f_d (desorption coefficient; Eq. 15), (b) f_v (decomposition rate coefficient; Eq. 16), (c) f_bp (partitioning coefficient from biochar carbon to SOC_p; Fig. 1), (d) f_ba (partitioning coefficient from biochar carbon to SOC_a; Fig. 1), and the biochar-related input variables of (e) Rate_BC and (f) Age_BC. The yellow lines and the dotted green line in the box plots are the respective median and mean values of the changes in model output (i.e., change in ΔSOC; Eq. 19). The mean values of change in output ΔSOC at calibrated sites are shown in panel (g).

4 Discussion

4.1 Cropland SOC-related processes

4.1.1 Missing processes in the MIMICS model

We presented a framework to quantify the impact of density-dependent microbial turnover, sorption, and soil moisture effects on cropland SOC dynamics. Regulatory mechanisms (e.g., competition) may limit microbial population sizes, and neglecting this process could lead to indefinite microbial biomass growth with increasing litter inputs (Georgiou et al., 2017). Our analysis demonstrates that restricting microbial biomass size through density-dependent microbial turnover (MIMICS-T) slightly improves the model performance (Fig. 4), but further including sorption and soil moisture effects (MIMICS-TS and MIMICS-TSM) has a negligible contribution to the model performance. One possible reason for this is that the inclusion of these new processes greatly increases model complexity; however, in the parameter calibrations, SOC is still the only observational variable to constrain all of the newly added processes. Therefore, the model parameters may not be fully constrained due to the lack of direct observations of these new processes.

Another reason is that some other possibly important processes are missing from the model. For example, the MIMICS-TS version considers the impacts of soil clay on the adsorption capacity of SOC_a, but soil pH, ionic strength, and mineral content are also found to be important to the sorption–desorption of SOC (Kothawala et al., 2009; Mayes et al., 2012). The metal ion Ca²⁺ can form bonds between negatively charged clay minerals and available SOC via cation bridging, enhancing the adsorption of available SOC by soil clay minerals (Roychand and Marschner, 2014; Setia et al., 2013). Soil pH can also impact SOC sorption by altering the ionization degree and the surface change in SOC molecules (Shen, 1999). Moreover, iron minerals can preferentially bind to lignin components through sorption and coprecipitation, protecting them from microbial degradation and, consequently, increasing SOC (Liao et al., 2022). In addition to the sorption–desorption process that is associated with microbial accessibility to SOC, other factors that influence microbial activity are also underrepresented in MIMICS, such as the soil nutrient availability (e.g., nitrogen), which greatly impacts microbial use efficiency and growth rate (Manzoni et al., 2017). These processes improve our understanding of the mechanisms of SOC dynamics and should be incorporated in the model to represent the microbe–mineral processes realistically and mechanistically.

However, as shown in our results, without further observational constraints on each respective process, the model performance only relying on calibrations against total SOC contents may not improve. Therefore, various versions of MIMICS, representing different levels of our understanding of microbe–mineral processes, are retained in our study for further calibrations when sufficient observations emerge.

4.1.2 Cropland management impacts

Cropland management disturbs soils frequently, and the assumed equilibrium state of SOC may not be realistic, which also partly explains the mismatch between simulated and observed SOC. Agricultural management (e.g., irrigation and tillage) is an important factor that affects SOC decomposition and accumulation in croplands. The poor performance of MIMICS for rice is probably due to the inability of MIMICS to simulate SOC dynamics under anaerobic conditions stemming from the irrigation practice (Figs. S6–S7). Tillage may disrupt soil aggregates and release physically protected SOC, which is more susceptible to decomposition than that protected by soil aggregates (Six et al., 1999). Juice et al. (2022) modeled tillage effects on SOC loss via the transfer of protected SOC into unprotected pools, i.e., from SOC_p to SOC_a in this study, and the model can well capture the historical SOC dynamics in the agricultural system with intensive managements, such as different crop types, tillage, or fertilization. Our results from variable perturbation suggest that BD is the key driver of SOC changes, followed by MAT and crop NPP (Fig. 5), suggesting that processes related to these variables have a great effect on the SOC. The soil BD was also found to be affected by tillage practices (Osunbitan et al., 2005), and crop NPP may vary due to crop rotation and fallow practice, which are missing in the model.

In addition, managements such as fertilizer application and possible residue retention can increase SOC stock. Previous evidence indicates that SOC is increased by 11.3 % with residue return (X. Wang et al., 2020) and by 13.3 % with straw return and balanced nitrogen–phosphorus–potassium (NPK) fertilizer (Islam et al., 2023) compared with residue removal. However, these management processes of fertilization and residue retention are not represented in the MIMICS model due to the absence of quantitative management data and the poor understanding of the mechanisms. This may explain the underestimation of SOC by the calibrated MIMICS models at sites with a high carbon density (Fig. 4). Therefore, field measurements of the effects of agricultural practices on SOC dynamic are urgently needed to further improve the model processes (Campbell et al., 2007; Congreves et al., 2015).

4.2 Biochar-related processes

4.2.1 Tested processes in the MIMICS-BC model

Biochar can absorb SOC due to its large specific surface area, high porosity, and further promotion of soil macroaggregate formation (Han et al., 2020; Huang et al., 2018). Consistently, the optimized desorption coefficient ($f_{\mathrm{d}}=-\mathrm{0.0121}$ and −0.0122 for the short term and long term, respectively; Table S3) in MIMICS_TSMb-BC_D is negative, indicating that the carbon desorption from SOC_p to SOC_a is reduced with biochar addition. Incorporating the biochar impacts on microbial decomposition velocity in MIMICS_TSMb-BC_DV further improved the model with biochar addition, but the correlations between model–observation biases and the input variables of biochar application rate (Rate_BC) and soil moisture (SM) are still significant in the long term (p<0.05; Fig. 8), implying that some processes related to these variables are not well represented in the model. For example, biochar addition could increase soil moisture (Razzaghi et al., 2020) and further alter the composition and activity of soil microbial communities (Lehmann et al., 2011). Moreover, the direction and magnitude of biochar effects on SOC are dependent on the biochar addition rate and Age_BC (Ding et al., 2017). Compared with MIMICS_T-BC_DV, MIMICS_T-BC_DV-SOCa performed better with respect to the short-term SOC response to biochar addition but worse in the long-term response (Fig. 6). When biochar is applied, the labile carbon fraction may be immediately utilized by microbes; thus, including the biochar effects on the SOC_a process is important for the short-term response. In the long term, the SOC mineralization may be gradually suppressed via physical protection (Zimmerman et al., 2011), but both MIMICS_T-BC_DV and MIMICS_T-BC_DV-SOCa do not include the adsorption process. In the MIMICS_TSMb-BC versions, which include the adsorption process, the available SOC may be partly adsorbed by soil minerals and become physically protected. This could lead to the positive PE of biochar on SOC being less evident (Fig. 7).

4.2.2 Missing processes in current MIMICS-BC model

The effects of biochar on SOC are controlled by various factors, such as soil physicochemical and biological properties (e.g., clay, pH, and microbial activity), biochar properties (e.g., feedstock and pyrolysis temperature), and incubation conditions (e.g., periods and crop types) (Ding et al., 2017; Han et al., 2020). Some of these effects are not explicitly considered in the MIMICS biochar version. Biochar addition may change the composition of microbial community, and a previous study reported increased copiotrophic bacteria with a higher growth rate and decreased oligotrophic bacteria in acidic soils with biochar addition (Sheng and Zhu, 2018). This is related to the competition between r- and k-strategy microbes in MIMICS. Microbial carbon use efficiency (CUE) determines the relational proportions of microbial carbon uptake between growth and respiration (H. Zhou et al., 2017), and an increased CUE and reduced turnover time ($\mathrm{1}/\mathit{\tau }$) of microbial biomass were found with biochar addition, although the changes depended on the soil texture (Pei et al., 2021). This is consistent with our results that the ΔSOC is more sensitive to changes in τ and soil clay than other parameters and variables (Fig. S16). Therefore, processes and parameters related to τ and soil clay need to be accounted for in future with additional evidence.

In the MIMICS-BC version, we assumed that biochar, with a longer turnover time (about 1000 years; Schmidt et al., 2002) than SOC, is evenly mixed with SOC and is treated as a homogenous pool without an explicit vertical profile, which may also introduce uncertainties. In addition, due to the lack of long-term biochar addition experiments, the extended long-term SOC concentrations with biochar addition are calculated as the sum of SOC at the control site without biochar addition and the remaining biochar carbon based on the biochar degradation curve (Fig. S4; Wang et al., 2016). Although they are not direct observations and may induce uncertainty, the long-term model validation is important to assess the model ability to simulate the SOC stability with biochar addition. Therefore, long-term and comprehensive field measurements of SOC and other soil and microbe properties after biochar addition are urgently needed to understand the underlying mechanisms of biochar impacts on SOC changes, all of which will help improve the model performance.

5 Conclusions

Our study shows that the updated MIMICS versions with new processes (e.g., adsorption and soil moisture) improves the model performance with respect to simulating SOC dynamics on croplands. The model versions implemented with biochar processes can generally capture the SOC changes after biochar application from observations. Biochar is believed to have a large CDR potential, and its application on soils would affect the soil carbon and nutrient cycles. These impacts need to be incorporated ESMs to accurately simulate the mitigation potential of biochar under future climate change.