Improved model simulation of soil carbon cycling by representing the microbially derived organic carbon pool

Abstract

During the decomposition process of soil organic carbon (SOC), microbial products such as microbial necromass and microbial metabolites may form an important stable carbon (C) pool, called microbially derived C, which has different decomposition patterns from plant-derived C. However, current Earth System Models do not simulate this microbially derived C pool separately. Here, we incorporated the microbial necromass pool to the first-order kinetic model and the Michaelis–Menten model, respectively, and validated model behaviors against previous observation data from the decomposition experiments of ¹³C-labeled necromass. Our models showed better performance than existing models and the Michaelis–Menten model was better than the first-order kinetic model. Microbial necromass C was estimated to be 10–27% of total SOC in the study soils by our models and therefore should not be ignored. This study provides a novel modification to process-based models for better simulation of soil organic C under the context of global changes.

Introduction

The pool size of global soil organic carbon (SOC) to 1 m depth is estimated to be 1417–1469.5 Pg C [1, 2], which is much larger than the amounts of carbon (C) in plants and atmosphere. The uncertainty of current process-based models in estimating global SOC size is 50% [3], meaning that about 630–735 Pg C still cannot be accurately simulated. Measures have been proposed to improve Earth System Models, such as consideration of nitrogen constraints [4, 5], or incorporation of density-dependent formulation of microbial turnover [6]. Although other sources besides our understanding of the potential processes may also contribute to the model uncertainty such as error in observation data given the low precision of soil carbon measurement, uncertainty in initial conditions, and uncertainty in parameters due to the limited available data for constraining model parameters, it is still important to incorporate the mechanisms already confirmed in observation experiments into process-based models for better simulation and prediction of C cycling and climate changes.

SOC pools in process-based models are usually divided into different parts based on their decomposition rates, such as the slow pool and the fast pool. However, SOC is composed of a very complex mixture of polymers of microbial and plant residues and degradation products [7], and different compounds not only differ in decomposition rate, but also differ in decomposition pattern. For example, plant-derived organic N and microbially derived organic N have been found to have different decomposition curves [8]. In addition, the microbial necromass pool (i.e., mass from the cell death and subsequent lysis and fragmentation of soil microbes) may interact with soil minerals more easily than the plant residues [9], making it necessary to simulate the microbial necromass pool separately from the plant residue pool. Although the living microbial biomass in the soil is a relatively small C pool, microbially derived SOC may accumulate to form a large proportion pool if microbial biomass turnover rate is higher than plant litter input rate [10], or the chemical nature of microbial necromass is not labile or microbially derived SOC is protected by the soil mineral matrix [11,12,13,14]. Recently, a variety of methods have been used to estimate the contribution of microbial necromass C to SOC in different soils, and a range of 24–80% was determined [11, 15,16,17]. Therefore, the role of microbial necromass in the formation of organic matter may have been seriously underestimated, and should be considered in models [18].

In most process-based models, the decomposition of SOC is expressed by the first-order kinetic equation [19]. The decomposition coefficient is generally determined by the environmental factors [20], and the roles of the microbial community play in the decomposition of organic matter have not been well represented in these models [5]. Recently, the Michaelis–Menten equation was used to describe the decomposition rate of organic matter in many models, such as the MEND [21,22,23], GER [24], AWB [25, 26], MIMICS [27], RESOM [28, 29], CORPSE [30], ORCHIMIC [5] models. These models consider the microbial biomass carbon (MBC) pool separately and have been found to provide better simulation results than the first-order kinetic models [3, 31]. However, the microbial necromass pool has not been considered as a separate pool with a different decomposition pattern from plant residue in these models. We propose a conceptual four-pool model including MBC and microbial necromass pools (Supplementary Fig. 1) because once MBC pool is separated, it is undoubtedly logical to separate the necromass pool to better represent the formation process of microbially derived SOC.

Recent findings have made it possible to realize the conceptual model because the decomposition pattern of microbial necromass has been observed [8] and the responses of the microbial necromass pool to environmental factors have been studied [32]. In addition, although fungi, gram-positive bacteria (including Actinomycetes), and gram-negative bacteria have different molecular structures, their turnover rates were found to be similar, and the quality of microbial necromass was found to have no significant effect on their decomposition rates [33], suggesting that the same decomposition rate can be used for different groups of microbial necromass in the model.

In this study, we established a new model accounting for the microbial necromass pools, and demonstrated its validity and feasibility using observation data from previous experiments of the decomposition of ¹³C-labeled microbial necromass [33,34,35]. Two archetypal models were selected [24,25,26] and the microbial necromass pools were added into these two models and we named these two models as the Michaelis–Menten necromass decomposition (MIND) model and the first-order necromass decomposition (FOND) model (Fig. 1). The aims of this study were: (1) to evaluate the necessity of incorporating the microbial necromass pool into ESMs and to compare the performance of MIND and FOND models; (2) to investigate the decomposition rate of microbial necromass C in different soils; and (3) to estimate the contribution of microbial necromass C to soil organic C under steady state. Our results suggest that microbial necromass C is an essential, but overlooked part of stable SOC, and model performance is improved by representing the microbially derived organic carbon pool in process-based models.

Fig. 1: The structure of the two newly proposed models with the microbial necromass pools.

A The Michaelis–Menten necromass decomposition model (MIND); B the first-order necromass decomposition model (FOND). The red dotted line represents the paths involved when ¹³C-labeled microbial necromass is added to the system. The valve symbol and black dashed line represent the regulation of the size of MBC pool on the processes according to Michaelis–Menten kinetics. C_P plant-derived carbon, C_B microbial biomass carbon, C_NF fast pool of microbial necromass carbon, C_N-MAOM mineral-associated pool of microbial necromass carbon, C_DP plant-derived dissolved organic carbon, C_DN microbially derived dissolved organic carbon (color figure online).

Method

Model structure

The model structures are shown in Fig. 1. In MIND model, the Michaelis–Menten equation was used to describe the nonlinear variation of substrate decomposition rate with microbial biomass. The FOND model assumed that the variation of the decomposition rate was controlled by environmental factors and in our simulations the decomposition rate did not change because we used incubation experiment data under controlled environment. In both models, after the death of soil microorganisms they are transferred to the fast pool and the mineral-associated pool of microbial necromass proportionally. We minimized the number of additional C pools to avoid introducing more parameters and overfitting.

MIND was expressed by the following equations:

$$\frac{{{d}{\mathrm{C}}_P}}{{{d}t}} = I - \left(\frac{{V_{\max ,P} \times {\mathrm{C}}_B \times {\mathrm{C}}_P}}{{K_{M,P} + {\mathrm{C}}_P}}\right)$$

(1)

$$\frac{{{\it{d}}{\mathrm{C}}_B}}{{{\it{d}}t}} = \, {\mathrm{CUE}}_P \times \left(\frac{{V_{\max ,P} \times {\mathrm{C}}_B \times {\mathrm{C}}_P}}{{K_{M,P} + {\mathrm{C}}_P}}\right) + {\mathrm{CUE}}_N \times \left(\frac{{V_{\max ,N} \times {\mathrm{C}}_B \times {\mathrm{C}}_{{\mathrm{NF}}}}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{NF}}}}}\right)\\ + {\mathrm{CUE}}_N \times R_{{\mathrm{MAOM}} - {\mathrm{F}}} \times \left(\frac{{V_{\max ,N} \times {\mathrm{C}}_B \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}}}\right) - k_B \times {\mathrm{C}}_B$$

(2)

$$\frac{{{d}{{{\rm{C}}_{{\rm{NF}}}}}}}{{{d}t}} = f_{{\mathrm{BNF}}} \times k_B \times {\mathrm{C}}_B - \left(\frac{{V_{\max ,N} \times {\mathrm{C}}_B \times {{{\rm{C}}_{{\rm{NF}}}}}}}{{K_{M,N} + {{{\rm{C}}_{{\rm{NF}}}}}}}\right)$$

(3)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}}}{{{d}t}} = (1 - f_{{\mathrm{BNF}}}) \times k_B \times {\mathrm{C}}_B - R_{{\mathrm{MAOM}} - {\mathrm{F}}} \\ \times \left(\frac{{V_{\max ,N} \times {\mathrm{C}}_B \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}}}\right)$$

(4)

$$\frac{{{d}{\mathrm{C}}_N}}{{{d}t}} = \frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}}}{{{d}t}} + \frac{{{d}{\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}}}{{{d}t}}$$

(5)

where C_P, C_B, C_NF, C_N-MAOM, and C_N are pool sizes (mg C g⁻¹ soil) of plant-derived organic carbon, microbial biomass carbon (MBC), microbial fast necromass carbon, microbial mineral-associated necromass carbon, total microbial necromass carbon (MNC), respectively; I is the annual average carbon input rate (mg C g−1 soil h⁻¹) to soil from plant-derived litter, which was calculated based on litterfall minus litter heterotrophic respiration (g C m⁻² s⁻¹) and then converted to mg C g⁻¹ soil h⁻¹ using bulk density (0–15 cm) data; V_max,P, V_max,N are maximum assimilation rate (mg C mg⁻¹ MBC h⁻¹) of C_P and C_N, respectively; K_M,P, K_M,N are half-saturation for assimilation (mg C g⁻¹ soil) of C_P and C_N, respectively; CUE_P, CUE_N are carbon use efficiency (unitless) of C_P and C_N, respectively; k_B is average mortality rate of the microbial community (h⁻¹); f_BNF is proportion of fast pool in MBC (unitless); R_MAOM-F is ratio of the decomposition rate of the mineral-associated pool to that of the fast pool of microbial necromass (unitless).

When the ¹³C-labeled microbial necromass C was added to the system, we can get the following equations:

$$\frac{{{d}{\mathrm{C}}_P^u}}{{{d}t}} = I - \left(\frac{{V_{\max ,P} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_P^u}}{{K_{M,P} + {\mathrm{C}}_P^u}}\right)$$

(6)

$$\frac{{{d}{\mathrm{C}}_B^u}}{{{d}t}} = \, {\mathrm{CUE}}_p \times \left(\frac{{V_{\max ,P} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_P^u}}{{K_{M,P} + {\mathrm{C}}_P^u}}\right) + {\mathrm{CUE}}_N\\ \,\times \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{NF}}}^u}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{NF}}}^u}}\right) + {\mathrm{CUE}}_N \times R_{{\mathrm{MAOM}} - {\mathrm{F}}}\\ \, \times \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^u}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^u}}\right) - k_B \times {\mathrm{C}}_B^{u}$$

(7)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}^u}}{{{d}t}} = f_{{\mathrm{BNF}}} \times k_B \times {\mathrm{C}}_B^u - \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{NF}}}^u}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{NF}}}^u}}\right)$$

(8)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^u}}{{{d}t}} = \, (1 - f_{{\mathrm{BNF}}}) \times k_B \times {\mathrm{C}}_B^u - R_{{\mathrm{MAOM}} - {\mathrm{F}}} \\ \, \times \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^u}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^u}}\right)$$

(9)

$$\frac{{{d}{\mathrm{C}}_N^u}}{{{d}t}} = \frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}^u}}{{{d}t}} + \frac{{{d}{\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^u}}{{{d}t}}$$

(10)

$$\frac{{{d}{\mathrm{C}}_B^l}}{{{d}t}} =\, {\mathrm{CUE}}_N \times \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{NF}}}^l}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{NF}}}^l}}\right) \\ + {\mathrm{CUE}}_N \times R_{{\mathrm{MAOM}} - {\mathrm{F}}} \\ \times \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l}}\right) - k_B \times {\mathrm{C}}_B^l$$

(11)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}^l}}{{{d}t}} = f_{{\mathrm{BNF}}} \times k_B \times {\mathrm{C}}_B^l - \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{NF}}}^l}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{NF}}}^l}}\right)$$

(12)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{N - MAOM}}}^l}}{{{d}t}} = (1 - f_{{\mathrm{BNF}}}) \times k_B \times {\mathrm{C}}_B^l - R_{{\mathrm{MAOM - F}}} \\ \times \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{N - MAOM}}}^l}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{N - MAOM}}}^l}}\right)$$

(13)

$$\frac{{{d}{\mathrm{C}}_N^l}}{{{d}t}} = \frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}^l}}{{{d}t}} + \frac{{{d}{\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l}}{{{d}t}}$$

(14)

$$\frac{{{d}{{{\rm{CO}}_2}}^l}}{{{d}t}} = \; (1 - {\mathrm{CUE}}_N) \times \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{NF}}}^l}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{NF}}}^l}}\right) \\ + (1 - {\mathrm{CUE}}_N) \\ \times R_{{\mathrm{MAOM}} - {\mathrm{F}}} \times \left(\frac{{V_{\max ,N} \times ({\mathrm{C}}_B^u + {\mathrm{C}}_B^l) \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l}}{{K_{M,N} + {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l}}\right)$$

(15)

where superscript u stands for “unlabeled” and l stands for “labeled” pool; the size of a pool was the sum of labeled proportion and unlabeled proportion. For example, total MBC was the sum of labeled MBC and unlabeled MBC.

The analytical steady-state solutions of MIND were as follows:

$${\rm{C}}_B = \frac{{{\rm{CUE}}_P \times I}}{{(1 - {\rm{CUE}}_N) \times k_B}}$$

(16)

$${\rm{C}}_P = \frac{{I \times K_{M,P}}}{{(V_{\max ,P} \times {\rm{C}}_B - I)}}$$

(17)

$${\rm{C}}_{{\rm{NF}}} = \frac{{f_{{\rm{BNF}}} \times k_B \times K_{M,N}}}{{(V_{\max ,N} - f_{{\rm{BNF}}} \times k_B)}}$$

(18)

$${\rm{C}}_{{\mathrm{N}} - {\rm{MAOM}}} = \frac{{(1 - f_{{\rm{BNF}}}) \times k_B \times K_{M,N}}}{{(R_{{\rm{MAOM}} - {\mathrm{F}}} \times V_{\max ,N} - (1 - f_{{\rm{BNF}}}) \times k_B)}}$$

(19)

$${\rm{C}}_N = {\mathrm{C}}_{{\rm{NF}}} + {\rm{C}}_{{\mathrm{N}} - {\rm{MAOM}}}$$

(20)

$${\rm{SOC}} = {\rm{C}}_P + {\rm{C}}_B + {\rm{C}}_N$$

(21)

FOND was expressed by the following equations:

$$\frac{{{d}{\mathrm{C}}_P}}{{{d}t}} = f \times I + f_{{\mathrm{DP}}} \times k_{{\mathrm{DP}}} \times {\mathrm{C}}_{{\mathrm{DP}}} - k_P \times {\mathrm{C}}_P$$

(22)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{DP}}}}}{{{d}t}} = (1 - f) \times I + f_P \times k_P \times {\mathrm{C}}_P - k_{{\mathrm{DP}}} \times {\mathrm{C}}_{{\mathrm{DP}}} - k_{{\mathrm{uptake}}} \times {\mathrm{C}}_{{\mathrm{DP}}}$$

(23)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{DN}}}}}{{{d}t}} \,=\, f_N \times (k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{NF}}} + R_{{\mathrm{MAOM}} - {\mathrm{F}}} \times k_{{\mathrm{NF}}} \\ \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}) - k_{{\mathrm{DN}}} \times {\mathrm{C}}_{{\mathrm{DN}}} - k_{{\mathrm{uptake}}} \times {\mathrm{C}}_{{\mathrm{DN}}}$$

(24)

$$\frac{{{d}{\mathrm{C}}_D}}{{{d}t}} = \frac{{{d}{\mathrm{C}}_{{\mathrm{DP}}}}}{{{d}t}} + \frac{{{d}{\mathrm{C}}_{{\mathrm{DN}}}}}{{{d}t}}$$

(25)

$$\frac{{{d}{\mathrm{C}}_B}}{{{d}t}} = {\mathrm{CUE}} \times k_{{\mathrm{uptake}}} \times ({\mathrm{C}}_{{\mathrm{DP}}} + {\mathrm{C}}_{{\mathrm{DN}}}) - k_B \times {\mathrm{C}}_B$$

(26)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}}}{{{d}t}} = f_{{\mathrm{BNF}}} \times k_B \times {\mathrm{C}}_B - k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{NF}}}$$

(27)

$$\frac{{{d}{\mathrm{C}}_{N - {\mathrm{MAOM}}}}}{{{d}t}} \,=\, (1 - f_{{\mathrm{BNF}}}) \times k_B \times {\mathrm{C}}_B + k_{{\mathrm{DN}}} \\ \times {\mathrm{C}}_{{\mathrm{DN}}} - R_{{\mathrm{MAOM}} - {\mathrm{F}}} \times k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}$$

(28)

$$\frac{{{d}{\mathrm{C}}_N}}{{{d}t}} = \frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}}}{{{d}t}} + \frac{{{d}{\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}}}{{{d}t}}$$

(29)

where C_D is dissolved organic carbon pool (DOC, mg C g⁻¹ soil), which is the sum of plant-derived DOC (C_DP) and microorganism-derived DOC (C_DN); f is the fraction of inputs into C_P (unitless); k_P, k_DP, k_DN, and k_NF are the decomposition rate (mg C mg⁻¹ C h⁻¹) of C_P entering C_DP, C_DP entering C_P, C_DN entering C_N-MAOM, C_NF entering C_DN, respectively; f_P, f_DP, and f_N are fraction (unitless) of C_P entering C_DP, C_DP entering C_P, C_N entering C_DN, respectively; k_B is turnover rate constant of C_B (h⁻¹); k_uptake is uptake rate constant of C_D (mg C g⁻¹ DOC h⁻¹); CUE is carbon use efficiency (unitless); f_BNF is proportion of fast pool in MBC (unitless); R_MAOM-F is the decomposition rate ratio of the mineral-associated pool to that of the fast pool (unitless).

When ¹³C-labeled microbial necromass C was added to the system, we can get the following equations:

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{DN}}}^l}}{{{d}t}} \,=\, f_N \times (k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{NF}}}^l + R_{{\mathrm{MAOM}} - {\mathrm{F}}} \times k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l) \\ - k_{{\mathrm{DN}}} \times C_{{\mathrm{DN}}}^l - k_{{\mathrm{uptake}}} \times {\mathrm{C}}_{{\mathrm{DN}}}^l$$

(30)

$$\frac{{{d}{\mathrm{C}}_B^l}}{{{d}t}} = {\mathrm{CUE}} \times k_{{\mathrm{uptake}}} \times {\mathrm{C}}_{{\mathrm{DN}}}^l - k_B \times {\mathrm{C}}_B^l$$

(31)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}^l}}{{{d}t}} = f_{{\mathrm{BNF}}} \times k_B \times {\mathrm{C}}_B^l - k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{NF}}}^l$$

(32)

$$\frac{{{d}{\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l}}{{{d}t}} \,=\, (1 - f_{{\mathrm{BNF}}}) \times k_B \times {\mathrm{C}}_B^l + k_{{\mathrm{DN}}} \times {\mathrm{C}}_{{\mathrm{DN}}}^l \\ - R_{{\mathrm{MAOM}} - {\mathrm{F}}} \times k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l$$

(33)

$$\frac{{{d}{\mathrm{C}}_N^l}}{{{d}t}} = \frac{{{d}{\mathrm{C}}_{{\mathrm{NF}}}^l}}{{{d}t}} + \frac{{{d}{\mathrm{C}}_{{\mathrm{N}} - {\mathrm{MAOM}}}^l}}{{{d}t}}$$

(34)

$$\frac{{{d}{\mathrm{CO}}_2^l}}{{{d}t}} \,=\, (1 - {\mathrm{CUE}}) \times k_{{\mathrm{uptake}}} \times {\mathrm{C}}_{{\mathrm{DN}}}^l + (1 - f_N) \\ \times (k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{NF}}}^l + R_{{\mathrm{MAOM - F}}} \times k_{{\mathrm{NF}}} \times {\mathrm{C}}_{{\mathrm{N - MAOM}}}^l)$$

(35)

The analytical steady-state solutions of FOND were as follows:

$${\rm{C}}_{{\rm{DP}}} = \frac{{(1 - f) \times I + f_P \times f \times I}}{{k_{{\rm{DP}}} + k_{{\rm{uptake}}} - f_P \times f_{{\rm{DP}}} \times k_{{\rm{DP}}}}}$$

(36)

$${{{\rm{C}}_{{\rm{DN}}}}} = \frac{{f_N \times {\rm{CUE}} \times k_{{\rm{uptake}}} \times {\rm{C}}_{{\rm{DP}}}}}{{k_{{\rm{DN}}} + k_{{\rm{uptake}}} - f_N \times k_{{\rm{DN}}} - f_N \times {\rm{CUE}} \times k_{{\rm{uptake}}}}}$$

(37)

$${\rm{C}}_P = \frac{{f \times I + f_{{\rm{DP}}} \times k_{{\rm{DP}}} \times {\rm{C}}_{{\rm{DP}}}}}{{k_P}}$$

(38)

$${\rm{C}}_B = \frac{{{\rm{CUE}} \times k_{{\rm{uptake}}} \times {\rm{C}}_D}}{{k_B}}$$

(39)

$${{{\rm{C}}_{{\rm{NF}}}}} = \frac{{f_{{\rm{BNF}}} \times k_B \times {\rm{C}}_B}}{{k_{{\rm{NF}}}}}$$

(40)

$${{{\rm{C}}_{{\rm{N}} - {\rm{MAOM}}}}} = \frac{{(1 - f_{{\rm{BNF}}}) \times k_B \times {\rm{C}}_B + k_{{\rm{DN}}} \times {{{\rm{C}}_{{\rm{DN}}}}}}}{{R_{{\rm{MAOM - F}}} \times k_{{\rm{NF}}}}}$$

(41)

$${\rm{C}}_D = {{{\rm{C}}_{{\rm{DP}}}}} + {{{\rm{C}}_{{\rm{DN}}}}}$$

(42)

$${\rm{C}}_N = {{{\rm{C}}_{{\rm{NF}}}}} + {{{\rm{C}}_{{\rm{N}} - {\rm{MAOM}}}}}$$

(43)

$${\rm{SOC}} = {\rm{C}}_P + {\rm{C}}_D + {\rm{C}}_B + {\rm{C}}_N$$

(44)

The mean residence time (MRT) for the necromass ¹³C in FOND model was calculated according to the following equation according to literatures [32, 36]:

$${\rm{MRT}} = f_{{\rm{BNF}}} \times \frac{1}{{k_{{\rm{NF}}}}} + (1 - f_{{\rm{BNF}}}) \times \frac{1}{{R_{{\rm{MAOM - F}}} \times k_{{\rm{NF}}}}}$$

(45)

Models MIND-O and FOND-O had the same structure as MIND and FOND except that the MNC pool was not considered (Supplementary Fig. 2). The equations for MIND-O and FOND-O models are listed in the figure legend of Supplementary Fig. 2.

Observation data

We searched the Web of Science for observation data on the decomposition of ¹³C-labeled microbial necromass and found three articles [33,34,35] and four data sets (named studies 1–4; Supplementary Tables 1 and 2 and Supplementary Fig. 3). In short, soil microorganisms were cultured in a matrix with resources enriched in ¹³C. The microorganisms used ¹³C as their carbon sources and were then killed to form ¹³C-enriched microbial necromass. The ¹³C-enriched microbial necromass was then put into soil for lab incubation in studies 1 and 2 and for in situ decomposition experiments in studies 3 and 4. Studies 1 and 2 used ¹³C-labeled E. coli necromass, and studies 3 and 4 used ¹³C-labeled microbial groups including fungi, gram-positive bacteria (including Actinomycetes), and gram-negative bacteria. The soils in studies 1 and 2 were from temperate farmland and forests ecosystems, respectively, and the soils in studies 3 and 4 were incubated in temperate and tropical forest ecosystems, respectively.

Parameter estimation

The model parameterization was based on the ¹³C recovery data from each observation study. Some parameters had fixed values from previous findings, while most parameters were given a range first, within which the parameter had uniform distribution (Supplementary Table 3). Then these parameters without reference values were estimated using the nonlinear least squares method. We first used the ODE45 solver to calculate the integration of the differential equations y = f(t,y) from t₀ to t_n, with an initial condition of y₀. The ODE45 solver was based on an explicit Runge–Kutta (4,5) equation (the Dormand–Prince pair). It is a single-step solver, which needs only the solution at the immediate preceding time, y(t_n–1) [37, 38]. Then, the lsqnonlin function was used to fit estimated ¹³C recovery to measured data to get the best combination of parameters.

The lsqnonlin function uses nonlinear least squares to fit the data and solve the problems with constraints of parameters within the given ranges. The lsqnonlin function uses the trust-region-reflective algorithm [39, 40] by default to adapt to the problem that Levenberg–Marquardt algorithm [41] cannot constrain the ranges of parameters.

Model validation

We used the leave-one-out cross-validation (LOOCV) [42], which is for model validation of small sample size like our data. The learning algorithm was applied once for each instance, using all other instances as a training set and then the selected instance was used as a single-item test set. This step continued until each sample was treated as a validation.

We also tested the model with the recovery of ¹³C in respired CO₂ values in each observation study, as well as SOC and MBC values under steady-state conditions from each observation studies. Linear regression was performed between estimated data and measured data.

In the four studies on the decomposition of ¹³C-labeled necromass, the units of the decomposition rates of microbial necromass in their calculations were different from those in our study and therefore cannot be compared to our values. Therefore, we compared the ratio of the decomposition rate of the mineral-associated pool to that of the fast pool (R_MAOM-F) of microbial necromass and the proportion of fast pool in MBC (f_BNF) of modeled results to those originally estimated in their studies because these are unitless values. The exponential model was used to calculate the decomposition rates of the fast and mineral-associated pools of necromass in studies 1 and 2 [34, 35]. In studies 3 and 4, the initial decomposition rate and the final decomposition rate were given, which approximately corresponded to the decomposition rate of fast and mineral-associated pools in our study, respectively [33]. Since f_BNF was not given in studies 3 and 4, we only compared it with the values in studies 1 and 2.

Slope, R², root mean squared error (RMSE), and mean absolute error (MAE) were calculated to quantify the accuracy of the model estimates, following equations:

$$R^{\mathrm{2}} = {\mathrm{1}} - \frac{{\mathop {\sum}\nolimits_{i = 1}^n {(\hat y_i - y_i)^2} }}{{\mathop {\sum}\nolimits_{i = 1}^n {(\bar y_i - y_i)^2} }}$$

(46)

$${\rm{RMSE}} = \sqrt {\frac{{\mathop {\sum}\nolimits_{i = 1}^n {(\hat y_i - y_i)^2} }}{n}}$$

(47)

$${\rm{MAE}} = \frac{{\mathop {\sum}\nolimits_{i = 1}^n {\left| {\hat y_i - y_i} \right|} }}{n}$$

(48)

where $\hat y$ is simulation value, y is observations value, $\bar y$ is average value, n is number of observations.

Uncertainty analysis

The Lsqnonlin function starts from the initial value y₀, and iteratively finds the minimum value of the sum of squares of the function between the measured value and the modeled value. Because the Lsqnonlin function depends on the settings of initial values and iteration intervals, it was a local optimal solution instead of global optimal solution. Thus, it is possible that the best solution was missed. We used the Monte Carlo random sampling method to extract subintervals and initial values randomly within the parameter range for 1000 times. The distribution of the parameter was selected when R² between modeled results and observation data was >0.9.

The Box–Cox transformation was used to normalize the data when mean of the parameter distribution was not in the center of the interval, and the data did not obey the normal distribution:

$$W_i = \; Y_i^{\lambda} ,\lambda \,\ne\, 0\\ W_i = \, \ln (Y_i),\lambda = 0$$

(49)

where Y_i is the initial value, λ is the transformation parameter.

The uncertainty of the parameter was presented as the 95% confidence interval of the normal distribution.

Sensitivity analysis

Sensitivity analysis was performed by changing from the lower end to the higher end of 95% confidence interval of the parameter distribution and analyzing the changes in the sizes of different carbon pools in the model. The sensitivity was expressed as follows Allison et al. [25]:

$$\frac{\left| {\rm{log}}_{10}\left| {\mathrm{High}}\,{\mathrm{output}} \right| - {\rm{log}}_{10}\left| {\mathrm{Low}}\,{\mathrm{output}} \right| \right|}{\left| {\rm{log}}_{10}\left| {\mathrm{High}}\,{\mathrm{parameter}} \right| - {\rm{log}}_{10}\left| {\mathrm{Low}}\,{\mathrm{parameter}} \right| \right|}$$

(50)

Global scale estimation

Global microbial necromass C under steady state was estimated using Eqs. (18)–(20). We used our data of the linear relationships between microbial carbon use efficiency (CUE) and mean annual temperature (MAT) and mean annual precipitation (MAP) to estimate microbial CUE of plant residue (CUE_p) in different areas (Eq. (51)): [43]

$${\rm{CUE}}_P = - 1.6 \times 10^{ - 2} \times {\rm{MAT}}({\,}^\circ {\rm{C}}) + 1.3 \times 10^{ - 5} \times {\rm{MAP}}({\rm{mm}}) + 0.41$$

(51)

Microbial CUE of microbial necromass (CUE_N) was assumed to be 0.1 higher than CUE_p because the C:N ratio of microbial necromass is lower than that of plant residue [44,45,46]. We established linear relationships between f_BNF, R_NSF, k_B and MAT or MAP based on the four observation studies (Eqs. (52)–(54)) and extrapolated these parameters to global estimates:

$$f_{{\mathrm{BNF}}} = {\mathrm{2}}{\mathrm{.74}} \times {\mathrm{10}}^{ - 2} \times {\mathrm{MAT}} + {\mathrm{0}}{\mathrm{.35}}$$

(52)

$$R_{{\mathrm{NSF}}} = - {\mathrm{3}}{\mathrm{.49}} \times {\mathrm{10}}^{ - 6} \times {\mathrm{MAP}} + 6.59 \times 10^{ - 3}$$

(53)

$$k_B = - 3.63 \times 10^{ - 4} \times {\mathrm{MAT}} + 1.06 \times 10^{ - 2}$$

(54)

V_max,N, K_M,N were estimated according to the Arrhenius equations (Eqs. (55)–(57)) using soil temperature (T) according to literatures [21, 25]:

$$V_{{\mathrm{max}} ,N} = V_{{\mathrm{max}} ,0} \times f(T)$$

(55)

$$K_{M,N} = K_{M,0} \times f(T)$$

(56)

$$f(T) = e^{\frac{{ - Ea}}{R}(\frac{1}{T} - \frac{1}{{{\mathrm{Tref}}}})}$$

(57)

The decomposition rate of microbial residues was estimated by soil moisture, soil pH, and soil clay content according to literatures [5, 20, 47, 48]:

$$\frac{{d{\mathrm{C}}_N}}{{dt}} = - f(\theta ) \times f({\mathrm{pH}}) \times f({\mathrm{Clay}}) \times \frac{{V_{\max ,N} \times {\mathrm{C}}_N \times {\mathrm{C}}_B}}{{K_{M,N} + {\mathrm{C}}_N}}$$

(58)

$$f(\theta ) = {\mathrm{max}} [0.25,{\mathrm{min}} (1, - 1.1 \times \theta ^2 + 2.4 \times \theta - 0.29)]$$

(59)

$$f({\mathrm{pH}}) = e^{\frac{{ - ({\mathrm{PH}} - {\mathrm{PH}}_{{\mathrm{opt}}}^2)}}{{{{{\rm{PH}}^2_{{\rm{sen}}}}}}}}$$

(60)

$$f({\mathrm{Clay}}) = 1-0.75 \times {\mathrm{Clay}}$$

(61)

where T is the soil temperature (K); T_ref is the reference temperature, set as 285.15 K; Ea is the activation energy of substrate decomposition (kJ mol⁻¹) and was 50 and 30 kJ mol⁻¹ for V_max,N and K_M,N, respectively [30, 49]. R is the ideal gas constant (0.008314 kJ mol⁻¹ K⁻¹); V_max,0 and K_M,0 are the pre-exponential coefficient, respectively, and were estimated to be 240 (mg C mg⁻¹ MBC day⁻¹) and 300 (mg C g⁻¹ soil), respectively, based on the four observation studies. θ is soil moisture (%); pH_opt is the optimal pH for substrate decomposition and was set as 6 [21]. pH_sen is the sensitivity parameter of substrate decomposition and was set as 1.66 according to the average value in a previous study [48]. Clay is the soil clay content (%). MAT is the mean annual temperature (°C). MAP is the mean annual precipitation (mm).

Data sources

Input rate I was from the community land model (CLM) spatial output data, and CLM5.0 data, which is from the National Center for Atmospheric Research (NCAR) (https://www.earthsystemgrid.org). The file of monthly litterfall rate was clm50_r270_1deg_GSWP3V1_iso_newpopd_hist.clm2.h0.LITFALL.185001-201412.nc and the file of monthly litter heterotrophic respiration was clm50_r270_1deg_GSWP3V1_iso_newpopd_hist.clm2.h0.LITTERC_HR.185001-201412.nc. We chose the subset between 2000 and 2014 to obtain the mean annual data.

Climate data MAT and MAP were from the worldclim database (http://www.worldclim.org) [50]. Soil moisture data were from the gap-free global annual soil moisture for 1991–2018 [51]. Soil temperature was derived from NCEP/NCAR 40 years reanalysis data [52]. Soil pH, soil clay content, and soil bulk density were obtained from harmonized data set of derived soil properties for the world (WISE30sec) [53]. SOC content was from a global 3D soil information system [54], and the average value of the top 0–15 cm was used. Global biomes were based on the classification by Olson and Dinerstein [55], with the associated spatial data available at http://maps.tnc.org/files/shp/terr-ecoregions-TNC.zip. All data were resampled to a spatial resolution of 0.125° × 0.125° for calculation.

Results

Model simulation

Generally, the decomposition curve of ¹³C-labeled microbial necromass C can be well simulated by both MIND and FOND models (Fig. 2 and Supplementary Fig. 4). The R² of the regression of the recovery in soil ¹³C of the MIND simulation results with observations was >0.9, the RMSEs were between 0.88 and 2.66, while the R² of the regression of FOND simulation results with measured data was >0.68, and the RMSEs were between 1.90 and 4.93 (Table 1).

Fig. 2: Comparison of modeled (lines) and observed (dots) recovery of ¹³C in soil organic C pool and recovery of ¹³C in respired CO₂ using data from four decomposition experiments of ¹³C-labeled microbial necromass.

A–D and E–H represent the simulation results of MIND and FOND using data sets 1–4, respectively. Shaded areas denote the simulated ranges from 2.5th to 97.5th percentiles (i.e., 95% range).

Table 1 Comparison of modeled recovery of ¹³C in soil with measured data in four studies using four different models.

MIND model showed that the recovery of ¹³C in MBC experienced a rapid growth stage at the beginning, and then gradually declined for the four data sets (Fig. 3 and Supplementary Fig. 5). Additionally, 10–25% of the microbial necromass C was quickly transferred to MBC (Fig. 3), and then was transferred to the mineral-associated pool of microbial necromass gradually as microbes died. In contrast, the FOND model showed that the recovery of ¹³C in MBC of four data sets increased only slightly and <5% of the MNC was transferred to MBC (Fig. 3).

Fig. 3: Modeled results of the recovery of ¹³C in microbial biomass carbon (MBC), dissolved organic carbon (DOC), the fast pool of microbial necromass carbon (C_NF) and the mineral-associated pool of microbial necromass carbon (C_N-MAOM) using data from four decomposition experiments of ¹³C-labeled microbial necromass.

A–D and E–H represent the simulation results of MIND and FOND using data sets 1–4, respectively. Shaded areas denote the simulated ranges from 2.5th to 97.5th percentiles (i.e., 95% range).

In all the four studies, the fast pool of MNC (C_NF) was simulated to be quickly decomposed within 100 days (Fig. 3). The difference among different studies was that the proportion of MBC transferred to the fast and mineral-associated pools of necromass was different (Table 2). The MIND model simulated that the recovery of ¹³C in C_NF still had a little at the later stage, while the FOND model simulated that the recovery of ¹³C in C_NF was almost 0 at 50–100 days of decomposition (Fig. 3).

Table 2 Model parameters estimated by the lsqnonlin method in MIND and FOND models.

The temporal variation of simulated mineral-associated pool of labeled MNC (C_N-MAOM) was relatively small in all four studies (Fig. 3). Because of a new round of microbial biomass turnover, MIND model simulated that the recovery of ¹³C in C_N-MAOM exhibited a slightly increasing trend, while FOND model simulated that the recovery of ¹³C in C_N-MAOM did not change much with decomposition time.

Both models estimated that C_N-MAOM had very slow decomposition rate, and the decomposition rate of the fast pool of microbial necromass C was at least 1000 times faster than that of the mineral-associated pool of microbial necromass C (Table 2, R_MAOM-F being 1.08 × 10⁻⁷ to 8.26 × 10⁻³ for MIND model, and 6.14 × 10⁻⁵ to 1.89 × 10⁻³ for FOND model), which were different from the original estimations in these four studies.

Because CORPSE model is the only model that simulates “dead microbes” in current process-based models, we compared our estimated parameters to theirs (Supplementary Table 4) and found that the mean residence time of the mineral-associated microbial necromass estimated by our models (24–67 years) were similar to theirs (45–75 years) except for study 4.

About 44–99% of the microbial necromass C was transferred to the fast necromass C pool simulated by MIND model (Table 2, f_BNF). About 35–88% of the microbial necromass C was transferred to the fast necromass C pool simulated by FOND model. Because a portion of microbial necromass C was transferred to DOC in FOND model, but remained in C_NF and C_N-MAOM in MIND model, it was reasonable that this f_BNF was less estimated by FOND than by MIND model.

Model validation

Based on LOOCV and observed recovery of ¹³C in respired CO₂, modeled and observed SOC pool and MBC pool at steady state in four studies, we found that both models had relatively high accuracy (Figs. 4 and 5 and Supplementary Figs. 6 and 7). The validation results of short-term studies studies 1 and 2 based on MIND and FOND were basically consistent. For the long-term studies study 3 and 4, their errors were greater than those of short-term ones, and the errors of FOND were larger than those of MIND. This result suggested that MIND had better simulation ability and our model had larger errors when simulating long-term observation data (Figs. 4 and 5). However, it is also possible that the observation data themselves had larger errors in the long-term experiments due to the changes in soil characteristics after long-term incubation. SOC was better simulated by both models than MBC (Supplementary Fig. 6).

Fig. 4: The leave-one-out cross-validation test results of MIND and FOND models.

The y-axis is the absolute value of the difference between simulation and observation. Box shows the upper and lower quartile with the line inside representing median. Whiskers show upper and lower extremes. Studies 1–4 are four decomposition experiments of ¹³C-labeled microbial necromass in different ecosystems and the details of the site information of studies 1–4 are in Supplementary Table 2.

Fig. 5: The regression analysis between modeled and observed recovery of ¹³C in respired CO₂ for the four decomposition experiments of ¹³C-labeled microbial necromass.

A–D represent the simulation results of studies 1–4, respectively.

In the four studies simulated by the two models, the percentage of MBC out of total SOC was between 1 and 1.6%, and the percentage of MNC out of total SOC was between 10 and 27% (Table 3 and Supplementary Table 5) at steady state. MIND and FOND models estimated different values for the proportion of MNC for study 1, but estimated similar values for the rest three studies.

Table 3 Proportion of microbial biomass carbon (MBC) and microbial necromass carbon (MNC) out of total soil organic carbon (SOC) under steady state simulated by MIND and FOND models using the four experiment data sets (studies 1–4) and the comparison to observation data.

Long-term and large-scale simulation results

Both models showed that the remaining microbial necromass ¹³C was within 0.15% after 1000 years of decomposition, and remaining microbial necromass ¹³C was within 20% after 100 years of decomposition (Fig. 6 and Supplementary Fig. 8). FOND model estimated that the averaged mean residence time (MRT) of necromass C was 11.1, 23.7, 21.6, and 64.5 years for studies 1–4, respectively.

Fig. 6: Temporal variations of necromass ¹³C recovery in soil after 1000 years of incubation (without new C inputs) for the four studies simulated by MIND and FOND models.

A, B are the recovery of ¹³C in soil simulated by MIND and FOND, respectively; C, D are the recovery of ¹³C in respired CO₂ simulated by MIND and FOND, respectively.

We estimated the average MNC content (0–15 cm) of 13 biomes and calculated the proportion of MNC out of total SOC content (Supplementary Table 6). “Montane grasslands and shrublands” had the highest MNC/SOC (51%) and “Tropical and subtropical broadleaf forests” had the lowest MNC/SOC (2%). Overall, the proportion of MNC in SOC decreased from the cold areas to the tropical areas (Supplementary Table 6).

Comparison of simulation results between models with and without the microbial necromass pool

In order to evaluate if the addition of the microbial necromass pool in the models improved models’ performance, we used another two models without the microbial necromass pool (MIND-O and FOND-O corresponding to MIND and FOND, respectively) to simulate recovery of ¹³C in soil using data of the four studies. Results showed that MIND-O and FOND-O had very poor performance and the uncertainties were extremely high which cannot be presented in the figure (Supplementary Fig. 9 and Supplementary Table 7). Based on the four indexes of R², RMSE and MAE, FOND and MIND models were generally better than the other two models (Table 1).

Sensitivity of model parameters

Sensitivity analysis showed that V_max,P, K_{M, P} were the two dominant parameters of the decomposition of organic matter in MIND model; C_NF was most sensitive to V_max,N, K_M,N, f_BNF, and k_B; C_N-MAOM was most sensitive to R_MAOM-F because it determines the amount of microbial carbon input to the mineral-associated pool (Supplementary Table 8).

In FOND model, C_P was most sensitive to k_P; C_D was most sensitive to k_uptake; C_NF was most sensitive to f_BNF and k_NF because they determine the proportion of microbial necromass input to C_NF and the decomposition rate of C_NF, respectively (Supplementary Table 8).

Discussion

Representation of the microbial necromass pool in process-based models

We found our newly proposed models representing the microbial necromass pool more effectively simulated the four experimental data set of the decomposition of ¹³C-labeled microbial necromass than models without the microbial necromass pool (Table 1). In models without the microbial necromass pool, microbial necromass was considered as a part of total soil organic matter and had the same decomposition pattern as plant-derived organic matter. Plant residues are mainly composed of polysaccharides and lignin, aliphatic biopolymers and tannins; while the cell walls of microbes are mainly composed of homopolysaccharides and heteropolysaccharides [56]. Some fungal cell walls also contain relatively high proportion of proteins, lipids and melanin [56]. Therefore, due to the difference in molecular structures of plant-derived and microbially derived organic matter, the decomposition pattern and rate are also expected to be different. Although the average decomposition rate of microbial necromass in soil was found to be faster than that of plant residues [8], some component of microbial necromass such as some glucosamine, galactosamine, or sorbic acid may accumulate for a long term in soil [56]. For example, melanin has been found to inhibit the decomposition of microbial cell wall [57], and more peptidoglycans in bacterial cell walls usually result in slower decomposition rates [56]. In addition, because microbes generally live on the surface of a particle in soil [11, 58, 59], microbial necromass is easily protected by minerals [11, 32, 60,61,62]. Once the necromass interacts with soil minerals, they can become inaccessible to microbes and have very slow decomposition rates [63, 64]. Microbial necromass has been found to bind to clay particles and promote the formation of micro-aggregates with <50 μm [65], and could be more tightly bound in soil matrix than plant debris [9, 66]. Therefore, the unprotected parts of microbial necromass (i.e., fast pool of microbial necromass) could be decomposed much more rapidly than the protected parts (i.e., mineral-associated necromass pool) [67, 68], showing different pattern from plant residues. The mineral-associated pool of microbial necromass had at least 1000 times slower decomposition rate than the fast pool in our simulation results, which supported this viewpoint. By adding the microbial necromass pool and dividing the pool into fast and mineral-associated components, process-based models can better simulate the decomposition pattern of SOC. Sulman et al. [30], in their model CORPSE, integrated two compartments of microbial necromass, called “unprotected dead microbe C” and “protected dead microbe C.” However, because the “protected dead microbe C” pool was not differentiated from “protected plant residue” and they had the same decomposition pattern and rates, we believe our model is the first attempt to realize this purpose of separately simulating microbial necromass C from plant-derived organic C decomposition processes. Between the two models representing the microbial necromass pool, the MIND model performed better than the FOND model. The reason was because microbial controls on SOC decomposition was better simulated in MIND model due to the consideration of microbial biomass effect on decomposition [3, 5, 21, 22, 69]. The Michaelis–Menten model has been found to better simulate the response of soil respiration to the drying-rewetting cycle than the first-order kinetic model because the consumption of SOC depends on microbial biomass [31]. We believe that by adding the microbial necromass pool to the Michaelis–Menten model, the paradoxical role of soil microbes in regulating SOC as both consumers of (via the Michaelis–Menten equation) and contributors to (via the necromass pool) SOC can be better represented and more accurately simulated. Furthermore, the Michaelis–Menten model may be more sensitive to disturbance such as warming and raising CO₂ [3, 30]. Therefore, under the context of global changes, we believe MIND model considering microbial controls on SOC decomposition could be more effective in the long term.

Temporal variations of ¹³C recovery in MBC were better simulated by MIND model than by FOND model (Fig. 3) because ¹³C recovery in MBC was estimated to be 6.8–7% at day 224 by MIND model in study 1, which is similar to the PLFA measurement results of 11 ± 9% in the original study [60]. The ¹³C recovery in MBC only increased a little in FOND simulation results and remained very low in models without the necromass pool (Fig. 3), both of which seemed unrealistic.

Analysis of model parameters

In order to incorporate our model in Earth System Models, we need to know more on the major model parameters newly introduced in our model such as the decomposition rate of microbial necromass and the microbial CUE of microbial necromass. The V_max and K_M of the Michaelis–Menten equation for SOM decomposition given in previous studies were 0.01 mg C mg⁻¹ MBC h⁻¹ and 250 mg C g⁻¹ soil, respectively [6, 26]. In MIND model, our results showed that V_max,P was generally lower than 0.01 and V_max,N was higher than 0.01 mg C mg⁻¹ MBC h⁻¹. Moreover, the K_M,P and K_M,N in MIND model were above and below previous estimation, respectively, except in study 4 (Table 2). We believe these results are reasonable because SOM is a mixture of plant-derived and microbially derived organic matter and the values of V_max and K_M for the whole SOM should be between the values for plant residue and those for microbial necromass. V_max is the maximal reaction rate and decreases with decreasing substrate quality and quantity [70]. K_M is an indicator of the affinity of an enzyme has for its substrate [71] and an increase in K_M indicates a decrease in overall enzyme function [24]. V_max and K_M for microbial necromass should be higher and lower than V_max and K_M for plant residue, respectively, because most components of microbial necromass are faster than plant residue. In the first-order kinetic model, the decomposition rate of SOC was reported to be 5.6 × 10⁻⁶ mg C mg⁻¹ C h⁻¹ in previous studies [6, 26]. In FOND model, we found the decomposition rate of plant residue was approximately equal to this number but the decomposition of the fast pool of microbial necromass was higher than this value and the decomposition of the mineral-associated pool of microbial necromass was much lower than this value (Table 2, k_P and k_NF). This result further agrees with the viewpoint discussed above that the fast part of microbial necromass can be quickly used while the mineral-associated part of microbial necromass may accumulate for a very long period. R_MAOM-F was the most sensitive parameter for the C_N-MAOM pool (Supplementary Table 8), suggesting the importance of this parameter and should be studied in different soils and under different conditions in the future.

Additionally, we also compared other major parameters existing in previous models. The decomposition rate of microbially derived DOC (k_DN) in our models was estimated to be 0.003–0.005 mg C mg⁻¹ C h⁻¹ (Table 2) and it is similar to the decomposition rate constant of DOC (k_D) in other models, which was estimated to be 0.001 mg C mg⁻¹ C h⁻¹ in previous studies [6, 26]. It was slightly larger in our study, probably because that DOC from microbial necromass is easier to be protected by soil minerals than plant-derived DOC [9]. The turnover rate constant of microbial biomass (k_B) was estimated to be 0.16–3.15 year⁻¹ in our study (Table 2) and it was 2.45 year⁻¹ in the previous modeling studies [6, 26] and 0.2–20 year⁻¹ in previous observation studies [72, 73]. The fraction of fast content in microbial biomass (f_BNF) was estimated to be 0.53 in FOND for study 1 and 0.52 in their original estimation [34]. Very few previous studies estimated this parameter, and the only available study estimated that the average f_BNF was 0.51 based on their modeling of temperature sensitivity of microbial respiration [74]. In a different study, the proportion of microbial necromass C transferred to DOC was estimated to be 0.172 [62], which is lower than our estimation because our fast pool of microbial necromass contains not only DOC, but also other compounds whose decomposition rates may be lower than DOC.

Our estimation of the mean residence time of mineral-associated microbial necromass was 24–67 years for studies 1–3, which is similar to CORPSE model estimates (45–75 years) (Supplementary Table 4), suggesting the stability of microbial necromass once it is protected by soil minerals. For study 4, our estimates of the mean residence time of mineral-associated microbial necromass was very high (533 years), which was because that the protection rate of MNC was very low (k_P in Supplementary Table 4, i.e., only a very small amount of microbial necromass entered the mineral-associated pool because the study was conducted in tropical areas). Therefore, although the decomposition rate of this pool was very low in study 4, the long-term accumulation (proportion) of MNC in SOC was comparable to other studies (Table 3). Overall, our estimated parameters were all within the reasonable range and future studies should be focused on how these parameters vary with temperature and moisture for model development.

Implication for long-term and large-scale modeling

After running the models for 1000 years, some ¹³C still remained in the soil (Fig. 6 and Supplementary Fig. 8) and the average MRT of microbial necromass was 11.1–64.5 years, suggesting the stabilization of microbial necromass. Although this result seems to be surprising, it is consistent with recent observation on microbial necromass decomposition, which found 33.1–39.5% of labeled necromass N remained in soil after 2.2 years of decomposition [8]. Under steady-state conditions, the estimated microbial necromass C accounted for 10–27% of the total SOC, which is at the lower end range of previous estimations, which were 24–80% [11, 15,16,17]. Some previous studies used the content of amino sugars in soil to estimate microbial necromass content and found microbial necromass C accounted for 23.6–58.3% of total SOC [16, 17]. However, the conversion factors used to convert amino sugar to necromass are still under debate. Other previous studies used the relative abundance of fatty acids and amino acids in microbial biomass and estimated the fraction of necromass C to SOC to be about 40% [11]. Fatty acids and amino acids exist not only in microbial necromass, but also in living microbial cells, plant cells, and root exudates [75, 76], and the fraction of necromass C to SOC was probably overestimated based on this method and therefore higher than our estimations. Sulman et al. [30] estimated that “dead microbe carbon” was 0.34% of total SOC based on CORPSE model, but their “dead microbe carbon” pool corresponds to our labile MNC pool because once microbial carbon enters the “protected carbon” pool (which is the sum of chemical resistant carbon, simple carbon, and dead microbe carbon), they are no longer counted as “dead microbe carbon” in their model. Our labile MNC pool was 0.01–0.06% of total SOC, which is comparable to their estimation.

Alternatively, it is possible that the different findings among different studies were due to site difference because the contribution of microbial necromass to soil organic C is affected by many factors, such as soil pH, salinity, clay content, fungal bacterial ratio, environmental stress, ecological type, land use change, soil C:N, microbial C:N, and nutrient status [16, 61, 77]. But it should be noted that of ¹³C-labeled necromass decomposition studies [33,34,35] used for our models only contain ~20 isolates, and therefore the necromass C decomposition rate for the whole microbial community should be examined in future studies.

Our estimation of the proportion of MNC out of total SOC content in 13 biomes suggested that “montane grasslands and shrublands” and “tundra” biomes had highest proportion of necromass C in SOC (31–51%, Supplementary Table 6), which is close to an estimation based on soil amino sugar content in the Tibet Plateau [78]. Generally, these areas have high altitude and low temperature, and the plant productivity is low [79, 80]. Due to the higher microbial CUE [81] and the better adaptability of microbes than plants under extreme conditions, it is reasonable that the proportion of microbial necromass out of total soil organic matter is high. In tropical areas, plant productivity is high and microbial CUE is low [81], microbial necromass C was estimated to be <14% of SOC by our models (Supplementary Table 6). Previous experiments also found that the contribution of microbial residues to SOC decreased from tropical to boreal forests [82]. Although these results at the large spatial scale are preliminary because the estimates of some parameters were based on four observation studies, it can still shed light on the importance of microbial necromass pool and the variations among different biomes.

A very interesting phenomenon is: century-scale projection using MIND showed some oscillation compared to FOND (Fig. 6). Usually more C pools in Michaelis–Menten-based microbial models would reduce such oscillation [26], and Georgiou et al. [6] reduced this oscillation by considering the turnover of microorganisms at the population level. We hope the addition of the necromass pool could also realize it and found that changes in carbon input, temperature and microbial CUE [6, 26] could still induce oscillation of SOC simulation results for both MIND and FOND models although the addition of necromass pool alleviated such oscillation. The amplitude and period of the oscillation depended on microbial CUE, V_max and K_M [83]. Future studies considering the turnover of microorganisms at the population level [6] may reduce the oscillation and further improve model performance.

Limitations and future studies

At present, only four sets of observation data based on the decomposition experiments of the ¹³C-labeled microbial necromass were available to calibrate our models. Although the four studies were conducted in different ecosystems and somewhat represent different conditions, more experiments are needed for better parameterization of the models. In addition, the four observation data were based on both lab incubation experiments (studies 1 and 2) and field experiments (studies 3 and 4) and whether or not the experiments were conducted in the field is important. In the future, more field-based microbial necromass decomposition studies are needed to better parameterize the model and better understand the behavior of microbial necromass decomposition. Especially, experiments under the changing environment are needed to test whether model behavior under climate change is consistent with observation and this will be the immediate next step of our work.

The number of observation on some parameters like microbial CUE has increased tremendously in recent 5 years [46, 84]. More direct measurements of the decomposition parameters for the microbial necromass pool (i.e., the decomposition rate of microbial necromass, the turnover rate of microbial biomass, etc.) will substantially increase the effort of parameterization. Previous studies suggested that the decomposition rates of microbial necromass nitrogen from different microbial groups were not significantly different [8, 32], suggesting that substrate quality may not influence the decomposition rates of microbial necromass. However, soil properties such as soil temperature, soil moisture, and soil clay content may influence these parameters and more experiments are needed to estimate these parameters under different conditions.

Nevertheless, our study provides an independent model estimation of microbial necromass C contribution to stable SOC, supplementary to previous model estimations [15, 85, 86]. Although slighter lower than previous estimations, this significant contribution ratio further proved that it is necessary to incorporate the microbial necromass pool in Earth System Models to better simulate carbon cycling under the context of global changes.