Corn stover harvest reduces soil CO2 fluxes but increases overall C losses

Harvesting corn (Zea mays L.) stover for use as bioenergy feedstock may provide short‐term economic benefits and perhaps improve grain yield in continuous corn systems, but excessive stover removal may lead to long‐term depletion of soil C stocks. To better quantify the impacts of stover harvest on the soil C balance, we investigated CO2 fluxes under three harvest treatments (none, moderate, and high) in three continuous corn systems: (a) no tillage and no biochar applications, (b) chisel plowing with biochar amendments, and (c) chisel plowing without biochar amendments were quantified. We sampled static chambers 14, 13, and 15 times in 2010, 2011, and 2012, respectively, and the measurements were used to calibrate the Agricultural Production Systems sIMulator model. Although CO2‐C emissions did not differ among the three systems, both moderate (approximately 30%) and high (approximately 60%) stover removal rates reduced simulated CO2‐C emissions by nearly 10% and more than 22%, respectively. Despite these reductions in CO2 flux, the sum of CO2‐C and stover C exceeded CO2‐C losses in those plots without stover removal. This finding suggested soil C depletion was occurring for the soil and crop management practices used, even though depletion was not evident after 3 years. Application of biochar increased soil C levels, suggesting it may be able to offset some C losses. We conclude that for sustainable, highly productive agroecosystems that include stover harvest, all agronomic practices must be optimized to minimize C losses and maintain soil quality.


| 895
O'BRIEN Et al. to the importance of SOC in delivery of agroecosystem services (Powlson et al., 2011), maintaining resilient production systems (Grandy, Fraterrigo, & Billings, 2012), and mitigating climate change effects (Lal, Follett, Steward, & Kimble, 2007). To stem or even reverse this trend, research quantifying C balances associated with current agronomic practices is essential.
Corn (Zea mays L.) stover harvest for biofuel and bio-product production, as well as for livestock feed and bedding, has increased significantly in recent decades Schmer, Brown, Jin, Mitchell, & Redfearn, 2017). Corn stover was identified as a potential U.S. biofuel feedstock because it is an abundant, renewable resource, and harvesting a sustainable amount may increase continuous corn grain yields in the upper Midwest (Rogovska, Laird, & Karlen, 2016;Sawyer, Woli, Barker, & Pantoja, 2017;Sindelar, Coulter, Lamb, & Vetsch, 2013). However, excessive long-term stover harvest has the potential to further decrease SOC stocks and adversely affect soil quality and soil health (Clay et al., 2015;Johnson et al., 2014;Laird & Chang, 2013;Ruis, Blanco-Canqui, Jasa, Ferguson, & Slater, 2017;Wilhelm, Johnson, Karlen, & Lightle, 2007). Greater C inputs from root biomass and exudates in high-yielding systems may partially offset C losses from stover removal (Clay et al., 2015;Johnson, Allmaras, & Reicosky, 2006), but stover removal will decrease total C input and may impact soil thermal and moisture regimes. This removal may alter both the quantity and composition of SOC, negatively affect soil health indicators, and further increase SOC losses via wind and water erosion (Newman, Kaleita, & Laflen, 2010;Thomas, Engel, & Chaubey, 2011). Furthermore, reducing C input may also affect SOC mineralization rate, as quantified by soil CO 2 fluxes (Jin et al., 2014).
Stover removal affects soil CO 2 fluxes in three ways. First, stover is a major microbial substrate, so removal not only reduces C input but also other nutrient inputs that influence microbial respiration (Lehman & Osborne, 2016). Second, stover removal alters the soil microclimate by exposing the soil surface and removing the layer of crop residue at the soil-atmosphere interface. These changes typically increase soil temperature and reduce soil water content (Duiker & Lal, 2000;Horton, Bristow, Kluitenberg, & Sauer, 1996;Kenney et al., 2015), two factors that have a strong influence on soil respiration rates (Borken & Matzner, 2009;Parkin & Kaspar, 2003). Finally, increases in grain yield following stover removal are typically accompanied by increased aboveground and belowground biomass production, and the quantity of root biomass and exudates also impacts total soil respiration (Johnson et al., 2006).
Interactions among the effects of stover removal, coupled with variability inherent in management and environmental factors (Jin et al., 2014), have led to inconsistent CO 2 flux response following stover harvest. For example, Baker, Fassbinder, and Lamb (2014) reported CO 2 -C emissions of 1,020 versus 1,103 g C/m 2 following 100% versus 0% removal when averaged for two seasons in MN. The difference (83 g C/m 2 ) was far less than total C removal via stover harvest, suggesting a negative C balance. Similarly, mean CO 2 fluxes were lower in a four year South Dakota study with 55% compared to 0% stover removal, which was attributed to reduced C inputs (Lehman & Osborne, 2016). Conversely, in another 3 year South Dakota study, removing crop residue increased soil CO 2 fluxes compared to the non-removal treatment. That finding was attributed to higher soil temperatures and increased soil microbial activity in plots where residue was removed (Wegner et al., 2018).
Application of biochar, a C rich byproduct of thermochemical bioenergy production systems, has been suggested as one means of compensating for the adverse effects of crop residue harvesting on SOC stocks and soil quality . Most of the C in biochar is highly recalcitrant to microbial degradation and may persist in soils for hundreds or even thousands of years; therefore, biochar applications directly increase soil C stocks (Laird et al., 2010;Zimmerman, Gao, & Ahn, 2011). The impact of biochar applications on soil respiration is variable in time and space. Biochar applications typically cause an initial pulse of CO 2 emissions (Jones et al., 2011), which is due to the combined effects of tillage used to incorporate the biochar, the hydrolysis of carbonates, and the mineralization of labile organic compounds that are added with the biochar (Fidel, Laird, & Parkin, 2017). The long-term effects of biochar applications on soil respiration are less consistent, as biochar may either enhance mineralization of biogenic soil organic matter (positive priming) or reduce the rate of biogenic SOC mineralization (negative priming) and enhance stabilization of C added to soils with crop reside inputs (Blanco-Canqui, Laird, Heaton, Rathke, & Acharya, 2020;Chang et al., 2016;Rogovska et al., 2016). Thus, interactions between biochar applications and stover removal and their impact on soil CO 2 fluxes require further investigation.
Individual stover harvest effects on grain yield, soil properties, and CO 2 fluxes, as well as the interactions among those effects, vary temporally and across locations due to different climates and management scenarios Raffa, Bogdanski, & Tittonell, 2015). Models integrating these multiple effects and interactions may therefore be useful tools to complement field measurements and better describe stover harvest effects on soil CO 2 fluxes (Liska et al., 2014;Lugato & Jones, 2015). For this study, we used the Agricultural Production Systems sIMulator (APSIM) model (Holzworth et al., 2014;Keating et al., 2003) to simulate daily CO 2 efflux for different tillage and stover harvest scenarios. The APSIM platform has been used to successfully simulate crop development and yields Togliatti, Archontoulis, Dietzel, Puntel, & VanLoocke, 2017), nitrogen dynamics and leaching Martinez-Feria et al., 2018;), as well as long-term soil organic matter dynamics (Aller et al., 2018;Basche et al., 2016) in corn-soybean (Glycine max (L.) Merr.) systems in Iowa. Modeling CO 2 fluxes using this predictive framework is useful because direct CO 2 flux measurements are often limited to a single weekly measurement, especially in complex or large-scale field experiments. Automated chambers that provide continuous data are available, but due to cost and/or logistical constraints, they are often not widely applied. Therefore, a predictive framework that accounts for major determinants of CO 2 fluxes, such as soil organic matter, soil temperature, soil water content, and plant root growth, can be valuable in expanding our understanding of season-long C dynamics.
Identifying appropriate stover management practices is essential for maximizing economic gain associated with increasing stover demand for livestock feed and/or biofuel feedstock while maintaining SOC stocks needed for long-term yield stability and soil health. Choosing these practices based on site-specific characteristics is crucial to avoid potential detrimental soil property effects (e.g., increased soil erosion or C stock depletion) associated with excessive stover harvest (Blanco-Canqui, 2013;Thomas et al., 2011;Wu & Liu, 2012). We evaluated soil CO 2 fluxes under three levels of stover harvest (none, moderate, and high) in continuous corn systems using (a) no tillage and no biochar applications or (b) chisel plowing with biochar amendments and (c) chisel plow without biochar applications. Our hypotheses were that: (a) stover removal would reduce the overall amount of C lost via soil CO 2 fluxes, and (b) the magnitude of response would vary based on tillage intensity and biochar application. By combining experimental field data with a comprehensive modeling framework, our goal was to quantify C budgets and better assess long-term impacts of tillage and stover removal on soil health for this U.S. Midwestern site.

| Study area, agronomic management, and experimental setup
This study was conducted at the Iowa State University Agricultural Engineering and Agronomy Research Farm (42.017584°, −93.76448°), with soils consisting primarily of Webster silty clay loam (fine-loamy, mixed, superactive, mesic Typic Endoaquolls) on 0% to 2% slopes and Clarion loam (fine-loamy, mixed, superactive, mesic Typic Hapludolls) on 2% to 5% slopes. Average growing season temperature (April through September) is 17.7°C, and typical annual precipitation is 990 mm (Iowa Environmental Mesonet, 2019); climate data for 2010-2012 are shown in Table 1. The research plots (85.3 m × 12.2 m each) had been under continuous corn production with stover removal since 2008. The plots were planted each spring (79,000 kernel/ha; 76 cm row spacing) and received annual nutrient applications of N (urea ammonium nitrate), P (triple super phosphate), and K (potassium chloride) according to either the amount of stover removed (N) or soil test levels (P and K; Table 1). Most N (80%-90%) was applied during sidedress applications, with smaller amounts applied as starter N at planting or with P and K fertilizers prior to fall tillage (Table 1). Corn grain and stover were harvested simultaneously using single-pass equipment .
We evaluated three management systems with three corn stover harvest rates and four replications (i.e., 36 plots) from 2010 to 2012. The three management systems were no tillage (NT), conventional tillage with fall chisel plowing and spring disking (CT), and CT plots amended with biochar (BC). Biochar produced by slow pyrolysis of mixed hardwoods was applied to select plots in 2007 at a rate of 18.4 Mg dry matter/ha (13.1 Mg C/ha) and was incorporated immediately following application. More information regarding biochar characteristics can be found in Rogovska et al. (2016). For the 3 years, mean (±standard deviation) corn stover harvest treatments in Mg/ha were none (0 ± 0), moderate (3.1 ± 1.1), or high (5.3 ± 1.3), respectively. These rates reflected 0%, 50% and 100% of what was mechanically feasible, with actual stover removal ranging from 20% to 40% and 50% to 70% of the aboveground biomass for the moderate and high treatments, respectively (Obrykci, Karlen, Cambardella, Kovar, & Birrell, 2018).

| CO 2 , soil temperature and water, and weather measurements
We installed PVC collars (20 cm diameter, 12 cm height) to a depth of 9 cm (i.e., 3 cm of collar exposed above soil surface) within each plot at both in-row and between-row locations. Twelve holes (2.54 cm diameter) centered at 4 cm below the soil surface allowed for lateral fluid movement and root growth inside the collars. Collars were initially installed after planting in 2010 and were removed and reinstalled to accommodate farm practices (e.g., tillage, planting, fertilizer application). For each sampling period, we sampled each collar using a static chamber (LI-8100-103, LI-COR Biosciences) connected to an infrared gas analyzer (LI-8100, LI-COR Biosciences), which has a sampling volume of 4,843 cm 3 and requires a measurement time of approximately 60 s. In addition to the gas measurement period, each sampling period included a prepurge and postpurge cycle of approximately 30 s to avoid errors introduced by prior measurements. We sampled each collar on 14, 13, and 15 different days in 2010, 2011, and 2012, respectively. Due to weather conditions and instrument constraints, the amount of data reported varies from four to eight measurements per treatment per sampling date. For each sampling interval, we averaged values from in-row and between-row locations to provide a mean value for each plot. All gas measurements were taken between 8:00 a.m. and 10:00 a.m. to best represent the mean daily flux (Parkin & Kaspar, 2003). In continuous chamber measurements from select plots at this experimental site, no daily means exceeded 8.5 μmol CO 2 /m s (data not presented). Since our objective for the chamber samples was to represent the daily mean flux, measurements exceeding 9 μmol CO 2 /m s (44 measurements, 1.5% of total dataset) were deemed not reflective of the daily mean and were therefore rejected from analysis to reduce model error. When the collars were present, soil temperature and volumetric water content were measured near each at a depth of 5 cm (5TM, METER Environment). Measurements were taken hourly and stored in data logger (EM50, METER Environment).

| Soil and plant sampling and analysis
At harvest each year, we collected a soil core (7.3 cm diameter) within each collar using a slide hammer to a depth of 15 cm. Cores were divided into two depth increments, 0-5 and 5-15 cm, and assessed for organic C, total N, and bulk density. After screening to ensure no inorganic C was present, soil C and N were quantified by dry combustion method (Flash 1112EA Nitrogen and Carbon analyzer; Thermo Fisher Scientific). Bulk density was determined as the ovendry soil mass in a known volume from the soil cores. Grain and harvested stover subsamples were also analyzed for total C by dry combustion.
We customized model input for a series of practices to match site-specific management, including tillage (depth and intensity), planting (variety, rate, depth, and spacing), fertilizer application (rate, type, and timing), and dates for grain harvest and residue removal. We modeled each plot individually because of variation in soil types, fertilizer application rates, biochar application rates (BC plots only), and measured CO 2 fluxes. When a plot consisted of more than one soil type, we created a simulation for each soil type. We calculated the final output for each plot based on weighted averages determined by the areal proportion of each soil type within each plot. We began all simulations on 20 September 2007 and ceased simulations on 31 December 2012. We obtained climate data for this period from the Iowa Environmental Mesonet (2019).
APSIM simulates CO 2 emissions as the sum of mineralization of soil organic C and the decomposition of fresh inputs (i.e., residue, root biomass). Soil organic C is partitioned into three pools: microbial biomass C (fast decomposing), humic soil C (moderately decomposing), an inert soil C (no decomposing), and a biochar C pool (slow decomposing; Archontoulis et al., 2016;Probert et al., 1998). The rate of decomposition is driven by a potential decomposition rate, which is mediated by C:N ratios and soil temperature and water content modifiers at each soil layer. Additionally, APSIM simulates a fourth pool of potential CO 2 emissions from the fresh organic matter, which includes crop residue and roots (Probert et al., 1998;Thorburn et al., 2001). A fifth pool of potential soil CO 2 emissions is the autotrophic root respiration, but APSIM does not account for this pool because it uses the radiation use efficiency approach to drive biomass production. Given that root respiration can comprise a significant portion of overall CO 2 fluxes (Werth & Kuzyakov, 2009), we estimated the daily contribution of root respiration (R root ) using a Q 10 equation (Equation 1; Teixeira, Moot, & Brown, 2009), where M root is the dry weight biomass of the root (simulated by APSIM), 0.48 is the conversion factor for biomass C, Q 10 is estimated at 2.2, and T soil is the simulated soil temperature at 5 cm depth, and T ref is the reference soil temperature (18°C; Teixeira et al., 2009).
Thus, the simulated CO 2 flux output was the sum of CO 2 emitted from these five pools: microbial biomass C, humic soil C, residual soil C, fresh organic inputs, and root respiration.
We calibrated the model based on the measured values collected for each plot by manipulating the water table depth to mimic known trends from nearby experimental locations  and the temperature dependence of the root respiration to improve model performance. We focused on calibration of the root respiration because our initial model runs showed good agreement between simulated values and measured values early in the growing season before root respiration became a factor ( Figure S1). Additionally, we refrained from manipulating the four pools of soil respiration simulated by APSIM because changing the soil C parameters is complex (see Archontoulis et al., 2016), and it is risky to alter soil C coefficients without robust and multi-faceted datasets. Thus, we maintained the documented robustness of APSIM in simulating long-term soil C dynamics that has been documented across a range of conditions (e.g., Luo et al., 2014;Luo, Wang, Sun, Smith, & Probert, 2011;Puntel et al., 2016).
A single corn cultivar of 105 day maturity was used for all years from the APSIM database. We assessed model performance using root mean squared error (RMSE; Equation 2) and relative RMSE (RRMSE; Equation 3), the absolute and relative error, respectively, where performance is improved as these values approach 0 (Archontoulis & Miguez, 2015), as well as model efficiency (ME; Equation 4).
where n is the number of observations, R i is the measured value for a given parameter, R avg is the average of measured values, and R * i is the modeled value for a given parameter.

| Statistical analysis
We analyzed soil and plant samples, as well as cumulative CO 2 fluxes, using three-way analysis of variance (ANOVA) using year, system, and stover removal rate as factors with the R 3.6.2 statistical environment (R Core Team, 2018). When results were significant at the α = 0.05 level, we conducted pairwise comparisons within each factor using a post hoc Tukey test with package "agricolae" (de Mendiburu, 2017). Each CO 2 sampling date was analyzed using a repeated measures three-way ANOVA within each system, and we used a post hoc Tukey test to identify differences between stover removal rates for each sampling date.

| CO 2 survey chamber measurements
The CO 2 measurements followed expected annual temperature trends, where respiration was highest during the summer and very low in early spring, late fall, and winter ( Figure 1). Stover removal rate rarely affected CO 2 fluxes for the CT or NT systems; fluxes were diminished under higher stover removal relative to no stover removal rates only three times throughout the 3 year study for the CT and NT systems. We identified consistent significant differences between stover removal rates for only one system, BC, in one year, 2010 ( Figure 1).

| Soil C and N, microclimate, and grain yields
Soil C and soil C:N ratios were similar for all stover removal rates in all three systems except for minor differences in 2011 (Table 2), likely due to inherent variability in sampling for soil C; no differences in soil C and soil C:N ratios were evident in either the NT or CT systems.
Differences among the removal rates and systems in soil temperature were highest (up to 5°C) early in the growing season ( Figure 2), with minimal differences as the season progressed. Similarly, the greatest differences in soil water content were observed early in the growing season (up to 0.10 m 3 /m 3 ; Figure 2). (1) .
( While precipitation in 2011 was close to the 30 year average (828 mm), 2010 had nearly twice as much precipitation as 2012, 1,227 mm and 637 mm, respectively (Table 1). These trends in precipitation were reflected in grain yields (Table 4). Overall, grain yields were not significantly different among either systems or removal rate, but the interaction of system and rate was significant (Table 2). In CT systems, stover removal increased yields compared to those with no stover removal in all 3 years, while it only increased yields in BC and NT in the wettest year, 2010 (Table 4).

| Model performance
The APSIM model simulated CO 2 fluxes for 2010, 2011, and 2012 with an RMSE of 15.96 (kg CO 2 -C/ha), RRMSE of 0.67, and ME of 0.35. The model performed consistently among all systems and stover removal rates in simulating CO 2 fluxes (Table 3), as well as in simulating soil temperature, volumetric water content, and amount of stover removed (Table S1). The model showed a slight tendency to underpredict high measured values (i.e., >50 kg CO 2 -C/ha), though these values accounted for only 10% of the measured dataset. These high values most often occurred during the middle of the growing season ( Figure 1) when contribution from root respiration was highest. Despite this tendency, the simulated values closely matched yearly trends of measured values for each treatment (Figure 5), where CO 2 fluxes increase with soil temperature and plant growth during the summer, while fluxes decrease as temperatures decline and root growth ceases.

| Cumulative CO 2 flux estimation
Simulated cumulative CO 2 -C emissions were different among years and stover removal rates, but they did not differ among systems (Table 2). Despite diminished yields and, likely, root production, CO 2 fluxes remained high in 2012. Notably, CO 2 -C emissions in CT did not decrease with moderate stover removal in any of the 3 years, which contrasts with decreases in both NT and BC in 2011 and 2012 (Table 5).

| CO 2 survey chamber measurements
The annual trends in measured CO 2 fluxes (Figure 1) are caused by the exponential relationship between soil heterotrophic respiration and temperature and the contribution of root respiration during the peak of the growing season (Lloyd F I G U R E 2 Soil temperature (above) and soil volumetric water content (below) in plots from three systems (conventional tillage, no tillage, and biochar-amended plots with conventional tillage) under different stover removal rates (none, moderate, high , 1994;Rochette, Desjardins, & Pattey, 1991;Xu & Baldocchi, 2004). The lack of consistent response in CO 2 fluxes to stover removal contrasts with other researchers who found that stover removal caused CO 2 fluxes to either decrease (Wegner et al., 2018) or increase (Jin et al., 2014). The lack of significant differences between either stover removal treatments or tillage systems in this study may be partially explained by high variability in all the plots, as every treatment had coefficients of variation exceeding 42% (error bars not shown in Figure 1 to improve readability). The only consistent response was observed in BC in 2010, and the lower fluxes in this system-year under both moderate and high stover removal may be partially explained by the potential of biochar to have a negative priming effect on the mineralization of soil C (Blanco-Canqui et al., 2020), such that the elevated CO 2 fluxes are caused predominantly by the increased residue inputs; however, this effect did not persist in 2011 or 2012.
Many of the soil properties that govern CO 2 fluxes were not strongly affected by the stover removal treatments in this study (Table 2), possibly explaining why our survey chamber measurements were not sensitive to stover removal rates on these Mollisols with relatively high levels of SOC. For example, decreases in CO 2 emissions following stover removal have been attributed to changes in soil C:N ratios Jin et al., 2014;Johnson et al., 2006). Soil C and C:N ratios were significantly higher in the BC plots relative to the NT and CT plots, and this difference is attributed to the presence of recalcitrant biochar C. Biochar applications often cause a short-term increase in CO 2 flux rates (Fidel et al., 2017;Jeffery, Verheijen, van der Velde, & Bastos, 2011). Since CO 2 fluxes in BC were not elevated compared to the other systems until the summer of 2010 (Figure 1), ~3 years after the Fall 2007 biochar application, any short-term pulse in CO 2 fluxes caused by the biochar application had ended.
In addition to decreasing the total amount and changing the composition of microbial substrate added to the soil, stover harvesting also alters the microclimate at the soil surface, which may increase CO 2 fluxes (Wegner et al., 2018). Removal of crop residue increases the amount of bare soil, which may increase soil temperatures and decrease soil water content, as well as reducing boundary resistance between the soil surface and the atmosphere (Swan, Kaspar, & Erbach, 1996). Increasing temperatures typically increases soil CO 2 fluxes (Parkin & Kaspar, 2003), while the influence of soil water content may vary widely (Castellano et al., 2011;Daigh, Sauer, Xiao, & Horton, 2014;Kaspar & Parkin, 2011). In our study, differences in soil temperature likely did not affect CO 2 fluxes during the chamber sampling because temperature differences among the stover harvesting treatments were only evident during the beginning of the growing season (Figure 2). Once the crop canopy closed, soil temperatures were similar among all treatments and systems for the remainder of the year. Soil CO 2 fluxes showed a strong relationship with soil temperature for all treatments (Figure 3). Since soil CO 2 fluxes are low when soil temperatures are low, the early season differences that may have been present in soil temperature (Figure 2) had little effect on the overall soil CO 2 emissions. Similarly, the weak relationship between soil water content and CO 2 fluxes indicates that differences in water content did not consistently affect soil CO 2 emissions.
No measures of root growth were collected in this study, but CO 2 from root respiration can exceed 40% of total CO 2 fluxes under corn during the growing season (Werth & Kuzyakov, 2009). Measurements from nearby fields (Nichols et al., 2019;Ordonez et al., 2018) indicated that corn root mass is about 1 Mg/ha and reaches a peak mass, and therefore peak respiration, around the R3 to R4 corn growth stages (typically between DOY 205-225). Root respiration was assumed to occur at a similar rate among all treatments, such that differences would likely only occur due to different root biomass and rhizomicrobial respiration. To assess differences in root biomass without a direct measurement, some research has linked belowground C dynamics with aboveground C dynamics using root:shoot ratios (Johnson et al., 2006). Even though grain yields were increased with stover removal in CT throughout the study (Table 4), as well as in NT and BC in 2010, the magnitude of increase suggests that additional root biomass C and rhizodeposited C in higher yielding plots are unlikely to offset losses from stover removal. This assumption is supported by research showing that increased root biomass does not always equate to increased soil respiration (Hirte, Leifeld, Abiven, Oberholzer, & Mayer, 2018;Nichols, Miguez, Sauer, & Dietzel, 2016).
Differences in three major drivers of soil CO 2 flux, soil microbial substrate availability, soil microclimate, and vegetation (root respiration), were either offsetting or not large enough to cause quantifiable differences in CO 2 flux in the chamber surveys, with the sole exception of BC in 2010. This finding may be partially the result of sampling schedule. For example, large pulses of CO 2 emissions may be expected immediately following tillage (Reicosky, Dugas, & Torbert, 1997) or biochar applications (Jones et al., 2011), but no chamber samples showed evidence of those pulses. These pulses were likely absent because they are typically short-lived, and no sampling events occurred within 7 days of either fertilizer application or tillage (Table 1; Figure 1). Additionally, the sampling schedule was based on technician availability, such that chambers were surveyed on a given day of the week rather than on days with favorable sampling conditions. This type of schedule removes bias that may occur when fluxes are measured only on sunny, warm days that may show the largest treatment differences.
Finally, CO 2 fluxes are often highly variable in space and time, and individual measurements may reflect only subtle differences. Rather, the true effect of these different management strategies may only be evident in the cumulative fluxes over the course of an entire growing season. Identifying treatment differences with direct measurements is difficult because the survey chambers capture only a small fraction of each research plot, and the time-intensive nature of sampling typically makes capturing temporal variability unfeasible. One path forward may be to incorporate more automatic chambers that can provide much greater temporal resolution, though they may be costly to F I G U R E 3 Relationship between measured CO 2 flux values and soil temperature (left) and soil volumetric water content (right) measured in survey chambers from 2010 to 2012 in plots from three systems (conventional tillage, no tillage, and biochar-amended plots with conventional tillage) under different stover removal rates (none, moderate, high). ***p < .001, ns, not significant For each year, different UPPERCASE letters indicate differences between stover removal rates (i.e., differences within each row), whereas different lowercase letter indicate differences between systems (i.e., differences within each column) implement and risk creating chamber-induced environmental effects that can affect CO 2 flux measurements (Gorres, Kammann, & Ceulemans, 2016). Since we did not have the ability to improve either spatial or temporal resolution of our CO 2 measurements, we chose to incorporate a modeling framework, APSIM, that could simulate CO 2 fluxes between sampling periods. While we did not observe consistent statistically significant differences in our survey chamber measurements, we were able to identify relationships between management, soil type, soil temperature, soil water content, and vegetation production that could be used to calibrate APSIM to improve model performance for these treatments.

| Model performance
The parameters indicating model performance, RMSE, RRMSE, and ME, have substantially higher values for CO 2 fluxes than those for soil temperature and grain yield (Table 3; Table S1), though the ME is similar to the simulations for soil volumetric water content (Figure 4). Similarly, these values exceed errors reported in other studies for crop biomass (RRMSE = 0.13; Dietzel et al., 2016) or soil temperature (RRMSE = 0.12; Basche et al., 2016), and the CO 2 flux RRMSE for each year is greater than 0.20, a threshold introduced in Ma et al. (2011) for what may been considered satisfactory for models of many agricultural variables. However, CO 2 fluxes have much greater spatiotemporal variability than yield or biomass production, such that these thresholds for model performance may be unsuitable. For example, Iqbal et al. (2018) used the DAYCENT model to predict CO 2 and N 2 O fluxes and had comparable RMSE values (ranging from 11 to 19.5 kg CO 2 -C/ha d) and ME (ranging from 0.25 to 0.37). Additional uncertainty is introduced to APSIM by the assumption that a single measurement of CO 2 flux is used to represent the daily flux for comparison with the model.
Nonetheless, we have confidence in the model performance because both the RMSE and cumulative fluxes are comparable to other CO 2 flux models used on similar soil types near these research plots (Daigh et al., 2014;Daigh, Sauer, Xiao, & Horton, 2015). Additionally, visual examination of temporal patterns ( Figure 5) may be more beneficial in gauging model performance than RMSE values, which are typically higher for complex soil processes. While these simulations are a rigorous effort at quantifying cumulative CO 2 fluxes, the higher uncertainty (e.g., RRMSE, ME) compared to other parameters (e.g., yield, soil temperature) reflect the need to continue improving both the direct measurement and estimation of CO 2 fluxes.

| Cumulative CO 2 flux estimation
The simulated cumulative CO 2 fluxes differed among years (Table 2), likely driven by an interaction of weather variability and plant inputs. The higher fluxes in 2012 were likely due to increased temperatures during the growing season compared to the other 2 years (Table 1), despite lower plant growth. Conversely, the simulated cumulative CO 2 fluxes did not differ among systems. This finding was unexpected, as tillage typically increases CO 2 flux compared to no till systems (McCourty, Gyawali, & Stewart, 2018), and biochar has been shown to reduce CO 2 fluxes by enhancing stabilization of crop residue C (Blanco-Canqui et al., 2020;Chang et al., 2016). However, the sustained levels of CO 2 flux in NT systems due to greater residual soil C and more stable pore networks appears to offset any large CO 2 burst associated with tillage in NT or BC that may have occurred.
The lack of differences between the systems may also be partially explained by the interaction effect of system and stover removal (Table 5). Increasing stover removal rate resulted in greater decreases of CO 2 fluxes in both NT and BC, but fluxes in CT were not different between no removal and moderate removal. Stover removal may be expected F I G U R E 4 Simulated daily values for soil temperature (T; panel a), soil volumetric water content (VWC; panel b), and grain yield (panel c) plotted against experimental values. These plots show all data combined from 2010 to 2012 for plots from three systems (conventional tillage, no tillage, and biochar-amended plots with conventional tillage) under different stover removal rates (none, moderate, high). ME, model efficiency; RMSE, root mean squared error; RRMSE, relative root mean squared error to have a greater impact in NT since it both alters gas exchange at the soil-atmosphere interface and eliminates substrate for mineralization. In tilled systems, conversely, the soil-atmosphere interface is altered by tillage and residue is incorporated into profile, such that the only difference is the amount of residue available for decomposition (Table 2). The moderate rate of removal, therefore, is not high enough to create a difference in the CT system, which may also be influenced by the increased crop production following stover removal in CT (Table 4). The BC system, although subjected to the same tillage management as CT, responded differently. This difference is likely the result of F I G U R E 5 Examples of yearly APSIM model output from each treatment for 2011, where the shaded area represents simulated fluxes and points (with standard error bars) represent measured data from three systems (conventional tillage, no tillage, and biochar-amended plots with conventional tillage) under different stover removal rates (none, moderate, high)

| 905
O'BRIEN Et al. biochar catalyzing the stabilization of crop residue and root C, an effect that is more apparent when stover is removed (Blanco-Canqui et al., 2020).

| Implications for long-term soil C balance
While CO 2 fluxes declined as stover was removed, they did not decline enough to offset the amount of C removed from the system via stover harvest. Thus, the combination of CO 2 efflux and stover removal results in increased C losses in both tilled systems ( Figure 6), with an overall average close to 1,000 kg C/ha lost each year. This annual loss of 1,000 kg C/ha is not an estimation of system-wide C losses because CO 2 -C and stover C are not the only components in the C cycle of these systems. Notably, this value does not account for increases in CO 2 -C that may result from increased root biomass and associated heterotrophic respiration; in these cases, increased CO 2 -C losses via emissions may be offset by increased C inputs from belowground biomass and likely do not represent a net decrease in soil C.
Nonetheless, these data indicate that stover removal contributes to long-term losses in soil C. These differences may not have been reflected in soil C at harvest sampling (Table 2) because residue decomposition can comprise widely variable fractions of soil CO 2 fluxes, ranging from 20% to 70% (Glenn et al., 2011), and soil C pools may be slow to respond to management factors, such that a longer study period is required to quantify soil C changes. However, just as Figure 6 does not account for belowground C inputs that would offset some CO 2 -C losses, nor does it account for other annual C losses associated with cropland production (Manies, Harden, Kramer, & Parton, 2001), especially those exacerbated by tillage (Blanco-Canqui, 2013;McCourty et al., 2018;Thomas et al., 2011).
On the other hand, the BC plots had significantly higher SOC levels (~3 g/kg) relative to SOC in the CT plots (Table 2). This finding suggests that applying biochar partially compensates for the loss of C due to residue harvesting, as the results show that CO 2 fluxes from the biochar plots were not elevated relative to CO 2 fluxes from the NT or CT plots. Thus, biochar application in stover removal systems is a promising strategy to mitigate declines soil C levels from stover removal and standard agricultural management, even with only periodic application (e.g., every 5 years).

F I G U R E 6
Combined losses of C via simulated soil CO 2 -C fluxes (filled portion) and Measured stover removal (open portion) from 2010 to 2012 in plots from three systems (conventional tillage; no tillage; biochar-amended plots with conventional tillage) under different stover removal rates (none, moderate, high). Error bars correspond to sum of soil CO 2 -C and stover C. Different letters indicate significant differences compared among all nine treatments at α = 0.05 level. The 3-year cumulative pane denotes the sum of daily CO 2 -C losses via these two processes from 1 January 2010 through 31 December 2012

| CONCLUSIONS
Corn stover removal decreased soil CO 2 -C losses each growing season, but not enough to offset C removed from the system from stover harvesting. Corn grain yield and presumably belowground biomass production in this study were not sufficient to support stover removal, thus increasing net C loss for all three systems. Application of biochar increased soil C levels and reduced C mineralization, but persistence of those benefits is uncertain without long-term studies. Furthermore, CO 2 emissions and stover removal are not the only sources of C loss, especially if soils are tilled. Therefore, without using cover crops or other conservation practices to significantly increase crop yields, the consistent response across time and the management systems we evaluated suggests that annual stover removal will ultimately reduce soil C levels and diminish potential productivity.