Responsiveness of miscanthus and switchgrass yields to stand age and nitrogen fertilization: A meta‐regression analysis

Optimal management of the perennial bioenergy crops, miscanthus and switchgrass, requires an understanding of their responsiveness to nitrogen (N) fertilizer at different maturity stages across locations and growing conditions. Earlier studies that have examined the yield response of these crops to N and stand age using field experiments or meta‐analysis techniques provide mixed evidence. We extend earlier studies by applying a multi‐level mixed‐effects (MLME) meta‐regression model to conduct a more extensive multivariate regression of yield response of these crops to N and stand age, while controlling for climate and location conditions and unobserved factors related to study design. Our findings are based on 1403 and 2811 yield observations for miscanthus and switchgrass, respectively, from experiments conducted between 2002 and 2019 across the rainfed region in the United States. We find statistically significant evidence that an additional year of maturity increases miscanthus and switchgrass yields but at a decreasing rate; yields peak at the 7th and 6th year respectively, for the observed range of applied N rates and stands. We also find that an increase in N application increases yield by a statistically significant level, but at a declining rate; the magnitude of the yield response to N is, however, small and varies with the age of the crop. The impact of N is larger on older compared to younger and middle‐aged stands of miscanthus. In contrast, the impact of N on switchgrass is larger on middle‐aged compared to younger and older stands of switchgrass. We do not find a statistically significant effect of soil productivity on yield for either crop. This analysis provides a basis for developing N application recommendations and optimal rotation age for miscanthus and switchgrass and shows that these energy crops can grow just as productively on low productivity land as on high productivity land.

model to conduct a more extensive multivariate regression of yield response of these crops to N and stand age, while controlling for climate and location conditions and unobserved factors related to study design. Our findings are based on 1403 and 2811 yield observations for miscanthus and switchgrass, respectively, from experiments conducted between 2002 and 2019 across the rainfed region in the United States. We find statistically significant evidence that an additional year of maturity increases miscanthus and switchgrass yields but at a decreasing rate; yields peak at the 7th and 6th year respectively, for the observed range of applied N rates and stands. We also find that an increase in N application increases yield by a statistically significant level, but at a declining rate; the magnitude of the yield response to N is, however, small and varies with the age of the crop. The impact of N is larger on older compared to younger and middle-aged stands of miscanthus. In contrast, the impact of N on switchgrass is larger on middle-aged compared to younger and older stands of switchgrass. We do not find a statistically significant effect of soil productivity on yield for either crop. This analysis provides a basis for developing N application recommendations and optimal rotation age for miscanthus and switchgrass and shows that these energy crops can grow just as productively on low productivity land as on high productivity land.

K E Y W O R D S
age, meta-analysis, meta-regression, miscanthus, nitrogen, switchgrass, yield 1 | INTRODUCTION Miscanthus and switchgrass are appealing feedstocks for biofuel production because, they are high yielding, can be grown productively in marginal growing conditions, have lower fertilization requirements, and have higher nitrogen (N) use efficiency compared to row crops (Christian et al., 2008;Fewell et al., 2016;Heaton et al., 2004;Jiang et al., 2018;Lewandowski et al., 2000;Stoof et al., 2015). There are also multiple environmental benefits of cultivating these crops including increased soil carbon, reduced nitrous oxide (N 2 O) emissions, and reduced nitrate (NO 3 ) leaching compared to row crops (Chamberlain & Miller, 2012;Davis et al., 2012;Hanssen et al., 2020;Hudiburg et al., 2015). These are perennial crops that have been shown to grow productively in rainfed conditions in the United States (US), but are yet to be grown at a commercial scale (Heaton et al., 2010;Lee et al., 2014;Parrish & Fike, 2005;Pyter et al., 2009). Though there is a large and growing experimental literature examining how these crops respond to various agronomic management practices, generalizable knowledge about the performance of these crops is still limited. Specifically, there is limited evidence on how these crops respond to N fertilization and whether this response varies linearly or not. There is also varying evidence of the effects of crop age on their productivity over time and its implications for their productive lifespan. More importantly, there is inconclusive evidence on whether yield declines with stand age because of depletion of soil nutrients and if N fertilization can replenish soil fertility and increase yields. Furthermore, systematic knowledge of the productivity response of these crops to differing locations, soil types, and climatic conditions is still lacking. As a result, it is challenging to develop any generalizable agronomic recommendations for managing these crops for optimal performance.
The purpose of this research is to quantify the responsiveness of yield to N and stand age and analyze the extent to which higher N can offset the effects of age. We also examine non-linearity in yield response to N and stand age while accounting for their interaction effect. We estimate a production relationship by applying a multi-level mixed-effects (MLME) meta-regression model. Our analysis controls for study and location-specific variation including differences in these across harvest years. Our meta-analytic approach also controls for cross-sectional (between locations) and temporal (between harvest years) dependencies among repeated observations from the same study. We undertake this analysis using data from field experiments conducted between 2002 and 2019 in the rainfed region of the US. These data include 1403 observations of miscanthus yield obtained from 27 locations in 14 states and 2811 observations of switchgrass yield obtained from 23 locations in 15 states.
Several studies have pooled data from multiple experiments to conduct a meta-analysis to examine the response of yield to N, while controlling for differences in stand age and diverse climatic conditions (Chen et al., 2019;Gunderson et al., 2008;Heaton et al., 2004;Jager et al., 2010;LeBauer et al., 2018;Miguez et al., 2008;Wang et al., 2010;Wullschleger et al., 2010). Some of these studies estimate a univariate response of yield to N fertilization (Chen et al., 2019;Gunderson et al., 2008), while others estimate multivariate functions that control for climate and other factors (Heaton et al., 2004;Jager et al., 2010;LeBauer et al., 2018;Miguez et al., 2008;Wang et al., 2010;Wullschleger et al., 2010). Many of these studies focus on estimating the directionality of the impact of N on crop yield, but not the magnitude of the impact or how it varies with stand age and changes in growing conditions (Heaton et al., 2004;Miguez et al., 2008;Wang et al., 2010). For example, Heaton et al. (2004) conducted a quantitative review of 21 articles representing 174 yield observations of 3 years and older miscanthus and switchgrass stands and reported a significant positive response to N. However, they use a mixed model analysis of variance (ANOVA) approach and did not control for other factors that can influence yield. Using an identical approach, Wang et al. (2010) report a significant positive response of N on switchgrass yield for monocultures compared to mixtures.
A few studies have examined the responsiveness of yield to N using meta-analysis, but do not control for stand age and various site-specific and study-specific factors (Chen et al., 2019;Gunderson et al., 2008;Jager et al., 2010;LeBauer et al., 2018;Wullschleger et al., 2010). For example, LeBauer et al. (2018) present a meta-analysis to estimate the effects of N, stand age, temperature, and precipitation on yields following Heaton et al. (2004) while including location and year within the location as random effects, following Wang et al. (2010). With a dataset of over 966 yield observations across 52 locations, their study confirms a positive response of yield to N for both crops. Similarly, Chen et al. (2019) constructed a dataset including 50 field studies across 65 locations, representing 919 observations with N rates ranging from 0 to 896 kg ha −1 up to 17 harvest years. Using linear and non-linear univariate regressions, they identify significant positive N effects on yields of both crops, which show a declining trend with the number of fertilization years. Their analysis shows the positive N effects on miscanthus yield increase with annual N rates of 60 to 300 kg ha −1 whereas switchgrass yield increases and peaks at N rates of 120 to 160 kg ha −1 for 5 to 6 fertilization years, but decreases for rates more than 300 kg ha −1 and beyond 7 fertilization years. These studies provide mixed evidence on the yield response of both energy crops to N-with some finding a statistically significant positive yield response (Heaton et al., 2004;LeBauer et al., 2018) and others finding a declining yield response after attaining peak yields (e.g., Chen et al., 2019). Specifically, LeBauer et al. (2018) find that a 1 kg ha −1 increase in N increases miscanthus and switchgrass yields by 0.028 and 0.013 Mg ha −1 , respectively. Investigating the impact of maturity using similar univariate analysis, Jager et al. (2010) report a decline in switchgrass yield after several years of harvest whereas LeBauer et al. (2018) report a linear decline in miscanthus and switchgrass yields following the third year of establishment, i.e., statistically significant reductions of 1.4 and 1.2 Mg ha −1 in miscanthus and switchgrass yields, respectively, for a 1-year increase in stand age. To investigate yield response to N and age interaction, Wang et al. (2010) report a significant positive response of N on switchgrass yield which is constant for newly established and matured stands but do not provide the magnitude of yield responsiveness. Although most of these studies examine directionality and/or magnitude of yield responsiveness to N and/or stand age, they do not account for possible non-linear effects including their interactions and do not control for site-specific and study-specific factors. Our study extends previous studies by estimating the magnitude of yield responsiveness to N and stand age including their non-linearities taking into account interaction effects while controlling for climate and location conditions and unobserved factors related to study design.
Beyond N fertilization and stand age, existing metaanalyses have examined the influence of environmental conditions on miscanthus and switchgrass yields. Although we expect all biological organisms to respond to temperature, previous studies conducting meta-analysis obtain mixed findings of this relationship. This could be due to differences in their construction of the variable that captures variation in the temperature as well as the other climate variables they include as additional controls and the extent of collinearity among these variables. For example, some studies have used mean seasonal temperature while others use mean annual temperature or seasonal GDD for capturing the impact of temperature on yields. Investigating the impacts of climatic factors, Heaton et al. (2004) show that miscanthus yield is strongly influenced by growing season precipitation, but not temperature (measured as growing degree days, GDD); Chen et al. (2019) do not find any significant relationship between mean annual temperature and miscanthus yield, but a significant positive impact of annual precipitation in their linear model. In contrast, LeBauer et al. (2018) report positive yield response to growing season precipitation and GDD immediately preceding harvest with miscanthus being more responsive to temperature than switchgrass. Studies have reported increased switchgrass yield with increased temperature to a point followed by a decrease whereas increased yield with increased growing season precipitation up to around 600 mm, which stabilizes thereafter (e.g., Gunderson et al., 2008;Wullschleger et al., 2010). Jager et al. (2010) also report the quadratic response of switchgrass yield to temperature along with positive yield response to increased precipitation and minimum winter temperature. Wang et al. (2010) report a significant positive response of switchgrass yield to precipitation, while GDD has no significant effect on the yield. In contrast, Chen et al. (2019) report a positive response of switchgrass yield to mean annual temperature but no significant relationship with annual precipitation. Like N and stand age, there is mixed evidence of the influence of climatic factors on miscanthus and switchgrass yields.
We consider the possibility of non-linearity in productivity response of miscanthus and switchgrass yields to N and stand age by including quadratic and interaction effects. However, generalizing yield responsiveness to N and stand age based on pooled data from multiple experiments requires caution because of differences in location and climatic features including variabilities due to experimental designs. We hypothesize that the inclusion of appropriate controls for location and climatic features as well as potential site-year-specific and study-year-specific unobserved variabilities will lead to more reliable estimates of yield responsiveness to N stand age. Our study goes beyond previous studies that have typically relied upon the mixed model ANOVA method to estimate the direction of yield response but do not estimate yield responsiveness to management and climate variables (Heaton et al., 2004;Wang et al., 2010). More importantly, our study improves upon those studies that estimate the magnitude of yield responsiveness to N and stand age, but do not account for their possible non-linear effects including N and age interactions or control for climatic variables and sitespecific features (Chen et al., 2019;Gunderson et al., 2008;Jager et al., 2010;LeBauer et al., 2018;Miguez et al., 2008;Wullschleger et al., 2010). We also add climatic factors as controls and investigate their linear and non-linear impacts on miscanthus and switchgrass yields.

| EMPIRICAL STRATEGY
The empirical analysis presented in this study examines the role of crop-specific factors, such as the age of the plant and genetics, environmental conditions, and management practices on the yield of miscanthus and switchgrass and allows for non-linear responses to these variables. We develop a generalized meta-regression model with yield as the dependent variable, and N and age as key explanatory variables including their quadratic terms, N × age and N × soil productivity interactions. We include temperature and precipitation including their quadratic terms and soil productivity as control variables. Our inclusion of soil productivity and N × soil productivity interaction is based on the experimental literature, which shows that the yield impact of N depends on the establishment of soil conditions Thomason et al., 2005;Zumpf et al., 2019). Since switchgrass has many cultivars and their yields differ considerably, we include genetic variability as the additional control variable. We include broadly classified regions (Great Plains, Midwest, Northeast, and Southeast) for controlling spatial variability in miscanthus yield at the regional level. Furthermore, we expect that there are unobserved factors pertaining to the design of the experiment, site location (latitude-longitude of the experimental site) and its management, and the conduct of the experimental work, which affect the yield. Since this information is not available, we use an empirical approach to control for this unobserved heterogeneity.
Our approach distinguishes between a deterministic effect related to the observed variables in the study and a random effect that is associated with the unobserved factors that could affect yield. We empirically estimate both the response and magnitude of responsiveness of miscanthus and switchgrass yields to N fertilization and stand age including possible non-linearities, while controlling for (i) growing season temperature (GDD) and precipitation including possible non-linearities, and cultivar and region characteristics (ii) harvest year, location, and study-specific unobserved heterogeneities. From the methodological perspective, we utilize the multi-level meta-regression, which has an advantage over standard regression that ignores hierarchy and dependency of data, i.e., yields measured at study, location, and harvest year levels and their potential correlations. The standard regression model estimated, for instance, by ordinary least squares (OLS), potentially results in biased estimates of standard errors thus leading to biased parameter estimates when hierarchy and dependency in data measured at different levels are ignored, i.e., the model estimated without multi-level random effects.
We hierarchically cluster the data by three different measurement levels, i.e., study, location, and year. Furthermore, each study and location have repeated observations for a specific year attributed to harvests from multiple replications (Figure 1). Since the data is clustered, there could be unobserved heterogeneity in yield, i.e., the conditional mean of the yield observations varies across studies, locations, and harvest years for reasons that are unobserved at those levels. Thus, we apply a MLME meta-regression model to estimate the crop production function.
where y itjk is the observed yield i from the crop harvested in year t at location j under study k.
The measure of yield is log-transformed to ensure nonnegative predictions, since the response variable only takes positive values and there is otherwise no way to prevent negative out-of-sample yield predictions for values of the explanatory variable outside the range of explanatory variables used in the meta-regression. The log-likelihood value for the log-transformed model is smaller than for the untransformed model-therefore, it is the preferred model. The estimated MLME meta-regression model in Equation (1) with log-transformed miscanthus and switchgrass yields results in the log-likelihood values of −1184 and −1149, respectively, compared to −4503 and −7934 from corresponding models, which have levels of yields as the response variable. The vectors of continuous F I G U R E 1 Hierarchical (multi-level) representation of observed yield by study, location, and year explanatory variables including their interactions are denoted by x itjk and corresponding parameters to be estimated are . The vectors of dummy variables for regions and corresponding parameters to be estimated (region fixed effects) are denoted by d r and r , respectively. The vectors of dummy variables for cultivars and corresponding parameters to be estimated (cultivar fixed effects) are denoted by d c and c , respectively. The heterogeneity in energy crop yield from different studies, i.e., mean and year-specific random intercepts for each study are captured by u k and r tk , respectively, whereas v jk and s tjk capture heterogeneity in the energy crop yield between locations from the same study, i.e., mean and year-specific random intercepts, respectively, for each location within the same study. The residual error term at the observation level is denoted by e itjk .
The specific form of x is as follows: where N is the annual nitrogen application rate, age is the number of years after planting at the time of harvest, GDD is growing degree days, precip is growing season cumulative precipitation, and NCCPI is National Commodity Crop Productivity Index (a measure of soil productivity). The total effect or observed yield is partitioned into deterministic effect and random effect. The deterministic effect represented by x , d r r , and d c c captures the effects of time-varying explanatory variables that include management and climate factors, time-invariant unobserved region, and cultivar fixed effects, respectively. The random effect or unexplained portion of the observed yield is a combination of factors that are unmeasured at harvest year, location, and study levels. A detailed explanation of the econometric procedure and interpretation is provided in SI Section A.
We estimate four different specifications for each of the energy crops, i.e., Model A includes deterministic effect without interaction terms; Model B includes deterministic effect with interaction terms (N 2 , age 2 , N × age, GDD 2 , precip 2 , N × NCCPI); Model C includes deterministic effect with interaction terms, study, and location random effect; Model D includes deterministic effect with interaction terms, harvest year augmented study, and location random effect. Model A ignores the possible unobserved heterogeneity in yield across studies, locations, and harvest years. Model B is the same as Model A, and it also includes interaction terms for N, age, GDD, precipitation, and NCCPI. In contrast, model C controls for the hierarchy and dependency of data across experiments at the same location or study and their potential correlations by including mean study and location random effects. Model D extends Model C by including harvest year-specific study and location random effects. All specifications are estimated using restricted maximum likelihood considering the degrees of freedom from the deterministic effect thus producing less biased estimates for random effects, i.e., variance components. The appropriate model is selected based on the lowest values for restricted log-likelihood (equivalently, deviance), Akaike's information criterion (AIC), and Bayesian information criterion (BIC).
There is some evidence that lower winter temperatures adversely affect miscanthus in the establishment year the severity of which varies over planting sites and times, but the matured stands of miscanthus generally tolerate cold winters (Dong et al., 2019;Maughan et al., 2012). As a sensitivity analysis, we also estimate a fully-fledged MLME meta-regression model, i.e., model D, with the inclusion of linear and quadratic terms of heating degree days (HDD) for the winter season in examining the impact of winter temperatures and the robustness of the findings.

| Variable construction
We compile field trial data spanning the past two decades from 10 published articles, 1 published dissertation, and 1 published dataset in the US (SI Section B). For each plot level observation, we obtain dry biomass yield (Mg ha −1 ), location, latitude, longitude, planting year, growing year, harvest year, nitrogen application rate (kg ha −1 ), and cultivar type (switchgrass). We extract relevant data from texts and tables of the publications or corresponding datasets obtained through personal communications with the authors. For switchgrass, we include dummy variables to distinguish the effects of three cultivars-Alamo, Blackwell, and Cave-In-Rock and group the other cultivars for which we have a few observations into the Other category (a composite category including Ceres 1102, Ceres EG1101, Liberty, Shawnee, SL 93 2001-1, Summer, Sunburst, and Timber) as the reference. We aggregate these cultivars since the number of observations for each of these cultivars-Ceres 1102, Ceres EG1101, Liberty, Shawnee, SL 93 2001-1, Summer, Sunburst, and Timber-is statistically too small to include them as independent categories. For Alamo, Blackwell, Cave-In-Rock, and Other, we have 2208, 84, 352, and 167 yield observations, respectively. For miscanthus, we include three regions-Great plains (NE, OK, SD), Midwest (IL, IA, MI, WI), and Southeast (AR, (2) GA, KY, LA, MS, TN, VA) with Northeast (NJ, NY) as the reference category. We compute stand age for both crops as the difference between the growing year (the same year as harvest year or 1 year before the harvest year) and planting year.
To evaluate soil characteristics in terms of productivity, weighted National Commodity Crop Productivity Index (NCCPI) values available at 150 m grid cells are obtained from the gridded Soil Survey Geographic (gSSURGO 2020) Database. The NCCPI is a rating index capturing soil capacity to produce commodity crops in rainfed (nonirrigated) conditions. We use the weighted average of NCCPI values for major commodity crops including corn, soybean, and cotton. These values are extracted for each location based on latitude-longitude using ArcGIS. Climate information is obtained from the parameter-elevation regressions on independent slopes model (PRISM 2020) dataset at each latitude-longitude corresponding to each location. The climate-related variables are calculated for the growing season (April-November for miscanthus and April-September for switchgrass). We compute growing season cumulative precipitation (precip) using daily precipitations (mm). We calculate GDD (°C) in a growing season to measure the thermal time using the following relationship (Arundale et al., 2014;: where d is a particular day with n denoting number of days in a growing season, T max,d and T min,d are maximum and minimum temperatures in any day d, respectively, and b is the base temperature below which growth is supposedly restricted. We use 6°C as the base temperature for miscanthus  and 10°C for switchgrass (Arundale et al., 2014). Similarly, we calculate HDD (°C) for the winter season (December-February preceding the harvest), which defines the cumulative winter temperature below the base temperature as follows:
The differences in the mean of the yield observations across sites for miscanthus ( Figure 2) and switchgrass ( Figure 3) reflect substantial spatial variation in the yield of these crops. In general, lower latitudes have lower miscanthus yields, whereas higher latitudes have lower switchgrass yields on average. The middle latitude has a higher average yield pattern for both crops. Table 1 describes the variables used in the estimated models and their summary statistics. The miscanthus dataset includes a wide range of N rates from 0 to 202 kg ha −1 with a mean and standard deviation of 69.1 and 59.5 kg ha −1 , respectively, across stand ages of 2-11 years with a mean age and standard deviation of the age of 3.9 and 1.9 years, respectively. We exclude miscanthus observations using fertilizer application rates of 224-448 kg N ha −1 . These higher N rate yield observations are from a single staggered start experiment  and exceed the range of N rates consistent across all other studies in the miscanthus dataset. The seasonal precipitation and GDD for the miscanthus dataset ranges from 388.6 to 1401.8 mm and from 2034.2 to 4434.7°C, respectively. The mean and standard deviation for precipitation are 852.7 and 202.9 mm, respectively, whereas the mean and standard deviation for GDD are 2786.9 and 366.9°C, respectively. For the switchgrass dataset, N rates range from 0 to 202 kg ha −1 with a mean and standard deviation of 89.2 and 72.3 kg ha −1 , respectively, across stand ages of 2-10 years with a mean and standard deviation of 4.8 and 1.9 years, respectively. The seasonal precipitation and GDD for the switchgrass dataset range from 314.9 to 1180.3 mm and from 1060.6 to 2984.5°C, respectively. The mean and standard deviation for precipitation are 709.1 and 172.2 mm, respectively, whereas the mean and standard deviation for GDD are 2207.9 and 352.8°C, respectively. The winter season HDD for the miscanthus and switchgrass dataset ranges from 80.3 to 1614.3°C and from 178.4 to 2138°C, respectively. The mean and standard deviation for HDD in the miscanthus dataset are 810.3 and 317.2°C, respectively, whereas the mean and standard deviation for HDD in the switchgrass dataset are 738.8 and 348.1°C, respectively.

| RESULTS AND DISCUSSIONS
Relevant model fitness and test statistics point toward the appropriateness of the multi-level mixed-effects model D with study × year and location × year random effect (see SI Sections C and D for miscanthus and switchgrass, respectively). We report parameter estimates and relevant statistics of four different specifications for miscanthus and switchgrass in Tables 2 and 3, respectively. We also examine the impact of winter temperatures as a sensitivity analysis by estimating a fully fledged MLME meta-regression model, i.e., model D, with the inclusion of linear and quadratic terms of winter season HDD (see parameter estimates in SI Section E). We neither find statistical significance of the included linear and quadratic terms for HDD nor any improvement in the model fitness.

| Production function estimation and yield prediction for miscanthus
The deterministic and predicted (total) yields from model D as well as observed yields are shown in Figure 4. In general, the deterministic portion of yield increases with age but at a decreasing rate. The deterministic yield peaks at 7th year and then again spikes at 9th year before declining. The deterministic yields in the 5th, 7th, and 9th years are 17.4, 19.2, and 20.4 Mg ha −1 , respectively, which are more than 25% to 75% higher than in the 2nd year (11.5 Mg ha −1 ). We find that the random factors affecting yields in the miscanthus dataset are more important in middle and later years than in earlier years, that is, in later years, the random effects explain a larger portion of the observed yield compared to the deterministic factors. Observed yields in the younger stages of maturity are identical to those

F I G U R E 3 Mean observed switchgrass yield by location
Note: Numbers in parentheses denote actual yield observations.

Variable Definition
Miscanthus ( expected based on stand age alone, but different in the middle and later stages of maturity than those expected based on stand age alone due to random factors. For example, deterministic effects explain 97% and 94% of the predicted yields in the 2nd and 3rd years, respectively. In contrast, deterministic effects in the 4th and 5th years explain 69% and 90% of the predicted yields, respectively. Deterministic effects in the 8th and 10th years are 86% and 133% of the predicted yields, respectively. Our predicted yields across maturity stages are in close agreement with the observed yields. The mean predicted yield peaks at the 4th year, remains constant until the 9th year, and declines thereafter following a similar pattern as the observed yield. The mean predicted yield at the 4th year is 21.3 Mg ha −1 (range of 9.9-32.1 Mg ha −1 ), which is close to the mean observed yield of 21.8 Mg ha −1 (range of 6.4-44.4 Mg ha −1 ). As the plant grows older, the range of the observed yield becomes smaller. Our prediction also indicates a similar pattern for the older stands. For example, the predicted yield falls between 12.6 and 30.6 Mg ha −1 (mean of 20.2 Mg ha −1 ) in the 9th year, which is close to the observed yield between 13.5 and 29.7 Mg ha −1 (mean of 20.5 Mg ha −1 ).
The estimated MLME meta-regression (model D) results in the lowest AIC and BIC values of 2405 and 2499, respectively, compared to 2903 and 2987 in standard regression, which ignores hierarchy and dependency of data (i.e., OLS in model B). The mean OLS yield at the 4th year is 18.4 Mg ha −1 with a range of 13.3-25.4 Mg ha −1 . Similarly, the mean OLS yield falls between 11.3 and 39.4 Mg ha −1 (mean of 18.9 Mg ha −1 ) in the 9th year ( Figure 4). These values deviate from the observed yields as opposed to our predicted yields. This shows potential site-year-specific and study-year-specific unobserved variabilities in a MLME model lead to better yield prediction and thus more reliable estimates of N and stand age compared to OLS. The mean random effect for each location from model D is shown in Figure 5. It represents the mean of the yield variation attributed to a particular harvest year within that particular location for a given experiment (study). This effect is an important component of the overall yield prediction, even though it is not explained by observable factors. It varies across regions and ages. For example, the mean estimated location random effect in IA ranges from 1.4 to 2.5 Mg ha −1 , whereas it ranges from 0.63 to 1.3 Mg ha −1 in IL.

| Response of miscanthus yield to age and nitrogen
The linear impact of age on the productivity of miscanthus is positive and statistically significant at 1% level. However, miscanthus has a diminishing marginal productivity response to age as shown by the negative coefficient estimate on the squared term, which is also statistically significant at 1% level. The coefficient of the N × age interaction has a positive sign and is statistically significant at 10% level, which shows that the diminishing marginal productivity to age can be partially offset by increasing N. Accounting for all the interaction effects, the marginal impact of age is statistically significant at 1% level; i.e., an additional year of maturity increases miscanthus yield on average by 13.9%. The yield increases at a decreasing rate till the 7th year, beyond which, yield starts decreasing at an increasing rate. For instance, for a 2nd, 6th, and 10th-year miscanthus stand, an additional year of maturity changes yield by on average by 2.6, 0.47, and −2.5 Mg ha −1 , respectively.

| Response of miscanthus yield to nitrogen and soil productivity (NCCPI)
The linear impact of N on the productivity of miscanthus is not statistically significant but miscanthus has a F I G U R E 4 Mean observed, OLS, deterministic, and predicted miscanthus yields Note: OLS yield is from standard pooled regression without multi-level random effects (model B); deterministic and predicted yields are from MLME meta-regression (model D); predicted yield includes deterministic yield and multi-level random effects. The figure compares the deterministic as well as (total) predicted yields using the MLME meta-regression (model D) with the predicted yields from the standard regression, which ignores hierarchy and dependency of data (i.e., OLS in model B). It shows the impact of the random factors in the observed yields for different stages of maturity, and how the potential site-year-specific and study-yearspecific (random) unobserved variabilities in a MLME model lead to better yield predictions as compared to the OLS model. diminishing marginal productivity response to N as depicted by the negative coefficient of the squared term, which is statistically significant at 5% level. The coefficients of NCCPI and N × NCCPI are statistically not significant. Accounting for the interaction effects, the marginal impact of N on miscanthus yield is statistically significant at 1% level, i.e., an additional 100 kg N ha −1 increases miscanthus yield on average by 11.2%. The yield increases at a decreasing rate till 140 kg N ha −1 , beyond which, yield starts decreasing at an increasing rate. For instance, considering N application rates of 0, 100, and 200 kg ha −1 , an additional kg of N ha −1 changes miscanthus yield on average by 0.03, 0.01, and −0.02 Mg ha −1 , respectively. The overall marginal impact of NCCPI is statistically not significant.

| Response of miscanthus yield to temperature and precipitation
The linear impact of growing season GDD on the productivity of miscanthus is positive and statistically significant at 1% level. However, miscanthus has a diminishing marginal productivity response to growing season GDD as shown by the negative coefficient estimate on the squared term, which is statistically significant at 1% level. Overall, the marginal impact of growing season GDD on productivity is statistically significant at 1% level; i.e., an additional 100°C of growing season GDD increases miscanthus yield on average by 5.4%. In contrast to growing season GDD, there is no statistically significant yield response to growing season cumulative precipitation.

| Production function estimation and yield prediction for switchgrass
We eliminate regional dummies in estimating switchgrass production function to avoid multicollinearity with cultivar dummies. For example, the cultivar Alamo has 2208 (out of 2811) observations in total from a total of three studies (Arundale, 2012;Boyer et al., 2013;Lee et al., 2018), which are all conducted in the Southeast region comprising 2220 (out of 2811) observations. In general, observed data follows regional growth patterns across individual switchgrass cultivars; i.e., Alamo as a lowland cultivar grown in the Southeast, Blackwell as an upland cultivar grown in the Great Plains, and Cave-In-Rock as an upland cultivar grown in the Midwest and Northeast. This is consistent with previous literature suggesting the

F I G U R E 5 Mean estimated miscanthus random effect (RE) by location and region
Note: The figure shows the expected impact of the random factors at the location level, i.e., the site-year-specific and study-year-specific (random) unobserved variabilities, in predicting the total (observed) yields. suitability of different growing regions for different cultivars (e.g., Parrish & Fike, 2005). The deterministic and predicted (total) yields from model D as well as observed yields are shown in Figure 6. In general, deterministic yield increases with age but at a decreasing rate. The deterministic yield peaks at the 6th year and then starts declining with a sharp drop at the 8th year. The deterministic yields in the 3rd, 5th, and 7th years are 8.8, 12.9, and 12.8 Mg ha −1 , respectively, which are more than 25% to 75% higher than in the 2nd year (5.7 Mg ha −1 ). We find that random factors affecting yield in the switchgrass dataset are more important in earlier and later years than in the middle years. Observed yields in the middle stages of maturity are identical to those expected based on stand age alone but different in the earlier and later stages of maturity than those expected based on stand age alone due to random factors. For example, deterministic effects in the 4th and 5th years are 105% and 100% of the predicted yields, respectively. In contrast, deterministic effects explain 64% and 62% of the predicted yields in the 2nd and 3rd years, respectively. Similarly, deterministic effects explain only 49% and 40% of the predicted yields in the 8th and 10th years, respectively.
Our predicted yields across maturity stages closely replicate the pattern of the observed yields. The mean predicted yield peaks at the 3rd year and remains constant until the 8th year, declining thereafter. The mean predicted yield at the 3rd year is 14.3 Mg ha −1 (range of 1.8-20.1 Mg ha −1 ) which is close to the mean observed yield of 15.4 Mg ha −1 (range of 1.8-39.7 Mg ha −1 ). Similar to the observed data in which the range of observed yield becomes smaller for older stands our prediction shows a similar pattern. For example, the predicted yield has a range of 2.8-20.1 Mg ha −1 (mean of 14.2 Mg ha −1 ) at the 8th year, which is close to the range of 2.5-27.4 Mg ha −1 for the observed yield (mean of 14.3 Mg ha −1 ).
The estimated MLME meta-regression (model D) results in the lowest AIC and BIC values of 2334 and 2441, respectively, compared to 3670 and 3765 in simple regression, which pools data ignoring hierarchy (i.e., OLS in model B). In contrast to our predicted yields, which are close to the observed yields, the mean OLS yield at the 3rd year is 10.8 Mg ha −1 with a range of 3.5-14.6 Mg ha −1 . Similarly, the mean OLS yield has a range of 4.1-16.5 Mg ha −1 (mean of 12.2 Mg ha −1 ) at the 8th year ( Figure 6). This again shows that the inclusion of potential site-year-specific and study-year-specific unobserved variabilities in a MLME model leads to better yield prediction and thus more reliable estimates of N and stand age compared to OLS. The mean random effect for each location from model D is shown in Figure 7. This effect is an important component of the overall yield prediction, even though it is not observed. For example, the mean estimated location random effect in SD ranges from 1.1 to 2.4 Mg ha −1 , whereas it ranges from 0.59 to 1.7 Mg ha −1 in IL.

| Response of switchgrass yield to age and nitrogen
The linear impact of age on the productivity of switchgrass is positive and statistically significant at 1% level. However, switchgrass has a diminishing marginal productivity response to age as shown by the negative coefficient estimate on the squared term, which is also statistically significant at 1% level. The N × age interaction shows that the diminishing marginal productivity to age can be partially offset by increasing N, which however is statistically not significant. Accounting for all the interaction effects, the marginal impact of age on switchgrass yield is statistically significant at 1% level; i.e., an additional year of maturity on average increases switchgrass yield by 11.7%. The yield increases at a decreasing rate till the 6th year, beyond which, it starts decreasing at an increasing rate. For instance, considering F I G U R E 6 Mean observed, OLS, deterministic, and predicted switchgrass yields Note: OLS yield is from standard pooled regression without multi-level random effects (model B); deterministic and predicted yields are from MLME meta-regression (model D); predicted yield includes deterministic yield and multi-level random effects. The figure compares the deterministic as well as (total) predicted yields using the MLME meta-regression (model D) with the predicted yields from the standard regression, which ignores hierarchy and dependency of data (i.e., OLS in model B). It shows the impact of the random factors in the observed yields for different stages of maturity, and how the potential site-year-specific and study-yearspecific (random) unobserved variabilities in a MLME model lead to better yield predictions as compared to the OLS model. 2nd, 6th, and 10th-year switchgrass stands, an additional year of maturity changes yield on average by 2.9, −0.97, and −2.1 Mg ha −1 , respectively.

| Response of switchgrass yield to nitrogen and soil productivity (NCCPI)
The linear impact of N on the productivity of switchgrass is positive and statistically significant at 1% level, but switchgrass has a diminishing marginal productivity response to N, as shown by the statistically significant negative coefficient estimate on the squared term. Even though the NCCPI impact is statistically not significant, the N × NCCPI interaction effect shows that there is a negative impact of adding N on land with higher NCCPI, which is statistically significant. Accounting for the interaction effects, the marginal impact of N on switchgrass yield is statistically significant at 1% level; i.e., an additional 100 kg N ha −1 on average increases switchgrass yield by 32.5%. The yield increases at a decreasing rate till 155 kg N ha −1 , beyond which, it starts decreasing at an increasing rate. For instance, considering N application rates of 0, 100, and 200 kg ha −1 , an additional kg of N ha −1 changes switchgrass yield on average by 0.05, 0.03, and −0.02 Mg ha −1 , respectively. The overall marginal impact of NCCPI is statistically not significant.

| Response of switchgrass yield to temperature and precipitation
The linear impact of growing season GDD on the productivity of switchgrass is positive and statistically significant at 1% level. However, switchgrass has a diminishing marginal productivity response to growing season GDD, as shown by the negative coefficient estimate on the squared term, which is statistically significant at 10% level. The linear impact of growing season cumulative precipitation on the productivity of switchgrass is positive and statistically significant at 1% level. Similar to growing season GDD, switchgrass has a diminishing marginal productivity response to growing season cumulative precipitation, as shown by the negative coefficient estimate on the squared term, which is statistically significant at 1% level. Overall, the marginal impact of growing season GDD is significant at 1% level; i.e., an additional 100°C of growing season GDD on average increases switchgrass yield by 8.3%. Similar to growing season GDD, there is an overall positive yield response to growing season cumulative F I G U R E 7 Mean estimated switchgrass random effect (RE) by location and region Note: The figure shows the expected impact of the random factors at the location level, i.e., the site-year-specific and study-year-specific (random) unobserved variabilities, in predicting the total (observed) yields. precipitation with a yield increase of 6.8% on average for an additional 100 mm of growing season cumulative precipitation. This effect is statistically significant at 5% level.
Since we are interested in yield effects due to climate, management, cultivar, and region that we can observe and explain, we provide a detailed discussion of the deterministic portion of the predicted yields in subsections that follow.

| Predicted deterministic yields by age and nitrogen interaction
To better understand the N and age interaction effect, we provide predicted yields at various N rates for selected plant maturity levels assuming other explanatory variables (controls) remain at their observed levels across the differing N rates. Using predicted values, we compute magnitudes of changes in miscanthus and switchgrass yields for various N levels across maturity stages. In general, both the crops at various stages of maturity show an increased response to N up to a specific level and then follow a declining trend with a further increase in N rate.
We predict that the deterministic portion of the yield of both crops varies over N rates with differences across younger (1-3 years), middle-aged (4-7 years), and older (8-10 years) stands (Figure 8). We find that the positive yield response of miscanthus to N is higher for the older stands compared to younger and middle-aged stands. In contrast, we find that the positive yield response of switchgrass to N is higher for the middle-aged compared to younger and older stands. For example, given yields of 10.8, 16.5, and 11.1 Mg ha −1 in the 2nd, 6th, and 10th-year miscanthus stands with 0 kg N ha −1 , respectively, the yield increases by 2.1 Mg ha −1 in a 10th year stand compared to a yield increase of 0.7 and 2.1 Mg ha −1 in the 2nd and 6th year stands, respectively, when the N rate increases from 0 to 50 kg ha −1 .
In contrast, given yields of 3.8, 9.6, and 2.1 Mg ha −1 in the 2nd, 6th, and 10th-year switchgrass stands with 0 kg N ha −1 , respectively, the yield increased by 3.1 Mg ha −1 in the 6th year stand compared to yield increase of 1.2 and 0.7 Mg ha −1 in the 2nd and 10th year stands, respectively, when the N rate is increased from 0 to 50 kg ha −1 (Figure 9, see the details in SI Section F).
The predicted deterministic yield of the 6th year miscanthus stand without N application (16.5 Mg ha −1 ) and 10th year stand with N rate of 150 kg ha −1 (16.6 Mg ha −1 ) are similar; the latter would have otherwise been 11.1 Mg ha −1 without N application. Similarly, the predicted yield of the 6th year switchgrass stand without N application (9.6 Mg ha −1 ) and 8th year stand with N rate of 100 kg ha −1 (9.7 Mg ha −1 ) are similar, which could have been 6.1 Mg ha −1 without N application (SI Section F). These results show that the decline in productivity in older stands can be offset with additional N.
Yields of various age miscanthus stands reach a maximum level with N rates between 100 and 250 kg N ha −1 (Figure 9). In contrast, yields of all age switchgrass stands reach a maximum at 150 kg N ha −1 . We find that the mean predicted yields of the younger and middleaged miscanthus stands decrease when the N rate is increased from 150 to 200 kg ha −1 . In contrast, we find that when the N rate is increased from 150 to 200 kg ha −1 , switchgrass shows yield decline for all age stands. These findings show the yield declining impact of excessive F I G U R E 8 Mean predicted deterministic yield across maturity stages at various N levels for miscanthus and switchgrass N on miscanthus and switchgrass irrespective of plant maturity but the predicted yields are still greater than without N fertilization.
To determine differential impact at lower N levels across maturity stages, we further investigate the change in mean predicted miscanthus and switchgrass yields between 0 and 50 kg ha −1 at an interval of 10 kg N ha −1 (SI Section G). We find variation in yield response across younger, middleaged, and older stands with both crops at various stages of maturity showing an increased response with an increase in N rate till 50 kg ha −1 . Even with lower N levels, we find a higher positive yield response for the older compared to younger and middle-aged stands in miscanthus, whereas a higher positive yield response for the middleaged compared to younger and older stands in switchgrass. Compared to LeBauer et al. (2018), who find an average yield increase of 0.28 Mg ha −1 per 10 kg ha −1 increase in N, we find that miscanthus yield increases by 0.14, 0.42, and 0.43 Mg ha −1 on average per 10 kg ha −1 increase in N rates between 0 and 50 kg ha −1 across 2nd, 6th, and 10th year stands, respectively. Similarly, we find that switchgrass yield increases by 0.25, 0.62, and 0.14 Mg ha −1 on average per 10 kg ha −1 increase in N rates between 0 and 50 kg ha −1 across 2nd, 6th, and 10th year stands, respectively, compared to LeBauer et al. (2018), who find an average yield increase of 0.13 Mg ha −1 per 10 kg ha −1 increase in N.

| Predicted deterministic yields across region and cultivar
We provide predicted yields at various maturity stages across regions assuming other explanatory variables (controls) remain at their observed levels across the differing ages (Figure 10). In general, the mean predicted miscanthus yields are the highest in the Great Plains followed by Midwest, Southeast, and Northeast, across all maturity stages with maximum yield achieved around 6-7 years. In contrast, the mean predicted switchgrass yields are the highest in the Southeast followed by Midwest, Great Plains, and Northeast, across all maturity stages with maximum yield achieved around 5-6 years.
We also provide predicted yields at various maturity stages for the switchgrass cultivars, assuming all other explanatory variables (controls) remain at their observed levels across the differing ages ( Figure 11). In general, the mean predicted yields are the highest for Alamo followed by Other, Cave-In-Rock, and Blackwell across all maturity stages with maximum yield achieved around 5-6 years.
Our findings on predicted yields across combinations of N rate and stand age should be interpreted with caution. First, both miscanthus and switchgrass respond positively to N with subsequently smaller increases in the yields with additional increases in N rate. After attaining peak yield, both crops respond negatively to N with subsequently larger decreases in the yields with additional increases in N rate. Second, yield increases at higher rates in the older compared to younger and middle-aged stands in miscanthus. In contrast, yield increases at higher rates in the middle-aged compared to younger and older stands in switchgrass. There is a positive but small impact of N on various age (2nd to 10th year) miscanthus stands, ranging from 0.17 to 0.50 Mg ha −1 and switchgrass stands between 0.14 to 0.64 Mg ha −1 when N rate is increased from 0 to 10 kg ha −1 . This positive impact of N depreciates quickly and becomes negligible as N is further increased. For example, the impact of N on various F I G U R E 9 Mean predicted deterministic yield across N levels at various maturity stages for miscanthus and switchgrass age (2nd to 10th year) miscanthus stands ranges from 0.11 to 0.45 Mg ha −1 and switchgrass stands between 0.13 to 0.59 Mg ha −1 when N rate is increased from 40 to 50 kg ha −1 . We also find that miscanthus stands reach maximum predicted yields ranging between 11.8 and 20.8 Mg ha −1 depending on crop maturity (2nd to 10th year) for N levels between 100 and 250 kg ha −1 , whereas switchgrass stands reach a maximum at 150 kg N ha −1 irrespective of crop maturity with yields ranging between 3.6 and 16.1 Mg ha −1 depending on crop maturity.

| CONCLUSIONS
Our meta-regression analysis applying multi-level mixed-effects modeling extends our understanding of the impacts of N and stand age interaction on productivity of miscanthus and switchgrass. We find that yield of both of these crops increases at declining rates as they mature till a peak, where the yield remains stagnant for a few years and then starts to decline. Both crops respond positively to N fertilization for the observed range of applied rates with the effect being more pronounced across older stands in miscanthus and middle-aged stands in switchgrass. However, the positive response to N fertilization is small in magnitude and quickly deteriorates. The crops at various stages of maturity show an increasing but negligible response to N fertilization up to a specific level of N and then a decline in yield with further increases in N rate. Even though both energy crops exhibit a positive yield response to growing season GDD, switchgrass has a higher yield response than miscanthus. Switchgrass also responds positively to growing season cumulative precipitation but this is not the case for miscanthus. We find that miscanthus stands reach maximum yields for N levels between 100 and 250 kg ha −1 depending on crop age, whereas switchgrass stands reach a maximum yield at 150 kg N ha −1 irrespective of crop age. Our findings also provide statistical evidence that yields of miscanthus and switchgrass do not differ across land with different biophysical properties and that they can be grown productively on low quality land.
Our parameter estimates can be used to lead to recommendations for optimal N application rates as well as the optimal age to harvest miscanthus and switchgrass that can maximize yield. They can also be used to show that N rates that maximize yield may not necessarily be the rates that will maximize profits, since these will depend on market factors such as biomass price and N price. We leave the application of our findings here to determine F I G U R E 1 0 Mean predicted deterministic yield across maturity stages and regions for miscanthus and switchgrass F I G U R E 1 1 Mean predicted deterministic switchgrass yield across maturity stages and cultivars recommendations for profit-maximizing N rates across maturity stages and market prices for future research. Nitrogen use is of much significance since excessive application not only leads to economic inefficiency but becomes counterproductive by releasing N 2 O emissions and NO 3 leaching which have well-known greenhouse gas (GHG) implications (Cadoux et al., 2012;Chen et al., 2019). The findings of our analysis can improve both the economic and environmental outcomes from lignocellulosic biofuel production in terms of farm profitability and GHG reduction. The estimates obtained here can also guide the validation of biophysical (mechanistic) crop growth models and inform regional assessments of bioenergy suitability and profitability relative to conventional crops, which can affect incentives for land-use change.