Boosting climate change mitigation potential of perennial lignocellulosic crops grown on marginal lands

Renewable lignocellulosic biofuels are one of the few mature technologies with significant potential to reduce fossil fuel consumption and greenhouse gas (GHG) emissions in the transportation sector (Gelfand et al 2020). Supplying sufficient lignocellulose feedstock to satisfy mandates for low-C biofuels will require massive land-use change, possibly involving the conversion of millions of hectares of existing cropland and marginally productive sites into perennial biomass crops (Gelfand et al 2013). Land-use change in such a scale, however, could have significant impacts on food security, biodiversity, and GHG emissions if bioenergy production directly or indirectly displaces food crops or sensitive habitats (Fargione et al 2008, Melillo et al 2009, El Akkari et al 2018). Furthermore, the potential carbon neutrality of lignocellulose-derived biofuel hinges on the ability of the production system to offset upstream and downstream supply-chain GHG emissions with the sequestration of atmospheric CO₂ in soil organic carbon (SOC) pools (Robertson et al 2011), which is a process highly dependent on the pedo-climatic conditions (Agostini et al 2015, Paustian et al 2016). Thus, careful selection of sites and management practices used for feedstock cultivation has become the core of strategies seeking the sustainable and efficient deployment of lignocellulose-based bioenergy systems (Somerville et al 2010, Gelfand et al 2013, 2020, Field et al 2018, O'Neill and Maravelias 2021).

Nitrogen fertilizer (N_F) stands as one of the major uncertainties in the production of perennial bioenergy crops. Application of N_F can increase plant growth and agronomic yields (Mclaughlin and Kszos 2005, Ruan et al 2016), and the additional biomass may contribute to climate mitigation by increasing availability of biofuel feedstocks and improving efficiency of biofuel supply chains (O'Neill and Maravelias 2021). However, the responses of perennial bioenergy crops to N_F additions are quite variable and poorly understood (Fike et al 2017, Monti et al 2019). Enhancing plant productivity with fertilizer additions also generally means increased C inputs to soils via plant litters (Miesel et al 2017), which may result in greater rates of SOC sequestration (Poffenbarger et al 2017). Yet, evidence suggests that in environments with low soil N availability to plants and microbes, N_F additions can accelerate the rates of SOC decomposition and reduce soil C stocks, particularly in the mineral-associated fractions (Valdez et al 2017), with large uncertainty of whether N_F will result in higher or lower SOC sequestration at a given candidate site. Moreover, N_F also promotes microbial nitrification and denitrification processes, resulting in the release of nitrous oxide (N₂O), a potent GHG with a global warming potential ∼265 times that of CO₂ (IPCC 2014). These emissions are on top of those incurred by the fossil fuel energy used for manufacturing, transportation, and application of the N_F itself, which on the aggregate could negate the GHG mitigation benefit from any additional feedstock produced (Crutzen et al 2008).

From a climate mitigation perspective, the question of whether bioenergy crops ought to receive N_F—and at what rate—needs to balance the benefits and costs of the practice in terms of GHG emissions. The objective here is to minimize the net GHG contribution of N fertilization to the life-cycle greenhouse warming intensity (GWI, Mg CO₂e ha⁻¹ yr⁻¹) of the producing lignocellulosic biofuel. That can be accounted as the differences in GWI between the fertilized (${\text{GW}}{{\text{I}}_{{{\text{N}}_{\text{F}}}}}$) and unfertilized (${\text{GW}}{{\text{I}}_{{{\text{N}}_{\text{0}}}}}$) systems, as show in equation (1):

$$\begin{align}{\text{GH}}{{\text{G}}_{{\text{Net}}}} = {\text{GW}}{{\text{I}}_{{{\text{N}}_{\text{F}}}}} - {\text{GW}}{{\text{I}}_{{{\text{N}}_{\text{0}}}}}.\end{align} \tag{ 1 }$$

Solving this expression (methods) reveals that ${\text{GH}}{{\text{G}}_{{\text{Net}}}}$ can be quantified if we know how much additional GHG is emitted due to N_F synthesis and distribution (${\text{GH}}{{\text{G}}_{{\text{FertSyn}}}}$), field application (${\text{GH}}{{\text{G}}_{{\text{FertApp}}}}$), the changes in soil N₂O emissions ($\Delta {\text{GH}}{{\text{G}}_{{{\text{N}}_{\text{2}}}{\text{O}}}}$), SOC sequestration ($\Delta {\text{GH}}{{\text{G}}_{\Delta {\text{SOC}}}}$), and the negative emissions due to the additional replacement of gasoline by the biofuel ($\Delta {\text{GH}}{{\text{G}}_{{\text{BioFuel}}}}$):

$$\begin{align}{\text{G}}{{\text{G}}_{{\text{Net}}}}& = {\text{ GH}}{{\text{G}}_{{\text{FertSyn}}}} + {\text{ GH}}{{\text{G}}_{{\text{FertApp}}}} + { }\Delta {\text{GH}}{{\text{G}}_{{{\text{N}}_{\text{2}}}{\text{O}}}}\nonumber\\ & \quad + \Delta {\text{GH}}{{\text{G}}_{\Delta {\text{SOC}}}} + \Delta {\text{GH}}{{\text{G}}_{{\text{BioFuel}}}}\end{align} \tag{ 2 }$$

here we examine the potential GHG mitigation benefit of fertilizing perennial bioenergy crops at the N_F rate that minimizes GHG_Net, at the US Midwest regional scale. Our analysis focuses on the production of upland ecotypes of switchgrass (Panicum virgatum), the leading candidate crop for deployment in the region (Mclaughlin and Kszos 2005), across 7.1 million ha of non-irrigated marginal lands identified in the Midwest. Understanding the magnitude of the potential benefits at each candidate production site, as well as their context-specific drivers, can aid in the design of targeted agronomic recommendations and policy incentives for the effective use of N_F in bioenergy feedstock production.

2.1. Predictive model development

Given that data on upland switchgrass productivity and ΔSOC is available only for a few experimental locations in the Midwest, we project its N_F response to heterogeneous soils, climate and management using a mixed-modeling approach that combines the strengths of process-based crop simulation models (CSMs) and machine learning (ML; figure 1(a)). To do this, we trained ML algorithms to replicate the ensemble mean of four CSMs using weather, soil, and management data summaries readily computed from available public databases. We selected the gradient boosting machine (GBM) among several algorithms examined because of its superior performance in prediction with this simulated dataset (supplementary figure 1 available online at stacks.iop.org/ERL/17/044004/mmedia). The CSM ensemble consisted of four process-based models (APSIM (Ojeda et al 2017)), DAYCENT (Field et al 2016), EPIC (Jones et al 2018) and SALUS (Martinez‐Feria and Basso 2020), all of which have been previously used to simulate upland ecotypes of switchgrass in the Midwest, and represent diverse model structures. We tested each of these CSMs and their ensemble against a dataset with plant and SOC measurements that encompass three decades of experimentation in the region and a range of N_F application rates (figures 1(b), (c) and supplementary table 1). Further details about the development of the predictive models for switchgrass yield and ΔSOC are included in the supplemental information.

Zoom In Zoom Out Reset image size

2.2.  In-silico nitrogen fertilizer rate experiments on marginal lands

2.3. Estimating net emissions associated with nitrogen fertilizer additions

The impact of fertilizer additions on the life-cycle emissions of biofuel production was calculated as the net difference between the total GWI under a given N_F rate minus the GWI under unfertilized management (N₀) as shown in equation (1). The GWI of switchgrass production can be accounted as the sum of the GHG emitted by fossil fuel used for field operations (GHG_Fuel),seed production (GHG_Seed), agrochemicals other than N_F (GHG_AgCh) and N_F synthesis (GHG_FertSyn), N_F application operations (GHG_FertApp), the biogenic direct and indirect N₂O emissions (GHG_N2O) and changes in soil C sequestration (GHG_ΔSOC) associated with fertilizer management, and downstream emissions related to conversion of the biomass feedstock into biofuel (GHG_Feedstock):

$$\begin{align}{\text{GWI}}\left( {{{\text{N}}_{\text{F}}}} \right)& = {\text{ GH}}{{\text{G}}_{{\text{Fuel}}}} + {\text{GH}}{{\text{G}}_{{\text{Seed}}}}\nonumber\\[1pt] & \quad + {\text{ GH}}{{\text{G}}_{{\text{AgCh}}}} + {\text{ GH}}{{\text{G}}_{{\text{FertSyn}}}}\left( {{{\text{N}}_{\text{F}}}} \right)\nonumber\\ & \quad + {\text{ GH}}{{\text{G}}_{{\text{FertApp}}}}\left( {{{\text{N}}_{\text{F}}} > 0} \right)\nonumber\\ & \quad + {\text{ GH}}{{\text{G}}_{{{\text{N}}_{\text{2}}}{\text{O}}}}\left( {{{\text{N}}_{\text{F}}}} \right)\nonumber\\[1.2pt] & \quad + {\text{GH}}{{\text{G}}_{\Delta {\text{SOC}}}}\left( {{{\text{N}}_{\text{F}}}} \right) + {\text{ GH}}{{\text{G}}_{{\text{BF}}}}\left( {{{\text{N}}_{\text{F}}}} \right){\text{.}}\end{align} \tag{ 3 }$$

Because GHG_Fuel, GHG_Seed and GHG_AgCh are not a function of N_F, these can be ignored in the calculation of GHG_Net. Additionally, GHG_FertSyn and GHG_FertApp equal zero under N₀. Therefore:

$$\begin{align}{\text{GH}}{{\text{G}}_{{\text{Net}}}}& = \left[ {\text{GH}}{{\text{G}}_{{\text{FertAp}}{{\text{p}}_{{{\text{N}}_{\text{F}}}}}}} + {\text{ GH}}{{\text{G}}_{{\text{FertSy}}{{\text{n}}_{{{\text{N}}_{\text{F}}}}}}} + {\text{GH}}{{\text{G}}_{{{\text{N}}_{\text{2}}}{{\text{O}}_{{{\text{N}}_{\text{F}}}}}}} + {\text{GH}}{{\text{G}}_{\Delta {\text{SOC}}}}_{{{\text{N}}_{\text{F}}}} + {\text{ GH}}{{\text{G}}_{{\text{BioFuel}}}}_{{{\text{N}}_{\text{F}}}}{ } \right]\nonumber\\ & \quad - \left[ {{\text{GH}}{{\text{G}}_{{{\text{N}}_{\text{2}}}{{\text{O}}_{{{\text{N}}_{\text{0}}}}}}} + {\text{GH}}{{\text{G}}_{\Delta {\text{SOC}}}}_{{{\text{N}}_{\text{0}}}} + {\text{GH}}{{\text{G}}_{{\text{BioFuel}}}}_{{{\text{N}}_{\text{0}}}}{ }} \right]\end{align} \tag{ 4a }$$

$$\begin{align}\phantom{0000000} & = {\text{ GH}}{{\text{G}}_{{\text{FertAp}}{{\text{p}}_{{{\text{N}}_{\text{R}}}}}}} + {\text{ GH}}{{\text{G}}_{{\text{FertSy}}{{\text{n}}_{{{\text{N}}_{\text{R}}}}}}} + \left[ {\text{GH}}{{\text{G}}_{{{\text{N}}_{\text{2}}}{{\text{O}}_{{{\text{N}}_{\text{R}}}}}}} - {\text{GH}}{{\text{G}}_{{{\text{N}}_{\text{2}}}{{\text{O}}_{{{\text{N}}_{\text{0}}}}}}} \right] + \left[ {{\text{GH}}{{\text{G}}_{\Delta {\text{SO}}{{\text{C}}_{{{\text{N}}_{\text{R}}}}}}} - {\text{GH}}{{\text{G}}_{\Delta {\text{SO}}{{\text{C}}_{{{\text{N}}_{\text{0}}}}}}}} \right]\nonumber\\ & \quad + \left[ {{\text{GH}}{{\text{G}}_{{\text{BioFuel}}}}_{{{\text{N}}_{\text{R}}}} - {\text{GH}}{{\text{G}}_{{\text{BioFuel}}}}_{{{\text{N}}_{\text{0}}}}{ }} \right].\end{align} \tag{ 4b }$$

GHG_FertSyn was assumed to be 4.5 kg CO₂ per kg of N, and GHG_FertApp as 26 kg CO₂e ha⁻¹. Direct and indirect N₂O emissions were estimated following the tier 1 method in the IPCC Guidelines for National Greenhouse Gas Inventories (IPCC 2019), which include N₂O from nitrogen fertilizer, dead biomass residues and SOC loss (supplemental table 3). GHG_ΔSOC were the predicted yearly changes in SOC times 44/12 to convert to CO₂e. Finally, for a given N_R, GHG_BioFuel is calculated as the product of the predicted agronomic biomass yield (Yield_Ag; Mg dm ha⁻¹) and the net GHG offset (GHG_offset; Mg CO₂e Mg⁻¹ dm), that is, the emissions avoided from gasoline combustion:

$$\begin{align}{\text{GH}}{{\text{G}}_{{\text{BioFuel}}}}\left( {{{\text{N}}_{\text{F}}}} \right) = {\text{Yiel}}{{\text{d}}_{{\text{Ag}}}}\left( {{{\text{N}}_{\text{F}}}} \right) \times {\text{GH}}{{\text{G}}_{{\text{offset}}}}.\end{align} \tag{ 5 }$$

The offset was estimated as −0.737 Mg CO₂e Mg⁻¹ of biomass. This value assumes a 1:1 replacement of gasoline by bioethanol on an energy basis and accounts for the supply-chain emissions associated with feedstock transportation (GHG_Transp) and storage (GHG_Storage), transportation and storage, conversion in biorefinery (GHG_Bioref), distribution of biofuel (GHG_Dist), and the CH₄ and N₂O emissions produced by biofuel combustion (GHG_Combust):

$$\begin{align}{\text{GH}}{{\text{G}}_{{\text{offset}}}}& = { }\frac{{{\text{GH}}{{\text{G}}_{{\text{Transp}}}}}}{{1 - {L_{\text{T}}}}}\nonumber\\ & \quad + \frac{{{\text{GH}}{{\text{G}}_{{\text{Storage}}}} + {\text{ Yiel}}{{\text{d}}_{{\text{BioFuel}}}} \times {\text{NC}}{{\text{V}}_{{\text{BF}}}} \times \left( {{\text{GH}}{{\text{G}}_{{\text{BioRef}}}} + {\text{ GH}}{{\text{G}}_{{\text{Dist}}}} + {\text{GH}}{{\text{G}}_{{\text{Combust}}}} - {\text{GW}}{{\text{I}}_{{\text{Gasoline}}}}} \right){ }}}{{\left( {1 - {L_{\text{T}}}} \right) \times { }\left( {1 - {L_{\text{S}}}} \right)}}\end{align} \tag{ 6 }$$

where L_T and L_S (unitless) are the portion of biomass lost during transportation and storage, respectively; Yield_BioFuel is the yield from biomass-to-ethanol conversion, NCV_BF is the net caloric value of bioethanol and GWI_Gasoline is the life-cycle GHG emissions intensity of gasoline. Values for each of the parameters in equation (6) were populated from values reported in the GREET (Greenhouse gases, Regulated Emissions, and Energy use in Transportation) model database (Argonne National Laboratory 2020) and studies shown on supplementary table 3.

3.1. Projected switchgrass response to nitrogen fertilizer additions on marginal lands

Consistent with previous CSM intercomparison studies (Rosenzweig et al 2013, Bassu et al 2014, Martre et al 2015, Basso et al 2018), the ensemble means of the CSMs proved superior to the prediction of any single model, reducing root mean square error of the prediction by 16% and 12% for yield and ΔSOC, respectively (supplementary figure 4). Loss of prediction accuracy by the GBM emulator was negligible (figure 1(f)), as demonstrated by the performance against the measured data which was similar to that of the CSM ensemble (supplementary figure 5).

Projected across all the marginal lands identified in the Midwest and 20 years of weather (figure 2), site-mean agronomic yields without N_F are expected to range between 2.1 and 11.9 Mg ha⁻¹, averaging at 5.6 Mg ha⁻¹ (figure 2(a)). Our predictions, however, show substantial spatial variability due to differences in soil characteristics and climates (figure 2(b)). Highest yields of unfertilized switchgrass are predicted mainly in the upper Midwest, eastern areas of the Northern plains, and scattered pockets of Indiana and Ohio. Unfertilized yields are projected to be relatively low in the drier parts of the central and northern plains (west of the 99th meridian) and western Illinois. Further, most marginal land cultivated with upland switchgrass (∼95%) is expected to accumulate SOC under unfertilized conditions, with ΔSOC projected across all marginal land averaging 0.58 Mg C ha⁻¹ yr⁻¹ (figure 2(d)). Negative ΔSOC under unfertilized conditions is predicted in some areas in North Dakota and South Dakota, and isolated areas in Minnesota, Wisconsin, Ohio and Michigan (figure 2(e)).

Zoom In Zoom Out Reset image size

Analysis of emulator predictions reveals that N_F rate is among the most important determinants of agronomic yields, ranking 3rd most important among all 40 variables included, and accounting for 17% of the total variance (supplementary figure 6). With annual applied N_F (up to 200 kg N ha⁻¹) upland switchgrass is modeled to respond positively in most marginal sites following diminishing-returns pattern. With full N_F, the projected long-term mean yield boost averages 3.1 Mg ha⁻¹ (or 54%), although yield increases are generally predicted to be small beyond 100 Kg N ha⁻¹ rate (12% on average, figure 2(c)). This result is supported by findings from literature syntheses that suggests that 100–180 kg N ha⁻¹ are generally needed to maximize yields (Parrish and Fike 2005, Chen et al 2019, Monti et al 2019).

Conversely, ΔSOC is projected to have a decreasing trend with increasing N_F rate. This is consistent with the theory that N_F in switchgrass has an overall dampening effect on the sequestration potential of in soil C pools via enhanced SOC decomposition (Valdez et al 2017). Yet, the mean effect is minimal, decreasing from 0.58 Mg C ha⁻¹ yr⁻¹ with no N_F to 0.46 Mg C ha⁻¹ yr⁻¹ with full N_F (figure 2(f)). In fact, N_F rate seems to have much less influence in the emulator for predicting ΔSOC, ranking 9th in variable importance and accounting for less than 2% of the variance (supplementary figure 6). Nevertheless, both predicted responses to N_F exhibit substantial spatial variation, explained in part by fundamental differences in soil characteristics and climates. For example, the yield response to N_F is expected to be much more marked in sites with high precipitation, coarser texture and low SOC, whereas the negative ΔSOC response is more pronounced with increasing stand age, in sites with low precipitation and in soils with mesic levels of SOC (∼1%–2% C; supplementary figure 7).

3.2. Optimized nitrogen fertilizer rate based on minimizing net GHG emissions

Next, we contextualize these projected site-specific N_F responses in terms of how they could impact the GHG emissions of the biofuel production system as defined in equation (2). In this calculation, we include the projected changes in yields and SOC sequestration, and the biogenic N₂O emissions due to N_F application using the IPCC tier 1 methodology (IPCC 2019). In addition, we apply a life cycle assessment (LCA) methodology to account for the upstream emissions due to fertilizer synthesis, transportation and application, and changes in the downstream emissions related to the manufacturing of the biofuel, its combustion, and the value replacement of gasoline fuel (methods).

A graphical example of the calculated net GHG balance as a function of N_F rates is shown on figure 3. Here we display the mean GHG balances across all marginal land sites. In figure 3(a) we show both GHG_FertSyn and GHG_N2O increasing monotonically with NF applied, whereas GHG_BioFuel (i.e. the emissions offset credit due to the additional biofuel produced) becomes more negative mirroring the yield response (figure 2(c)). Furthermore, GHG_ΔSOC is counted as net positive emissions, given the small dampening effect of N_F on the SOC sequestration rate predicted by the model (figure 2(f)). Including all GHG sources and credits indicates net GHG emissions mitigation (i.e. negative GHG_Net) for a wide range of N_F levels up to 140 kg N ha⁻¹.

Zoom In Zoom Out Reset image size

Focusing solely on the average response, however, would masks important variation. Indeed, further examination of a gradient of SOC content levels reveals that yield increases in sites with greater than 2% SOC in the top 25 cm are, on average, predicted to be insufficient to overcome the additional life-cycle emissions incurred by applying fertilizer (figure 3(b)). Thus, from a GHG mitigation perspective, these C-rich soils should not generally receive N_F. Yet, these sites represent a small fraction of marginal lands in our dataset (<5%). Optimal N rates are projected to increase monotonically with decreasing topsoil SOC content, up to ∼120 kg N ha⁻¹ in soils with very low topsoil SOC content (figure 3(b)).

Across all marginal lands, the modeled optimum N_F rate averages at 74.8 kg N ha⁻¹, although its distribution is quite wide and spatially variable as shown in figure 4(a). Areas with low optimal rates (<50 kg N ha⁻¹) are mainly projected in western states, and large areas in Minnesota, southern Iowa and northern Missouri, and isolated areas in the eastern states. These low optimal rates occur mostly in sites with low precipitation, SOC-rich soils, or a combination of both as is the case of areas of North and South Dakota. Indeed, these conditions in many marginal lands in the Great Plains, and to a lesser extend areas of Missouri and northern Wisconsin, lead to environments where the optimum N_F rate is much less likely to be greater than zero (figure 4(b)). On the other hand, high optimal N_F rate (>100 kg N ha⁻¹) is projected mainly for the eastern states, particularly in the southern tip of Indiana, Illinois and Missouri (figure 4(a)).

Zoom In Zoom Out Reset image size

With optimal N_F rate, agronomic yields across all marginal lands are expected to increase by an average of 1.99 Mg ha⁻¹ (figure 4(c)), which aggregated at the regional scale would amount to 14.1 million tons of additional lignocellulosic feedstock. Accounting for transportation and storage biomass losses and assuming a conversion yield of 355 l per Mg of biomass at the biorefinery (Argonne National Laboratory 2020) (supplementary table 3), the boost in productivity would amount to nearly 4.6 million m³ of additional bioethanol. This represents about 12% of the total statutory cellulosic biofuel production goal of 39 million m³ for 2020 in the US (National Research Council 2012). The optimal N_F is also projected to have a more mixed effect on SOC sequestration, with only a few marginal lands seeing increase SOC sequestration, particularly in the Great Lakes states (Michigan, Minnesota, Wisconsin), areas of Iowa, and the southern tip of Indiana, Illinois and Missouri. Overall, however, the effect for the entire region is projected to be neutral, with a mean change in soil C sequestration rate of 0.006 Mg C ha⁻¹ yr⁻¹.

Given these responses, optimum N_F additions would mitigate emissions by 0.66 Mg CO₂e ha⁻¹ yr⁻¹ on average across all marginal lands (figure 5(a)). Areas surrounding the Great Lakes are projected to have the greatest N_F-related GHG mitigation potential (up to 2.8 Mg CO₂e ha⁻¹ yr⁻¹), likely a result of concomitant gains in SOC sequestration and a large yield response to N_F (figure 4). On the aggregate, optimal NF additions could provide up to 4.63 Tg CO₂e yr⁻¹ of GHG mitigation compared to unfertilized switchgrass production across all marginal lands examined.

Zoom In Zoom Out Reset image size

Overall, these results reveal the significant contributions that optimal N_F could provide to the GHG mitigation of future lignocellulosic biofuel production in the Midwest. However, these projections should be considered as a hypothetical ceiling given that on-the-ground implementation of site-specific optimal N_F management could face important challenges due to the uncertainty in production costs (Fike et al 2017). Therefore, we further examine how optimum N_F would compare to a more simplistic blanket recommendation of 75 kg N ha⁻¹. This value not only represents the average optimum N_F rate across all marginal lands, but also approximates the long-term average N removal of the harvested biomass (Monti et al 2019).

The blanket-rate strategy is projected to result in a net GHG mitigation of 0.35 CO₂e ha⁻¹ yr⁻¹ on average, or 2.48 Tg CO₂e yr⁻¹ on the aggregate, when compared with the unfertilized scenario (figures 5(a) and (b)). This GHG mitigation is about 54% of that projected to be achievable with the optimal rate. Both of these scenarios compare favorably to a situation where N_F is applied at the rate that maximizes long-term yields (∼160–180 kg N ha⁻¹; supplementary figure 8 and Monti et al 2019), resulting in a slight increase in the GHG emissions over the unfertilized scenario (0.036 CO₂e ha⁻¹ yr⁻¹ on average or 0.26 Tg CO₂e yr⁻¹ aggregated; figures 5(a) and (b)). It should be noted that this latter scenario exhibits the strongest spatial divergence, with areas of large positive net GHG in the Great Plains states and southern Iowa, and large negative net GHG in the Great Lakes states and areas of Illinois and Indiana (figure 5(b)).

Our findings seem to challenge the prevalent belief that N_F application is generally detrimental to bioenergy GHG balances and that its use should be avoided in the production of lignocellulosic feedstock (Erisman et al 2009, Smith et al 2010). On the contrary, here projections suggest that the long-term yield advantages of relatively low amounts of N_F in marginal lands would surpass the GHG emissions stemming from its use. Even a blanket recommendation of 75 kg N ha⁻¹ that would forego nearly half of the potential GHG mitigation benefit, is predicted to be a net positive on the aggregate for Midwestern marginal lands (figure 5(c)). The empirical evidence, though limited, seems to concur. Indeed, a field study in Michigan (Ruan et al 2016) that examined field-level GHG measurements found that N_F additions could support GHG mitigation of switchgrass production up to 140 kg N ha⁻¹, although GHG_Net was minimized with a N_F application of 56 kg N ha⁻¹. Importantly, our projections seem to imply that the benefit could be much greater in marginal lands that are limited by poor N fertility as is the case of much of Michigan, Wisconsin, and isolated areas of the eastern Corn Belt (figure 4). Meanwhile, areas with SOC-rich soils (Iowa, Minnesota) and/or limited by low precipitation (Nebraska, Kansas) or short growing seasons (North and South Dakota) are predicted to be less likely to respond to N_F. Thus, the potential GHG mitigation benefit in these marginal lands is minimal. Lastly, it should be noted that the objectives of maximizing yield and minimizing GHG_Net do not seem to overlap in all but a few environments (figure 5(c)). As such, widespread implementation of fertilizing strategies aiming to maximize yield will likely result in the overall reduction of the climate benefit of cellulosic biofuels (figures 5(b) and (c)). (We include a map classifying Midwestern marginal lands in four levels of optimal N rates in the supplementary figure 9 to facilitate discussions with stakeholders).

Of course, it is likely that the optimal N_F rates estimated here, which are based on minimizing GHG emissions, are different from those N_F rates that are most economically efficient (i.e. where N_F marginal revenue equals marginal cost). Indeed, a recent study (Fike et al 2017) that examined production sites in five states concluded that typical N_F recommendations (56–112 kg N ha⁻¹) added between $37 and $74 ha⁻¹ to production costs and required switchgrass biomass prices in excess of $70 Mg⁻¹ to be economically justifiable under the observed yield responses. Whether the bioenergy industry or ancillary carbon markets can support such pricing is yet to be determined and lies beyond the scope of our analysis.

We advocate caution when interpreting these results given the various sources of uncertainty contained in our projections. Foremost are the uncertainties associated with our characterization of the responses to N_F additions across the pedo-climatic conditions of Midwestern marginal lands. Although the process-based CSMs here used are reasonably capable of capturing the variation in biomass production and SOC changes and their dynamic interaction with agricultural management (supplementary figure 4 (Surendran Nair et al 2012, Field et al 2016, Jiang et al 2017, Liu and Basso 2017, Ojeda et al 2017, Jones et al 2018, Jarecki et al 2020, Martinez‐Feria and Basso 2020)), their complex structure, large sets of parameters, and statistical misspecification (Wallach 2011) generally means that prediction uncertainty is large when CSMs are applied outside the calibrated domain. These uncertainties are partially addressed by using the ensemble mean of independent CSMs, which has been shown to provide more accurate predictions outside of sample populations (Rosenzweig et al 2013, Bassu et al 2014, Martre et al 2015, Basso et al 2018). The emulation step unlikely added significant uncertainty at this scale because of the independence between CSM and ML errors, which would tend to cancel out when aggregating across soil types and geographical locations.

Nevertheless, there are a few known weaknesses shared by all CSMs that could bias our projections. Chief among them is the representation of associative N fixation (ANF) by free-living diazotrophs, which is increasingly recognized as a potentially important source of N to switchgrass (40–60 kg ha⁻¹ yr⁻¹), especially in N-limiting environments (Roley et al 2018, Bahulikar et al 2021). Yet, in contrast to biological N fixation (BNF) in legumes, ANF is largely episodic and irregularly detectable (Roley et al 2019), and our current understanding of the environmental factors controlling ANF rates is rather incomplete (Bahulikar et al 2021). Thus, current switchgrass CSMs, either do not account for this pathway (EPIC) or use instead BNF routines as a surrogate (APSIM, SALUS). Only the DAYCENT model has an explicit routine that, in the absence of soil phosphorus data, calculates free-living N fixation rates as a linear function of annual precipitation (Necpálová et al 2015). We are not aware of any work verifying the robustness of these routines for capturing ANF environmental variation in switchgrass, or whether it is even possible at this stage given the paucity of data available. The likely effect of underestimating non-N_F sources under certain environments would be the overestimation of the yield response, which would lead to an overestimation of the optimal N rate. Although we found no widespread evidence of overestimation of the site-specific yield response in the experiments (supplementary figure 10), a conservative interpretation of our data may still be warranted until further examination of this mechanism within the models can be performed.

An additional source of potential bias in our optimal N_F rate estimates is the use of the offline calculation of direct and indirect soil N₂O emissions following the IPCC tier 1 methodology. We opted not to use simulated N₂O because not all the CSMs we selected contain an explicit nitrification-denitrification routine, or have not been comprehensively validated in production of perennial bioenergy grasses. On the other hand, the uncertainties of the IPCC methodologies are more thoroughly understood. Namely, that the response of annual N₂O fluxes to N_F applied tends to be exponential rather than linear (Philibert et al 2012). Thus, the use of a constant emission factor as required by IPCC tier 1, tends to overestimate N₂O fluxes at lower N_F rates, and underestimate them at higher rates. This pattern seems to hold in switchgrass production (Ruan et al 2016). Using a non-constant emissions factor that increases with N_F rates to calculate N₂O fluxes (as shown by Ruan et al 2016) results in overall lower optimum N_F rates with a much narrower distribution (supplementary figure 11). Nevertheless, it is not known how well this exponential relationship, which is derived from a few experiments, would represent emissions at the regional scale. Thus, more work is needed to improve the characterization of soil N₂O emissions when estimating GHG balances of bioenergy.

A final, possibly large source of uncertainty are the assumptions made for calculating the biofuel GHG offset, which essentially scales the projected yield and ΔSOC responses into GHG emissions (supplementary table 3). Here, we used a purely static LCA approach, which assumes constant values derived from the GREET model (methods) notwithstanding the spatial distribution of future infrastructure (e.g. depots, biorefineries, etc) or biomass end use (e.g. pelletization, direct combustion, energy cogeneration). These would result in changes in the efficiency of utilization of the feedstock supply (O'Neill and Maravelias 2021), particularly because lower switchgrass yields would require larger collection radius, and therefore greater GHG emissions of switchgrass transport (Gelfand et al 2013). Yet, the uncertainty associated with these parameters may be substantial. A sensitivity analysis where each of these values are randomly varied within ±20% of the original value reveals a potentially large departure from the mean GHG response (supplementary figure 12), though the computed optima tend to cluster around the mean found in our main analysis (75 kg N ha⁻¹) suggesting low bias. Yet, two parameters the yield from biomass-to-ethanol conversion (Yield_BioFuel) and the life-cycle GHG emissions intensity of gasoline (GWI_Gasoline), resulted to have the greatest impact in the calculation of the biofuel offset (supplementary figure 10). Thus, much of the uncertainty can be reduced by finding best values for those parameters.

Despite the uncertainties, our study represents an important step forward, as it is the first to examine prospective responses of switchgrass at such scale and resolution, and the emulator model developed could provide the groundwork for building dynamic tools to guide future decisions in agronomic management, technoeconomic analysis, and policy formulation (Ruan et al 2016). Only through this kind of high-resolution predictive tools, biofuel supply-chain network design can move toward integrating temporally and spatially explicit data to reflect the reality of heterogeneity in biomass productivity and SOC sequestration capacity across the landscape, and to allow consideration of land allocation and crop management decisions in addition to infrastructure and operational decisions (O'Neill and Maravelias 2021). This is becoming more critical given the interest in targeted integration of perennial bioenergy crops within the margins of row-crop fields or other marginal subfield areas (Brandes et al 2018, Basso and Antle 2020), which could not only avoid GHG through the production of biofuels, but also by replacing low yielding, highly polluting croplands (Basso et al 2019). Therefore, further research is needed that comprehensively evaluates the modeling uncertainties here outlined to advance prospects for efficient site-specific N_F management in perennial bioenergy crops, and to maximize the GHG mitigation potential of this technology.

This research was supported in part by the Great Lakes Bioenergy Research Center, US Department of Energy, Office of Science, Office of Biological and Environmental Research (Awards DESC0018409 and DE-FC02-07ER64494), the USDA National Institute of Food and Agriculture (Awards: 2015-68007-23133, 2018-67003-27406, 2019-67012-29595) and Michigan State University AgBioResearch. Authors would like to thank Max de Antoni Migliorati (DAYCENT model support), Sasmita Sahoo (EPIC model support), the APSIM initiative, and Lydia Price (spatial data support).

The data that support the findings are within this study in the cited references. R code and the database needed to run the developed GBM emulators has been made public at the following repository: https://doi.org/10.5281/zenodo.5116736.

The data that support the findings of this study are available upon reasonable request from the authors.

The authors declare no competing interests.

R M F and B B conceived and designed the research; R M F assembled the datasets and developed predictive models. S K defined LCA values. R M F analyzed the predicted datasets and visualization of the results. All authors contributed to interpretating the results and writing the manuscript.