Ecosystem model parameterization and adaptation for sustainable cellulosic biofuel landscape design

Renewable fuel standards in the US and elsewhere mandate the production of large quantities of cellulosic biofuels with low greenhouse gas (GHG) footprints, a requirement which will likely entail extensive cultivation of dedicated bioenergy feedstock crops on marginal agricultural lands. Performance data for such systems is sparse, and non‐linear interactions between the feedstock species, agronomic management intensity, and underlying soil and land characteristics complicate the development of sustainable landscape design strategies for low‐impact commercial‐scale feedstock production. Process‐based ecosystem models are valuable for extrapolating field trial results and making predictions of productivity and associated environmental impacts that integrate the effects of spatially variable environmental factors across diverse production landscapes. However, there are few examples of ecosystem model parameterization against field trials on both prime and marginal lands or of conducting landscape‐scale analyses at sufficient resolution to capture interactions between soil type, land use, and management intensity. In this work we used a data‐diverse, multi‐criteria approach to parameterize and validate the DayCent biogeochemistry model for upland and lowland switchgrass using data on yields, soil carbon changes, and soil nitrous oxide emissions from US field trials spanning a range of climates, soil types, and management conditions. We then conducted a high‐resolution case study analysis of a real‐world cellulosic biofuel landscape in Kansas in order to estimate feedstock production potential and associated direct biogenic GHG emissions footprint. Our results suggest that switchgrass yields and emissions balance can vary greatly across a landscape large enough to supply a biorefinery in response to variations in soil type and land‐use history, but that within a given land base both of these performance factors can be widely modulated by changing management intensity. This in turn implies a large sustainable cellulosic biofuel landscape design space within which a system can be optimized to meet economic or environmental objectives.


Introduction
The Energy Independence and Security Act of 2007 expanded the US renewable fuel standard to require the use of large quantities of 'advanced' and 'cellulosic' biofuels with lifecycle greenhouse gas (GHG) emissions reductions of 50% and 60%, respectively, relative to a conventional gasoline baseline (110th Congress of the United States, 2007). This mandate is predicated on the wide availability of biomass feedstocks with low direct environmental impacts and causing minimal disruption to agricultural commodity markets, which could lead to indirect leakage effects (Searchinger et al., 2008). A variety of low-impact feedstock provisioning strategies have been envisioned including the collection of agricultural residues, forestry residues, and municipal wastes, as well as the cultivation of dedicated woody and herbaceous crops on marginal or non-agricultural lands (Campbell et al., 2008;Robertson et al., 2008;Tilman et al., 2009). In particular, perennial grasses such as yield potential, relatively low input requirements, high nitrogen use efficiency, and ability to sequester carbon through soil organic matter formation (Vogel et al., 2002;Tilman et al., 2006;Heaton et al., 2008;Walter et al., 2014). A comprehensive review of US feedstock availability by the Department of Energy suggests that widespread cultivation of high-yielding dedicated perennial grasses will be necessary to achieve the target of displacing 30% of U.S. petroleum consumption with biofuels (U.S. Department of Energy, 2011), and it has been estimated that biofuel supply chains based on such feedstocks will have highly favorable lifecycle GHG impacts as compared to first-generation biofuel technologies Davis et al., 2012;Wang et al., 2012).
Despite their promise, agronomic experience with perennial grass feedstock crops is still relatively limited, and questions around the best management practices to balance the often-competing goals of maintaining high yields while minimizing environmental impacts have not yet been resolved. The debate between the relative merits of low-intensity cultivation over large areas (often referred to as 'land-sharing') vs. more intensive production on a more limited land base ('land-sparing') is still being waged (Anderson-Teixeira et al., 2012). Process-based ecosystem modeling can play an important role in extrapolating limited existing perennial grass field trial results to make more general estimates of productivity, environmental impacts, and optimal management and landscape design strategies.

Modeling cellulosic feedstock yields & environmental impacts
The use of crop models for assessing managementenvironment interactions and predicting bioenergy feedstock crop productivity was thoroughly reviewed by Nair et al. (2012). Crop models such as APSIM, BioCro, and ALMANAC have been applied at regional or national scales to assess the productivity of first-and second-generation feedstock crops as affected by broadscale climate-soil associations (Bryan et al., 2010;Miguez et al., 2012;Behrman et al., 2014). Bioenergy system design is not just a question of yield, however, and understanding the biogeochemical cycling of carbon, nitrogen, and water through these systems is essential for quantitatively evaluating their sustainability . Perennial grass feedstock crops are often associated with a high potential for soil carbon sequestration. One recent meta-analysis suggests that switchgrass increases soil organic carbon (SOC) levels at a median rate of ca. 0.7 t C ha À1 yr À1 when cultivated on carbon-depleted agricultural soils, though performance is more neutral on pastureland or areas that were not previously cropped (Qin et al., 2015a). Application of nitrogen fertilizers is typically required to replace losses during harvest and maintain yield levels, but emissions of nitrous oxide (N 2 O, a byproduct of soil microbial activity and a potent GHG) increase non-linearly with increasing N rate (Hoben et al., 2011;Shcherbak et al., 2014) and can threaten the overall lifecycle GHG footprint of any bioenergy system based on feedstocks with inefficient nitrogen cycling (Crutzen et al., 2008). The biogeochemical cycling of C, N, and H 2 O and associated fluxes of biogenic GHGs are tightly linked in all agroecosystems by fundamental mechanisms including plant tissue stoichiometry, photosynthetic pathway (C 3 vs. C 4 vs. CAM), stomatal conductance, and microbial mineralization/immobilization, and are sensitive to local climate, soil type, and land-use history.
Detailed reviews of process-based biogeochemical model use to capture these interactions in the context of bioenergy system sustainability assessment are provided by Thomas et al. (2013) and Robertson et al. (2015). The CENTURY model was among the first to be applied to bioenergy sustainability assessment, and it and its derivative DayCent model have been widely used to evaluate corn grain production, corn stover removal, and the dedicated cultivation of switchgrass and Miscanthus from the level of individual sites to national scales (Sheehan et al., 2003;Kim & Dale, 2005;Chamberlain et al., 2011;Davis et al., 2012;Lee et al., 2012;Duval et al., 2013). The Environmental Policy Integrated Climate model has also been applied extensively to bioenergy feedstocks in the context of economic analyses (Jain et al., 2010;Egbendewe-Mondzozo et al., 2011) and environmental sustainability assessments (Gelfand et al., 2013) at scales from regional (Zhang et al., 2010) to global (Kang et al., 2014).
Experience with using models to assess the productivity and biogeochemical implications of fine-scale variations in environmental factors such as soil type or topography across landscapes or individual farms is more limited. Studies of switchgrass cultivation in the southeastern US using crop production models have come to different conclusions as to whether yield is sensitive (Woli, 2012) or insensitive (Persson et al., 2011) to underlying soil type. At the finest spatial scale, a group at Idaho National Laboratory has used corn yield data from a precision agriculture system in Iowa to drive DayCent and models of soil erosion and identify areas of low nitrogen use efficiency to target for conversion to switchgrass (Abodeely et al., 2013). Many studies conduct large numbers of fine-scale simulations to make regional-scale estimates of feedstock productivity and environmental performance, though they typically do not report on soil-climate-management interactions explicitly but rather emphasize more aggregate descriptions of landscape performance (e.g. Davis et al., 2012;Gelfand et al., 2013).

Challenges in bioenergy ecosystem modeling
The use of ecosystem models to assess different bioenergy landscape design strategies is complicated by challenges around the representation of marginal lands, adequate bioenergy crop parameterization, and selection of the most appropriate spatial resolution and agronomic management practices for simulation. While biogeochemical process models are increasingly used to simulate conversion of marginal agricultural lands to bioenergy feedstock cultivation (Bandaru et al., 2013;Qin et al., 2015b), these scenarios are particularly challenging from a modeling perspective. The definition of 'marginal' land itself is not straightforward or consistent across studies, and the term is applied to both land with low productivity potential and to land with ample productivity but vulnerable to long-term degradation (erosion, loss of soil organic matter, etc.) under conventional cropping systems. In many bioenergy assessment studies the marginal designation has been based on unfavorable biophysical properties as judged using land suitability ratings (Gelfand et al., 2013) or remote sensing techniques (Cai et al., 2011). Another more direct basis for the designation considers past transitions in an out of agricultural production as inferred from land-use datasets (Campbell et al., 2008), remote sensing (Wright & Wimberly, 2013), or sector-level economic modeling . When considering lands designated as marginal based on their low productivity potential, bioenergy feedstock crops are often not immune to the factors that make such lands challenging for conventional crops. Recent perennial grass field trials purposely conducted on marginal sites indicate reduced productivity relative to performance on the prime lands typically encountered at many agricultural field stations (Mooney et al., 2009;Shield et al., 2012;Boyer et al., 2013). From an ecosystem modeling perspective, accurate assessment is only possible to the extent that underlying biophysical limitations on productivity (e.g. unfavorable climates, soil texture extremes, shallow soils, low soil organic matter levels, site drainage problems, slope, vulnerability to erosion, etc.) are represented directly or indirectly in model data inputs and the processes simulated.
It is well understood that process-based ecosystem models require proper parameterization specific to the agroecosystems being simulated in order to achieve good performance, and that such models have limited predictive power when extrapolated far beyond their parameterization scope (Thomas et al., 2013). However, many bioenergy assessment studies are based on models parameterized under prime conditions and then extrapolated to highly marginal sites, or lacking an explicit independent validation of performance (e.g. Gelfand et al., 2013;Kang et al., 2014). In the case of the DayCent model, parameterization typically involves adjusting study site and crop parameters by hand in order to match observed real-world performance for a small number of field trial cases for which extensive data are available . While this approach can often yield a high degree of fidelity across a range of performance criteria for the sites in question (Hudiburg et al., 2015), the very large number of empirically-determined crop and site parameters in the model makes the process vulnerable to over-parameterization. In such cases, model fit to the training dataset is improved via mechanisms that lack broader underlying ecological significance, reducing the generality of the resulting model for other geographic areas, environmental conditions, or management regimes (Necp alov a et al., 2015). It is also possible to introduce bias with the selection of the parameterization cases themselves, if the researcher gravitates to focus on studies that confirm a priori assumptions of how a system 'should' perform.
There are additional challenges around the spatial resolution of landscape modeling and assumptions about crop agronomic management. Management factors including tillage intensity (Adler et al., 2007), fertilizer application rate , and rotation length (Py€ or€ al€ a et al., 2014) can potentially change the lifecycle GHG performance of a bioenergy system from positive to negative, an effect termed the 'management swing potential' . Management recommendations for bioenergy crops are not always well defined. For example, nitrogen fertilizer recommendations for Miscanthus have been widely debated (Arundale et al., 2014a), with important implications for the overall GHG footprint of production (Roth et al., 2015). To the extent that there are interactions between best management practices and sitelevel ecosystem properties (soil texture, land-use history, etc.), assumptions about management used within an assessment study should ideally be implemented at the level of management decision-making, that is, the field-scale. Some landscape modeling studies have started to endogenize management intensity questions by simulating productivity and environmental impacts at different rates of N application (Gelfand et al., 2013) or different levels of tillage (Zhang et al., 2010) for a given simulation run, for example, accounting for realistic variations in best management practices across the study landscape.

Study goals
Our study used the DayCent biogeochemical process model to assess perennial grass productivity and associated biogenic GHG emissions as a function of land quality and management intensity. Implications for bioenergy landscape design were explored through a case study of switchgrass production around a newly constructed commercial-scale cellulosic biorefinery in an area with substantial heterogeneity in soils and land use. This investigation expands on previous work in two main ways: 1. We conducted an extensive model parameterization and validation effort based on a data-diverse, multicriteria approach, using a large parameterization dataset collected from the literature spanning wide gradients of climate, soil texture, and management intensity. 2. We evaluated the impacts of management intensification at the full spatial resolution of the assessment, estimating optimal levels of nitrogen fertilizer application for each point on the landscape in order to either maximize yield or minimize biogenic GHG emissions.
Our objective was to develop a rigorous, well-validated spatially explicit biogeochemical modeling capability that can serve as the basis for future integrated assessment and landscape optimization efforts.

Case study introduction
We performed a landscape assessment case study simulating the cultivation of switchgrass (Panicum virgatum) to supply biomass to the Abengoa cellulosic biorefinery located outside the town of Hugoton in southwestern Kansas (Peplow, 2014), which began operations in fall 2014. While the plant will initially produce 25 million gallons of ethanol per year using corn stover as the primary lignocellulosic feedstock, switchgrass has been mentioned as an advanced cellulosic feedstock of interest and Biomass Crop Assistance Program Project Area 7 is sponsored by the company and targets the establishment of 20 000 acres of switchgrass production in the area (U.S. Department of Agriculture FSA, 2011). The case study region has long been at the center of agricultural sustainability and energy issues, having been deeply affected by the Dust Bowl in the 1930s (Kansas Historical Society, 2014) and being the site of the earliest hydraulic fracturing trials in the U.S. (Borowski, 2012).
Today, the surrounding Stevens County is a highly diverse and productive agricultural area. In 2012, 21.4% of the county was dedicated to highly-productive irrigated corn cultivation (average yield of 12.1 Mg ha À1 , or 192 bushels acre À1 ) and 11% to dryland wheat (1.1 Mg ha À1 , or 18 bu acre À1 ), with smaller fractions devoted to other crops and pasture/rangeland, supporting an inventory of 45 500 head of cattle including calves (U.S. Department of Agriculture NASS, 2014). In this study we investigated the biogeochemical implications of converting non-irrigated cropland and rangeland in Stevens County and its six neighboring counties in southwestern Kansas and the Oklahoma panhandle to switchgrass cultivation (see Fig. 1), examining tradeoffs between productivity and associated biogenic GHG emissions as a function of underlying soil type and management intensity, specifically nitrogen fertilizer application rate. Issues of land ownership, conservation easement status, and other land-use policy limitations are excluded here, but explored in a subsequent publication dedicated to bioenergy landscape optimization.

The DayCent model
Productivity and net fluxes of biogenic CO 2 and N 2 O from soils under switchgrass cultivation were modeled with the DayCent biogeochemistry model (Parton et al., 1998;Del Grosso et al., 2011). DayCent is a semi-empirical process-based model that simulates cycling of C, N, and water in natural and agroecosystems based on site-specific biophysical factors, land-use history, and management practices (e.g. tillage, fertilizer application, irrigation, etc.). The spatial and temporal scope of the model lies in between that of dedicated crop growth models (Miguez et al., 2012) and generalized earth climate system models Hallgren et al., 2013). DayCent has been used extensively to predict yields and environmental impacts of switchgrass cultivation (Adler et al., 2007;Chamberlain et al., 2011;Davis et al., 2012;Lee et al., 2012) and is also used to predict agricultural soil GHG emissions for the annual Inventory of U.S. Greenhouse Gas Emissions and Sinks (U.S. Environmental Protection Agency, 2014a).
DayCent computes soil temperature and moisture for different layers of the soil profile (resolved separately based on soil texture, bulk density, and pH) using daily climate data inputs. Crop growth (net primary productivity, or NPP) is simulated using species-specific parameters describing photosynthetic efficiency, tissue C/N ratio limits, above-and below-ground C partitioning, and phenology, many of which are determined empirically using model parameterization datasets as described previously. Daily biomass growth potential (NPP pot ) is derived from top-of-atmosphere radiation (srad), corrected for atmospheric transmission losses and multiplied by a series of 0-1 factors that represent deviations from ideal temperature (f temp ) or soil moisture (f H2O ) or limitations due to canopy immaturity or self-shading (f canopy ). This potential growth is then adjusted down if the available soil mineral nitrogen supply is limiting, as determined based on a maximum incremental biomass C : N ratio at the given level of plant maturity: In addition to plant productivity DayCent also estimates soil carbon and nitrogen cycling including net changes in SOC levels and N 2 O emissions, the major constituents of the agricultural soil GHG balance. Carbon dynamics are simulated for soil surface pools and the top soil layer (0-20 cm by default), with organic matter represented by two litter pools (metabolic and structural) and three SOC pools-an 'active' pool representing microbial biomass and associated microbial products with a rapid turnover rates, and two others representing chemically/ physically stabilized carbon with decadal-('slow' pool) and century-scale ('passive' pool) turnover times. Nitrogen mineralization and immobilization rates for each pool are controlled by the maximum and minimum permissible C : N ratios for each pool, soil temperature, soil moisture, and microbial efficiency as a function of soil texture (Parton et al., 1987). The soil nitrogen balance considers synthetic N fertilizer addition, manure and other organic N amendments, atmospheric deposition, volatilization, leaching, plant uptake, and N mineralization and immobilization associated with soil organic matter transformations. The model simulates nitrification of ammonium (NH 4 + ) to nitrate (NO 3 À ) including NO X and N 2 O by-products, as well as denitrification of nitrate to gaseous products (N 2 O, NO X , N 2 ). Nitrification is simulated by multiplying available soil ammonium by a maximum potential nitrification rate adjusted based on soil temperature, water-filled pore space (WFPS), and pH limitations (Del Grosso et al., 2000;Parton et al., 2001). The overall rate of denitrification and the N 2 O/N 2 ratio of its products are modeled based on the availability of nitrate and organic matter substrates (as inferred from heterotrophic respiration rate), and local soil micropore redox state and gas diffusion rates as inferred from WFPS and heterotrophic respiration rate.

Ecosystem model parameterization & validation
We undertook an extensive parameterization and validation of the DayCent switchgrass growth model to improve overall simulation accuracy and verify the model response across gradients of land quality and management intensity. A large dataset of switchgrass field trials for the continental United States was assembled from the peer-reviewed literature. Studies were included in the dataset provided they (a) specified the underlying soil in sufficient detail that a corresponding Natural Resource Conservation Service Soil Survey Geographic database (SSURGO; Ernstrom & Lytle, 1993) map unit could be identified using the Web Soil Survey (U.S. Department of Agriculture NRCS); (b) specified key management variables (switchgrass ecotype planted, N fertilizer application rate, harvest date) in sufficient detail to define key DayCent simulation parameters; and (c) reported disaggregated annual yield and GHG results (studies reporting results averaged across multiple sites or N treatments were excluded). All modeled-vs.measured performance comparisons were made on the basis of time-averaged results across the duration of the field trial, so studies that report either annual yield data or treatmentaveraged data were included. If multiple cultivars of the same ecotype were included in a study, results were averaged across cultivars into a single representative value for the ecotype for simplicity. Additional details on the parameterization and validation dataset are provided in the Appendix S1. Fig. 1 Map of all switchgrass field trial sites included in the full model parameterization and validation dataset prior to ecotype/latitude filtering. Ring size indicates the number of experimental treatments (e.g. different ecotypes or nitrogen fertilizer application rates) conducted at that site, and color represents the types of data available. The 7-county case study area is highlighted in pink, with a star designating the biorefinery location.
A total of 24 appropriate studies were identified, covering 573 annual biomass yield points across 147 unique combinations of site (soil, climate) and management treatments (N rate, harvest date, etc.). Initial exploratory analysis confirmed the need to exclude the first two seasons of yield data before switchgrass yields stabilize (Lesur et al., 2013;Arundale et al., 2014b), and to filter out treatments for ecotypes grown outside their typical latitude range (up to 54°N for lowland varieties and down to 34°N for upland, as per Casler, 2012). The remaining dataset was then randomly split at the level of individual studies 70 : 30 into a parameterization dataset and an independent validation dataset (Table 1). Several of these field trials included more detailed information on biomass partitioning or nitrogen content (Frank et al., 2004;Dohleman et al., 2012;Anderson-Teixeira et al., 2013), long-term changes in SOC based on either repeat measurements or paired sites Liebig et al., 2008;Follett et al., 2012;Anderson-Teixeira et al., 2013;Bonin & Lal, 2014), and/or time-resolved N 2 O emissions (Niki ema et al., 2011;Hong et al., 2012;Schmer et al., 2012;Smith et al., 2013). There was insufficient data available to perform an independent validation of these model performance criteria, so all of these studies were included in the parameterization dataset. While the raw SOC dataset was very noisy, eliminating individual data points that were not reported as statistically significant yielded a final reduced modeling dataset that behaved more predictably and was consistent with recent switchgrass SOC meta-analysis results (i.e. Qin et al., 2015a).
The combined parameterization and validation dataset covers a wide range of latitudes and longitudes (Fig. 1) and temperature/precipitation regimes (Fig. S1), as well as a wide range of soil textures and Natural Resource Conservation Service land capability class (LCC) ratings (Helms, 1992; Fig. S2). Note that LCC ratings are reflective of a variety of land-use limitations, some of which are explicitly simulated in the Day-Cent model (e.g. dry climates, extreme textures, shallow soils) and some of which are not (e.g. erosion susceptibility, drainage class). The switchgrass crop parameterization was further informed with data from greenhouse or growth chamber experiments looking at productivity response across gradients of temperature (Balasko & Smith, 1971;Hsu et al., 1985;Reddy et al., 2008;Kandel et al., 2013;Wagle & Kakani, 2014) or moisture (Xu et al., 2006).
The parameterization process started with a default switchgrass crop parameter set based on previous work (Adler et al., 2007;Davis et al., 2012) and focused on refining parameters relating to productivity, temperature and moisture stress response, nitrogen management, shoot vs. root partitioning, and tissue death and turnover rates, with separate parameterizations  (2002) for both upland and lowland ecotypes as appropriate. Initial exploratory analysis suggested that capturing differences in phenology between upland and lowland ecotypes was essential for accurate yield simulation, consistent with the current understanding of maturation based on photoperiod being a strong determinant of yield differences between different cultivars grown at a given latitude (Casler et al., 2004). We set green-up dates uniformly for both ecotypes as function of latitude based on a variety of literature sources (Sanderson, 1992;Sanderson et al., 1997;Wang et al., 2013) as illustrated in Fig. S5. Peak biomass dates were predicted as a function of both latitude and ecotype as inferred from a variety of sources (Sanderson, 1992;Hopkins et al., 1995;Sanderson et al., 1997;Vogel et al., 2002;Frank et al., 2004;Berdahl et al., 2005;Casler et al., 2007;Anderson-Teixeira et al., 2013;see Figs S6 and S7) and used to trigger plant senescence events within DayCent.
After crop phenology was set other parameter adjustments were implemented manually, using the model automation routine described in the next section to rapidly evaluate parameter changes against the 67-point parameterization dataset. The parameterization process focused on maximizing modeled-vs.measured fidelity (based on visual inspection of the plotted data and calculation of modeled-vs-measured regression parameters, Pearson correlation coefficient, and root mean squared error; e.g. Smith et al., 1996) for upland and lowland ecotype yields, changes in SOC, and growing season N 2 O emissions, but also took into account time-resolved shoot : root and C : N ratio data where available for a multi-criteria evaluation of parameter set performance (e.g. Fig. S13). Once the parameterization process was complete, independent validation of upland and lowland yield performance was conducted based on the data held in reserve (i.e. using the holdout validation method).

Spatial data inputs, model initialization & automation
A variety of spatially explicit data inputs are necessary to initialize and run the DayCent model for a large-scale parameterization or landscape analysis, including data on climate, soil type, and land-use history. Data sources used in this analysis are summarized in Table 2. Soil texture, rock fraction, and pH for different soil profile layers of the dominant soil component for each map unit were taken directly from the SSURGO database (Ernstrom & Lytle, 1993), and bulk density, field capacity, wilting point, and saturated hydraulic conductivity were computed using the Saxton equations (Saxton et al., 1986). Climate data on a 32 km grid was derived from the North American Regional Reanalysis database (NARR; Mesinger et al., 2006).
Land-use history and current land management practices were compiled from a variety of sources. Current land use was determined from the National Land Cover Database 2006 (NLCD; Wickham et al., 2013), re-sampled from the native 30 m resolution to 240 m for ease of use, and re-classified into the simplified categories of annual agricultural lands ('cultivated crops', 'pasture/hay'), rangeland ('dwarf shrub', 'shrub/ scrub', 'grassland/herbaceous', 'sedge/herbaceous'), and excluded (all other categories including forested and developed lands). Irrigated areas were identified using the MIrAD-US database (Pervez & Brown, 2010), and federally-owned lands were identified using the USGS Federal Lands of the United States data layer (U.S. Geological Survey, 2015) and excluded from further analysis (3.4% of the landscape, most part of Cimarron National Grassland). These five GIS layers were then intersected and small slivers were eliminated by merging all polygons smaller than 1 ha into the neighbor with which they shared the longest border. This yielded 39 320 polygons of a variety of sizes across the seven-county Hugoton case study area (Fig. S16), representing 3779 unique combinations of model inputs requiring individual simulation, which we refer to as DayCent modeling 'strata'.
For each strata, the DayCent model was pre-initialized using the same pre-settlement and historical land use assumptions as used in the EPA Inventory of U.S. Greenhouse Gas Emissions and Sinks (U.S. Environmental Protection Agency, 2014a) and described in detail by Ogle et al. (2010). Model initialization included an equilibration run of several thousand years duration reflecting the natural state of the land prior to conversion to agriculture in order for all soil C and N pools to achieve steady-state values. Historical management between initial plow-out and the modern period was simulated with crop rotations and management practices compiled at regional scale from a variety of historical and modern sources (Ogle et al., 2010). The future switchgrass simulations were then executed across part of a 29-node, 288-processor cluster computing system at the Colorado State University Natural Resource Ecology Laboratory. Parallel execution was implemented in Python (http://www.python.org/) using forking operations to take advantage of multiple cores within a given node.

Landscape design analysis scenarios, results processing & sensitivity analysis
For the landscape analysis case study we simulated conversion of all non-irrigated, non-federally owned polygons within the seven-county case study area to rain-fed lowland switchgrass cultivation. We conducted 30-year forward simulations to assess long-term productivity and trends in soil C and N cycling, recycling the full range of the NARR historic weather record to represent future weather conditions. In order to assess the response of crop productivity and GHG performance to management intensity, seven different rates of nitrogen fertilizer application were simulated for each strata (0-150 kgN ha À1 in 25 kgN ha À1 increments). We assumed that switchgrass would be replanted every 10 years after field preparation consisting of chisel plow and field cultivator operations, and that the crop would be neither fertilized nor harvested the year of establishment in order to limit competition from weeds and ensure robust crop establishment, as per local extension recommendations. These assumptions are highly conservative as switchgrass is often established in this region without tillage, the need for periodic replanting is widely debated, and first-year harvest can be possible if the crop achieves sufficient first-year productivity. Switchgrass harvest yields, changes in SOC levels, and annual N 2 O emissions were then averaged over the 30-year simulation period for each strata. Average annual N 2 O emission values were converted into CO 2 equivalents using a 100-year global warming potential (GWP 100 ) value of 298 (Forster et al., 2007) then added to the CO 2 flux values associated with average annual net SOC changes for an estimate of total direct biogenic emissions. To determine the biogenic GHG intensity of production (Mg CO 2 eq Mg À1 biomass harvested) total emissions per hectare were divided by the associated simulated switchgrass yield. Continuous functions of biomass yield and biogenic GHG intensity vs. nitrogen fertilizer application rate were developed for each strata by applying a cubic regression to the 30-year averaged simulation results for the different N rate simulations, and the yield-maximizing and GHG balance-minimizing N rates were then interpolated. The simulated yield and GHG intensity associated with these optimized N rates was then reassociated with the appropriate landscape polygons, and aggregated across the full landscape in order to develop curves illustrating total potential landscape productivity and biogenic GHG emissions balance when strata are managed for these different objectives. The sensitivity of these landscape results to key crop parameters, landscape characterization, and switchgrass cultivation scenario assumptions was assessed as detailed in Table 4 in order to determine the overall robustness of our conclusions. All results analysis routines were automated in Python through a combination of SQLite database operations (http://sqlite.org/ ) via the sqlite3 module (http://docs.python.org/2/library/sqli-te3.html), data manipulation in the native Python list data type, and figure generation using the matplotlib.pyplot module (http://matplotlib.org/api/pyplot_api.html).

Ecosystem model parameterization & validation
A total of 79 different switchgrass parameter set iterations were ultimately developed and tested. Final upland and lowland ecotype parameter values that differ from the default DayCent switchgrass crop parameterization are detailed in Table 3. The most significant changes were: • Increased plant potential NPP rate -The most recent version of DayCent explicitly models solar radiation atmospheric transmission losses, and the revised potential NPP parameter PRDX(1) must be adjusted higher relative to previously-published versions to reflect the new growth calculation on a canopy photosynthetically active radiation (PAR) basis rather than a top-of-atmosphere PAR basis. Further fine adjustments were made to PRDX(1) to optimize the observed yield difference between the different ecotypes and to offset slightly increased belowground C partitioning.
• Adjusted temperature and moisture stress response curves -Crop temperature (Fig. S8) and moisture stress (Fig. S9) response curves were set based on the greenhouse and growth chamber studies listed in Table 1. In the case of temperature response, fine adjustments to the edges of the curve where direct empirical data were lacking were implemented to improve overall modeled-vs.-measured yield performance across the full parameterization dataset. The same temperature and moisture response curves were used for both ecotypes. Comparison of measured and modeled yield ranges binned by site average growing degree day accumulation (Fig. S10) or annual precipitation (Fig. S11) verifies that the model accurately captures increasing switchgrass productivity at warmer, wetter sites.
• Increased belowground partitioning and root turnover -The default parameterization slightly underestimated belowground biomass, significantly underestimated observed SOC increases, and over-predicted N 2 O emissions. A small increase in belowground partitioning coupled with a large increase in root turnover rates resulted in more carbon being cycled into the soil and more mineral N being taken up by the plant, improving model performance on all three criteria.
Model parameterization and validation results for yields and soil GHG fluxes are shown in Fig. 2. Sufficient data was available to perform holdout method independent validation of yield predictions for both the upland and lowland ecotypes (Fig. 2a). The out-of-sample validation root mean square errors (RMSE) are 3.7 and 4.1 Mg ha À1 for the upland and lowland ecotypes, respectively, with minimal bias (mean difference of À0.2 and +2.3 Mg ha À1 , respectively). When these data are binned by nitrogen fertilizer application rate (Fig. 2c) we see that the model is able to capture the general trend toward increased switchgrass productivity with increasing management intensity. In contrast, yield response to land quality was more ambiguous. Neither measured nor modeled yields exhibited a clear relationship with soil texture across the full parameterization and validation dataset (Fig. 2d), and the weak trend towards lower yields at sites with higher LCC rating (more marginal) was not replicated in our model (Fig. S12).
Field trial data on soil GHG balance were sparser, and the results presented in Fig. 2b reflect within-sample model performance against the parameterization dataset rather than an independent validation. Observed annual changes in SOC under switchgrass were much larger than measured growing season N 2 O emissions when compared in CO 2 -equivalent terms. The within-sample RMSE value for the combined SOC and N 2 O dataset is 0.69 MgCO 2 eq ha À1 yr À1 across the wide range of sites, climates, and nitrogen application rates represented in the underlying studies (n = 19, Fig. 1). Bias calculations show that the model tends to err in the direction of overestimating N 2 O and underestimating SOC accumulation, so the resulting predictions of switchgrass GHG balance are somewhat conservative. Additional detail on SOC and N 2 O performance is available in Figs S14 and S15.

Landscape design case study
Simulated lowland switchgrass yields as a function of nitrogen fertilizer application rate across the 3779 unique DayCent strata in the Hugoton case study area are shown in Fig. 3 for both crop land and rangeland conversion, color-coded to the underlying soil texture of each strata. The maximum attainable yield under arbitrarily well-fertilized conditions shows significant variation with soil texture, ranging from more than 10 Mg ha À1 in certain clay and sandy soils down to about 6 Mg ha À1 in the more moderate-textured silty soils. In this semi-arid climate of intermittent precipitation events soil moisture levels are often near wilting point, and simulated average yields reflect a tension between the greater total water holding capacity of finer-textured soils (an advantage during relatively wet years) vs. more effective infiltration and less surface soil evaporation in coarser soils (the so-called 'inverse texture effect', beneficial during dry years; Noy-Meir, 1973;Epstein et al., 1997;Lane et al., 1998). Switchgrass yields on converted rangeland are generally less sensitive to fertilizer application rates than those on converted cropland, in some cases showing no response to increasing N rate. This is due to higher background mineral nitrogen levels from transient soil organic matter turnover following conversion in these areas. However, for most strata, full switchgrass yield potential is realized at N rates from 60 to 100 kgN ha À1 , as indicated with circular markers in the upper panels of Fig. 3.

Clayey
Loamy Sandy Silty When examining landscape assessment results under different levels of management intensity it is important to note a discrepancy between the metrics of total area GHG emissions or mitigation (MgCO 2 eq per hectare per year) and the GHG intensity of biomass production (MgCO 2 eq per Mg biomass grown). Biomass production GHG intensity results are shown as a function of nitrogen application rate across all simulation strata in the lower panels of Fig. 3. For most strata at most nitrogen fertilizer rates, 30 year-average soil carbon sequestration outweighs nitrous oxide emissions on a CO 2 -equivalent basis. Biomass has the lowest (most negative) direct biogenic GHG footprint when cultivated with no nitrogen fertilizer on previously cropped fine-textured soils due to their high potential for SOC accumulation. This sequestration value is increasingly offset at higher N rates due to marginal soil N 2 O emissions outpacing corresponding yield and SOC gains. In contrast, rangeland has higher initial SOC levels and thus less capacity for carbon sequestration after conversion to bioenergy feedstock cropping, which results in an overall net GHG intensity much closer to zero.
Total area biogenic GHG emissions are explored in relation to productivity in Fig. 4. Since switchgrass can be managed either to maximize yields or to optimize soil GHG balance (i.e. to maximize soil carbon sequestration while minimizing N 2 O emissions), we selected a random 10% of the landscape polygons, determined the nitrogen application rate that maximizes yield and the rate that maximizes GHG mitigation for each polygon, and aggregated the results for these different management strategies across our landscape sub-sample. The difference between the two curves is representative of the degree to which system productivity and GHG performance can be modulated by adjusting management on the same limited land base. Converting 10% of the landscape to switchgrass managed under either strategy results in >0.75 Mt biomass feedstock production annually, enough material to supply approximately two facilities the size of the Abengoa biorefinery, with the associated net GHG impact of sequestering >0.06 Mt (60 000 metric tons) of CO 2 . Figure 5 illustrates that the optimal switchgrass management intensity for maximum GHG mitigation is related to the distribution of soil types and current land use across the case study area, factors that are also correlated with one another as conventional cropping tends to be concentrated on the moderate textured soils of this landscape. Fig. 3 Simulated yield and GHG intensity response to increasing nitrogen fertilizer rate for switchgrass production on former crop land and former rangeland across the 3779 distinct DayCent strata in the case study landscape. Interpolated yield maxima and GHG intensity minima are marked with solid markers for cropland conversion and open markers for rangeland conversion. The color of the lines and markers indicates the soil surface texture of the strata, with yellow = sand, red = clay, green = silt, and brown = loam.

Sensitivity analysis
Our modeled landscape productivity and GHG balance are highly sensitive to certain key crop model parameters, particularly optimal growth temperature, potential NPP, and response to moisture stress (Table 4), highlighting the importance of careful parameterization and validation. While the first two parameters were reasonably well constrained as a result of the switchgrass parameterization process, empirical data on switchgrass moisture stress in the literature is rare, and uncertainty around those parameters in the model should be considered large. Interestingly, the aggregate landscape simulation results show little sensitivity to the length of the switchgrass growing season in early spring and late fall, even though accurate characterization of the phenology proved essential for determining crop parameter sets with adequate out-of-sample validation performance. Accurate characterization of the landscape in terms of climate and land-use history is important to the integrity of simulation results as well, with GHG mitigation particularly sensitive to past land-use history. Landscape results were generally less sensitive to crop management scenario assumptions around tillage intensity and fertilizer timing.

Challenges of model development and data-diverse parameterization
This study was grounded in an extensive model parameterization and validation effort using a data-diverse multi-criteria approach enabled by automation of all simulation runs and results analyses. The diversity of studies included in our switchgrass field trial dataset was intended to provide a highly general test of model performance independent of a single study, environment, or management practice, and the simultaneous consideration of different types of data (Table 1) was designed to ensure that accuracy of one model performance criteria was never improved at the expense of others. This approach proved useful for sorting out interrelated model responses, for example, the strong effects of belowground partitioning and root turnover rate parameters on harvested biomass yield, SOC changes, and N 2 O rates that must be balanced in a systematic manner during model parameterization.
Our resulting model explains approximately half of the observed variability in yield and GHG observations in our dataset, and realistically captures climate effects and responses to management intensity. While soil type has been observed to have an effect on bioenergy grass productivity in semi-arid climates outside the U.S. (Di Virgilio et al., 2007;o Di Nasso et al., 2015;Roncucci et al., 2015), we did not observe a strong texture signal in our domestic switchgrass productivity dataset, consistent with a similar previous large-scale model parameterization effort (Wullschleger et al., 2010). Future modeling work on soil-climate interactions in the U.S. would greatly benefit from additional field trails like Wilson et al. (2014) that include paired trails across multiple soil types on a landscape, in order to tease apart soil-precipitation interactions without the other larger sources of variability present in composite datasets.
While we believe the approach presented here represents an improvement over more limited model parameterization efforts, it is still possible to overparameterize to a dataset of this size and achieve good fits via unrealistic mechanisms that do not translate well out-of-sample. In our experience, reliable yield performance was only achieved once crop phenology was adequately captured. This parameterization dependency is somewhat challenging, as there are a limited number of studies in the literature reporting detailed phenology in the form required for this type of generalized modeling effort (i.e. dates for green-up and peak biomass or senescence as a function of crop cultivar and site latitude). Additionally, since maturation dates can vary significantly even for cultivars within a single ecotype (Frank et al., 2004), focusing on only two ecotype groupings introduces additional errors. Accurate representation of phenology can be superseded to a certain extent by narrowing crop temperature response curves to truncate early and late-season productivity, but our initial parameterization efforts in this direction performed very poorly on validation.
After several dozen parameter set iterations we were approaching the limit of what could be accomplished with setting parameter values manually, as model Fig. 4 Cumulative total soil GHG mitigation vs. cumulative switchgrass biomass production for a random 10% of the case study landscape under different management goals: managing each land parcel to maximize switchgrass yields (dashed line, highest-productivity sites aggregated first), or managing each land parcel to maximize ecosystem GHG mitigation (solid line, strongest mitigation sites aggregated first). responses became more subtle and antagonistic across the different performance criteria. This, coupled with the relatively large number of crop parameters in the DayCent model, suggests that future efforts could benefit greatly from systematic parameter optimization techniques such as inverse modeling and/or j-fold cross validation. Such techniques provide a more systematic and transparent approach, facilitate maximum extraction of information from a field trial dataset, and help to identify and avoid over-parameterization issues. While inverse modeling techniques have been demonstrated to improve model performance against small datasets (Necp alov a et al., 2015), they are not commonly applied to crop growth and biogeochemistry models used in bioenergy sustainability analysis at the scale of a full parameterization and validation effort such as the one presented in this study. More broadly, model performance is dependent on underlying model structure in addition to parameterization, and DayCent model development efforts to improve its performance in perennial grass systems are ongoing. For example, accurate yield estimates depend on capturing crop phenology and canopy dynamics, and efforts are underway to represent early-season leaf canopy closure and late-season senescence dynamics more realistically. Similarly, while past model development efforts have largely focused on the upper 0-20 cm soil layer most relevant for annual crops, representation of soil carbon changes under perennial grasses may be improved through more explicit simulation of deep rooting, dissolved organic carbon movement, and other aspects of soil organic matter dynamics deeper in the soil profile (Campbell, 2015).

Landscape design case study interpretation and climate accounting issues
The Hugoton case study site was selected because it featured the most heterogeneity in soils and current land use among the first three commercial-scale cellulosic biorefinery sites in the U.S. (the other two being located in the Iowa corn belt; Peplow, 2014). However, the dry climate at this site proved challenging from a modeling perspective, with landscape simulation results showing high sensitivity to interactions between crop moisture stress parameterization and soil texture, and with N 2 O response to increasing N rate lower than what might be expected in a wetter climate. Overall, our landscape simulation results are similar to others in the literature that find the GHG balance of perennial grasses dominated by soil carbon sequestration immediately postestablishment, with only modest N 2 O emissions (e.g. Gelfand et al., 2013). In field trials where baseline N 2 O emission rates are reported and an Intergovernmental Panel on Climate Change (IPCC) Tier 1 estimate can be made, observed switchgrass N 2 O emissions rates are often near or below the lower bound of the IPCC N 2 O emissions factor range (0.3-3% fertilizer N emitted as N 2 O-N; Fig. S15). This suggests that switchgrass is an efficient nitrogen cycler, and that the critique of biofuel (a) (b) (c) Fig. 5 Map illustrating (a) soil texture distribution, (b) current land use, and (c) the switchgrass nitrogen fertilizer application rates associated with maximum area GHG mitigation (as opposed to biomass GHG intensity explored in Fig. 3), for all non-federal non-irrigated cropland and rangeland across the Hugoton case study area. Aerial imagery is included as a background layer for scale, and the boundary of the 7-county case study area is shown in white with the biorefinery location marked with a yellow star.
GHG mitigation benefits being outweighed by direct and indirect N 2 O emissions (Crutzen et al., 2008) is likely overstated for second-generation perennial grass feedstocks. Our predicted SOC sequestration rates after cropland conversion are very similar to median metaanalysis results in Qin et al. (2015a), while the finding of a small net positive sequestration with the conversion of grassland is consistent with, but on the optimistic edge of, their observed range. Full quantification of the confidence interval around these landscape assessment results is highly challenging due to multiple levels of uncertainty in problem specification, landscape input data, model structure, and model parameterization (Walker et al., 2003). While previous work has shown that uncertainty in DayCent-estimated changes in SOC is dominated by model structural uncertainty (Ogle et al., 2010), repeating this empirical uncertainty assessment approach here is impractical as only one of the 145 treatments in our parameterization and validation dataset includes data for all three of the factors that determine biomass GHG intensity (yield, SOC change, and N 2 O; Fig. 1). Here we rely on sensitivity analysis to identify priorities for future model improvement efforts. When discussing model sensitivity and uncertainty issues it is also important to be cognizant of issues in the underlying field trial datasets. Field measurements of N 2 O are very challenging due to emissions variability on extremely fine spatial (Li et al., 2013) and temporal (Jørgensen et al., 2012;van der Weerden et al., 2013) scales, and studies based on sampling at weekly or every-other-week frequency (as was the case for all N 2 O studies in our dataset) are vulnerable to systematic biases of up to 20% and 60%, respectively (Parkin, 2008;Smemo et al., 2011). Additionally, small-scale field trials do not necessarily reflect imperfections in agronomic management (e.g. uneven fertilizer application) and real-world harvest losses, possibly introducing a systematic over-estimation of switchgrass productivity (Searle & Malins, 2014).
Finally, there are several issues around climate impact accounting relevant to this landscape study. While soil carbon sequestration will eventually attenuate even though N 2 O emissions will persist for as long as N fertilizer is being applied (Sheehan et al., 2003;Adler et al., 2007), our assessment for this semi-arid system suggests it will take 60-80 years for annual sequestration and N 2 O rates to reach parity (Fig. S17). A more dynamic climate impact accounting approach that takes transient forcing benefits into account (e.g. Holtsmark, 2015) would tend to further weight near-term SOC benefits against future N 2 O emissions, though no standardized accounting approach is yet widely accepted (U.S. Envi- Tillage intensity Switch to no-till crop establishment +0.13% +0.97% GHG, greenhouse gas; NARR, North American Regional Reanalysis. ronmental Protection Agency, 2014b). The current study is also somewhat limited in focusing on the climate impacts of biomass feedstock production from a purely biogeochemical perspective, ignoring potential biophysical impacts such as changes in surface albedo or evapotranspiration and water dynamics that are significant in some bioenergy production scenarios (Muñoz et al., 2010;Georgescu et al., 2011;Cherubini et al., 2012;Caiazzo et al., 2014), and ignoring any potential broader impacts on other ecosystem services (Chamberlain & Miller, 2012). Additionally, future changes to atmospheric CO 2 concentrations and climate are not considered here, though they may have large repercussions for landscape design (Bryan et al., 2010). The ability of current ecosystem models to accurately extrapolate to such future conditions is an active area of investigation (De Kauwe et al., 2013). Additional sensitivity analysis on these points would improve our understanding of bioenergy landscape performance potential, and highlight priorities for future research efforts.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Appendix S1. Switchgrass parameterization and validation dataset development. Appendix S2. Parameterization and validation detail. Appendix S3. Landscape analysis detail. Appendix S4. Dynamics of SOC changes and N 2 O emissions. Figure S1. Climate range covered in the full switchgrass calibration & validation dataset as per the NARR database. Figure S2. Classification of calibration and validation dataset field trial sites by soil surface texture and NRCS land capability class (LCC) rating. Figure S3. Scatter matrix of site location and climate parameters in the switchgrass parameterization & validation dataset. Figure S4. Scatter matrix of soil parameters in the switchgrass parameterization & validation dataset. Figure S5. Switchgrass green-up as a function of latitude only. Figure S6. Heading dates estimated as a function of ecotype and site latitude. Figure S7. Peak biomass as estimated to occur three weeks after heading, as compared to empirical observations for upland switchgrass at multiple sites. Figure S8. Temperature stress response curve based on experimental data from a variety of sources. Figure S9. Moisture stress response curve compared to normalized experimental data from Xu et al. (2006) and model defaults. Figure S10. Modeled and measured yield ranges binned by site average annual growing degree day accumulation (calculated for the range 12-30°C). Figure S11. Modeled and measured yield ranges binned by site average annual precipitation. Figure S12. Modeled and measured yield ranges binned by site NLCD land capability classification. Figure S13. Example of an observed-modeled comparison for a study  where time-resolved data is available for more detailed comparison, in this case for total aboveground carbon and belowground live carbon. Figure S14. Modeled vs. measured changes in soil organic carbon by study for all studies in the parameterization & validation dataset that include usable SOC measurements. Figure S15. Modeled vs. measured cumulative annual emissions of nitrous oxide by study for all studies in the parameterization & validation dataset that include usable N 2 O measurements. Figure S16. Histogram showing the area distribution of the polygons generated during the GIS intersect operation. Figure S17. Annual fluxes, assuming moderate tilling during field preparation for initialswitchgrass crop establishment and every 8 years when the crop is replanted.