Improving a Biogeochemical Model to Simulate Microbial-Mediated Carbon Dynamics in Agricultural Ecosystems.

Soil microbes drive decomposition of soil organic matter (SOM) and regulate soil carbon (C) dynamics. Process-based models have been developed to quantify changes in soil organic carbon (SOC) and carbon dioxide (CO2) fluxes in agricultural ecosystems. However, microbial processes related to SOM decomposition have not been, or are inadequately, represented in these models, limiting predictions of SOC responses to changes in microbial activities. In this study, we developed a microbial-mediated decomposition model based on a widely used biogeochemical model, DeNitrification-DeComposition (DNDC), to simulate C dynamics in agricultural ecosystems. The model simulates organic matter decomposition, soil respiration, and SOC formation by simulating microbial and enzyme dynamics and their controls on decomposition, and considering impacts of climate, soil, crop, and farming management practices (FMPs) on C dynamics. When evaluated against field observations of net ecosystem CO2 exchange (NEE) and SOC change in two winter wheat systems, the model successfully captured both NEE and SOC changes under different FMPs. Inclusion of microbial processes improved the model's performance in simulating peak CO2 fluxes induced by residue return, primarily by capturing priming effects of residue inputs. We also investigated impacts of microbial physiology, SOM, and FMPs on soil C dynamics. Our results demonstrated that residue or manure input drove microbial activity and predominantly regulated the CO2 fluxes, and manure amendment largely regulated long-term SOC change. The microbial physiology had considerable impacts on the microbial activities and soil C dynamics, emphasizing the necessity of considering microbial physiology and activities when assessing soil C dynamics in agricultural ecosystems.

Introduction

Soils contain more carbon (C) than plant biomass and the atmosphere combined at the global scale (IPCC, 2013; Lehmann & Kleber, 2015). Through soil C sequestration or converting and releasing soil organic carbon (SOC) into atmospheric carbon dioxide (CO2) or methane (CH4), soils exert large impacts on global climate change (IPCC, 2013; Lal, 2004). Moreover, soil organic matter (SOM) contains nutrients that are essential for plant growth. Therefore, SOM dynamics influence the availability of soil nutrients, thereby impacting plant productivity in terrestrial ecosystems (Lal, 2008). Despite considerable research, large uncertainty remains over the spatiotemporal variability and magnitude of the SOC changes and the associated CO2 fluxes in agricultural ecosystems due to complex processes and interactions involved (e.g., Dignac et al., 2017; Lal, 2018). At the same time, there is a considerable interest in promoting soil C sequestration in agricultural soils (e.g., Amelung et al., 2020; Smith, 2004), which will require advances in methods for measurement and monitoring, increased process‐level understanding, improved process‐based modeling of agricultural soils and impacts of farming activities, and expanded and enhanced geospatial databases of soil properties and farming management (e.g., Paustian et al., 2016).

Changes of SOC in agricultural ecosystems ultimately depend on the balance between C inputs (e.g., plant litter, organic manure) and outputs (e.g., SOC losses due to decomposition or erosion), which are affected by farming management practices (FMPs), plant biomass production, and SOM formation and decomposition. Agricultural management practices, such as fertilization, input of crop straw, and organic manure amendment, have been regarded as important regulators of the SOC changes, although the trend and magnitude of the SOC changes in response to agricultural activities are highly variable (e.g., Liu et al., 2014; Maillard & Angers, 2014; Powlson et al., 2012). In addition, SOC changes and CO2 fluxes from decomposition are largely controlled by soil microbial physiology and activity, given that soil microbes play a central role in SOM decomposition, formation, and depletion (Kallenbach et al., 2016; Miltner et al., 2012). However, the responses of microbial‐mediated decomposition processes to agricultural activities and the impacts of microbial physiology on SOC changes remain unclear (Kallenbach et al., 2016; Prommer et al., 2020).

Process‐based biogeochemical SOC models such as DeNitrification‐DeComposition (DNDC), DayCent, and Roth‐C have been developed and applied to quantify SOC dynamics and CO2 fluxes in agricultural ecosystems as well as evaluate their responses to variations in FMPs and environmental factors (e.g., Parton et al., 2015; Smith et al., 1997; Van Wesemael et al., 2010). While many improvements have been made in modeling agricultural SOC dynamics by incorporating farming activities, vegetation, and biogeochemical processes related to C cycling into the framework of biogeochemical models (e.g., Li et al., 1994), key limitations still exist. For example, in most biogeochemical models SOC decomposition is simulated as a first‐order process and is directly proportional to the size of SOC pool. Microbial processes important to SOC decomposition and their associated microbial communities often have not been, or are inadequately, represented in these models, which limit the model's capacity of predicting SOC responses to changes in microbial activities (Georgiou et al., 2017; Sulman et al., 2014). In recent years, several microbial models of SOC decomposition have been proposed or developed to predict soil C dynamics by explicitly representing microbial or enzymatic decomposition of SOC (e.g., Abramoff et al., 2017; Allison et al., 2010; Gao et al., 2020; Guo et al., 2020; Wang et al., 2013; Wieder et al., 2015). Microbial SOC models have incorporated processes to simulate microbial dynamics and feedbacks between SOC changes and microbial activities. However, they usually lack detailed processes to explicitly simulate plant growth and litter production, and have not incorporated FMPs (e.g., harvest, tillage, fertilization, irrigation, organic manure amendment) that influence soil C dynamics (e.g., Sulman et al., 2018), limiting their applications in agricultural ecosystems. Furthermore, there are large uncertainties regarding their behaviors in modeling SOC changes and CO2 fluxes due to very limited model testing against field‐scale observations.

The process‐based biogeochemical model, DNDC, has incorporated a relatively complete suite of biophysical and biogeochemical processes and FMPs, which enables it to simulate vegetation growth, transport and transformations of C and nitrogen (N), SOC changes, and greenhouse gas (GHG) fluxes in terrestrial ecosystems (e.g., Deng et al., 2018, 2020; Gilhespy et al., 2014; Giltrap et al., 2010; Li, 2000; Li et al., 1992a, 2012). However, DNDC does not explicitly simulate microbial or enzymatic decomposition of SOM, and therefore is unable to reproduce microbial regulations on soil C cycle. In this study, our objectives were: (a) to develop a microbe‐driven SOM decomposition model, based on the DNDC model, to predict microbial activities and C dynamics, and evaluate the model's performance in quantifying CO2 fluxes and SOC changes from agricultural ecosystems and (b) to explore the impacts and evaluate the relative importance of FMPs, SOM, and microbial physiology on soil microbial activity, SOC change, and CO2 flux.

Materials and Methods

The DNDC Model

DNDC is a process‐based biogeochemical model developed for quantifying ecosystem C sequestration and the exchange of C and N gases between terrestrial ecosystems and the atmosphere (Li, 2000; Li et al., 1992a, 1992b). The model has been extensively evaluated against datasets of SOC changes and GHG fluxes measured worldwide (Gilhespy et al., 2014; Giltrap et al., 2010). DNDC is comprised of six interacting sub‐models: soil climate, plant growth, decomposition, nitrification, denitrification, and methanogenesis. The soil climate, plant growth, and decomposition sub‐models convert the primary drivers, such as climate, soil properties, vegetation, and anthropogenic activity, into soil environmental factors, such as soil temperature and moisture, pH, redox potential, and concentration of substrates for the simulated biogeochemical processes. The nitrification, denitrification, and methanogenesis sub‐models simulate C and N transformations that are mediated by soil microbes and controlled by soil environmental factors and substrate concentrations (Li, 2000; Li et al., 2012).

In DNDC, net ecosystem CO2 exchange (NEE) is calculated as the difference between net primary production (NPP) and soil microbial heterotrophic respiration. NPP is simulated at a daily time step by considering impacts of several environmental factors (e.g., solar radiation, air temperature, soil moisture, and N availability) on plant growth (Deng et al., 2014). The impacts of solar radiation, soil moisture, or N availability on NPP are calculated based on the potential photosynthetically active radiation (PAR), water, or N required by optimum plant growth and the availability of these factors. The model simulates the production of plant biomass and litter, and incorporates the residue litter into the SOM pools. Soil heterotrophic respiration is simulated by tracking the decomposition of SOM. The model divides SOM into four organic pools: litter, living microbes, humads (i.e., active humus), and passive humus. Each pool is further divided into labile and resistant sub‐pools with specific C:N ratios and decomposition rates (Figure 1). In DNDC, the size of the living microbes pool does not affect the decomposition rates of other SOM pools. Decomposition of each SOM pool depends on its specific decomposition rate and pool size, as well as soil thermal, moisture, and mineral N conditions (Li et al., 1994, 2012). When organic matter, such as litter or organic manure, has been applied into the soil, the model partitions it into different litter and/or humads pools based on its C:N ratio (Figure 1; Li et al., 1994, 2012). Organic C is converted to CO2, microbial biomass, humads, and finally passive humus through decomposition. Labile C is gradually lost and resistant C becomes relatively more abundant in the soil during decomposition of exogenous organic matter (Figure 1). The model also simulates soil N dynamics by tracking a series of biogeochemical reactions: decomposition (mineralization), microbial immobilization, plant uptake, ammonia volatilization, ammonium adsorption, nitrification, denitrification, nitrate leaching, and fluxes of N gases (i.e., NH3, NO, N2O, and N2). Soil mineral N concentration is calculated at a daily time step and is a factor regulating SOM decomposition (Li et al., 1992a). Further details regarding the DNDC structure, inputs, and outputs, as well as the physical, chemical, and biogeochemical processes incorporated into the model's framework, are available in Gilhespy et al. (2014), Li (2000), and Li et al. (2012).

Figure 1

Structure of the soil organic matter (SOM) decomposition module in the modified DNDC model. Boxes represent simulated organic matter pools, solid arrows represent fluxes among the pools, and P1 to P10 represent simulated processes. The model simulates decomposition of six SOM pools (green boxes). SOM pool decomposition rates are directly controlled by the soil enzyme concentration, which is regulated by dynamics of soil living microbes, dissolved organic carbon, and dissolved organic nitrogen (collectively DOM in the figure).

Modification of DNDC

In this study, we improved the DNDC model by explicitly simulating the microbe‐driven SOM decomposition. The improved model simulates cycling and storage of multiple SOM pools by explicitly simulating microbial and enzyme dynamics and their controls on SOM decomposition, in addition to simulating impacts on C dynamics of crop growth and common FMPs in agricultural ecosystems (Figure 1).

SOM Pools and Decomposition

The model primarily simulates decomposition of six SOM pools (Figure 1). The very labile, labile, and resistant litter pools are organic matter from plant residue or organic manure amendments, and have different C:N ratios and relatively high, moderate, and low potential decomposition rates, respectively. The other three SOM pools are labile humads, resistant humads, and passive humus. The humads pools receive products from SOM decomposition or microbial death, and have different potential decomposition rates as well. The model simulates decomposition of these six pools, with the decomposition rate of each pool directly regulated by soil enzyme concentration (Equation 1, below, and Equations (A1), (A2), (A3), (A4), (A5), (A6), (A7) in Appendix A). During decomposition, decomposed litter and SOM are primarily converted into dissolved organic carbon (DOC) and dissolved organic nitrogen (DON). The DOC and DON can be assimilated by soil microbes to form microbial biomass, released as CO2 and mineralized into ammonium (NH4‐N), respectively, and/or leached out of the soil. The model also explicitly simulates dynamics of soil microbial and enzyme pools, and dead microbes and decayed enzyme are major components of humads and an important source of DOC and DON (Figure 1). In addition, the model considers C flows from humads to humus to represent humification of the humads pools (Li et al., 1992a; Molina et al., 1983).

The basic concept underlying a microbial‐driven decomposition model is that the decomposition of SOM is controlled by the activity of exoenzymes, instead of simply concentration of SOM (Schimel & Weintraub, 2003). Specifically, the rate of decomposition of each SOM pool is simulated using the Arrhenius and Michaelis‐Menten equations (e.g., Allison et al., 2010; Wang et al., 2013); with the rate directly controlled by concentrations of soil enzyme and each SOM pool, soil temperature and moisture, and clay content (Equation 1). where D SOM is the decomposition rate (mg C h−1) of each SOM pool, VmaxSOM is the maximum rate of SOM decomposition at reference temperature (20°C) when SOM concentration is not limiting, Enzyme is the aggregate concentration of soil exoenzymes (mg C g−1) that is simulated based on soil microbial activities (Section 2.2.3), SOM is the concentration of each SOM pool (mg C g−1), K SOM is the Michaelis half‐saturation constant for decomposition of each SOM pool, and f(W) and f(ST) are functions calculating the impacts on SOM decomposition of soil water content and clay content, respectively (Equations A9 and A10). The VmaxSOM values are different for each SOM pool (Table 1), and the maximum rate of SOM decomposition under a specific temperature is calculated using the Arrhenius equation (f(T); Equation A8).

Table 1

Model Parameters for Simulating Microbe and Enzyme Dynamics, and SOM Decomposition

Parameter	Description	Values	Sources
Parameters selected for sensitivity analysis
VmaxDOCuptake	Maximum uptake rate of DOC at a reference temperature of 20°C, mg C (mg microbe C)−1 hr−1	0.24	Huang et al. (2018)
R MicrobeMaintenance	Rate of microbial maintenance respiration, h−1	0.002	Calibrated a
R MicrobeDeath	Microbial turnover rate, h−1	0.002	Allison et al. (2010)
CAE	Potential microbial carbon assimilation efficiency	0.6	Sinsabaugh et al. (2013)
F MICtoDOC	Fraction of dead microbial biomass that allocated to DOC	0.5	Allison et al. (2010)
			Wang et al. (2013)
R EnzymeProduction	Enzyme production rate, h−1	10–5	Calibrated a
R EnzymeDecay	Enzyme decay rate, h−1	10–3	Abramoff et al. (2017)
VmaxLitter_vl	Maximum decomposition rate of very labile litter at reference temperature, mg C (mg Enzyme C)−1 hr−1	81	Wang et al. (2012)
VmaxLitter_l	Maximum decomposition rate of labile litter at reference temperature, mg C (mg Enzyme C)−1 hr−1	81	Wang et al. (2012)
VmaxLitter_r	Maximum decomposition rate of resistant litter at reference temperature, mg C (mg Enzyme C)−1 hr−1	5.6	Calibrated b
VmaxHumads_l	Maximum decomposition rate of labile humads at reference temperature, mg C (mg Enzyme C)−1 hr−1	13	Calibrated b
VmaxHumads_r	Maximum decomposition rate of resistant humads at reference temperature, mg C (mg Enzyme C)−1 hr−1	1.3	Calibrated b
VmaxHumus	Maximum decomposition rate of humus at reference temperature, mg C (mg Enzyme C)−1 hr−1	0.43	Calibrated b
Other parameters
Ea Litter_vl	Activation energy for decomposition of very labile litter in the Arrhenius equation, KJ mol−1	37	Wang et al. (2012, 2013)
Ea Litter_l	Activation energy for decomposition of labile litter in the Arrhenius equation, KJ mol−1	37	Wang et al. (2012, 2013)
Ea Litter_r	Activation energy for decomposition of resistant litter in the Arrhenius equation, KJ mol−1	53	Wang et al. (2012, 2013)
Ea Humads_l	Activation energy for decomposition of labile humads in the Arrhenius equation, KJ mol−1	47	Wang et al. (2012, 2013)
Ea Humads_r	Activation energy for decomposition of resistant humads in the Arrhenius equation, KJ mol−1	47	Wang et al. (2012, 2013)
Ea Humus	Activation energy for decomposition of humus in the Arrhenius equation, KJ mol−1	53	Wang et al. (2012, 2013)
K Litter_vl	Michaelis half‐saturation constant for decomposition of very labile litter, mg C g−1 soil	1.20	Calibrated c
K Litter_l	Michaelis half‐saturation constant for decomposition of labile litter, mg C g−1 soil	1.20	Calibrated c
K Litter_r	Michaelis half‐saturation constant for decomposition of resistant litter, mg C g−1 soil	12.0	Calibrated c
K Humads_l	Michaelis half‐saturation constant for decomposition of labile humads, mg C g−1 soil	12.0	Calibrated c
K Humads_r	Michaelis half‐saturation constant for decomposition of resistant humads, mg C g−1 soil	60.0	Calibrated c
K Humus	Michaelis half‐saturation constant for decomposition of humus, mg C g−1 soil	100.0	Calibrated c
K DOC	Michaelis DOC half‐saturation constant for DOC uptake, mg C g−1 soil	0.30	Calibrated d
K_O_2	Michaelis O2 half‐saturation constant for DOC uptake, mmol cm−3	0.00015	Li (2016)

R MicrobeMaintenance and R EnzymeProduction were calibrated within their uncertainty ranges of 10−4 to 8 × 10−3 hr−1 and 10−5 to 8 × 10−5 hr−1, respectively (He et al., 2015).

VmaxLitter_ was calibrated within the uncertainty range of 0.2–33.0 mg C (mg Enzyme C)−1 hr−1, and VmaxHumads_, VmaxHumads_, and VmaxHumus were calibrated within the uncertainty range of 0.05–22.0 mg C (mg Enzyme C)−1 hr−1 (Wang et al., 2013).

K Litter_, K Litter_, K Litter_, KmaxHumads_, KmaxHumads_, and KmaxHumus were calibrated within the uncertainty ranges of 0.01–100 mg C g−1 soil for decomposition of litter pools and 0.01–500 mg C g−1 soil for decomposition of humads or humus pool (Huang et al., 2018; Wang et al., 2013).

K DOC was calibrated within the uncertainty range of 0.14–0.38 mg C g−1 soil (Wang et al., 2013).

DOC and DON Dynamics and Soil Microbial C and N Uptakes

Soil DOC and DON are produced through decomposition of each SOM pool, from death of microbes and decay of enzymes, added through root exudation, and/or through amendments of organic manure (e.g., slurry manure). They can be taken up and used by soil microbes for biomass production and respiration (Schimel & Weintraub, 2003), or leached out of soils. In addition, the model simulates DOC consumption by denitrification (Li et al., 1992a) and methanogenesis (Deng et al., 2017) as well as DON mineralization. Specifically, DOC and DON dynamics are simulated as (Equations 2 and 3): where DOCManure and DONManure are DOC and DON added through organic manure and are calculated based on model inputs of manure amendment rate, manure type, and C:N ratio of manure (Li, 2016; Li et al., 2012); DOCRoot is the DOC added through root exudation, which is simulated as 45% of the C transferred to roots from photosynthetic production (Zhang et al., 2002), DONRoot is the DON added through root exudation, which is calculated as the DOCRoot divided by the C:N ratio of root exudation. D Litter , D Humads, and D Humus are decomposition of litters, humads, and humus, respectively (Equations (A1), (A2), (A3), (A4), (A5), (A6), (A7)), DeathMicrobe and D Enzyme are the microbe death and enzyme decay at each time step, CNLitter, CNHumads, CNHumus, CNMicrobe, and CNEnzyme are the C:N ratios of litters, humads, humus, soil microbes, and enzyme respectively, F HumadsToDOC and F MICToDOC are the fractions of decomposed humads and dead microbes that allocated to DOC, respectively, DOCUptake and DONUptake are the DOC and DON uptakes at each time step, DOC and DOC are the DOC consumed through denitrification (Li et al., 1992a) and methanogenesis (Deng et al., 2017), respectively, and DOCLeaching and DONLeaching are the DOC and DON leached out of soils, which are simulated based on subsurface drainage flows and DOC and DON concentrations (Deng et al., 2011; Li et al., 2006). Note that DOC and DON flux rates in Equations 2 and 3 are regulated by soil environmental factors and are zero under some conditions (e.g., DOCM is zero in uplands).

DOC uptake by soil microbes is modeled based on the Arrhenius and Michaelis‐Menten equations (Abramoff et al., 2017; Allison et al., 2010) with the potential uptake rate (DOCPUptake) controlled by the soil temperature, soil living microbes, and soil DOC and oxygen concentrations (Equation 4). where VmaxDOCUptake is the maximum rate of DOC uptake at reference temperature (20°C) when DOC and oxygen concentration are not limiting, DOC and O2 are the soil DOC and oxygen concentrations, respectively, and K DOC and are their corresponding Michaelis half‐saturation constants. The maximum rate of DOC uptake under a specific temperature is calculated using the Arrhenius equation (f(T); Equation A8).

The actual DOC uptake is constrained by soil DOC and DON availability, and is calculated as (Equation 5):

After determining the DOCUptake, the actual DON uptake is calculated as DOCUptake/CNDOC, where CNDOC is the C:N ratio of DOC. Please note that the actual DON uptake is constrained by soil DON availability because the CNDOC is calculated using the DOC and DON concentrations at each time step.

Dynamics of Microbes and Enzyme

The model explicitly simulates soil microbial dynamics by considering growth and mortality of microbes, and production of enzymes (Equation 6): where GrowthMicrobe is the growth of microbial biomass, DeathMicrobe is the microbe mortality, P Enzyme is the enzyme production.

Microbial growth is controlled by DOC uptake, microbial maintenance respiration, microbial carbon assimilation efficiency (CAE; defined as the maximum fraction of DOC uptake that can be allocated to microbial biomass), and availability of total dissolved nitrogen (DN) in soils (Schimel & Weintraub, 2003). If DOC uptake exceeds maintenance respiration, the excess C can be assimilated by soil microbes. After meeting the demand of microbial growth, the element in excess is released, by growth or overflow respiration for C or by N mineralization (Manzoni et al., 2012; Schimel & Weintraub, 2003). The actual growth of microbial biomass is also constrained by the DN availability. Specifically, the model simulates the microbial growth as: where M Maintenance is the microbial maintenance respiration, DN is total dissolved nitrogen, including both DON and mineral N from either SOM mineralization or inputs through FMPs (e.g., fertilization), is simulated through the N biogeochemical reactions in DNDC, and N Demand is the potential N required for microbe growth.

The microbial maintenance respiration is simulated as a constant fraction (R MicrobeMaintenance) of living microbial biomass (Equation 8):

The potential N required for microbial growth is determined by potential microbial growth without DN limitation and the microbial C:N ratio (Equation 9):

The mortality rate of soil microbes is modeled as a density‐dependent process of living biomass (Abramoff et al., 2017) (Equation 10):

The enzyme concentration directly controls the decomposition of each SOM pool, and is determined by the enzyme production (P Enzyme) and decay (D Enzyme) at each time step (Equation 11).

By following Allison et al. (2010), we modeled enzyme production as a constant fraction (R EnzymeProduction) of living microbial biomass (Equation 12) and modeled enzyme decay as a first‐order process with a rate constant R EnzymeDecay (Equation 13).

Field Data

SOC Data From the Rothamsted Experiment

To evaluate the performance of the modified DNDC in simulating long‐term SOC changes, we used SOC observations from the Rothamsted classical broadbalk winter wheat experiment (Powlson et al., 2012). The experiment started in 1843, and field observations at Harpenden in Hertfordshire in Southern England (N51°48′36″, W0°22′30″) included SOC changes of winter wheat croplands during 1843–2010. The soil is classified as Chromic Luvisol in the FAO Soil Classification, or Aquic Paleudalf in the U.S. Soil Taxonomy. The soil (0–23 cm depth) has a pH of 7.0–7.5, and its texture is 25% sand, 50% silt, and 25% clay (www.era.rothamsted.ac.uk/Broadbalk/bbksoils).

Winter wheat has been grown continuously since 1843, except for occasional fallows during the period from 1926 to 1966 to control weeds. Wheat was sown in autumn and harvested in summer of the next year. We evaluated the model using field data of SOC changes under three treatments: control, conventional fertilizer application (NPK), and farmyard manure amendments (FYM). The control treatment has received no fertilizer or organic manure since 1843. The NPK is local conventional practice of applying 144 kg N ha−1 yr−1, 35 kg P ha−1 yr−1, and 90 kg K ha−1 yr−1. The FYM treatment has received farmyard manure at a rate of 35 metric ton ha−1 yr−1 (about 2,800 kg C ha−1 yr−1) since 1843 (Table 2). Changes in SOC (0–23 cm), climate, soil properties, and FMPs were recorded during the experimental period. Long‐term SOC changes under the three treatments indicated that both application of synthetic fertilizers and amendment of organic manure caused increasing SOC. After about 150 years, the SOC stock (0–23 cm) at the NPK plot was about 0.25 times (25%) higher than the control plot, and the SOC stock (0–23 cm) at the FYM plot was about 2.5 times higher than the control plot (Johnston et al., 2009). The long‐term SOC data provided a unique data set for assessing the impacts of fertilization and manure management on SOC changes and evaluating the model's capability in predicting the SOC responses to these different FMPs.

Table 2

Primary Model Input Parameters for the Rothamsted and Oklahoma Winter Wheat Sites

Model parameters a	Rothamsted	Oklahoma
Soil bulk density, g cm−3	1.25	1.3
Clay content, %	25	21
Initial SOC content d , %	1.0	1.6
pH	7.7	5.9
Fraction of crop residue return d	0.5	1.0
N fertilization b ,  d , kg N ha−1	0 or 144	72 or 62
Organic manure application c ,  d , kg C ha−1	0 or 2,800	0
C:N ratio of manure d	12.5	None

Model input parameters were from field records.

Amount of synthetic N fertilizer applied into fields was 0 kg N ha‐1 under the control and FYM treatments and 144 kg N ha−1 under the NPK treatment at the Rothamsted site, and was 72 kg N ha−1 in the 2014 to 2015 crop season and 62 kg N ha−1 in the 2015 to 2016 crop season at the Oklahoma site.

Amount of organic manure applied on Rothamsted fields was 0 kg C ha−1 under the control and the NPK treatment and 2,800 kg C ha−1 under the FYM treatment.

The parameters selected for sensitivity analysis.

NEE Data From Winter Wheat Eddy Flux Tower Sites in Oklahoma

Field data used to evaluate the simulations of NEE from the modified DNDC model were measured from a winter wheat cropland at the US Department of Agriculture‐Agricultural Research Service (USDA‐ARS), Grazinglands Research Laboratory (GRL, N35°34′7, W98°3′21″), in El Reno, Oklahoma. This area has a temperate continental climate with an average air temperature of 14.9°C and an average annual rainfall of 860 mm from 1971 to 2000 (Fischer et al., 2012). Field data included daily NEE measured during October 2014 to September 2016 (Bajgain et al., 2018). The soil of the study field is characterized as deep, well‐drained, loam with clayey or loamy subsoil. During the study period, wheat was sown in September and was terminated in June of the next year, and the field was kept fallow during the summer seasons with weed control by tillage and herbicide application. The field received C inputs through wheat residue return and N inputs through synthetic N fertilizations during the study period (Table 2). Continuous CO2 fluxes were measured using an eddy covariance (EC) tower, and daily and annual values of NEE (positive values represent net CO2 fluxes into the atmosphere and negative fluxes represent net CO2 uptake by the wheat field) were used for evaluating the model. In addition to NEE measurements, climate, soil properties, and FMPs for the wheat field were also recorded.

Model Application

Model Evaluation

The modified DNDC was run for the Rothamsted and Oklahoma winter wheat sites to evaluate the model's performance in simulating long‐term SOC changes and seasonal to annual CO2 fluxes, respectively. The Rothamsted site was simulated from 1840 to 2010, and the Oklahoma site was simulated from 2014 to 2016. Daily meteorological data required for driving the model (i.e., air temperature and precipitation for Rothamsted and air temperature, precipitation, and solar radiation for Oklahoma) were obtained from on‐site measurements for both sites. The primary soil input parameters, including soil texture, clay fraction, bulk density, pH, initial SOC content, were also determined using on‐site records (Table 2). The input parameters of FMPs, including planting and harvest dates, tillage, fertilization, wheat residue return, and amendments of organic manure, were estimated from field records (Table 2).

New model parameters are needed to simulate microbial and enzyme dynamics as well as decomposition of each SOM pool. These parameters were set either to literature values or calibrated against the field observations of SOC changes under NPK in the Rothamsted classical broadbalk winter wheat experiment through the trial and error method (Table 1). The calibrated parameters included rates of microbial maintenance respiration (R MicrobeMaintenance) and enzyme production (R EnzymeProduction), maximum decomposition rates of resistant litter (VmaxLitter), labile humads (VmaxHumads), resistant humads (VmaxHumads), and humus (VmaxHumus) at reference temperature, and Michaelis half‐saturation constants for decomposition of each SOM pool and DOC uptake.

R MicrobeMaintenance and R EnzymeProduction were estimated by calibrating simulated SOC changes with observed SOC changes under the NPK treatment at the Rothamsted site. The maximum decomposition rates at 20°C were estimated as 5.6, 13, 1.3, and 0.43 g C (g Enzyme C)−1 hr−1 for the resistant litter, labile humads, resistant humads, and humus pools, respectively, by calibrating against observed SOC changes under NPK. The values of these parameters were within their reported uncertainties (Wang et al., 2013) and generally reflected the trend of decreased potential decomposition rate from labile humads to resistant litter and humads, and to humus in DNDC (Li et al., 1992a).

The Michaelis half‐saturation constants for decomposition of each SOM pool were estimated after fixing their corresponding maximum decomposition rates. We calibrated the Michaelis half‐saturation constants (Table 1) by considering the reported uncertainties of these parameters (Huang et al., 2018; Wang et al., 2013) and the steady‐state concentrations of each SOM pool by referring to the analysis by Wang et al. (2013). In general, the input parameters described above were primarily determined by calibrating against observed SOC changes under NPK at the Rothamsted site, and the calibrated model was then validated by comparing simulations against observations of SOC changes under the control and FYM treatments at Rothamsted, and NEE at the Oklahoma site. The new model parameters for simulating microbe and enzyme dynamics, and SOM decomposition (Table 1) were identical between the Rothamsted and Oklahoma sites.

Two statistical measures, the relative root mean squared error (RRMSE) and the correlation coefficient (R), were used to quantify the accordance and correlation between model predictions and field observations (Moriasi et al., 2007).

Impacts of FMPs, SOM, and Microbial Physiology on Soil Carbon Dynamics

To investigate the impacts of FMPs, SOM, and microbial physiology on soil microbial activity, CO2 fluxes, and SOC changes, we conducted a series of simulations with the modified DNDC by varying 18 parameters (Tables 1 and 2). These parameters are relevant to FMPs, initial SOC content, maximum SOM decomposition rate, or microbial physiology, and can impact SOC change and CO2 flux from soil heterotrophic respiration (e.g., Six et al., 2006). The baseline was set by referring to the environmental conditions and FMPs under the control treatment of the Rothamsted classical broadbalk winter wheat experiment. We conducted simulations for the years from 1840 to 2010 (around 170 years), with continuous cultivation of winter wheat during these years, which was slightly different with the setting of occasional fallow years under the control treatment. In the baseline simulation, there was no fertilization and manure amendment, and the fraction of returned crop residue was set to 0.5. In order to fully represent potential combinations of the selected parameters, we created alternative scenarios by varying these parameters simultaneously. The variations for the parameters related to FMPs were set generally by referring to the farming management in the Rothamsted experiment, and the ranges were from 0.0 to 1.0 for the fraction of the returned crop residue, 0–144 kg N ha−1 year−1 for N fertilizer application rate, 0–2,800 kg C ha−1 year−1 for organic manure input, and 10 to 15 for manure C/N ratio. For the other parameters, we set two ranges of variation (±20% and ±50% of the baseline value) to represent relatively small and large variations for parameters related to SOM or microbial physiology. These two sets of scenarios were denoted as V20% and V50%, respectively. The changes of all these parameters were randomly picked from their corresponding ranges, assuming uniform frequency distributions. Two thousand scenarios were run for each range of variations by using the Latin Hypercube Sampling strategy (Helton & Davis, 2003). The modified DNDC model was run from 1840 to 2010 with the varied parameters for both V20% and V50% (in total 4,000 runs). Simulated mean annual SOC changes, CO2 fluxes from soil heterotrophic respiration, soil living microbe biomass, and enzyme pool size from 1840 to 2010 under different scenarios were analyzed.

In addition, we evaluated relative importance of the selected parameters to SOC change, CO2 flux from soil heterotrophic respiration, and soil microbial activity by using a global sensitivity analysis, the Smirnov test (Saltelli et al., 2008; Tarantola & Becker, 2016). The Smirnov test can calculate a sensitivity index that quantifies the relative importance of the selected parameters (Tables 1 and 2) to outputs (i.e., SOC change, CO2 flux from soil heterotrophic respiration, soil living microbes, and enzyme pool in this study) by fully exploring the variance space of the parameters. A parameter with higher sensitivity index is more sensitive (or important) than another parameter with lower sensitive index under their corresponding variations (Saltelli et al., 2008).

Results

Model Evaluation

SOC Change

At the Rothamsted site, changes of SOC storage were clearly different among the three treatments with different FMPs (Figure 2). During 1840 to 2010, the observed SOC storage decreased by about 15% under control, slightly increased by about 4% under NPK, and substantially increased by 172% under FYM. The amendment of FYM was identified as an important factor regulating the SOC changes through the Rothamsted experiment (Johnston et al., 2009; Powlson et al., 2012). Current SOC storage under the FYM treatment is substantially higher than under control and NPK, by 220% (78.4 vs. 24.5 Mg C ha−1) and 162% (78.4 vs. 30.0 Mg C ha−1), respectively. SOC storage was also influenced by the N fertilization (Figure 2), which was attributed to the higher organic matter returns through wheat residue in the NPK plot (Johnston et al., 2009). Simulations of SOC change were similar to corresponding field observations for both the NPK treatment used for model calibration and control and FYM treatments used for model validation (Figure 2), with the model capturing the differences of SOC change among the different treatments (Figure 2). From 1840 to 2010, simulated SOC storage decreased by 20% under control, increased by about 9% under NPK, and increased by 161% under FYM. These simulated SOC changes were in agreement with observations for both model calibration and validation treatments. The RRMSE values were 7% for NPK (calibration), and 8% and 5%, respectively, for control and FYM (validation).

Figure 2

Simulated (lines) and observed (circles) SOC (0–23 cm) from 1840 to 2010 under the control, fertilizer (NPK), and farmyard manure (FYM) treatments at the winter wheat field in the Rothamsted Agricultural Station, Harpenden, UK.

NEE

Simulated and observed daily NEE showed similar seasonal patterns at the OK winter wheat site, with net CO2 uptake increasing following winter wheat planting, strong net CO2 uptake for most days during wheat growing seasons from January to May, and net CO2 release following wheat harvest (Figure 3). During summer fallow seasons, net CO2 release occurred on most days, and peaks of CO2 release were simulated and observed following the events of wheat root or residue incorporation and/or tillage. The simulated and observed daily NEE values were significantly correlated in both rotational years from 2014 to 2015 and 2015 to 2016 (P < 0.0001), with R values of 0.78 and 0.82 in the first and second rotational years, respectively. These results suggest that the modified DNDC model captured the seasonal patterns of daily NEE at the Oklahoma wheat field. The model also captured the magnitudes of the observed NEE. Annual total simulated NEE values were 495 and 389 kg CO2‐C ha−1, respectively, in the first and second rotational years, which were close to the corresponding observations of 559 and 397 kg CO2‐C ha−1. The RRMSE values were 13% in the first rotational year and 2% in the second rotational year.

Figure 3

Simulated (lines) and observed (circles) daily net ecosystem exchange (NEE) of CO2 during (a) October 2014 through September 2015 and (b) October 2015 through September 2016, at a winter wheat field in El Reno, OK, USA. The black arrows indicate the dates of roots and residues return and the blue arrows indicate the dates of tillage events. The correlations between the simulated and observed daily NEE were significant for all cases (P < 0.001). The gray lines are simulations of daily NEE from the original DNDC model that does not explicitly simulate microbial decomposition of soil organic matter. The R values between the original DNDC's simulations and observations of daily NEE were 0.70 and 0.77 in the rotational years from 2014 to 2015 and 2015 to 2016, respectively.

Impacts of Microbial Physiology, SOM, and FMPs on Carbon Dynamics

Figure 4a illustrates the simulated average annual SOC changes across the scenarios with different FMPs, SOM property, and microbial physiology. The SOC stock was predicted to be relatively stable under the baseline scenario. The simulated average annual SOC change rate was −42 kg C ha−1 year−1 (i.e., net C loss), and the simulated SOC stock was 14% lower than the initial stock after about 170 years under the baseline scenario (Figure 4a). The average annual change in simulated SOC varied from a net decrease of −119 kg C ha−1 year−1 to a net increase of 782 kg C ha−1 year−1 across the scenarios in V20% due to the variations in FMPs, SOM property, and microbial physiology (Figure 4a). Increasing the variations in initial SOM content, maximum SOM decomposition rate, and microbial physiology from ±20% to ±50% substantially increased the variations of the SOC change, with the simulated average SOC changes ranged between −264 kg C ha−1 year−1 (net decrease) and 1,582 kg C ha−1 year−1 (net increase) across the scenarios in V50% (Figure 4a). Simulated SOC stocks tend to be increasing along with the increases of wheat residue input, organic manure amendment, carbon assimilation efficiency (CAE), microbial turnover rate (R MicrobeDeath), and enzyme decay rate (R EnzymeDecay) (P < 0.01; Table 3), and the decreases of initial SOC stock, maximum decomposition rate of humus (VmaxHumus), fraction of dead microbial biomass that allocated to DOC (F MICtoDOC), and enzyme production rate (R EnzymeProduction) (P < 0.01; Table 3). However, relatively strong correlations (|r| > 0.3) were calculated only for ±50% variations in CAE, VmaxHumus, F MICtoDOC, R EnzymeProduction, and R EnzymeDecay as well as amount of organic manure amendment.

Figure 4

Ranges of simulated (a) average SOC change, (b) average annual soil heterotrophic respiration, (c) average annual mean living microbe mass, and (d) average annual mean enzyme pool size over a 170‐year period under variations in selected parameters representing farming management practices, soil organic matter (SOM) properties, and microbe physiology. Parameters are defined in Tables 1 and 2. For parameters related to SOM or microbial physiology, variation ranges are within ±20% (V20%), or within ±50% (V50%). Bars show min and max values, asterisks indicate the 10th and 90th percentiles, boxes represent the bounds of 25th, 50th (median), and 75th percentiles, squares represent the average values, and blue triangles represent the baseline simulations.

Table 3

Correlation Coefficients Between Model Outputs, Including SOC Change, CO2 Flux, Living Microbes, and Enzyme, and Selected Model Parameters That are Relevant to FMPs, SOM, or Microbial Physiology

	V20% b				V50% b
	SOC change	CO2 flux	Living microbes	Enzyme	SOC change	CO2 flux	Living microbes	Enzyme
VmaxDOCuptake	−0.052	−0.015	−0.009	0.007	−0.051	0.002	0.001	0.012
R MicrobeMaintenance	−0.057	−0.016	−0.343*	−0.241*	−0.069	−0.012	−0.492*	−0.277*
R MicrobeDeath	0.157*	0.017	−0.110*	−0.066*	0.217*	−0.014	−0.209*	−0.102*
CAE	0.245*	−0.078	0.424*	0.323*	0.365*	−0.167*	0.649*	0.390*
F MICtoDOC	−0.242*	0.022	−0.004	0.001	−0.348*	0.083*	0.025	0.020
VmaxLitter_vl	0.017	0.001	−0.004	−0.030	0.030	−0.001	−0.009	−0.033
VmaxLitter_l	−0.001	−0.009	−0.013	0.016	−0.007	−0.012	−0.008	0.021
VmaxLitter_r	0.039	0.021	0.018	0.025	0.037	0.014	0.012	0.020
VmaxHumads_l	0.017	0.032	0.011	0.021	−0.008	0.037	0.002	0.016
VmaxHumads_r	0.014	0.020	0.004	0.016	0.010	0.023	0.004	0.011
VmaxHumus	−0.252*	0.051	0.038	−0.006	−0.356*	0.112*	0.045	−0.018
SOCInitial	−0.115*	0.050	0.056*	0.069*	−0.201*	0.075*	0.061*	0.067*
F Residue	0.198*	0.571*	0.594*	0.443*	0.108*	0.549*	0.317*	0.199*
N fertilizer	0.023	−0.011	−0.012	−0.016	0.030	−0.010	−0.011	−0.016
Manure	0.758*	0.785*	0.486*	0.334*	0.425*	0.756*	0.283*	0.139*
CNManure	0.020	0.020	0.026	0.024	0.022	0.018	0.038	0.024
R EnzymeProduction	−0.224*	0.046	−0.014	0.439*	−0.322*	0.097*	−0.024	0.517*
R EnzymeDecay	0.244*	0.002	−0.032	−0.478*	0.345*	−0.050	−0.063*	−0.577*

The "*" represents significant correlations (P < 0.01) between model outputs and selected parameters.

V20% and V50% denote variation range ±20% and ±50%, respectively, for parameters related to SOM or microbial physiology.

CO2 flux

The simulated average annual CO2 flux from soil heterotrophic respiration was 1,919 kg CO2‐C ha−1 yr−1 across the roughly 170 simulation years in the baseline scenario. The simulations of average annual CO2 flux varied from 1,497 to 5,707 kg CO2‐C ha−1 yr−1 across the scenarios in V20%. Increasing the variations in initial SOM stock, maximum SOM decomposition rate, and microbial physiology from ±20% to ±50% slightly increased the variation range of the average annual CO2 flux by 8%–1,431–5,992 kg CO2‐C ha−1 yr−1 (Figure 4b), implying that the variation of annual CO2 flux from soil heterotrophic respiration was primarily due to variations in FMPs. For example, increasing the organic manure input from 0 to 2,800 kg C ha−1 yr−1 and the fraction of wheat residue return from 0.5 to 1.0 nearly tripled the baseline CO2 flux (5,730 vs. 1,919 kg CO2‐C ha−1 yr−1; Figure 5). Changing other parameters, such as CAE and F MICtoDOC, resulted in smaller variations in the simulated average annual CO2 fluxes in comparisons with the changes induced by modifying C inputs (Figure 5). In addition, simulated annual CO2 fluxes were strongly correlated with the fraction of wheat residue return and organic manure amendments (r > 0.5, P < 0.01; Table 3).

Figure 5

Simulated average SOC change (kg C ha−1 yr−1; yellow boxes) and annual CO2 flux (kg C ha−1 yr−1; blue boxes) from soil heterotrophic respiration across the 170 simulation years under selected scenarios with changes in fraction of wheat residue return and manure C input (different rows), and potential microbial carbon assimilation efficiency (CAE), and fraction of dead microbial biomass that allocated to dissolved organic carbon (F MICtoDOC) (different columns). Values in the white boxes are litter (upper box) and manure (lower box) inputs (kg C ha−1 yr−1) for each case. Simulated change in SOC sequestration was primarily driven by C input, particularly by C input from manure amendments, and regulated by microbe physiology. Simulated SOC stocks tend to increase with increasing CAE and decreasing F MICtoDOC. Note that litter inputs include C from root and are not equal under the scenarios with identical fraction (50% or 100%) of wheat residue return due to differences in wheat productivity.

Microbial and Enzyme Dynamics

The simulated average annual mean living microbial and enzyme pools were 69.0 and 0.69 kg C ha−1, respectively, across the 170 simulation years under the baseline scenario. The average annual mean living microbial and enzyme pools varied from 50.4 to 163.2 kg C ha−1 and 0.41–2.07 kg C ha−1, respectively, across all the scenarios in V20% (Figures 4c and 4d). Increasing the variations in SOM properties and microbial physiology from ±20% to ±50% substantially increased the variations of the simulated living microbial and enzyme pools. The simulated average annual mean living microbial and enzyme pools varied from 35.8 to 286.1 kg C ha−1 and 0.19–4.97 kg C ha−1across the scenarios in V50% (Figures 4c and 4d). Simulated living microbial pools were positively correlated with wheat residue return, amount of organic manure amendments, CAE, and initial SOC stock (P < 0.01; Table 3), and negatively correlated with rates of microbial maintenance respiration and turnover (P < 0.01; Table 3). However, relatively strong correlations (|r| > 0.3) appeared only for wheat residue return, amount of organic manure amendments, CAE, and rate of microbial maintenance respiration. Responses of simulated soil enzyme pools were similar to the responses of the simulated soil living microbes because the enzyme production was modeled as a constant faction of living microbial biomass (Equation 12). In addition, simulated soil enzyme pools were positively correlated with enzyme production rate and negatively correlated with enzyme decay rate (P < 0.01; Table 3).

Sensitivity of SOC Change, CO2 flux, and Soil Microbial and Enzyme Pools

The most sensitive factor controlling the average annual SOC change was the amount of manure amendments, although the simulated SOC changes were also strongly regulated by the maximum decomposition rate of humus and microbial physiological parameters, primarily CAE, F MICtoDOC, and enzyme production and decay rates under the typical changes in FMPs and ±20% variations in SOM properties and microbial physiology (Figure 6a). The simulated CO2 fluxes from soil heterotrophic respiration were predominantly controlled by the rate of C input from manure amendments or wheat residue return at the study site (as demonstrated by high sensitivity index values for manure amendments and fraction of residue return, Figure 6c). The simulated soil living microbial and enzyme pools were largely regulated by both the C input rate from manure amendments or residue return and microbial physiology (Figures 6e and 6g). The microbial physiological parameters exerting relatively large impacts were CAE and R MicrobeMaintenance for the simulated living microbial pool and R EnzymeProduction, R EnzymeDecay, and CAE for the simulated enzyme pool.

Figure 6

Sensitivity of SOC change, soil heterotrophic respiration (CO2 flux), living microbes, and enzyme pool to changes in selected parameters, with (left column) ±20% variation in initial soil organic matter (SOM) stock, maximum SOM decomposition rate, and microbial physiology parameters (V20%), and (right column) ±50% variation in initial SOM stock, maximum SOM decomposition rate, and microbial physiology parameters (V50%). Parameters are defined in Tables 1 and 2. The variable outcomes are more sensitive to parameters with higher sensitivity index values (i.e., they are more important) than to parameters with lower sensitive index values.

Increasing the variations of the initial SOM stock, maximum SOM decomposition rate, and microbial physiology noticeably increased their influences on the simulated SOC changes as well as soil living microbial and enzyme pools. In comparison with the C input rate, CAE and R MicrobeMaintenance exerted larger impacts on the simulated soil living microbial pools and CAE, R MicrobeMaintenance, R EnzymeProduction, and R EnzymeDecay exerted larger impacts on the simulated soil enzyme pools (Figures 6f and 6h). However, their influences on the simulated SOC changes were similar to or smaller than the amount of manure amendments (Figure 6b). In addition, increasing the variations of these parameters only slightly increased their influence on simulated CO2 fluxes from soil heterotrophic respiration (Figures 6c and 6d). Therefore, as expected, C input from organic manure or crop residue played an important role in regulating SOC change and CO2 flux, even under the conditions with relatively large variations in SOM properties and microbial physiology.

Discussion

Model Performance

In this study, we developed a microbe‐driven SOM decomposition model based on a traditional biogeochemical model, DNDC, and then investigated the impacts of FMPs, SOM, and microbial physiology on SOC changes and CO2 fluxes from soil heterotrophic respiration. The developed model provides a framework to predict SOC changes and CO2 fluxes by explicitly simulating both microbial regulations, impacts of plants, and common FMPs on soil C cycling.

The model was evaluated against the observed NEE and long‐term SOC changes from the two winter wheat fields. In comparisons with the field observations, the model generally captured the magnitudes and seasonal dynamics of the daily NEE at the OK wheat field (Figure 3). In particular, the incorporation of the microbial processes on SOM decomposition and DOC uptake appreciably improved DNDC's simulation of the peak CO2 fluxes induced by the wheat residue return (Figure 3) primarily because the modified model can capture priming effects of wheat residue inputs by explicitly simulating microbial activity (Blagodatsky et al., 2010), while the original DNDC could not simulate the priming effects due to missing relevant microbial dynamics. The R values between the improved DNDC's simulations and observations of daily NEE values were 0.78 and 0.82 in the first and second rotational years, respectively, and were higher than the corresponding values of 0.70 and 0.77 calculated based on simulations from the original DNDC. The improvements in the model's performance in simulating CO2 fluxes highlight the importance of representing key microbial mechanisms for simulating seasonal dynamics of CO2 fluxes from agricultural ecosystems. However, we note that adding the new processes also introduced new parameters that are not well‐enough constrained by literature values and require calibration. Our parameter calibration was based only on long‐term SOC simulations from the Rothamsted winter wheat NPK treatment, not the OK winter wheat field.

The modified model also captured the differences of the SOC changes under the three treatments with different C and N inputs at the Rothamsted site (Figure 2), suggesting that the microbe‐driven SOM decomposition model can be used to quantify the long‐term SOC changes for the studied winter wheat fields. However, the performance of simulating long‐term SOC dynamics by the modified model is similar to that of the original DNDC in simulating SOC changes of the same fields (Li et al., 1994; Smith et al., 1997), even though the original version is equipped with first‐order kinetics to simulate SOM decomposition. Both the original and modified models captured the differences of SOC change among the different treatments. The RRMSE values (5%–8% among the treatments) calculated using the simulations from the modified DNDC were comparable with the RRMSE (7%–8%) based on the original DNDC simulations. Therefore, the incorporation of the microbial processes and activities on SOM decomposition in this study did not considerably improve the model's performance in capturing the observed SOC changes over decadal to century time scales. This may be because any missing microbial processes and activities of long‐term SOM decomposition could be compensated by calibration of the SOM decomposition rates.

We note some discrepancies between the modeled results and field measurements. The modified DNDC model over‐predicted CO2 release rates following the wheat residue incorporation and slightly under‐predicted CO2 release rates in August and September in 2016 at the OK site. The over‐predictions of the CO2 release rates may have resulted from over‐predictions of either litter availability or microbial growth and activity following the wheat residue incorporation, causing DNDC to predict relatively higher litter decomposition rates, DOC availability, and subsequently CO2 production and release. Over‐predicting initial litter decomposition could reduce the litter availability in the following August and September, and therefore be partially responsible for the under‐predicted CO2 release rates. The model also under‐estimated CO2 release rates in May and June in 2015. The under‐estimations of the CO2 fluxes could be due to under‐estimation of wheat autotrophic respiration considering that soil heterotrophic respiration would be low because there was no wheat residue incorporation or tillage during this period. Further studies, and additional years of field observations, are needed to better explain the discrepancies between the simulations and observations. In addition, it should be noted that there are uncertainties in the eddy covariance measurements of NEE due to uncertainties associated with the instrument, source heterogeneity, gap‐filling, and the turbulent nature of the transport process (Mauder et al., 2013; Richardson et al., 2006). These measurement uncertainties can also contribute to discrepancies between simulated and measured NEE. For example, random measurement errors were estimated from around 20%–80% of the NEE values at several croplands (Mauder et al., 2013; Richardson et al., 2006). These measurement uncertainties are larger than the discrepancies between the simulated and measured annual total NEE in this study.

Impacts on Soil Carbon Dynamics

The model predictions under different scenarios quantified the sensitivity of SOC change and CO2 flux from soil heterotrophic respiration to changing FMPs, microbial physiology, and SOM property. In general, the simulations demonstrated that SOC changes were largely regulated by the C input from manure amendments and the soil CO2 fluxes were predominantly controlled by the C inputs through manure amendments and wheat residue return (Figure 6). Simulated SOC changes were strongly correlated with the manure amendments, and simulated soil CO2 fluxes were strongly correlated with C input rates through wheat residue return across the simulated scenarios (Table 3). These results are consistent with previous studies that have highlighted the importance of FMPs in controlling SOC change and CO2 fluxes in agricultural ecosystems (e.g., McLauchlan, 2006).

Furthermore, in addition to the microbial physiology parameters, manure amendments and wheat residue return rates exerted large impacts on microbial activities (Figure 6; Blagodatskaya & Kuzyakov, 2013; Six et al., 2006). For example, the model predicted relatively low living microbial biomass under scenarios with low C inputs even though the microbial physiology parameters (e.g., high V MaxDOCuptake or CAE) were favorable for microbial growth. Although SOC changes were also regulated by parameters other than C inputs (e.g., CAE, F MICtoDOC, and VmaxHumus), particularly for the variations of ±50% of the baseline for these parameters, we note that the sensitivities depended on the variation ranges of the parameters (Figure 6). While the variations of FMPs were set to represent common changes in agricultural management, the ±50% variations for the microbial physiology might be overestimated for a specific cropland (i.e., the Rothamsted winter wheat field in this study). For example, the CAE was the most important microbial parameter regulating the SOC change and CO2 flux (Figure 6) and the ±50% variation of the baseline value was generally larger than the CAE modifications investigated in the previous studies (e.g., Allison et al., 2010; Wang et al., 2013). Therefore, the importance of CAE in regulating the SOC change and CO2 flux under the ±50% variations might be overestimated. Based on the large impacts of FMPs in governing SOC change, CO2 flux, and microbial activities as predicted in this study and reported in other studies (e.g., Blagodatskaya & Kuzyakov, 2013; Böhme et al., 2005), we conclude that SOC dynamics might be largely regulated by FMPs in agricultural ecosystems, although microbes played a central role in decomposing SOM.

The model also predicted differential impacts on SOC changes and CO2 fluxes of different FMPs through incorporating processes to characterize these FMPs. For example, different responses of SOC changes and CO2 fluxes were predicted between the applications of FYM and wheat residue. Consistent with the field measurements at the Rothamsted site, the model simulated that most C inputs from wheat straw were released as CO2 fluxes and a larger fraction of added manure C was retained in the soil compared to C inputs from wheat straw (Figure 5; Johnston et al., 2009). As a result, the impact on simulated SOC changes of manure amendment was larger than the impact of wheat straw return (Figure 6). The predicted higher efficiencies of organic carbon retention by applying manure were primarily attributed to the different qualities of organic matter. In comparison with the wheat straw, which had a C/N ratio of about 65, the organic manure had a lower C/N ratio of 12.5 (Table 2). Therefore, for a given amount of organic matter input, the model predicted more manure C converted into microbial biomass during decomposition due to lower N limitations for microbial growth, and thereby predicted more C retained in the soil derived from microbial death. In addition, the model partitions a fraction of manure into the humads pool (Li et al., 2012) that could be directly decomposed into the passive humus pool with low decomposition rates.

We note that microbial physiology may exert considerable impacts on SOC changes and soil CO2 fluxes. Simulated SOC increase rates tended to increase with increasing CAE and decreasing F MICtoDOC if the FMPs, SOM, and other microbial physiology parameters were not changed (Figure 5). Simulated SOC changes were positively correlated with CAE (P < 0.01) and negatively correlated with F MICtoDOC (P < 0.01) even under conditions with relatively large variations in C input (0.0–1.0 for fraction of returned crop residue, and 0–2,800 kg C ha−1 year−1 for organic manure input) (Table 3). In addition, the sensitivity of simulated SOC changes to changes in CAE or F MICtoDOC were comparable with that to changes in wheat residue return and manure application rate under the V50 scenarios (Figure 6b). Simulated annual CO2 fluxes tended to be lower in scenarios with higher F MICtoDOC and lower CAE if the other factors were not changed (Figure 5). Annual CO2 fluxes also significantly correlated with these two parameters under the V50 scenarios (Table 3), indicating appreciable impacts of microbial physiology on soil heterotrophic respiration in agricultural ecosystems. These results are consistent with previous studies that have suggested the importance of microbial physiology in controlling SOC change and CO2 fluxes (e.g., Liang et al., 2019; Six et al., 2006; Wang et al., 2021). Furthermore, the enzyme production and decay rates exerted comparable impacts with manure amendments in regulating the SOC changes under ±50% variations in microbial physiology parameters (Figure 6). Therefore, the function of enzymes in catalyzing SOC decomposition is considerable, although the enzyme concentrations are very small compared to other C pools. These results highlight the central role of soil microbes in converting the C inputs to resistant SOC, and suggest that microbial physiology and activities need to be considered when assessing SOC changes even in agricultural ecosystems where FMPs usually play a major role in controlling SOC dynamics.

Integrating Soil Microbial Dynamics Into Agricultural Ecosystem Models

Modeling microbial‐mediated SOM dynamics in agricultural ecosystems is in an early stage, and much additional work is required for a more complete framework of all of the processes involved. As shown in this and other studies (e.g., Bastian et al., 2009; Böhme et al., 2005; Bossio et al., 1998), microbial biomass and activity are strongly regulated by crop residue incorporation and other FMPs, which stresses the importance of correct simulations of crop dynamics and parameterization of FMPs when predicting responses of microbial activities to environmental changes and their influences on SOM decomposition and CO2 fluxes in agricultural ecosystems. However, most microbial models that investigate impacts of soil microbes on SOM turnover are based on static or prescribed C input (e.g., Abramoff et al., 2017; Schimel & Weintraub, 2003; Sulman et al., 2018). Therefore, biases and uncertainty in predicting SOM dynamics may result from neglecting crop dynamics and associated changes in crop residue return. Furthermore, the impacts of FMPs on microbial and enzyme activities are large but highly uncertain (Burns et al., 2013). Therefore, further studies need to be performed to quantify how FMPs affect microbial dynamics, how this in turn affects vegetation and SOM dynamics, and how the interactions among FMPs, soil microbes, crop growth, and SOM turnover are best incorporated into process‐based biogeochemical models.

Another uncertainty in the DNDC model framework is the representation of SOM pools. The multiple SOM pools in DNDC were designed to empirically represent decay rates of a wide variety of organic matter compounds, with the multiple SOM pool dynamics simulating overall actual SOM dynamics (Li et al., 1992a). The approach of multiple SOM pools is applicable for predicting SOC dynamics under varied FMPs (e.g., Figure 2; Li et al., 1994, 1997), and may be necessary for simulating SOM dynamics in agricultural ecosystems where organic amendments and SOM consists of materials at various stages of decomposition and that have different characteristics (De Graaff et al., 2010; Lal, 2018). However, the pools based on decay rates are conceptual and difficult to relate to measurements (Schmidt et al., 2011), although this approach may be more informative than the model with a single SOC pool (e.g., Allison et al., 2010). In addition, the modified DNDC model presented here simplifies the processes of microbial and enzyme dynamics. For example, the model does not differentiate microbial communities responsible for SOM decomposition (i.e., it simulates a single microbial pool) in order to minimize the number of poorly constrained microbial parameters and to maintain the applicability of DNDC to agricultural ecosystems with a wide range of management practices. Further steps to incorporate microbial communities and responses of microbial and enzyme activities to substrate chemistry and microbial stoichiometry (Sinsabaugh et al., 2008) into mathematical models will be required for a complete model framework.

High uncertainty in parameters related to soil microbes and SOM decomposition were also detected. Several new parameters were estimated by calibrating the simulations against the field observations of SOC changes under the NPK treatment due to a lack of specific parameter values (Table 1). Most of the calibrated parameters were comparable with published values. For example, the maximum rates of decomposition of various SOM pools and DOC uptake were within the reported variations (e.g., Wang et al., 2013). However, little information can be found for some parameters, such as the half‐saturation constants in the Michaelis‐Menten equations for calculating the decomposition of the SOM pools, because these parameters could be variable across different environmental conditions and are poorly constrained by observations (Moorhead & Sinsabaugh, 2006). Evaluating the model against more observations of microbial activities and C dynamics from different locations and under different environmental conditions should help to constrain uncertain parameters. Furthermore, it is hard to estimate the plausible variations of the parameters related to soil microbes and SOM decomposition for a specific ecosystem or site (e.g., the winter wheat system investigated in this study). Insufficient constraints on the parameters suggest that further studies are required to identify plausible ranges for parameters regarding microbial dynamics and SOM decomposition, and how these parameters depend on environmental conditions. Our efforts incorporating microbial regulations into a traditional biogeochemical model and investigating responses of SOC change and CO2 flux to varied FMPs, SOM, and microbial physiology provide critical guidance on identifying processes/parameters for which soil C stocks and CO2 fluxes are strongly sensitive.

Conclusions

We developed a microbial‐mediated SOM decomposition model, based on a biogeochemical model (DNDC), to simulate C dynamics in agricultural ecosystems. The model simulates decomposition of SOM and external organic matter, soil respiration, and SOC formation by explicitly simulating microbial and enzyme dynamics and their controls on SOM decomposition, in addition to simulating impacts on C dynamics of climate, soil properties, crop, and common FMPs. The model was evaluated against field observations of NEE and SOC change in two winter wheat systems to assess its performance in predicting both CO2 fluxes and SOC changes. The model successfully captured seasonal variations of NEE and SOC changes under different FMPs. The incorporation of microbial processes improved the model's performance in simulating peak CO2 fluxes induced by the wheat residue return primarily because the improved model captured priming effects of wheat residue inputs. We applied the model to investigate impacts of microbial physiology, SOM, and FMPs on soil microbial activity, SOC change, and CO2 fluxes. The results show that C input through crop residue or manure drove microbial activity and predominantly regulated the soil CO2 fluxes, and that manure amendment largely regulated SOC change. Microbial physiology also exerted considerable impacts on the microbial activities, SOC changes, and soil CO2 fluxes, emphasizing the necessity of considering microbial physiology and activities when assessing soil C dynamics in agricultural ecosystems where FMPs usually play a major role in controlling soil C.

Variable	Definition, unit
CAE	Maximum microbial carbon assimilation efficiency
Clay	Soil clay content, %
CNLitter , CNHumads , CNHumus , CNMicrobe , CNEnzyme	C:N ratios of litters, humads, humus, soil microbes, and enzyme respectively
D Litter , D Litter_vl , D Litter_r , D Humads_l , D Humads_r , D Humus	Decomposition rates of litters, very liable litter, liable litter, resistant litter, liable humads, resistant humads, and humus, respectively
D Enzyme	Enzyme decay rate, h−1
DeathMicrobe	Mortality of microbial biomass
DN	Dissolved nitrogen, including both DON and mineral N
DOC	Concentration of DOC, mg C g−1 soil
DOCLeaching	DOC transferred out of soils through water leaching
DOC M	DOC consumed through methanogenesis
DOCManure	DOC input from organic manure
DOC N	DOC consumed through denitrification
DOCRoot	DOC input from root exudation
DOCUptake	DOC uptake by microbes
DON	Concentration of DON, mg N g−1 soil
DONLeaching	DON transferred out of soils through water leaching
DONManure	DON input through organic manure
DONRoot	DON input through root exudation
DONUptake	DON uptake by microbes
Ea Litter_vl , Ea Litter_l , Ea Litter_r , Ea Humads_l , Ea Humads_r , Ea Humus	Activation energies in the Arrhenius equation for decomposition of very labile litter, labile litter, resistant litter, labile humads, resistant humads, and humus, respectively, KJ mol−1
F HumadsToDOC	Fractions of decomposed humads that allocated to DOC
F MICtoDOC	Fraction of dead microbial biomass that allocated to DOC
GrowthMicrobe	Growth of microbial biomass
Humads_l	Concentration of labile humads, mg C g−1 soil
Humads_r	Concentration of resistant humads, mg C g−1 soil
Humus	Concentration of humus, mg C g−1 soil
InputLitter	Litter input from crop residue
K Litter_vl , K Litter_l , K Litter_r , K Humads_l , K Humads_r , K Humus	Michaelis half‐saturation constant for decomposition of very labile litter, labile litter, resistant litter, labile humads, resistant humads, and humus, respectively, mg C g−1 soil
K DOC	Michaelis DOC half‐saturation constant for DOC uptake, mg C g−1 soil
K_O_2	Michaelis O2 half‐saturation constant for DOC uptake, mmol cm−3
Litter_vl	Concentration of very labile litter, mg C g−1 soil
Litter_l	Concentration of labile litter, mg C g−1 soil
Litter_r	Concentration of resistant litter, mg C g−1 soil
Microbe	Concentration of soil living microbes, mg C g−1 soil
MN	Soil mineral nitrogen
N Demand	Potential N required for microbe growth
O2	Concentration of soil oxygen, mmol cm−3.
P Enzyme	Enzyme production
R	The gas constant, 8.314 J K−1 mol−1
R MicrobeMaintenance	Rate of microbial maintenance respiration, h−1
R MicrobeDeath	Microbial turnover rate, h−1
R EnzymeProduction	Enzyme production rate, h−1
R EnzymeDecay	Enzyme decay rate, h−1
SW	Soil moisture in water‐ filled fraction of total porosity
T	Soil temperature, °C
VmaxLitter_vl , VmaxLitter_l , VmaxLitter_r , VmaxHumads_l , VmaxHumads_r , VmaxHumus	Maximum decomposition rate of very labile litter, labile litter, resistant litter, labile humads, resistant humads, and humus, respectively, at reference temperature, mg C (mg Enzyme C)−1 hr−1
VmaxDOCuptake	Maximum uptake rate of DOC at a reference temperature of 20°C, mg C (mg microbe C)−1 hr−1