Introduction

Sustainability denotes the capacity to meet current needs while also maintaining options for future generations to meet their needs ([7], volume 2). For bioenergy systems, sustainability refers to maintaining high crop yield while promoting positive environmental, economic, and social effects while minimizing potential negative impacts of bioenergy production [13]. In particular, sustainable bioenergy crop production often refers to maintaining high levels of plant productivity with minimal fertilizer and irrigation under a wide range of environmental conditions, including on marginal or degraded lands [6, 21, 41, 43, 58]. In addition to productivity, other indicators may be used to evaluate environmental impacts of bioenergy crop production, such as soil quality, water quality and quantity, greenhouse gases (GHG) emissions, biodiversity, and air quality ([7, 40], volume 2). Comprehensive field experiments that measure several environmental impacts require detailed data collection over many years and are costly. For these reasons, crop modeling has become an important tool for assessing the sustainability of biofuel production systems in a timely manner and over a wide range of environmental conditions [12, 34, 36].

Crop models simulate several dynamic processes such as water balance, soil nutrient cycling, and crop growth under site-specific weather, soil, and management conditions [14]. Crop models have been used to estimate yield for a variety of soils, response to climate change, best management practices, trade-offs among different cultivars, and gene by environment interactions for quantitative traits with benefits for plant breeding [4, 6, 24, 53]. Recently, crop models have also been used to assess sustainability of complex agricultural systems or ecosystems as they can predict long-term impacts on GHG emissions, nitrogen loss, and soil carbon sequestration [21, 26]. For instance, Izaurralde et al. [26] used the Erosion/Productivity Impact Calculator (EPIC) model to determine that decreasing tillage can change a watershed from acting as a carbon source to a carbon sink. Gelfand et al. [21] also used the EPIC model to determine that biomass production on marginal lands in the Midwestern USA can reduce GHG emissions.

Over time, many crop models have become more complex and some of the equations and parameters have become model-specific, resulting in a large number of biophysical and physiological parameters (see Di Paola et al. [14] for a review of the available crop models for studies and assessments in agriculture). Crop models are diverse, and some parameters in specific models are unknown or difficult to experimentally measure. Therefore, these parameters may be estimated from empirical models or borrowed from a similar crop, resulting in a potentially large uncertainty for these parameters [25, 57, 61]. Understanding how model predictions may change due to uncertainty in model parameters is important for evaluating and improving model performance. Additionally, identifying the influential parameters for a specific process may also help guide experimental research to most efficiently parameterize these complex models.

Sensitivity analysis is a commonly used method to assess parameter uncertainties and dependencies [55]. Although sensitivity analysis is not an obligatory element within the modeling process, it is recommended as a tool to ensure the quality of models and model prediction [16]. The aim of sensitivity analysis is to determine how sensitive the model output is, with respect to the parameters of the model, which are subject to uncertainty or variability. Sensitivity analysis may help guide experimental data collection by highlighting pathways or parameters with large impact. For example, sensitivity analysis for wheat yield showed complex interactions between phenological and physiological traits in high and low yielding environments; the rate and duration of grain filling were identified as the most influential parameters for wheat yield [3]. Wang et al. [62] conducted a sensitivity and uncertainty analysis with EPIC and found that radiation use efficiency, harvest index, heat units during growing season, and soil water content were the influential parameters for corn yield and soil organic carbon. To our knowledge, there are no systematic sensitivity analyses of crop models performed on a wide range of parameters for perennials. However, Kiniry et al. [33] and Kim et al. [29] have changed the values of a few selected crop parameters and concluded that the yield of switchgrass was sensitive to radiation use efficiency, water use efficiency, and nitrogen application.

Switchgrass (Panicum virgatum L.) is one of the leading second generation perennial crops capable of maintaining or enhancing soil quality through impacts on features such as soil carbon [1, 56, 60]. Previous studies of switchgrass have largely focused on modeling biomass production by parameterizing aboveground growth, because the most commonly measured parameters are those governing leaf area development and aboveground biomass accumulation [6, 11, 23, 32]. For example, Behrman et al. [6] modeled the differential growth for two lowland cultivars (Alamo and Kanlow) and two upland cultivars (Cave-in-Rock and Blackwell) using ALMANAC (Agricultural Land Management and Numerical Assessment Criteria) by parameterizing potential heat units during the growing season, maximum leaf area index, and light extinction coefficient while all belowground parameters remained constant for all the cultivars. This is largely due to the difficulty of measuring root growth and very few studies have attempted to quantify the belowground growth parameters necessary for crop simulation models [39]. However, these data might become more abundant in the future due to the development of advanced phenotyping techniques for root growth [9, 27]. Until this data is available, assessing the impact of belowground parameters together with aboveground parameters is a necessary first step to further evaluate model performance, identify critical model parameters for future data collection, and identify interesting allometric and tradeoff relationships that contribute to plant performance.

Recent studies have examined several environmental aspects of sustainable bioenergy production of switchgrass, including long-term GHG impacts and changes in soil organic carbon [12, 34, 36]. For instance, Lee et al. [36] modeled the long-term (100-year) production of six switchgrass cultivars along with soil carbon sequestration and N2O emission in a Mediterranean climate using the DAYCENT model. Their results suggested that nitrogen fertilizer and irrigation are necessary for switchgrass cultivation to maintain biomass production and act as a soil carbon sink. Gelfand et al. [21] evaluated the GHG emission in ten Midwest US states using EPIC and found that successional herbaceous vegetation including switchgrass had a direct GHG emission mitigation capacity. However, identifying which model parameters are influential for switchgrass productivity and how these parameters impact switchgrass sustainability remains unclear, as many of the plant growth parameters are interrelated and can influence more than one core process.

The objective of this study is to identify influential aboveground and belowground parameters on three metrics of sustainability of switchgrass using the APEX (Agricultural Policy/Environmental eXtender, [66]) model, based on two experiments conducted at USDA-ARS Grassland, Soil and Water Research Laboratory in Temple, TX (TMP), and Lady Bird Johnson Wildflower Center in Austin, TX (WFC), for three lowland switchgrass genotypes. First, the APEX model (v. 1501) was parameterized to simulate annual aboveground growth dynamics of switchgrass based on measured data collected at TMP and WFC, respectively, for each genotype. Next, sensitivity analyses of model parameters were conducted for three sustainability-related variables (productivity, nitrogen loss, and soil carbon change at the top 30-cm soil on an annual basis) using the enhanced Morris method. Lastly, the individual relationship and two-way interactions between each influential parameter identified with the three response variables were analyzed.

Materials and Methods

APEX Model Description

The APEX model was developed as an extended tool of EPIC to simulate management on whole farms or small watersheds to assess sustainable production efficiencies and maintenance of environmental quality [66]. The APEX uses data on weather, soil properties, and management practices to characterize the system under study and to track water, nutrients, carbon, and plant growth on a daily basis. The APEX is based on hydrology, water movement, and can be used to evaluate the interactions between fields within a small watershed with respect to nutrient and pesticide transport [18, 46, 50, 59]. The APEX can also be used to assess the impact of land management practices, conservation practices, and global climate change. Furthermore, the APEX has the ability to assess the sustainability of agricultural system on a whole farm for a long term (1–4000 years) [17, 20].

The APEX simulates hydrology, soil erosion, nutrient cycling, and plant growth processes, which are all inherently interrelated. Water enters the system as precipitation or irrigation and infiltrates the soil, then is lost through surface runoff, evapotranspiration, percolation, and lateral subsurface flow. Nutrients are added as fertilizer or manure, which are either lost through leaching, runoff, denitrification, and volatilization, or taken up by plants and later removed by harvest.

The APEX model contains three main categories of input data. The first category is site-specific input data, which includes daily weather (maximum and minimum temperatures, precipitation, etc.), soil physical properties (soil hydrologic group, texture class, etc.) and chemical properties (soil pH, organic C content, etc.), topographic properties (field slope, field size, etc.), and crop management practices (planting date, fertilizer date and amounts, harvest date, etc.). The second category is definitional input data, such as crop parameters (potential radiation use efficiency, harvest index, etc. contained within the CROP file), tillage practices (tillage depth, harvest efficiency, etc.), fertilizer composition (organic N and P content, etc.), and pesticides (pesticide solubility, half-life in soil, etc.). Lastly, there is a general parameter file (PARM), which contains many equation coefficients, S-shape curve coefficients, and other general parameters used for model calibration. Detailed descriptions of the model components and processes can be found in the theoretical documentation of APEX [59].

APEX Model Setup and Parameterization

APEX model simulations were set up to represent two experiments conducted at USDA-ARS Grassland, Soil and Water Research Laboratory in Temple, TX (TMP, 31.04 N, 97.35 W, elevation 199 m above sea level), and Lady Bird Johnson Wildflower Center (WFC) in Austin, TX (WFC, 30.18 N, 97.87 W, elevation 149 m). The two sites differed largely in soil texture and depth. The switchgrass genotypes were transplanted from gallon pots and well-watered to facilitate establishment under a rainout shelter in 2011. A number of precipitation treatments were implemented in subsequent years based on historic climatic conditions for Central Texas. However, here, we explore switchgrass growth across ambient precipitation conditions implemented as an irrigation treatment from 2012 to 2015 (see Aspinwall et al. [2] for more details). Data collected on the green-up day, flowering day, and leaf area index (LAI) provided the basis for parameterizing the aboveground growth parameters of the APEX model. Green-up day was recorded and identified visually as the day of year when 50% of crown tillers were green. Flowering day was recorded as the day of year when 50% of plant tillers reached full anthesis. LAI was measured using the AccuPAR LP-80 ceptometer (Decagon Devices, Inc., Pullman, WA, USA). Two perpendicular light measurements were taken through the center of the plant at a height of 10 cm off the ground [2].

For this modeling exercise, data collected on the three lowland genotypes (WBC, AP13, and KAN) grown under ambient water conditions in 2012 was used for model calibration. These three lowland genotypes represent cultivated switchgrass genotypes and genotypes found currently in a natural population. AP13 and KAN are clones of individual plants originating from cultivars ‘Alamo’ and ‘Kanlow’, respectively, both of which were bred for their high yield and adaptability. ‘Alamo’ was originally from Live Oak County, TX, and is adapted throughout the majority of the USA. ‘Kanlow’ was originally from Wetumka, OK, and was selected for use in low wetland areas and land areas subject to excessive moisture. WBC (West Bee Cave) is a clone of an individual plant collected in West Bee Cave, TX (30.32 N, 97.98 W), and is presumably locally adapted to central Texas [2, 42]. Similar to other lowland genotypes, WBC exhibits delayed flowering, higher tiller number and higher tiller mass than upland genotypes [42].

As mentioned earlier, the APEX requires three main categories of input data. For the site-specific input data, weather station data spanning 20 years (1996–2015) including daily solar radiation (MJ m−2), maximum temperatures (°C), minimum temperature (°C), precipitation (mm), for TMP and WFC, were obtained from USDA-ARS (https://www.ars.usda.gov/, accessed in 2016) and NOAA (https://www.ncdc.noaa.gov, accessed in 2017), respectively. A 20-year simulation was used to evaluate relatively long-term environmental effects of switchgrass growth and ensure a wide variety of year-to-year climate variation that switchgrass may encounter for biomass production. The annual precipitation ranged from 364 to 1288 mm, and 408 to 1523 mm for these 20 years at TMP and WFC, respectively (Supplemental Fig. S1). The average monthly maximum temperature was from 10.7 to 39.0 °C for TMP and 12.3 to 40.5 °C for WFC, and the average monthly minimum temperature ranged from 1.1 to 25.1 °C for TMP and 2.6 to 25.8 °C for WFC (Supplemental Fig. S2). Soil data including physical, chemical, and textural properties of the Austin silt clay soil series for TMP and the Speck clay loam soil series for WFC, were obtained from USDA-NRCS Web Soil Survey (https://websoilsurvey.sc.egov.usda.gov/, accessed in 2016). TMP soils are fine-textured with medium to rapid runoff and moderate to low permeability. WFC soils are coarse-textured with low runoff and slow permeability. A 1-ha field with field slope of zero was assumed adequate to represent the topographic properties for this study. To be consistent with the management practices from 2012 to 2015, the seedlings were transplanted on March 10 of the first simulation year and plants were hand-harvested every subsequent year on November 10 with a cutting height of 25 cm and aboveground harvest efficiency of 0.95. Automatic fertilization was turned on, assuming no fertilizer stress for the 20-year simulation and reflecting the field situation from 2012 to 2015. The methods selected for hydrology and other components of the model are listed in Supplemental Table S1 and follow the suggestions made by Wang et al. [63].

One type of definitional input required by the APEX model is the specification of the crop growth parameters in the CROP file. For this study, the crop parameters were calibrated to the data collected from the 2012 season. Specifically, switchgrass parameters were derived by assuming green-up was the day when 10% of the maximum leaf area was reached, flowering was assumed to be 80% of maximum leaf area was reached, and leaf area was assumed to start declining after flowering. Based on green-up and flowering assumptions, DLAP1 and DLAP2 (two points representing the leaf area development curve, see Table 1) were parameterized to 10.35 and 35.99, respectively, for all three genotypes. Heat units accumulated during the growing season (PHU) from 1 week before green-up to harvest on November 10. Heat units were calculated as the total sum of the daily average temperature (average of the daily maximum temperature and minimum temperature) minus the base temperature for switchgrass growth (TBS). In this study, TBS and TOP (optimum temperature for switchgrass growth) were set to 12 and 25 °C, respectively [32]. Thus, PHU was calculated for each year with an average of approximately 3500 °C, which was different from the PHU values measured to maturity that are commonly reported in the literature [31]. Other CROP parameters were specified based on values from literatures or left as a default value from the APEX model [59] as exemplified for WBC in Table 1. For the PARM file, parameter P6, which causes dormancy in plants over winter was set to 0.5 h, the average value between the default ranges, to represent a perennial grass with winter dormancy. Other PARM parameters used the default setups in the APEX model due to lack of information for parameterization. The growth pattern and biomass production for the other two lowland genotypes, AP13 and KAN, were very similar (supplemental Fig. S5), and so we utilized the same parameter settings for these genotypes. The parameterization for AP13/KAN was also similar to that of WBC (Table 1) except for lower radiation use efficiency (WA, 25 kg ha−1/MJ m−2) and maximum leaf area index (DMLA, 8) values.

Table 1 Biologically relevant CROP and PARM parameters selected from the APEX for the sensitivity analysis
Full size table

The parameterized APEX model produced reasonable intra-annual growth accumulation compared to phenology data and measured LAI accumulation as exemplified for the lowland genotype WBC at TMP in 2012 (Fig. 1a). Modeled LAI tracked the predicted green-up day and flowering day, falling into the ranges of observed green-up days (83 ± 5 days) and flowering days (203 ± 27 days), respectively. The model also tracked the measured LAI accumulation well, except for one measurement in the middle of growing season, which might be caused by high temperatures during that period, thus the decline was not included in the parameterization to avoid overfitting. The validation of the model with WBC biomass production from 2012 to 2015 showed good parameter calibration since the simulated biomass fell nicely between the upper and lower quantiles of the box plot for the observed biomass (Fig. 1b). Note that the biomass was measured on per plant basis, causing large variation from plant to plant. The 20-year simulation also displayed relatively stable biomass production across years, ranging from around 20–40 Mg ha−1 except the low biomass production in 2010 (Supplemental Fig. S3). The low biomass may be due to the high aeration stress from flood at Temple, TX, in 2010 that was captured by the model.

Fig. 1
figure 1

a Modeled LAI versus measured LAI in reference to phenology observations recorded on green-up date and flowering date in 2012 for lowland genotype WBC, along with daily maximum temperature (TMAX) and daily average temperature (T ave) at Grassland, Soil and Water Research Laboratory, Temple, TX. b Boxplot of the observed biomass versus simulated biomass for lowland genotype WBC from 2012 to 2015

Full size image

Additionally, the WBC crop parameters (Table 1) calibrated for TMP were used to independently validate the LAI development at WFC by changing the soil series and weather station. The simulated LAI curve of WBC at WFC was similar to that from TMP and was comparable to the observed LAI (supplemental Fig. S4), indicating that the parameters were realistically calibrated to represent WBC at locations differing in soils and weather. Similarly, calibration for AP13/KAN was based on LAI data collected in 2012 at TMP. The calibrated model captured the LAI development of these two lowland genotypes in 2012 and was validated for the biomass production from 2012 to 2015 at TMP (supplemental Fig. S5). The model was also well validated for the LAI development at WFC for AP13/KAN in 2012 (supplemental Fig. S6).

Sensitivity Analysis

There are approximately 60 parameters in the CROP file and 130 in the PARM file, which contains 30 two-part S-shaped parameters and 100 other parameters (coefficients, limits and fractions, etc.) in the APEX model. For this study, 25 CROP and 20 PARM parameters (shown in Table 1) were selected for sensitivity analysis of average annual total biomass (sum of above- and belowground biomass), total nitrogen loss (sum of soluble nitrogen in surface runoff, subsurface runoff, and percolation; sediment-bound nitrogen from water and wind erosion; nitrogen lost through denitrification and volatilization), and soil carbon change (difference between ending and beginning soil carbon at the top 30-cm soil) over the 20-year simulation time period. The CROP and PARM parameters were selected after going through all the CROP and PARM parameters and were considered as biologically relevant parameters to any one of the three metrics. In other words, these parameters were considered to be related to plant growth and plant nutrient cycling and were relevant to sustainability of the agroecosystem.

Sensitivity analysis identifies the governing parameters in complex simulation models by calculating how much a model output changes due to a minor change in its input parameter value [35, 47, 52]. For this study, a quantitative measure of sensitivity was conducted using the enhanced version of Morris method, the measure of which is the mean of the absolute values of Morris local measures and is believed to be good for ranking parameters in order of importance [8]. The Morris method, also called one-parameter-at-a-time method, deals efficiently with models containing a large number of parameters without relying on strict assumptions [45]. Specifically, the method starts by sampling a set of starting values within a defined interval for all parameters based on sample statistics and the associated probability distribution (normal, uniform, etc.) of each parameter. Then, focal parameters are changed one at a time, while all other parameters are held constant, to obtain the corresponding model outcomes across all study parameters. The procedure can be repeated with different sets of starting values (r) within the specific range for each parameter, which leads to r*(k + 1) model runs, where k is the number of parameters being evaluated in our study.

For this study, each of the selected parameters was subdivided into 10 equal intervals (assuming uniform distribution) and sensitivity analyses were conducted for CROP and PARM parameters separately due to their different magnitudes and ranges in parameter distribution. CROP parameters were set to plus and minus 25% of the respective starting values used for modeling each lowland genotype (Table 1). PARM parameter ranges were set to the default limits of each value ([66]; Table 1). The number of sets of starting values (r) was determined to be 3000 because the sensitivity indices had reached to the convergence at 3000 after testing different r values (1000; 3000; and 5000).

In the APEX, there are some two-part parameters that set the inflection points for the S-shaped curves. These are DLAP1, DLAP2, FRST1, and FRST2 (see variable definitions in Table 1). These two-part parameters are separated by a decimal. In addition to being two-part parameters, parameters with the same core names (DLAP and FRST) are also sequentially dependent on each other and therefore were treated differently. For example, DLAP1 and DLAP2 are the two points on the optimal leaf area development curve. The value before the decimal point is the percentage of the growing season PHUs and value after the decimal is the percentage of maximum potential LAI reached. For our model setup, DLAP1 is set to 10.39, which means that at 10% of the growing season PHUs, LAI can increase to 39% of the maximum LAI. Because DLAP1 and DLAP2, and FRST1 and FRST2 are related, an additional constraint was imposed on these two-part parameters such that both parts of each parameter in DLAP1 and FRST1 had to be less than DLAP2 and FRST2.

Model Run and Statistics

Parameter sets were generated using the Morris method in R package “sensitivity” [48]. When conducting sensitivity analyses for CROP parameters, the PARM parameters remained the same for all simulations. The same procedure was used for sensitivity of the PARM parameters, except the CROP parameters remained set for all simulations. Two indices from the enhanced Morris method were calculated for assessing the influence of each parameter on the three sustainability metrics. The first is the mean of the absolute values of Morris local measures (μ *), which represents the overall importance of a specific parameter. The second is the sigma square (δ 2), which is used to detect factors involved in interaction with other factors or whose effect is non-linear [45]. The rank of each influential parameter and the response value range (the minimum value, the maximum value) were summarized. Furthermore, the percentage change of the three sustainability metrics was also calculated for each of the influential parameter changes within its range (Table 1). The percentage change was calculated as the difference between the maximum and the minimum value of the response variable divided by the minimum value (Table 2). Subsequently, when investigating the functional form of the relationship between an influential parameter and a response variable (Fig. 3), all the other parameters remained constant and only the influential parameter was changed within the specified range (Table 1). A similar procedure was adopted when exploring the two-way interactions of the most influential parameters for three response variables (as shown in Figs. 4 and 5). The two parameters were changed jointly within their ranges while all other parameter values remained set.

Table 2 The rank of each influential parameter for CROP and PARM parameters calibrated for genotype WBC and genotypes AP13/KAN in Temple, TX, the response range, and the percentage changes of each sustainable metric (total biomass, Mg ha−1, nitrogen loss, kg ha−1 and soil carbon change, Mg ha−1) with respect to the changes of each of the influential parameters
Full size table

Results

Sensitivity Analysis Across Genotypes

Results of the sensitivity analysis for the respective sets of WBC and AP13/KAN parameters at TMP and WFC highlight the four most influential CROP and PARM parameters impacting total biomass, nitrogen loss, and soil carbon change based on the two Morris sensitivity indices, μ * and δ 2 (with criteria of μ * value larger than 1). Our sensitivity results showed that the importance of parameters were not dependent on the genotypes (Fig. 2 for WBC and supplemental Fig. S7 for AP13/KAN at TMP; supplemental Table S2 at WFC), but differed slightly depending on the location for each sustainable metric. In other words, the same set of influential parameters was identified across genotypes at TMP and another same set of influential parameters was identified across genotypes at WFC. Additionally, TOP (optimum temperature for crop growth) and P2 (the coefficient describing the impact of bulk density on root growth) had the largest δ 2 (Fig. 2, and supplemental Fig. S7) across genotypes, showing high interaction with other parameters or non-linearities in affecting total biomass, nitrogen loss, and soil carbon change in the model.

Fig. 2
figure 2

Influential CROP and PARM parameters from the APEX model for total biomass, nitrogen loss, and soil carbon change in the top 30-cm soil of lowland genotype WBC. a–c CROP parameters in assessing total biomass, nitrogen loss, and soil carbon change, respectively; d–f PARM parameters in assessing total biomass, nitrogen loss, and soil carbon change, respectively

Full size image

Sensitivity Analysis Across Locations

A different set of influential parameters was identified for each location. However, the top three CROP and PARM parameters at both locations remain consistent for each sustainability metric, despite small changes in the ranking order (Table 2 and supplemental Table S2). Specifically, total biomass was primarily impacted by TOP, TBS (base temperature for crop growth), and WA (potential radiation use efficiency). However, DLAI (the fraction of growing season in heat units when LAI begins to decrease due to senescence) was identified as the fourth most influential CROP parameters at TMP (Fig. 2a and supplemental Fig. S7A) while RDMX (maximum root depth in the soil) was identified at WFC (supplemental Table S2). The top three influential PARM parameters for biomass production were P2, P38 (the coefficient for plant water stress), and P6 (the photoperiod parameter for winter dormancy). The fourth most influential parameter identified at TMP was P69 (the coefficient for microbial activity in the top soil layer) (Fig. 2d) while P5 (soil water lower limit) was identified at WFC (Table S2). Additionally, RDMX was identified as an influential parameter for edge-of-field nitrogen loss at both locations; however, BN3 (fraction of nitrogen in biomass at maturity) was identified at TMP and DLAI was identified at WFC (Table 2 and S2). Not surprisingly, PARM parameters P72 (partitioning of volatilization/nitrification), P80 (upper limit of nitrification-volatilization), P4 (water storage nitrogen leaching), and P2 were identified as influential for nitrogen loss (Fig. 2e) at TMP. However, parameters P36 (upper limit of daily denitrification rate) and P6, instead of P80 and P72 at TMP, showed a big impact on nitrogen loss at WFC (supplemental Table S2), indicating that sensitivity analysis results may differ due to different soil series. In total, we found that 6 out of 25 CROP parameters and 9 out of 20 PARM parameters were very influential across these two locations, and other parameters showed relatively little impact on the three metrics of lowland switchgrass sustainability.

Sensitivity Analysis for each Sustainable Metric

Some parameters were influential across more than one sustainability metric. For example, the CROP parameters, TBS and TOP, were always the top two parameters for each sustainable metric across the three lowland genotypes at the two locations (Fig. 2a, b, c), thus highlighting the importance of interrelated parameters and processes in the APEX model. Additionally, both CROP and PARM parameters that had a large influence on total biomass also had a large influence on soil carbon change (Fig. 2a, c; Table 2 and S2). Since our analysis revealed that the sensitivity results are not highly dependent on genotype and the top three parameters are always the same for each sustainability metric across the two locations, for simplicity, the following results are only presented for the WBC genotype at the TMP location.

Relationships Between Influential Parameters and Sustainability Metrics

The four parameters identified to have the largest impact on total biomass, nitrogen loss, and soil carbon change for genotype WBC were further investigated by analyzing the functional form of the relationship between each parameter and each of the three sustainability metrics (Fig. 3). The response value range and percentage change for the three sustainability metrics with respect to each of the influential parameters are presented in Table 2. Among the influential CROP parameters, both TOP and TBS had a non-linear relationship with all three response variables (Fig. 3a, b, i, j, q, r). In contrast, WA and DLAI had linear or nearly linear positive relationship with total biomass and soil carbon change (Fig. 3c, d). It is not surprising that higher radiation use efficiency and later leaf area decline resulted in higher biomass production. RDMX and BN3 influenced nitrogen loss in a step-like and a linear function, respectively (Fig. 3k, l). Among PARM parameters, P2 had a non-linear relationship with total biomass and soil carbon change (Fig. 3e, u). Parameters P38 and P69 had a linear or nearly linear influence on total biomass and soil carbon change (Fig. 3f, g, h, w, x) while parameter P6 had a linear influence on total biomass but non-linear influence on soil carbon change (Fig. 3h, u). Parameters P2, P4, and P80 affected nitrogen loss approximately linearly when their values were at the lower end of the ranges assayed; beyond these ranges, these parameters had no impact on nitrogen loss.

Fig. 3
figure 3

The relationship between the top four individual influential CROP and PARM parameters identified for each sustainability metric: total biomass (a–h), total nitrogen loss (i–p), and change in soil carbon (q–x)

Full size image

For total biomass production, the most influential CROP and PARM parameters are TOP and P2, which resulted in 218 and 28% change in biomass production, respectively (Table 2). Total biomass increased linearly from 10.6 to 33.6 Mg ha−1 as TOP increased from about 19 to 25 °C (all the other parameters remained set to the values for WBC in Table 1). Total biomass peaked at 33.6 Mg ha−1 when TOP equaled 25 °C, and total biomass decreased to 27.0 Mg ha−1 as TOP increased beyond 25 to 31 °C (Fig. 3a). Biomass increased as P2 increased until P2 reached 1.5, at which point total biomass plateaued at 34.1 Mg ha−1 (Fig. 3e). On average, the four influential CROP parameters for biomass resulted in 80% change of biomass, while the four influential PARM parameters resulted in 13% change of biomass (Table 2). Therefore, switchgrass biomass production is relatively more sensitive to changes in the CROP parameters than that of the PARM parameters.

TOP and P2 are the most influential CROP and PARM parameters for nitrogen loss, resulting in 150 and 319% change in nitrogen loss, respectively (Table 2). Nitrogen loss declined from 49.7 to 19.9 Mg ha−1 as TOP increased from 19 to 21 °C and nitrogen loss showed little change as TOP continued to increase from 21 to 31 °C (Fig. 3i). Nitrogen loss increased linearly from 3.1 to 13.0 kg ha−1 as P72, the partitioning coefficient between volatilization and nitrification, increased from 0.05 to 0.5 (Fig. 3m). Additionally, RDMX, the belowground rooting depth parameter, was identified as influential for nitrogen loss, resulting in 51% change in nitrogen loss (Table 2). Nitrogen loss remained constant at 31.1 Mg ha−1 as RDMX increased from 1.5 to 1.9 m and decreased to 20.6 kg ha−1 for all RDMX greater than 1.9 m (Fig. 3k). On average, the four influential CROP parameters caused 56% changes in nitrogen loss, while PARM parameters caused 152% changes in nitrogen loss, indicating nitrogen loss in the APEX model is more sensitive to changes in the PARM parameters than that of the CROP parameters.

TOP, TBS, WA, and DLAI, which were identified as influential for modeling total biomass, were also identified as influential for soil carbon change in the top 30-cm soil (Fig. 2, Table 2), indicating there is a high correlation between total biomass and soil carbon change in the model. The most influential CROP and PARM parameters are TOP and P6, resulting in 278 and 32% change in soil carbon change, respectively (Table 2). The relationship of TOP with soil carbon change (Fig. 3q) is similar to the relationship of TOP with total biomass (Fig. 3a), implying that there are strong interactions between two processes. Soil carbon remained constant at 19.4 kg ha−1 as P6 increased from 0.0 to 0.5, but soil carbon change increased to 25.5 kg ha−1 as P6 increased to 1 (Fig. 3u). Similar to biomass production, soil carbon is more sensitive to changes in the CROP parameters (averaged 189% change in soil carbon) than that of the PARM parameters (averaged 26% change in soil carbon) (Table 2).

Two-Way Interactions of Influential Parameters

The two-way interactions between the four most influential parameters for total biomass, nitrogen loss, and soil carbon change are shown in Fig. 4 for CROP parameters and Fig. 5 for PARM parameters. The response range and percent change of each sustainability metric to changes in two influential parameters are summarized in Table 3. For biomass production, the interaction between TOP and TBS from CROP and the interaction between P2 and P38 from PARM are most influential, resulting in 983 and 145% change in biomass production, respectively. The interaction of TOP with TBS revealed that low biomass occurs when TBS is high and TOP is low (Fig. 4a). Total biomass declined as P38 increased from 0.5 to 1.0 while P2 was between 1.0 and 1.2 (Fig. 5a). On average, biomass production is more sensitive to the interaction of CROP parameters (517% change) compared to that of the interaction from PARM parameters (73%) (Table 3).

Fig. 4
figure 4

Two-way interactions between the most influential CROP parameter (TOP) with other top three influential CROP parameters for each sustainability metric: total biomass (a–c), nitrogen loss (d–f), and change in soil carbon (g–i)

Full size image
Fig. 5
figure 5

Two-way interactions between the most influential PARM parameter (P2) with other top three influential CROP parameters for each sustainability metric: total biomass (a–c), nitrogen loss (d–f), and change in soil carbon (g–i)

Full size image
Table 3 The sensitivity ranges and the percentage changes of each response variable (total biomass, Mg ha−1, nitrogen loss, kg ha−1, and soil carbon change, Mg ha−1) with respect to the changes of two influential parameters from CROP and PARM parameters
Full size table

The average of the four two-way interactions showed that nitrogen loss is more sensitive to the interactions from PARM parameters (average 304% change) than that from CROP parameters (average 268% change) (Table 3). However, the response range for the CROP parameters (18.9–110.8 kg ha−1) is larger than that for the PARM parameters (2.9–36.1 kg ha−1). The strongest response is the two-way interaction between P2 and P4, which resulted in 806% change in nitrogen loss as P2 and P4 changed jointly within their respective ranges defined in Table 1. When P2 and P4 were at their smallest values (1.0 and 0.02, respectively), relatively high nitrogen loss (36.1 kg ha−1) occurred and as values for both parameters increased, nitrogen loss declined (Fig. 5e). The step function response identified for RDMX is also carried over in the interaction between TOP and RDMX with an abrupt vertical decline in nitrogen loss when RDMX reached 1.9 m (Fig. 4b).

The two-way interactions of CROP parameters for change in soil carbon (Fig. 4g–i) were similar to total biomass (Fig. 4a–c), with the strongest interaction between TOP and WA. As TOP and WA increased, soil carbon increased and became positive (Fig. 4h). Generally, as TOP increased, soil carbon increased from negative values to positive values (Fig. 4g–i), indicating that the soil was switching from a carbon source to a carbon sink. For the two-way interactions among PARM parameters, change in soil carbon was always positive and soil was serving as a carbon sink. The average of all two-way interactions of CROP parameter resulted in 272% change in soil carbon compared to the 88% change from PARM interactions (Table 3), indicating that interactions between CROP parameter are more influential than PARM interactions in affecting soil carbon change.

Discussion

Quantifying model uncertainty is important for understanding complex model behavior and model parameterization [44]. To date, no formal sensitivity analyses have been performed on perennials grown for biofuel and most sensitivity analyses performed for annual crops have mainly focused on aboveground growth parameters [35, 52]. This study investigates 45 parameters relevant to growth and nutrients cycling of switchgrass to analyze the impact from both above and belowground parameters for perennial grass. In addition to the five aboveground CROP parameters, we identified one influential belowground CROP parameter, RDMX, and several influential PARM parameters related to belowground nutrient cycling in the APEX model (Fig. 2), emphasizing the importance of understanding influential belowground parameters in crop modeling efforts. Furthermore, we analyzed the relationships and interactions for the response variables with respect to the influential parameters, providing information on the model functionality and the interrelated processes.

We identified 15 parameters in the APEX model that are highly influential for estimating total biomass, nitrogen loss, and soil carbon change of three lowland switchgrass genotypes across two locations and a majority of the remaining parameters (30 parameters) had relatively little impact. Similar results were found by Richter et al. [52] when conducting a sensitivity analysis of Durum wheat (Triticum durum). They identified 10 and 12 parameters out of 46 crop and site-specific parameters for evaluating photosynthesis and transpiration, respectively, and they suggested that their model could be simplified. Liu et al. [37] identified five influential crop parameters out of 22 selected parameters for evaluating winter wheat yield at a regional scale. Confalonieri et al. [10] concluded that rice yield was mostly influenced by the parameters relating to energy to biomass conversion and partitioning. These studies suggested that some crop models could be simplified in terms of number of parameters depending on the response variables of interest. However, caution should be taken when simplifying the complex structure of crop models because some model parameters may have little effect on certain response variables, but they may play key roles in other processes [65]. Additionally, the dynamic characteristics of parameter sensitivity need to be considered because parameter sensitivity may evolve over time [35].

Understanding the sensitivity of model outputs to changes in crop parameters for switchgrass is important because different values for two of the most influential crop parameters, TBS and TOP, are used in the literature [5, 22, 32, 64]. For switchgrass, three different values of TBS (13, 12, and 10 °C) have been reported by Grassini et al. [22], Kiniry et al. [32], and Wang et al. [64], respectively. Similarly, different values of TOP (33 and 25 °C) have been reported by Grassini et al. [22] and Kiniry et al. [32], respectively. TBS is used to calculate the growing season heat units in the model, so it has several impacts on crop growth including growth onset and heat stress. The WA values reported for switchgrass range from 22 to 52 kg ha−1 MJ−1 m−2 depending on cultivar and soil [30, 38]. Therefore, selection of proper values for each genotype is necessary to accurately estimate switchgrass sustainability.

For switchgrass, biomass production and soil carbon change are more sensitive to changes in CROP parameters than those of PARM parameters. Switchgrass alternated between serving as a carbon sink and a carbon source depending on model parameterization (Figs. 3 and 4). Therefore, determining the appropriate value for a specific genotype is necessary to determine if soil is maintaining or increasing carbon storage and may help provide recommendations for alternate management practices to increase carbon storage. On the other hand, nitrogen loss at the edge-of-field is more sensitive to adjustments in PARM parameters than that of CROP parameters, indicating these influential PARM parameters are just as important as the CROP parameters in estimating nutrient cycling and associated nitrogen loss. Additionally, the identification of belowground processes-related CROP and PARM parameters indicated that belowground parameters are just as important as aboveground parameters. Field experiments that can estimate the values of these influential belowground processes-related parameters will help improve the model calibration and performance.

Crop models, like the APEX, are complicated and the equations describing the interactions between model parameters are often not easy to decipher or pull part. Here, we investigated the relationships of the influential parameters identified with respect to each measure of sustainability. Understanding the two-way interactions helps provide more information on the complexity of parameter interactions and model structure for estimating total biomass, nitrogen loss, and soil carbon change. Some of the linear or non-linear relationships are intuitively results of the equations embedded within the APEX model. For example, total biomass is a linear function of WA and WA is not interacting with other parameters in determining biomass production. On the other hand, non-linear relationships may be due to the interactions among parameters affecting an intermediate process. For instance, biomass production is simulated as a function of many intermediate processes, such as leaf area index, heat unit accumulation, crop stress factors, and crop developmental stages. Each of these intermediate processes is influenced by changes in the crop parameters, TBS, and TOP, as well as those governing the shape of potential leaf area development curve over time, such as DLAI. These complex interactions among many processes and parameters lead to non-linear interactions between the parameters and biomass.

The APEX model can simulate aboveground leaf development fairly well under dryland conditions based on the data collected on three lowland genotypes in Temple, TX, and Austin, TX; however, the model did not capture a measured decline in LAI during the middle of the growing season for all the genotypes at both locations. This decline in LAI was most likely due to the relatively high maximum temperature (close to 40 °C) or average temperature (close to 36 °C) at around the end of June (Fig. 1 and Supplemental Fig. S4-S6). Another study documented similar leaf area decline under high temperatures for switchgrass [28]. Additionally, during the high temperature period, we observed that switchgrass leaves rolled causing a possible under-measurement of LAI in the experimental data collected by a ceptometer. Therefore, the decline was not included in the parameterization to avoid overfitting. Meanwhile, more data samplings or multiple year data collection may help further test the model.

The sensitivity analyses presented in this study were conducted at two locations which differed in soils and climate. Consequently, a different set of influential parameters was identified. However, the CROP parameters, TBS and TOP, were always the top two parameters for each sustainable metric across the three lowland genotypes at the two locations. The different influential parameters identified at the two locations may be largely due to differences in texture and depth of soil profile. Therefore, some caution should be taken when transferring the sensitivity analyses results to other environments. This is also consistent with the results from other studies in which the importance of each parameter may change across large environmental gradients, for different crops, or for different response variables of interest [15, 54]. For example, Ruget et al. [54] found that the influence of the model parameters in STICS (Simulateur mulTIdisciplinaire pour les Cultures Standard) depended on growing conditions and the output variables under consideration. Dzotsi et al. [15] also found the ranking of influential crop parameters in the SALUS (System Approach to Land Use Sustainability) crop model was strongly related to production level, location, and crop.

Further evaluation of model responses across diverse environmental conditions (temperature, precipitation, soils, etc.) with different switchgrass genotypes will improve our understanding of important crop (or cultivar) by environment interactions. This may also be important for understanding the impact of future climate change. If these influential model parameters can be linked to plant traits, then identification of influential plant traits in a specific environment could help identify the conditions for sustainable switchgrass production or help inform future switchgrass breeding. In addition, this study only evaluated three indicators of sustainability. There are many other important environmental variables to consider such as GHG emission and water quality. Data collection and model evaluation of these additional processes in the future will provide a more complete picture of the environmental impacts of switchgrass grown for biofuel and provide a better understanding of these interrelated processes.

Conclusion

Fifteen out of the 45 selected above and belowground APEX model parameters were identified as influential in estimating total biomass, nitrogen loss, and soil carbon change for three lowland switchgrass genotypes at two locations in Central Texas. The most influential CROP parameter for all three measures of sustainability is TOP (optimum temperature), which is capable of causing a 218, 150, and 278% change in total biomass, total nitrogen loss, and soil carbon change, respectively, when changing from its smallest value (18.8 °C) to its highest value (31.3 °C) assessed for the natural lowland genotype WBC. The most influential PARM parameters are P2 (the coefficient describing the impact of bulk density on root growth), P72 (the partitioning coefficient between volatilization and nitrification), and P6 (the photoperiod parameter for winter dormancy) for total biomass, nitrogen loss, and soil carbon change, respectively. Belowground processes-related parameters are just as important as commonly known CROP parameters. Measurements on these influential PARM parameters will help targeting important plant traits and better understanding various interrelated processes, especially the impact of variation in belowground processes. Designing experiments that can determine influential parameters such as TBS and TOP may help characterize the non-linear relationships between influential parameters and response variables, and help improve APEX model performance and our understanding of the possible underlying physiological mechanisms.