A new bioenergy model that simulates the impacts of plant‐microbial interactions, soil carbon protection, and mechanistic tillage on soil carbon cycling

Advancing our predictive understanding of bioenergy systems is critical to design decision tools that can inform which feedstock to plant, where to plant it, and how to manage its production to provide both energy and ecosystem carbon (C) benefits. Here, we lay the foundation for that advancement by integrating recent developments in the science of belowground processes in shaping the C cycle into a new bioenergy model, FUN‐BioCROP (Fixation and Uptake of Nitrogen‐Bioenergy Carbon, Rhizosphere, Organisms, and Protection). We show that FUN‐BioCROP can approximate the historical trajectory of soil C dynamics as natural ecosystems were successively converted into intensive agriculture and bioenergy systems. This ability relies in part on a novel tillage representation that mechanistically models tillage as a process that increases microbial access to C. Importantly, the impacts of tillage and feedstock choice also influence FUN‐BioCROP simulations of warming responses with no‐till perennial feedstocks, miscanthus, and switchgrass, having more C that is unprotected and susceptible to warming than tilled annual feedstocks like corn–corn–soybean. However, this susceptibility to warming is balanced by a greater potential for increases in belowground C allocation to enhance soil C stocks in perennial systems. Collectively, our model results highlight the importance of belowground processes in evaluating the ecosystem C benefits of bioenergy production.


| INTRODUCTION
Bioenergy production represents a significant tool in society's portfolio to slow the pace of carbon dioxide (CO 2 ) emissions, particularly in sectors that are difficult to decarbonize (Reid et al., 2020). However, the extent to which bioenergy production can power certain economic sectors and at the same time maintain net zero or net negative carbon (C) emissions remains uncertain (Searchinger et al., 2009). One way to reduce this uncertainty is to choose the optimal combination of feedstock, management practices, and location (i.e., Who? How? & Where?) that maximizes yields aboveground and minimizes soil C losses or even leads to soil C gains. To optimize this choice, the bioenergy community has developed or adapted existing ecosystem models (e.g., Agro-IBIS, CENTURY, DayCent, DNDC) that expand the predictive power of limited empirical datasets that exist at the ecosystem level (Adler et al., 2007;Hudiburg et al., 2015;Li et al., 1994;Qin et al., 2016;VanLoocke et al., 2017). While these models have proven integral to our predictive understanding of the C economy of bioenergy production, the models have not kept pace with recent advances in the empirical understanding of belowground processes that govern the C cycle.
Empirical understanding of the importance of belowground processes in shaping the C cycle has been transformed over the past decade (Blankinship et al., 2018;Schmidt et al., 2011). This transformation is not limited to natural ecosystems, but instead is directly applicable to bioenergy feedstocks where the interaction between roots, symbiotic microbes, and saprotrophic microbes directly impacts the availability of nutrients to support yields and the ability of different feedstocks to increase soil C stocks (Briones et al., 2019;Zhu et al., 2018). Despite this growing recognition of the importance of belowground processes, most models employed in the bioenergy sphere do not explicitly represent microbial function or plant-microbial interactions. Using first-order decomposition functions and treating microbial biomass as a soil organic matter (SOM) pool that does not have a functional impact on the cycling of other organic matter (OM) pools, the majority of models do not represent the following critical processes: (1) the ability of plants to send C to symbiotic and free-living microbes in the rhizosphere to prime decomposition and enhance nutrient availability (Brzostek et al., 2013;Cheng et al., 2014;Meier et al., 2017), (2) the functional capabilities of soil microbes in producing enzymes to decompose SOM (Schimel & Weintraub, 2003), and (3) the role of key microbial physiological parameters like C use efficiency (CUE, the efficiency of the conversion of OM into microbial biomass) and microbial biomass turnover rate in controlling the production of microbial necromass products that are preferentially protected and stabilized in soils (Buckeridge et al., 2020;Cotrufo et al., 2013;Liang et al., 2019). Given that these processes have been shown to govern the response of a variety of ecosystems to elevated CO 2 (Terrer et al., 2018), shifts in nitrogen (N) availability (Carrara et al., 2018), and warmer temperatures (Allison et al., 2010;Qin et al., 2019), there is a clear need to integrate these processes into the foundational models that have guided our predictive understanding of bioenergy ecosystems.
Predictions of the consequences of different management practices including bioenergy feedstock choice, residue decisions, and tillage are also likely to be improved by integrating the critical processes highlighted above. For example, foundational models typically represent tillage impacts on decomposition by ramping up the decay rate (i.e., k) in first-order decomposition equations (Maharjan et al., 2018). While this framework can accurately capture soil C losses, it does not capture the mechanism, whereby tillage decreases the ratio of physically protected to unprotected C and increases the ability of soil microbes to access and decompose the formerly protected C (Jastrow & Miller, 1997;Reicosky et al., 1995). In addition, the distribution of unprotected and protected soil C can feed back on the response of bioenergy systems to warming because of their differential temperature sensitivity which varies across ecosystems (Lugato et al., 2021). Finally, there are distinct differences among bioenergy feedstocks in a number of traits that impact microbial physiology and function. Differences in root and aboveground litter chemistry, N acquisition strategies, and rhizosphere priming among feedstocks could lead to shifts in microbial physiology, the production of microbial necromass products, and decomposition rates. Collectively, advancing the representation of belowground processes in bioenergy models has the potential to predict the optimal combination of feedstock, management, and location, and lead to decisions that enhance bioenergy's ability to reduce CO 2 emissions.
Here, we present a new model framework that advances the representation of belowground processes in bioenergy models, FUN-BioCROP (Fixation and Uptake of Nitrogen-Bioenergy Carbon, Rhizosphere, Organisms, and Protection). FUN-BioCROP is a novel bioenergy model that accounts for rhizosphere priming, agricultural management, and plant traits in a framework that couples plant C investment in N acquisition with soil microbial decomposition and SOM stabilization. We developed the model to meet three objectives: (1) Formulate and test new representations of plantmicrobe interactions, microbial physiology, and mechanisms of stable soil organic carbon (SOC) creation and loss during tillage that can be incorporated into more established model frameworks (e.g., DayCent, Agro-IBIS; see Lawrence et al., 2019;Shi et al., 2016;Sulman et al., 2019); (2) assess the impact of these representations on model estimates of soil C across different bioenergy feedstocks; (3) test key hypotheses on how the incorporation of these representations impacts the response of bioenergy systems to changes in plant traits and environmental conditions. As such, our modeling efforts here represent a bridge that links the use of models to test and generate new hypotheses as well as assess the impacts of new theories or paradigms with the more traditional use of models to generate projections that can directly inform management and policy decisions (e.g., DayCent, Agro-IBIS, DNDC).

| MATERIALS AND METHODS
To develop FUN-BioCROP, we modified the coupled FUN-CORPSE model (Fixation and Uptake of Nitrogen-Carbon, Organisms, Rhizosphere, and Protection in the Soil Environment, Sulman et al., 2017) for bioenergy systems by incorporating a new tillage mechanism, cropspecific parameters, OM addition, fertilization, and harvest. Below we describe how we developed FUN-BioCROP and how we formulated model simulations to answer the following questions: (1) How do commonly used model representations of tillage vary in soil C projections as compared to the FUN-BioCROP tillage mechanism? (2) How do differences between annual and perennial bioenergy feedstocks in their integrated plant-microbial traits and tillage regimes impact model estimates of soil C? (3) To what extent do these model differences in the representation of feedstock traits alter the response of soil C estimates to warming temperatures? (4) Does the outcome of altering these plant traits (i.e., increased rhizosphere C exudation) on model estimates of soil C vary across different feedstocks?

| FUN-CORPSE model: Nitrogen acquisition and SOM dynamics
We developed the FUN-BioCROP model by modifying FUN-CORPSE (Sulman et al., 2017), which couples a plant N acquisition model ("FUN";Brzostek et al., 2014;Fisher et al., 2010;Shi et al., 2016) to a microbial SOM model ("CORPSE"; Sulman et al., 2014a). As an ecosystem-scale tool, FUN-CORPSE has previously demonstrated the importance of including rhizosphere processes, SOM protection, and plant-microbial interactions in models to accurately project soil C dynamics in forests (Sulman et al., 2017). A detailed description of the core model equations is available in the supplemental material of Sulman et al. (2017).
Briefly (see Brzostek et al., 2014 for full model description), FUN dynamically calculates the cost of plant N acquisition to optimize C allocation and maximize net primary productivity (NPP). After determining plant N demand to meet growth requirements, FUN's resistance framework allows simultaneous N uptake across 5 possible pathways. Each pathway's resistance (i.e., C cost) determines its uptake rate, thereby achieving the lowest possible C cost to the plant. The pathways include: (i) biological N fixation; (ii) retranslocation; (iii) nonmycorrhizal root uptake; (iv) mycorrhizal root uptake; (v) N from storage. Biological N fixation is carried out by both symbiotic bacteria and free-living N fixers, and its cost is modeled as a function of temperature. FUN uses the same N fixation representation for all the plants. It becomes a significant N source when soil N levels are very low. Retranslocation occurs prior to autumn leaf senescence, and its cost increases as leaf N decreases. The cost of root facilitated uptake (i.e., mycorrhizal, or non-mycorrhizal) is modeled as a function of soil N levels and root biomass, with their decrease leading to increased uptake cost.
We modeled soil C and N cycling with a modified version of CORPSE (described in detail in Sulman et al., 2014aSulman et al., , 2017, which is a microbially-explicit SOM decomposition model that includes physical protection. CORPSE has litter, rhizosphere, and bulk soil compartments that each contain three SOM chemical types representing fast (labile), slow (resistant), and microbial necromass C and N. The litter compartment has only unprotected SOM, while the rhizosphere and bulk soil compartments additionally have physically protected fast, slow, and microbial necromass SOM that is inaccessible to decomposition. Protection rates vary by SOM chemical type. Microbial biomass, temperature, and soil moisture determine decomposition rates, and microbial biomass grows through SOM decomposition. Leaf litter enters the litter layer, and root litter is divided between the rhizosphere and bulk soil compartments according to rhizosphere size. Immobilization or N mineralization occurs depending on the balance between the SOM C:N ratio and the microbial biomass C:N ratio, with overflow respiration occurring when N is insufficient to support microbial growth. Inorganic N is modeled as one pool without distinguishing between N species (e.g., ammonium and nitrate). Across all of these processes, a sensitivity analysis for CORPSE tested all major model parameters ±20% of the value and found the highest sensitivity of C turnover rate to microbial turnover rate and the Michaelis-Menten constant (Sulman et al., 2014b).
In the coupled FUN-CORPSE model (Sulman et al., 2017), FUN models plant root inputs of C into the soil through root exudation (i.e., C invested in non-mycorrhizal root uptake) and C transfer to mycorrhizae, and CORPSE simulates mineralized N available to plants and the fate of C spent belowground. Carbon spent on non-mycorrhizal root uptake in FUN is considered root exudation and enters the CORPSE rhizosphere fast C pool. To reflect the range of mycorrhizal hyphae throughout the soil profile, C spent on mycorrhizal uptake in FUN enters both rhizosphere and bulk soil fast C and N pools according to the exudate C:N and rhizosphere size. The percent total soil volume composed of the rhizosphere was calculated at each time step with live root biomass simulated by DayCent (see details below; Berardi et al., 2020;Moore et al., 2020;Parton et al., 1998) and plant-specific root morphology data (diameter and specific root length), assuming that root volume was cylindrical and the rhizosphere extended 1 mm from the root surface. Most model parameters were the same as previously published, which were calibrated using mean SOM pools and N mineralization rates measured across 45 deciduous forested plots in southern Indiana (Sulman et al., 2017). Because FUN-CORPSE relies on the chemistry of inputs and microbial physiology and stoichiometry, it is broadly applicable to various ecosystems as evidenced by the use of its framework at the global scale to explore the effects of increasing CO 2 on soil C (Sulman et al., 2019). Due to this broad applicability, we were able to adapt the model to bioenergy systems by including new parameterizations of plant traits including whole plant stoichiometry, leaf and root litter chemistry, and root morphology.
Tillage in FUN-BioCROP simulates the mechanism of actual C loss in tilled soils. Tilling soils disrupts soil aggregates by both increasing their turnover and reducing their formation (Six et al., 1999). This disruption releases physically-protected particulate organic matter (POM), which becomes more susceptible to decomposition (Jastrow & Miller, 1997) as opposed to that contained within aggregates (Besnard et al., 1996;Six et al., 1999). Because it includes SOC protection, FUN-BioCROP can realistically model SOC loss following tillage by transferring protected SOC into unprotected pools, making it available for microbial decomposition. To achieve this, we incorporated new tilled rhizosphere and tilled bulk soil compartments, resulting in five soil compartments: litter, tilled rhizosphere, tilled bulk, rhizosphere, and bulk ( Figure S1). The size of each compartment depends on two factors: (i) a tillage disturbance parameter (T) that assigns a fraction of the soil column to be tilled and the remainder to be untilled and (ii) the rhizosphere size which is determined by root morphology and biomass (see FUN-CORPSE description above). We assigned T a value of 0.15 to represent the proportion of the soil disturbed by tillage events. When tillage occurs in FUN-BioCROP, a fraction of the litter (L) and tilled rhizosphere (R till ) compartments mixes into the tilled bulk soil compartment (B till ), with the fraction (f R and f L for rhizosphere and litter, respectively) varying from 0.2 to 0.92 as determined by the type of tillage in the DayCent management schedules. In the following equations, the 0 and 1 subscripts indicate preand post-tillage, respectively.
Next, tillage events cause the conversion of a set percentage (P, set to 30% to approximate C and N loss following tillage) of protected C and N (SOM P ) into unprotected C and N (SOM U ) within the tilled rhizosphere and tilled bulk soil (Figure 1a), allowing for microbial decomposition to occur on previously inaccessible SOM.
Finally, all inputs to the rhizosphere and bulk, tilled and not tilled compartments are divided according to rhizosphere size and the tillage disturbance parameter.
Organic matter amendments enter the litter compartment C and N pools (Figure 1b), with the OM lignin fraction entering the slow C and N pools and the remainder entering the fast pools. Chemical N fertilizer applications are added to the shared inorganic N pool ( Figure 1c). Above-and belowground litter inputs are calculated as the daily change in C and N pools. Aboveground litter inputs are added to the fast and slow C and N litter pools according to the fast-decomposing fraction of the tissue (see below). Belowground litter inputs are fractionated the same way and are additionally divided according to the size of the soil compartments (e.g., soil enters tilled versus untilled compartments based on T (the tillage disturbance parameter), and rhizosphere vs. bulk depending on the rhizosphere size, Figure S1). Harvests reduce C and N inputs to the litter layer according to the percentage of aboveground biomass removed ( Figure 1d).
Plant-specific attributes in FUN-BioCROP include leaf and root litter chemistry, and root morphology (Table S1). The fast-decomposing fraction of the leaves and roots was determined by the average lignin:N ratio of the tissues used in DayCent (Hudiburg et al., 2015). The average lignin:N of young and old leaves determined the fast-decomposing fraction of leaf litter, and fine root lignin:N determined the root fast fraction. The root diameter and specific root length of each plant (Table S1) were either obtained from published values or measured in the field experiment at the University of Illinois at Urbana-Champaign (UIUC) Energy Farm (site description below) and used to calculate rhizosphere size. Soil samples for root analysis at the UIUC Energy Farm were taken to 30-cm depth, with 32 cores per crop (4 plots per crop × 8 cores per plot) taken during summers 2018 and 2019. Soils were frozen until analysis then serially sieved to 4 mm with any remaining roots manually removed from sieved soil. Roots were soaked for 15 min in sodium hexametaphosphate solution and washed over a 500-micron sieve. Roots were scanned at 600 dpi on an Epson V850 photo scanner (Epson America, Inc.) then dried at 65°C to constant mass. Root length and diameter (Adam von Haden, unpublished data) were measured using IJ_Rhizo software (Pierret et al., 2013) in ImageJ (Schneider et al., 2012). Specific root length was calculated as meters per gram root.
FUN-BioCROP input data are generated by DayCent, which has been calibrated for bioenergy crops at the UIUC Energy Farm and accurately simulates NPP and C and N in biomass pools (Hudiburg et al., 2015). Specifically, FUN-BioCROP uses above-and belowground NPP, and C and N in aboveground live vegetation, juvenile, and F I G U R E 1 Diagram of the FUN-BioCROP model (fixation and uptake of nitrogen-bioenergy carbon, rhizosphere, organisms, and protection). FUN-BioCROP is a version of FUN-CORPSE (fixation and uptake of nitrogen-carbon, organisms, rhizosphere, and protection in the soil environment, Sulman et al., 2017) that has been modified for use in bioenergy systems to simulate plant carbon (C) allocation to nitrogen (N) acquisition and microbial soil organic matter decomposition. Yellow and purple represent C and N, respectively. Letters indicate modifications to FUN-CORPSE to develop FUN-BioCROP. Specific soil compartments (i.e., litter, tilled and untilled rhizosphere soil, and tilled and untilled bulk soil) are not shown (see Figure S1). Storage N in the FUN model is not shown. Plant parameters include tissue chemistry and root morphology. NPP, net primary productivity; OM, organic matter mature live roots generated by DayCent (Figure 1e). We used the most recent version of the DayCent-CABBI biogeochemical model, developed for the process-based representation of perennial bioenergy crops Moore et al., 2020). DayCent-CABBI includes separate pools for perennial plant leaves and stems, allowing separate parameterization and temporal variation of leaf and stem chemical composition (lignin and C:N ratios) and biomass.
Historic and future DayCent simulations used daily climate data from 1980 to 2019 downloaded from the Daymet database (Thornton et al., 2016). Historic simulations were calibrated from 1847 to 2007 using agricultural history and site-specific management for the UIUC Energy Farm (Hudiburg et al., 2015). Modern simulations (2008-2019) of switchgrass, miscanthus, and corn-corn-soybean were conducted using site-specific management (see Moore et al., 2020 for corn-corn-soybean and switchgrass; Table  S2 for miscanthus). We parameterized crops for DayCent using the observed maximum leaf area index (Joo et al., 2016) and minimum and maximum C:N in leaves, stems, and roots (Smith et al., 2013) at the UIUC Energy Farm (Table S3). Switchgrass and corn-corn-soybean simulations were calibrated previously (Moore et al., 2020) using observed annual above-and belowground NPP and soil C, while we calibrated miscanthus specifically for this study using the same methods ( Figure S2). Previous site management was continued 80 years into the future for all crops. Because of age-related yield decline, miscanthus was replanted every 13 years. Switchgrass did not exhibit age-related yield decline at this site to the extent that would necessitate replanting.

| Spin up and model validation
We spun up FUN-BioCROP for 3,000 years with a tallgrass prairie plant community based on the Konza Prairie Biological Station (Bark, 1987) composed of big bluestem (Andropogon geradii), Indiangrass (Sorghastrum nutans), and little bluestem (Andropogon scoparius). All plant parameters were based on the average value of their tissue chemistry and root morphology (Table S1). Following spin up, we conducted a historical simulation from 1847 to 2007 using the agricultural schedule of events, plant biomass, and productivity data from DayCent, which has comprehensive and detailed agricultural land use history of cropping in the great plains. This includes historical changes in crop rotations and varieties, tillage methods, nutrient additions, and crop yields (Hartman et al., 2011).
We also validated the model with bioenergy crop field data from the UIUC Energy Farm located in Urbana, IL, USA (40° 3ʹ 46.209ʺ N, 88° 11ʹ 46.0212ʺ W, 220 m above sea level). In spring 2008, maize (Zea mays), miscanthus (Miscanthus x giganteus), and switchgrass (Panicum virgatum) were planted in 4 ha (200 × 200 m) plots that had been prepared with diammonium phosphates, potash, and lime amendments via variable-rate technology to achieve uniform soil nutrient distribution (see full details in Zeri et al., 2011). Corn (an annual crop) planting and management occurred according to standard regional agricultural practices, and the perennials (miscanthus and switchgrass) were managed according to known best practices (Zeri et al., 2011). Corn plots were planted with soybean (Glycine max) every third year in a typical corncorn-soybean rotation. Soil temperature and moisture were measured at 10 cm depth (Hydra Probe II, Stevens Water Monitoring Systems, Inc.) in the center of each plot (Moore et al., 2020;Zeri et al., 2011). To hold soil temperature and moisture constant, we averaged the values from 2008 to 2016 to have 1 year of mean data for all simulations, including the spin up and 1847 to 2007 historical period. This allowed us to explore the effect of bioenergy crop differences, rather than the effects of micrometeorological differences, on SOC.
We also validated FUN-BioCROP using SOC data and protection ratios (the ratio of the total soil protected C to unprotected C) measured at the UIUC Energy Farm. Soil carbon was measured to 30-cm depth in each plot (corncorn-soybean, miscanthus, and switchgrass) once per year from 2008 to 2016 (Ilsa Kantola, unpublished data). Protected and unprotected SOC was also measured in the top 20 cm of soil in the corn and miscanthus plots at the UIUC Energy Farm in 2019 (Joanna Ridgeway, unpublished data), and used for validation of protected to unprotected SOC ratios. Protected C was defined as the silt/ clay fraction that was denser than 1.85 g/ml and smaller than 53 μm. Unprotected C included the light fraction that floats with 1.85 g/ml density and the POM fraction that is larger than 53 μm (Lavallee et al., 2020).

| Model experiments: Tillage, warming, increased rhizosphere C
We compared the effects of the novel tillage mechanism in FUN-BioCROP to typical model representations of tillage. Many models simulate the effects of tillage using scalars to increase decomposition rates following a tillage event. These models parameterize the impact of tillage on decomposition rates to achieve accurate estimates of soil C. In contrast, after transferring a fraction of the tilled rhizosphere and litter compartments into the tilled bulk compartment, FUN-BioCROP transfers protected SOC (fast, slow, and necromass) in the tilled compartments into unprotected pools, making it accessible to microbial decomposition without scaling the decomposition rate. To explore the effects of different tillage mechanisms on SOC projections, we ran FUN-BioCROP using two different tillage schemes: (i) "physical protection loss" refers to our novel tillage mechanism that transfers protected SOM into unprotected pools in the tilled rhizosphere and bulk soil following tillage, thereby increasing microbiallyaccessible SOM without altering rate parameter values and (ii) "decomposition scalars" in which we turned off the conversion of protected C to unprotected C in FUN-BioCROP and ran it instead with scalars that altered the decomposition rate parameters of unprotected SOC pools for 30 days following a till event. Depending on the type of tillage, the scalars accelerated the decomposition in the fast and slow OM to different degrees across soil compartments (e.g., litter, bulk, or rhizosphere) and varied in depth and intensity. We ran FUN-BioCROP with both tillage mechanisms (physical protection loss and decomposition scalars) under historical conditions (1847-2007), and from 2008 to 2127 with a corn-corn-soybean rotation to investigate the long-term outcomes of different tillage simulations. Because we ran both tillage schemes using FUN-BioCROP, we were able to analyze the size of the protected and unprotected C pools resulting from each representation of tillage, which would not be possible running a model that does not include SOC protection.
As soils warm under climate change, SOC storage will depend on the balance between litter inputs and decomposition, with evidence indicating that warming can decrease soil C pools by altering microbial physiology and community structure (Melillo et al., 2017). Given that unprotected SOC can be more vulnerable to loss in warmer soils (Lugato et al., 2021), and plant litter quality determines protected versus unprotected SOC resulting from decomposition (Cotrufo et al., 2013), the variation in the C:N of bioenergy feedstocks could result in different levels of protected SOC by crop, with differing vulnerabilities to soil C loss in a warmer future climate. We, therefore, simulated a soil warming experiment that elevated soil temperatures by 5°C (as in Melillo et al., 2002) in each of the UIUC Energy Farm feedstocks to explore how soil C under different feedstocks changes in warmer conditions. We also modeled a soil warming experiment (elevating soil temperatures by 5°C) with both tillage schemes (decomposition scalars and physical protection loss) to examine the extent to which differences between the tillage regimes in estimates of the ratio of protected to unprotected SOC impact the response of the whole SOC pool to warming.
Similarly, as the climate changes, C exudation into the rhizosphere could increase as plants experience increased growth and shunt more C belowground (Phillips et al., 2006). With important consequences for ecosystem C balance, root C exudation can prime microbial decomposition, reducing soil C and altering nutrient cycling (Bengtson et al., 2012;Brzostek et al., 2013). Alternatively, rhizosphere C exudation can promote the formation of persistent mineral-associated organic matter (MAOM;  and soil macroaggregates (Baumert et al., 2018) which can protect C from microbial decomposition (Adu & Oades, 1978). To simulate the impact of enhanced rhizosphere C exudation in bioenergy systems, we increased C spent on both non-mycorrhizal and mycorrhizal root N uptake by 5% of daily average NPP in the UIUC Energy Farm feedstocks (10% daily increase total; corn-corn-soybean: 0.20 g C m −2 day −1 , fertilized miscanthus: 0.29 g C m −2 day −1 , unfertilized miscanthus: 0.22 g C m −2 day −1 , switchgrass: 0.21 g C m −2 day −1 ). The C spent via these root-facilitated pathways entered the rhizosphere fast C pool (nonmycorrhizal C), and the rhizosphere and bulk soil fast C pools (mycorrhizal C).

| Spin up
FUN-BioCROP spin up achieved stable SOC pools that equilibrated at 10.02-10.24 kg C m −2 for the tallgrass prairie prior to cultivation in 1847. The ratio of protected to unprotected SOC stabilized at 12.10 (Table 1), varying annually from 10.73 to 13.98. When stabilized, the C spent by plants on all N acquisition pathways averaged 6.1% of NPP, with root facilitated uptake (mycorrhizal and nonmycorrhizal) accounting for the majority (57%) of N acquired, followed by resorption (18%) and N fixation (8% ,  Table S4).

| Historical validation
The model estimates of SOC declined over the 1847-2007 period (Figure 2). The onset of cultivation resulted in immediate and steep SOC declines as tillage exposed new SOM to microbial decomposition. Over the 160-year period, the model projected a 36% decline in total estimated SOC stocks, and in 2007, it simulated SOC within the range of values measured at the UIUC Energy Farm prior to the initiation of the field experiment ( Figure 2). The ratio of protected to unprotected SOC at the end of the simulation was 14.63 (Table 1), ranging from 5.27 to 23.40 throughout the 160-years under different crop rotations and management regimes. During this period, the average yearly C spent on N acquisition by plants accounted for 9.4% of NPP, and 84.6% of N was acquired via root-facilitated uptake, followed by N fixation (11%), and finally very low resorption (1%) due to the dominance of annual crops (Table S4).
Examination of the relative effects of the new model modifications using the Energy Farm corn-corn-soybean data revealed that leaf chemistry, as depicted by the parameter for the fast-decomposing fraction of the tissue (hereafter, fast fraction), had the largest impact on FUN-BioCROP estimates of total soil C stocks (5.8 g C/m 2 / percent change in the leaf fast fraction by the end of the simulation). Next in the magnitude of the effect was the root fast fraction (2.8 g C/m 2 /percent change in the root fast fraction by the end of the simulation). These plantspecific parameters partition litter inputs into fast or slow decomposing material, thereby altering the rate at which microbes decompose and incorporate them into the soil. The size of the tillage disturbance parameter, which determines the amount of soil C affected by tillage, had a smaller impact on total soil C stocks by the end of the simulation (1.3 g C/m 2 /percent change in the till disturbance parameter by the end of the simulation). Finally, we found that the model was largely insensitive to different values of P, the percent of protected SOM transferred to unprotected pools following tillage (0.2 g C/m 2 /percent change in P by the end of the simulation) because the repeated disturbance resulting from multiple tillage events converted most of the protected SOM to unprotected SOM and precluded significant accumulation of protected SOM in tilled soils.

| Field data validation
FUN-BioCROP also approximated the magnitude of SOC stocks at the UIUC Energy Farm, projecting SOC values within the range of empirical observations during the first two decades of the 21st century (Figure 3). Percent difference for all simulated crop means fell within 5% of measured values (Figure 3a, corn-corn-soybean: 3.5%; unfertilized miscanthus: 0.3%, switchgrass: 5.1%). In model simulations extending to the year 2100, cultivation of perennial bioenergy feedstocks increased SOC pools over time, while the annual corn-corn-soybean rotation decreased SOC pools (Figure 3b). At the end of T A B L E 1 Model estimates of the total protected soil carbon (C) pool, unprotected soil C pool, and the ratio of protected to unprotected soil C as projected by FUN-BioCROP. Values represent the sum of all the C pools (fast, slow, and necromass) in the soil compartments (litter, tilled, and untilled rhizosphere and bulk) at the end of each simulation. Protection loss-physical protection loss tillage routine, scalars-decomposition scalars tillage routine the simulation, SOC had increased by 31% in fertilized miscanthus, 17% in unfertilized miscanthus, and 11% in switchgrass while the corn-corn-soybean rotation led to a 20% decrease in SOC. The 13-year planting rotation of miscanthus created a saw-tooth pattern in its modeled soil C trajectory caused by a large influx of belowground litter associated with replanting ( Figure 3b). The model also captured relative differences in the protection ratios (the ratio of the total soil protected C to unprotected C) between corn-corn-soybean and unfertilized miscanthus. As determined in the UIUC Energy Farm field study plots, the ratio of protected to unprotected SOC of corn-corn-soybean and unfertilized miscanthus were 5.46 and 3.61, respectively. The model simulated the corn-corn-soybean protection ratio as 16.59, and unfertilized miscanthus as 11.88 (Table 1). Although the magnitude of the measurements is different, in each case the corn-corn-soybean ratio is roughly 1.5 times that of unfertilized miscanthus. Over the course of the simulation, the protection ratios changed differently according to annual versus perennial plants. The modeled perennial plant SOC protection ratio decreased by 35.1% in fertilized miscanthus, 18.9% in unfertilized miscanthus, and 48.4% in switchgrass, while the corn-corn-soybean protection ratio increased by 13.3% ( Figure S3).
The bioenergy feedstocks varied in the proportion of NPP dedicated to N acquisition, as well as the dominant N acquisition pathway. Overall, corn-corn-soybean dedicated the highest amount of NPP to N acquisition (3.8%), followed by unfertilized miscanthus (3.2%), fertilized miscanthus (2.4%), and finally switchgrass (1.8%). The corncorn-soybean rotation also obtained substantially more N via root facilitated uptake (90.3%) than the perennials (fertilized miscanthus: 22.9%, unfertilized miscanthus: 22.0%, switchgrass: 27.3%; Table S4). In general, the perennials acquired a large proportion of N through retranslocation (~35%). Both fertilized and unfertilized miscanthus also acquired more N through fixation (7.5% and 7.0%, respectively) than did switchgrass and corn-corn-soybean (3% and 6%, respectively, Table S4). The comparatively low fixation in corn-corn-soybean reflects the dominance of corn, as opposed to soybean which forms a symbiotic relationship with N-fixing bacteria. Finally, corn-cornsoybean acquired more N through non-mycorrhizal root uptake than the perennials (22% vs. ~3.5%).

| Tillage experiment
The different tillage mechanisms (decomposition scalars vs. physical protection loss) resulted in divergent estimates of SOC and the ratio of protected to unprotected C, but similar percentages of NPP spent on N acquisition and proportions of N acquired through each pathway. Decomposition scalars and physical protection loss, respectively, resulted in 33% and 48% loss of SOC from 1847 to 2127, encompassing both the historical period (1847-2007) and 120 years of corn-corn-soybean rotation (Figure 4a). The SOC values projected by physical protection loss in 2008 approximated empirical SOC measurements that year taken at the UIUC Energy Farm, as opposed to the overestimation that resulted from using decomposition scalars (Figure 4a). At the end of the simulation, decomposition scalars resulted in a 98% increase in the ratio of protected to unprotected C as compared to physical protection loss (Figure 4b). The two tillage mechanisms resulted in a similar proportion of NPP being spent on N acquisition (decomposition scalars: 7.1%, physical protection loss: 7.4%), as well as very close estimations of N acquired through each pathway (Table S4). In general, root uptake accounted for the most N uptake (~88%), followed by N fixation (~11%), and finally resorption (~0.5%).
Warming soils by 5°C under the physical protection loss tillage scheme resulted in 3.3% additional loss of SOC as compared to control conditions over the 106 years experiencing warming (Figure 4a), and a 23% increase in the SOC protection ratio (Table 1; Figure 4b). Conversely, warming soils with tillage via the decomposition scalars mechanism resulted in a 2.1% gain in total SOC stocks by the end of the simulation as compared to the control conditions (Figure 4a), and a 34% increase in the protection ratio of the SOC (Table 1; Figure 4b).

| Warming experiment
Warming soils by 5°C resulted in decreased SOC pools, with relatively higher protection ratios, for all crops as compared to control conditions at the UIUC Energy Farm ( Figure 5; Figure S3). Perennial crop SOC decreased the most (3.7%-4.9%), and corn-corn-soybean experienced the smallest reduction (2.4%) due to soil warming. The ratio of protected to unprotected SOC also increased for all crops, with the percent increase relatively consistent across crops. Here, the annual average protection ratio of switchgrass increased the most compared to control conditions (27%), followed by miscanthus (unfertilized: 26%, unfertilized: 24%), and finally corn-corn-soybean (23%) experiencing the smallest increase (Table 1; Figure S3). Compared to control conditions, soil warming had a negligible effect on both the percent of NPP dedicated to N acquisition as well as the amount of N acquired via each possible pathway (Table S4).

| Rhizosphere C experiment
Increasing the mycorrhizal and non-mycorrhizal C flux into the rhizosphere increased estimates of SOC pools for all crops as compared to control ( Figure 5). Corn-cornsoybean experienced the largest SOC increase compared to control conditions (18%), followed by switchgrass (14.7%), and miscanthus (unfertilized: 13.2%, fertilized: 13.1%). The ratio of protected to unprotected C also increased for all crops with increased rhizosphere C flux, with perennials experiencing the largest increase (switchgrass: 16.6%, fertilized miscanthus: 10.1%, unfertilized miscanthus: 8.2%). Although the protection ratio of corn-cornsoybean soils under warming increased 3.1% by 2100 compared to controls, for much of the simulation increased rhizosphere C flux resulted in less protected C relative to unprotected C than in the control simulation ( Figure S3). Increased rhizosphere C inputs resulted in a greater percentage of NPP being spent on N acquisition by all crops as compared to control conditions. Fertilized miscanthus experienced the largest percent increase in NPP spent (16.8%), followed by corn-corn-soybean (15.8%), unfertilized miscanthus (11.8%), and finally switchgrass (10.0%). Although the magnitude of total N acquired via each pathway did not change greatly, in general, increased rhizosphere C inputs slightly decreased both mycorrhizal and non-mycorrhizal N uptake, and slightly increased N uptake through both resorption and N fixation (Table S4).

| DISCUSSION
FUN-BioCROP simulations of soil C stocks generally approximated multiple empirical data sources on soil C stocks at the UIUC Energy Farm from 2008 to 2016. The model estimated that 36% of soil C was lost upon prairie conversion to intensive modern agriculture ( Figure  2). This soil C loss is within empirical estimates which range from ~0% to 65% SOC loss following conversion of grasslands to cropping systems Guo & Gifford, 2002;Sanderman et al., 2017), with a median of 31% SOC lost from the upper 30 cm of the soil profile in a meta-analysis of 29 such studies (Sanderman et al., 2017). FUN-BioCROP also modeled soil C within 5% empirical error across 3 feedstocks and 8 years of measurements (Figure 3). While the model estimated trajectories in soil C that are not present in the empirical data, this result is not unexpected given difficulties in detecting and measuring SOC trends at sub-decadal scales (Conant et al., 2011). Overall, we acknowledge that the ability of the model to simulate soil C would be improved by validation across a larger dataset. However, given data limitations and the targeted scope of our modeling exercise in developing new representations that can be incorporated into larger ecosystem models, this general agreement between data and model estimates allows us to meet our objectives.
Comparing the modeled ratio of protected-tounprotected soil C against independent empirical measurements revealed that the model estimated the relative differences between the feedstocks, but the modeled value was ~threefold higher than that of the observed data (model/empirical protection ratio for corn-cornsoybean: 16.59/5.46, unfertilized miscanthus: 11.88/3.61). Some of this difference between simulations and observations may result from a mismatch in soil depths. The observed values were measured from the top 20 cm where this ratio would be expected to be somewhat lower than the 30-cm soil column simulated by FUN-BioCROP. More empirical data are needed to constrain this estimate, but this disparity in depths could partially explain the difference between the modeled and empirical ratios for these two crops. However, the magnitude of the mismatch in F I G U R E 4 (a) Estimated mean annual soil carbon (C) and (b) mean annual ratio of protected to unprotected soil C (i.e., protection ratio) as simulated by FUN-BioCROP under control (solid lines, 1847-2127) and 5°C soil warming conditions (dashed lines, 2021-2127) using two different representations of tillage: (i) the FUN-BioCROP tillage mechanism (physical protection loss, black lines) utilizes tilled rhizosphere and bulk soil compartments within which a till event causes a set percentage of the protected soil C to convert into unprotected soil C, making it accessible to microbial decomposition; and (ii) a tillage mechanism that scaled up the decomposition rate parameters of unprotected soil organic C pools for 30 days following a till event (decomposition scalars, red lines). Points in panel (a) show 2008 soil C field data measured at the University of Illinois at Urbana-Champaign Energy Farm. The dotted vertical line marking the year 2008 indicates the transition from the historical period (see Figure 2) to a 120-year corn-corn-soybean crop rotation and associated management schedule. FUN-BioCROP-fixation and uptake of nitrogen-bioenergy carbon, rhizosphere, organisms, and protection observed versus simulated values also reflects uncertainty in the parameter that controls the protection rate of new SOM inputs in the model, and in the definition of fast versus slow C fractions of leaf and root litter. Future empirical research that quantifies the protection rate by following the fate of labeled feedstock litter into protected and unprotected SOM pools has the potential to significantly improve model outcomes. Finally, relating modeled SOM fractions to empirically measurable pools is a welldocumented challenge stemming from uncertainty in both our conceptual understanding and ability to measure the pools (Abramoff et al., 2018;Campbell & Paustian, 2015;Dungait et al., 2012).
Differences among feedstocks in litter chemistry, belowground production, and management practices led to divergent estimates of the long-term trajectories between the perennial and annual systems in soil C (Figure 3). The three perennial systems all showed increasing soil C stocks due to no-till, establishment practices that limited soil disturbance, and low-quality, high C:N ratio litter that slowed the microbial decomposition of SOM leading to the accumulation of physically unprotected C. Among the perennial feedstocks, fertilized miscanthus had the highest soil C accumulation rate owing to greater root inputs (Table S5) and a lower C cost of N acquisition (Table S4) which limited rhizosphere decomposition. In contrast to the perennial systems, tillage and low root inputs led to long-term declines in soil C in the corn-corn-soybean simulation. These declines in soil C were counterbalanced by a greater ratio of protected-to-unprotected soil C in corn-corn-soybean than the perennial systems. This difference reflects the rapid decomposition and conversion of the high quality, low C:N ratio corn and soybean litter into microbial necromass which has the highest protection rate in the model. Overall, the model differences in soil C trajectories between perennial and annual bioenergy systems mirror previous modeling efforts (Davis et al., 2012;Duval et al., 2015;Hudiburg et al., 2015;Martinez-Feria & Basso, 2020) and reinforce the importance of perennial feedstocks to the sustainability and C balance of bioenergy production (Ledo et al., 2020;Whitaker et al., 2018).
Across all the feedstocks, the amount of NPP spent on N acquisition was relatively uniform ranging from 1.8% to 3.8%, values that are comparable to other ecosystems (Brzostek et al., 2013;Phillips & Fahey, 2005;Yin et al., 2014). However, feedstocks differed in the C cost of N acquisition and the dominant mode of N uptake (Table S4). These differences primarily resulted from varying growth rates and management among the crops. For example, switchgrass had the lowest C cost of N acquisition due to it having three times more belowground NPP than the corncorn-soybean rotation to scavenge nutrients, and it being able to retranslocate nearly 35% of its N from senescing biomass to support growth. In contrast, unfertilized miscanthus had the highest C cost of N acquisition. To support rapid growth, unfertilized miscanthus needed to spend F I G U R E 5 Total soil carbon (C) simulated by FUN-BioCROP at the University of Illinois at Urbana-Champaign Energy Farm under three sets of experimental conditions: control (black lines; same data as shown in Figure 3), 5°C soil warming (red lines), and increased rhizosphere C flux ("Rhizo C;" blue lines; augmented mycorrhizal C flux and non-mycorrhizal root C exudates by 5% of average daily net primary productivity, 10% total increase). Values on the figure indicate the difference between each experimental (i.e., warming, or increased rhizosphere C) final C pool value as compared to control, with the percent difference in parentheses. For example, increased rhizosphere C in corn-corn-soybean led to 18% higher soil organic C at the end of the simulation. Fert, fertilized; FUN-BioCROP, fixation and uptake of nitrogen-bioenergy carbon, rhizosphere, organisms, and protection; Unfert, unfertilized additional C on N fixation. Importantly, this reliance of miscanthus on N fixation to support growth supports hypotheses that N fixation is required to fill the "missing" N budget of miscanthus (Davis et al., 2010) and associated evidence of substantial N fixation in miscanthus even when soil N is non-limiting (Keymer & Kent, 2014). We observed this in our model output, as N fertilization in miscanthus led to greater N fixation, contrary to expectations. When miscanthus was fertilized, there was greater N fixation because of reduced C spent per unit N acquired through N fixation that corresponded to an overall reduction in the C cost of N acquisition (Table S4). Finally, corn-corn-soybean differed from the perennials because of its lack of retranslocation as well as the high N fertilizer additions allowing it to take up the largest amount of N directly through roots without the aid of mycorrhizae. These model estimates and the consequent impact of this C spent belowground on modeled soil C show that differences in nutrient acquisition strategies have the potential to be of nearly equal importance in driving biogeochemical cycling in bioenergy systems as what has been observed in natural ecosystems (Phillips & Fahey, 2007;Sulman et al., 2019;Terrer et al., 2016Terrer et al., , 2018. As such, our modeling exercise highlights an important avenue for future empirical research. The development of a mechanistic representation of tillage in FUN-BioCROP ("physical protection loss") represents an important advance in the modeling of bioenergy production that could enhance predictions once incorporated into more traditional models. Previous model efforts simulated the impacts of tillage on soil C by ramping up decomposition rate parameter values following a tillage event based on empirical observations ("decomposition scalars;" Maharjan et al., 2018). While this approach can be parameterized to accurately capture soil C losses due to tillage, it neglects the mechanism by which tillage enhances decomposition: the destabilization of protected soil C that increases microbial access to C (Jastrow & Miller, 1997). By imitating the real-world mechanism of soil C loss after tillage, FUN-BioCROP captured the magnitudes of SOC stocks without requiring site-specific parameterization (Figure 4a), making it more broadly applicable in larger ecosystem models. Furthermore, when we compared FUN-BioCROP results between routines using decomposition scalars and our novel mechanistic routine, we found that there was a large difference in the distribution of soil C in protected versus unprotected pools (Table 1; Figure 4b). The difference in the distribution of protected versus unprotected soil C between tillage routines is likely more important to model performance than the mismatch between the observed soil C and the simulated values estimated through decomposition scalars. Although most models that use decomposition scalars can match observed soil C values by parameterizing tillage impacts on decomposition rates, they do not have the capability to simulate how tillage impacts the distribution of protected versus unprotected soil C. Given that this distribution has been shown to control soil C response to global change (Lugato et al., 2021), we examined how differences between the tillage routines in estimates of the ratio of protected to unprotected soil C impacted the response of total soil C to warming. Elevated temperatures combined with the decomposition scalars tillage routine resulted in higher estimates of protected soil C (Table 1; Figure 4b) that were less vulnerable to loss with warming in the model than in the physical protection loss routine (Table 1; Figure 4a). These results suggest that mechanistic representations of tillage like the one we examined here have the potential to improve projections of bioenergy ecosystem responses to global change.
Simulations of feedstock responses to soil warming were highly dependent upon the distribution of protected versus unprotected soil C. We found that soil C in perennial systems was estimated to be lost at a higher rate than in the annual corn-corn-soybean system when soils were warmed ( Figure 5). This greater loss was caused by a large fraction of unprotected C in perennial systems that was preferentially lost as temperature increased the V max of microbial decomposition and uptake. In contrast, soil C in corn-corn-soybean was relatively stable with warming owing to limited unprotected soil C pools. These results are consequential because they show that while perennial feedstocks can enhance soil C stocks under current climate conditions (Anderson-Teixeira et al., 2009;Whitaker et al., 2018), there is the potential for the accumulated soil C to be more susceptible to global change. While few field studies have examined the differential sensitivity of unprotected (i.e., POM) and protected SOC (i.e., mineral associated organic matter, MAOM) to temperature (Lavallee et al., 2020), our results are in line with the findings from an incubation study which found POM decomposition was more sensitive to temperature than that of MAOM (Benbi et al., 2014). The few field studies that have examined this question produced mixed results, either showing no difference in the response of POM and MAOM (Schnecker et al., 2016) or greater loss of POM from forests but unexpectedly greater response to warming by MAOM in croplands (Lugato et al., 2021). This uncertainty in the empirical results as well as our modeling estimates highlight the need to empirically test whether the sensitivity of soil C to warming varies within ecosystems with significant gradients in plant traits.
Our model estimates showed that increasing the influx of rhizosphere C for all feedstocks increased both soil C pools and the amount of protected C relative to unprotected C ( Figure 5; Figure S3). This pattern estimated by the model coincides with empirical findings that plant rhizodeposition leads to greater formation of mineral-associated organic C, which persists relatively longer in the soil than particulate organic C Villarino et al., 2021). The similarity in modeled values and empirical observations reflects the mechanistic rhizosphere and microbial physiological process representations in FUN-BioCROP. In the model, when rhizosphere C enters the unprotected soil pools, it stimulates microbial growth and production of necromass, which has the highest protection rate. Increased decomposition, therefore, stimulated the production of relatively more protected soil C. Thus, selecting feedstocks or engineering new feedstocks with greater belowground C allocation could represent an efficient way to enhance SOC and improve the biofuel C balance.
FUN-BioCROP captured declines in soil C stocks following cultivation of native prairie as well as relative patterns in protected and unprotected soil C fractions, but there are numerous possibilities in parameterization and model development worth exploring to improve the applicability of the model framework we developed. First, although FUN-BioCROP simulated the relative proportion of protected to unprotected soil C for different feedstocks, the model overestimated the absolute values. This overestimation reflects the lack of empirical data to parameterize modeled protection rates. Experiments that follow the fate of isotopically labeled leaf litter into protected and unprotected soil C pools during decomposition are needed to reduce parameter uncertainty and improve model outcomes. Second, FUN-BioCROP currently models inorganic N as a single pool without differentiating between ammonium and nitrate. Plants and microbes can prefer different N forms (Fraterrigo et al., 2011;Harrison et al., 2007;Weigelt et al., 2005), with potential impacts on modeled N dynamics. Therefore, it is worth exploring the effect of simulating distinct ammonium and nitrate pools and fluxes on modeled N retention, an important ecosystem service. Third, while FUN-BioCROP advances our prognostic understanding of bioenergy systems by including microbial physiology, it only models one microbial pool that does not vary microbial traits or physiology. Recent empirical research shows that soil microbes differ in CUE, death rate, and their ability to use different substrate classes (Raczka et al., 2021). The challenge for the modeling community lies in distilling this communitylevel data into discrete values that can inform multiple microbial groups. Finally, across all these areas of future model development, there is also a clear need to balance model complexity with actual reductions in the uncertainty of model simulations. Given that FUN-BioCROP is significantly less complex than larger ecosystem models, it represents an ideal testbed to optimize complexity versus uncertainty tradeoffs. Thus, new theories and empirical findings can be distilled into simplified model representations that can be integrated into larger ecosystem models.

| CONCLUSION
FUN-BioCROP moves the field of bioenergy modeling forward by integrating new empirical paradigms of the role of belowground processes in shaping coupled C and N cycles into modeling representations that larger ecosystem models can ingest. Using a more mechanistic approach, FUN-BioCROP's representation of tillage as a process that enhances microbial access to soil C has important implications for modeling bioenergy sustainability. By altering the ratio of protected to unprotected C, FUN-BioCROP's simulation of tillage not only captures well-established differences in soil C stocks between annual and perennial feedstocks, but also raises new hypotheses that there could be divergent trajectories in the response of their soil C stocks to warming. Moving forward, the uniform increase in projected soil C stocks across all feedstocks when we increased the rhizosphere C flux highlights this trait as an important priority for feedstock development. Overall, FUN-BioCROP represents a critical step in updating the belowground processes in bioenergy models and moves the community closer to the goal of designing a modeling platform that can optimize management decisions (i.e., Which feedstock? How? & Where?) to ensure the sustainability of bioenergy production.