ORCHIDEE-MICT-BIOENERGY: an attempt to represent the production of lignocellulosic crops for bioenergy in a global vegetation model

Bioenergy crop cultivation for lignocellulosic biomass is increasingly important for future climate mitigation, and it is assumed on large scales in integrated assessment models (IAMs) that develop future land use change scenarios consistent with the dual constraint of sufficient food production and deep decarbonization for low climate-warming targets. In most global vegetation models, there is no specific representation of crops producing lignocellulosic biomass, resulting in simulation biases of biomass yields and other carbon outputs, and in turn of future bioenergy production. Here, we introduced four new plant functional types (PFTs) to represent four major lignocellulosic bioenergy crops, eucalypt, poplar and willow, Miscanthus, and switchgrass, in the global process-based vegetation model ORCHIDEE. New parameterizations of photosynthesis, carbon allocation, and phenology are proposed based on a compilation of field measurements. A specific harvest module is further added to the model to simulate the rotation of bioenergy tree PFTs based on their age dynamics. The resulting ORCHIDEE-MICT-BIOENERGY model is applied at 296 locations where field measurements of harvested biomass are available for different bioenergy crops. The new model can generally reproduce the global bioenergy crop yield observations. Biases in the model results related to grid-based simulations versus the point-scale measurements and the lack of fertilization and fertilization management practices in the model are discussed. This study sheds light on the importance of properly representing bioenergy crops for simulating their yields. The parameterizations of bioenergy crops presented here are generic enough to be applicable in other global vegetation models.

1 Introduction

Biomass-derived fuels serve as an alternative energy source to substitute fossil fuel and are used by many countries to meet renewable energy and climate targets (Karp and Shield, 2008; Meier et al., 2015; Robertson et al., 2017). Expanding bioenergy crop plantation is considered in future scenarios for energy security and climate change mitigation (Karp and Shield, 2008; Smith et al., 2016; Robertson et al., 2017). For bioenergy production to provide economic and climate benefits, cultivated plants must have a high productivity and a high yield of harvestable biomass (Karp and Shield, 2008; Whitaker et al., 2010; Robertson et al., 2017). The first generation of bioenergy crops usually refers to grain and high-sugar crops like maize and sugarcane (Karp and Shield, 2008). These crops have high nutrient requirements that demand fertilizer additions, causing high N₂O emissions to the atmosphere to achieve a high productivity (Searchinger et al., 2008; Melillo et al., 2009). These grain and high-sugar crops are unlikely to be planted on a large scale for the purpose of bioenergy production because of the food demand pressure for fertile land and fertilizer (Alexandratos and Bruinsma, 2012; Gerland et al., 2014; United Nations, 2017). Compared to the first generation, the second generation of bioenergy crops, known as lignocellulosic energy crops like giant miscanthus, switchgrass, and short-rotation trees, are adapted to a wider range of climatic and soil conditions and require less nitrogen fertilizer (Miguez et al., 2008; Cadoux et al., 2012). Those second-generation bioenergy crops have potentials to be deployed on marginal lands to avoid direct and indirect land use change (LUC) carbon emissions and damage to ecosystem services (Robertson et al., 2017). They also appear to have less greenhouse gas (GHG) emissions and higher energy efficiency than the first-generation bioenergy crops (Whitaker et al., 2010).

Bioenergy with carbon capture and storage (BECCS) is the main class of future negative emission technologies expected to result in net removal of atmospheric CO₂ (Smith et al., 2016). BECCS has been extensively assumed in integrated assessment models (IAMs) to develop land-based mitigation scenarios for low warming levels (Fuss et al., 2014; Popp et al., 2014). In most IAMs like IMAGE (Bouwman et al., 2006; Stehfest et al., 2014) and MAgPIE (Popp et al., 2011; Klein et al., 2014), second-generation bioenergy crops are used as primary energy carriers (Popp et al., 2014). One output from IAMs is future land use maps based on different environmental, socioeconomic, and policy constraints. These land use maps, after being translated into plant functional type (PFT) maps, can be used in grid-based dynamic global vegetation models (DGVMs) to simulate the terrestrial carbon dynamics, biogeochemical (e.g., LUC carbon emissions), and biophysical (e.g., albedo and transpiration changes) effects of land use processes (Brovkin et al., 2013; Wilkenskjeld et al., 2014). Global vegetation models can provide in return to IAMs some valuable information like spatially explicit biomass density, crop yield, and water availability (Stehfest et al., 2014; Bonsch et al., 2015, 2016). For example, dedicated bioenergy crop modeling has been implemented in a global vegetation model (LPJml; Beringer et al., 2011; Heck et al., 2016) to simulate biophysical yields and water availability as input data for MAgPIE (Bonsch et al., 2016).

In most global grid-based vegetation models, there are no dedicated PFTs to represent second-generation bioenergy crops. Instead, these plants are often represented by a generic crop PFT. Biases in simulated biomass production and the resulting carbon and energy balance thus arise when ignoring differences in carbon assimilation and phenology between generic crops and lignocellulosic bioenergy crops. Moreover, lignocellulosic woody bioenergy crops like eucalypt, poplar, and willow cannot be properly represented by a herbaceous crop PFT. For example, eucalypt has a high maximum rate of carboxylation (V_cmax) but relatively low leaf area index (LAI; Stape et al., 2004; Whitehead and Beadle, 2004). Miscanthus, on the contrary, has a relatively lower V_cmax (Wang et al., 2012; Yan et al., 2015) but a higher LAI (Heaton et al., 2008; Zub and Brancourt-Hulmel, 2010) than eucalypt (Whitehead and Beadle, 2004). Even if both Miscanthus and switchgrass are C₄ crops, Miscanthus can achieve a significantly higher yield than switchgrass because of a higher efficiency of converting intercepted radiation into aboveground biomass than switchgrass (Heaton et al., 2008). The water, nitrogen and light use efficiencies are also higher for Miscanthus than for switchgrass, resulting in a higher rate of leaf photosynthesis in the former (Dohleman et al., 2009). All these important differences between lignocellulosic bioenergy crops need to be considered, which calls for a dedicated new model PFT for each species.

Similarly, the way that harvest is implemented for generic crops in global models (usually removing a fixed fraction of biomass, typically on the order of 50 %) cannot be used for bioenergy crops. The harvest index (HI, harvested biomass as a fraction of aboveground biomass) is very different for grain crops and herbaceous bioenergy crops. In addition, most vegetation models currently do not account for realistic rotations of ligneous bioenergy plants (e.g., poplar and eucalypt). Modeling the harvest of woody bioenergy crops should be based on rotation practices of typically a few years rather than on assuming annual full harvest like for herbaceous crops. This requires simulating forest age dynamics (Yue et al., 2018) to accurately represent the ligneous biomass harvest.

In this study, we aim to model the biomass yields of four major lignocellulosic bioenergy crops in the global dynamic vegetation model ORCHIDEE. We introduce the new bioenergy crop PFTs, adjust the parameters relevant to the physiology, phenology, and harvest process of bioenergy crops based on observations, and evaluate the simulated biomass yields using a new global dataset of field measurements.

Table 1 Plant functional types (PFTs) in ORCHIDEE. The newly added bioenergy PFTs (PFT14 to PFT17) use the default setting of the original PFTs (all processes except harvest, see Sect. 2.2) but with new parameterizations (see Sect. 2.3).

2 Model development and parameterization

2.1 Model description

The proposed parameterizations of lignocellulosic bioenergy crops are based on an extended version of ORCHIDEE (Krinner et al., 2005), ORCHIDEE-MICT (Guimberteau et al., 2018), which contains relevant features of gross land use change, wood harvest, and forest age classes dynamics (Yue et al., 2018; Fig. S1). The model simulates energy exchange, water balance, and vegetation carbon processes in the ecosystem and is the land surface component of the French Earth system model (ESM) IPSL-CM (Krinner et al., 2005). The principal processes related to carbon cycling comprise photosynthesis, vegetation carbon allocation, autotrophic and heterotrophic respiration, plant phenology (e.g., leaf onset and senescence), and litter and soil carbon dynamics (Krinner et al., 2005). ORCHIDEE-MICT further includes high-latitude-related processes with new parameterizations of soil carbon vertical discretization, snow processes, and the SPITFIRE fire module (Guimberteau et al., 2018). Importantly, the representation of forest age dynamics in this version (Yue et al., 2018) allows us to simulate wood harvest based on rotation length practices, a prerequisite for simulating the woody yields.

There is another ORCHIDEE version including short-rotation coppice poplar plantations (ORCHIDEE-SRC; Fig. S1 in the Supplement, De Groote et al., 2015) based on the forest management module (Bellassen et al., 2010), but ORCHIDEE-SRC is more designed for studying specific coppicing processes and is evaluated using only two coppicing sites in Belgium. Although detailed forest management processes are not included in ORCHIDEE-MICT, this version includes explicit gross land use changes, i.e., the rotational transitions from other vegetation types to woody bioenergy crops and periodic clear-cut harvest of forests. These features are important to study the carbon emissions from bioenergy crops when their areas expand by converting other land use types in future BECCS scenarios. In addition, ORCHIDEE-MICT contains a bookkeeping system to track different forest age classes as separate land cohorts at a sub-grid scale (Yue et al., 2018). This functionality allows for the simulation of woody harvest based on rotation length tracking the carbon stock dynamics of different age classes of forests individually. In addition to the poplar plantation in Europe in ORCHIDEE-SRC (De Groote et al., 2015), we aimed to include herbaceous bioenergy crops like Miscanthus and switchgrass as well as other woody crops like eucalypt and willow in a more systematic way on the global scale.

Originally, there are 13 plant functional types (PFTs) in ORCHIDEE (Table 1; Krinner et al., 2005). In order to represent the bioenergy crops, we introduced four new PFTs (Table 1). PFT14 is a tropical tree representing eucalypt (Eucalyptus spp.); PFT15 is a temperate tree representing poplar (Populus spp.) and willow (Salix spp.); PFT16 and PFT17 are treated as crops representing Miscanthus and switchgrass (Panicum spp.), respectively. The reason for separating Miscanthus and switchgrass into two PFTs is that they are significantly different in biomass yields and resource use efficiency (Heaton et al., 2008; Dohleman et al., 2009). The default model equations of the four new bioenergy crop PFTs follow the ones of similar PFTs already defined in the model (Table 1), i.e., tropical broad-leaved evergreen (PFT2) for PFT14, temperate broad-leaved summergreen (PFT6) for PFT15, and C₄ crop (PFT13) for PFT16 and PFT17. Some parameters were, however, adjusted specifically for their corresponding bioenergy crops based on field experiment or measurement data in Sect. 2.3.

2.2 Bioenergy biomass harvest module

The new module represents the periodical harvest of bioenergy crops, consisting of two subroutines differentiating woody and herbaceous crops. For woody types, harvest is based on simulated forest age classes (see details in Yue et al., 2018). Briefly, each woody PFT is subdivided into six cohort functional types (CFTs) corresponding to different age classes. The boundary of age classes is set as PFT specific and defined based on maximum woody biomass (total of the sapwood and heartwood biomass). When the biomass of a young woody CFT reaches the upper boundary defining its age class, it is moved to the next older CFT sequentially until it reaches the oldest CFT (mature). The fractional harvested area of a woody crop PFT in each grid cell is externally prescribed. Then, the harvest algorithm starts from the second youngest CFT, continues with the next older CFT, and eventually reverts to the youngest CFT until the prescribed harvested area is met. For woody bioenergy crops, we adjusted the fraction of aboveground biomass that is harvested (the harvest index denoted HI) and put harvested biomass into a separate bioenergy harvest pool rather than mixing it with the modeled wood product pools existing for forest management harvest (Yue et al., 2018) or with an agricultural product pool for the two crop PFTs (PFT13 and 14; Table 1) as defined by Piao et al. (2009). The non-harvested biomass goes to litter. For herbaceous types, only the HI fraction of aboveground biomass is harvested (Sect. 2.3.4) after leaf senescence either at the end of growing season or if climate conditions like drought and low temperature trigger canopy senescence in the model. The remaining part of aboveground and belowground biomass goes to litter pools. Carbon in the bioenergy harvest pool is released to the atmosphere directly.

2.3 Parameterization of bioenergy crops

Most parameters in ORCHIDEE are PFT specific (Krinner et al., 2005). Since we aim to improve the biomass production performance of the four bioenergy crop PFTs, we adjusted parameters controlling carbon assimilation (Sect. 2.3.1), allocation (Sect. 2.3.2), phenology (Sect. 2.3.3), and harvest processes (Sect. 2.3.4) based on observed values at ecosystem or leaf scale (Table 2). The number of observations for each parameter varied due to the availability of data, and the sample may also be biased in terms of different species or climate conditions. For each parameter, we collected observational values by a detailed literature survey and used the observational medians first. We then evaluated the model predictions of biomass yields using yield observations. If there is a bias, we adjusted the parameter value within the observational range to reduce the misfit between predicted and observed yields.

2.3.1 Photosynthesis parameters

The photosynthesis process at leaf level for C₃ and C₄ plants in ORCHIDEE-MICT is based on the extended version (Yin and Struik, 2009) of the Farquhar, von Caemmerer, and Berry model (FvCB model; Farquhar et al., 1980). The related parameters generally follow Yin and Struik (2009) except for the maximum rate of Rubisco activity (V_cmax) and maximum rate of electron transport under saturated light (J_max). The setting of V_cmax and J_max for C₃ plants is based on Medlyn et al. (2002) and Kattge and Knorr (2007) in order to account for the acclimation of V_cmax and J_max to temperature. In ORCHIDEE, V_cmax 25 (V_cmax at 25 ^∘C) is prescribed for each PFT, and J_max is calculated from the ratio (r_JV) between J_max and V_cmax as

$$\begin{array}{}\text{(1)}& {\displaystyle }{J}_{\max }=V_{\mathrm{cmax}}\times r_{JV}\left(\mathrm{1}\right),\end{array}$$

where r_JV is a function of growth temperature (T_growth) (Kattge and Knorr, 2007):

$$\begin{array}{}\text{(2)}& {\displaystyle }{r}_{JV}=a_{rJV}+b_{rJV}\times T_{\mathrm{growth}},\end{array}$$

where a_rJV and b_rJV are the acclimation parameters derived by fitting data from 36 plant species (Kattge and Knorr, 2007). For C₄ plants, no acclimation is considered for V_cmax and J_max, and thus b_rJV= 0, and a_rJV is a fixed value (Table 2).

Because values of V_cmax and J_max are critical for determining the carbon assimilation by bioenergy PFTs, we searched for published experimental data on these parameters for eucalypt, poplar, willow, Miscanthus, and switchgrass and found 26 observation-based publications with 127 entries for V_cmax and 69 entries for J_max (Table S1 in the Supplement).

Some observations of V_cmax and J_max were derived at temperatures (Table S1) other than 25 ^∘C, and we thus normalized these two temperature-dependent variables to V_cmax 25 and J_max 25 (J_max at 25 ^∘C) using a modified Arrhenius function from Medlyn et al. (2002) and parameters for C₃ and C₄ plants from Yin and Struik (2009). The ranges of V_cmax 25, J_max 25, and r_JV 25 (r_JV at 25 ^∘C, only for the studies reporting both V_cmax and J_max) are shown in Fig. 1a–c. V_cmax 25 values generally decrease from eucalypt > poplar and willow > Miscanthus ≥ switchgrass. The interquartile range of V_cmax 25 is large for eucalypt (N= 42) from 75 to 126 µmol m⁻² s⁻¹ and for poplar and willow (N= 30) from 57 to 165 µmol CO₂ m⁻² s⁻¹. Miscanthus and switchgrass have a relatively smaller interquartile range of V_cmax 25 (17 to 32, N= 38 and 12 to 26, N= 17, respectively). We adjusted the prescribed parameters V_cmax 25 and a_rJV (Table 2) for each bioenergy crop PFT using a value close to the median value in the observation dataset (Fig. 1a, c, within a range of 10 % of the median values). We also verified that J_max 25 from Eq. (1) is in the range of independent J_max 25 observations (Fig. 1b). Importantly, the observation-based estimates of V_cmax 25 and J_max 25 for Miscanthus are significantly larger than for switchgrass (p= 0.02 and 0.09, respectively; Fig. 1a, b). Note that the ranges shown in Fig. 1 could be influenced by the sample size and number of studies.

Figure 1 V_cmax 25 (a), J_max 25 (b), r_JV 25 (V_cmax 25∕J_max 25, c), and specific leaf area (SLA, d) collected from measurements. The box plot indicates the interquartile range of measurements. The data size of measurements is shown below the box. The default values (open circles) and adjusted values (filled circles) for bioenergy crops are also shown. Because the model does not prescribe J_max 25 but rather calculates it from V_cmax 25, the J_max 25 values for ORCHIDEE shown here (circles in b) are calculated by V_cmax 25 × r_JV (circles in a and b).

We also adjusted other parameters including θ (the convexity factor of the response of the rate of electron transport to irradiance), α_(LL) (conversion efficiency of absorbed light into e-transport rate at strictly limiting light), and g₀ (residual stomatal conductance when irradiance approaches zero) in the leaf-level photosynthesis equations of ORCHIDEE to match higher productivity based on field measurements or empirical data (Table 2). The detailed effects of these parameters on photosynthesis in the FvCB model can be found in Yin and Struik (2009). Briefly, θ and α_(LL) are used in the calculation of J (photosynthesis rate limited by electron transport):

$$\begin{array}{}\text{(3)}& {\displaystyle }J={\displaystyle \frac{{\mathit{\alpha }}_{\left(\mathrm{LL}\right)}I+J_{\max }-{\sqrt{({\mathit{\alpha }}_{\left(\mathrm{LL}\right)}\mathrm{I}+{\mathrm{J}}_{\max })}}^{\mathrm{2}}-\mathrm{4}\mathit{\theta }{J}_{\max }{\mathit{\alpha }}_{\left(\mathrm{LL}\right)}I}{\mathrm{2}\mathit{\theta }}},\end{array}$$

where I is the photon flux density absorbed by leaf photosynthetic pigments. g₀ is an intercept related to the estimation of g_s (stomatal conductance):

$$\begin{array}{}\text{(4)}& g_{\mathrm{s}}=g_{\mathrm{0}}+{\displaystyle \frac{A+R_{\mathrm{d}}}{C_{\mathrm{i}}-C_{{\mathrm{i}}^{\ast }}}}{f}_{\mathrm{VPD}},\end{array}$$

where A is the net photosynthesis rate, R_d is the day respiration, and C_i and $C_{{\mathrm{i}}^{\ast }}$ are the intercellular CO₂ partial pressure and C_i-based CO₂ compensation point in the absence of R_d, respectively. f_VPD is factor of the effect of the leaf-to-air vapor pressure difference (Yin and Struik, 2009).

Table 2 Parameters adjusted for bioenergy crop PFTs in ORCHIDEE-MICT. Values of original PFTs are also shown for comparison.

Specifically for bioenergy crop PFTs, we increased θ to 0.8 for PFT14 (eucalypt) based on Yin and Struik (2017) and to 0.84 for PFT16 (Miscanthus) based on field measurements from Dohleman and Long (2009). Light use efficiency and productivity are high for bioenergy crops (see, e.g., reviews by Heilman et al., 1996; Lewandowski et al., 2003; Whitehead and Beadle, 2004; Karp and Shield, 2008; Zub and Brancourt-Hulmel, 2010; Forrester, 2013; Laurent et al., 2015; McCalmont et al., 2017), and we thus set α_(LL) and g₀ to the maximum boundary in their ranges from Yin and Struik (2009) to favor the high light use efficiency and productivity characteristic of bioenergy cultivars (Table 2). Morphological plant traits are also of key importance to the canopy-level productivity (Chang et al., 2015). The specific leaf area (SLA) in ORCHIDEE is a PFT-specific constant (Krinner et al., 2005). SLA for different bioenergy crops from our data compilation (164 entries in Table S1) is shown in Fig. 1d. A factor of 2 is used to convert the SLA unit from m² g⁻¹ dry matter to m² g⁻¹ C. Observation-derived SLA for eucalypt is lower than for the other bioenergy crops, and SLA for switchgrass is relatively larger. SLA is set to the median value of observations for PFT16 (Miscanthus) and PFT17 (switchgrass), and it is close to the 75th percentile value of the data we compiled for PFT14 (eucalypt) and PFT15 (poplar and willow; Fig. 1d and Table 2).

Another important plant trait for photosynthesis is the leaf orientation, which determines the radiation extinction in the canopy. Although the LAI of eucalypts is generally moderate (Anderson, 1981; Stape et al., 2004; Whitehead and Beadle, 2004), leaf angles are nearly close to vertical in mature eucalypt forest (Anderson, 1981; King, 1997), leading to a good distribution of radiation to the lower canopy layers. The light extinction coefficient (k) for PFT14 (eucalypt) is therefore set to 0.36 (Table 2) according to the measurement-based estimate by Stape et al. (2004). Similarly, a field study shows the seasonal average k ranging from 0.23 to 0.37 for poplars (Ceulemans et al., 1992; Heilman et al., 1996), and a median value of 0.3 was used for PFT15 (Table 2).

2.3.2 Carbon allocation parameters

The maximum carbon allocation to leaf biomass is controlled in ORCHIDEE by a predefined maximum LAI value (LAI_max) beyond which no carbon will be allocated to leaf (Krinner et al., 2005). We adjusted this parameter to match the observed maximum LAI in the field for the four selected bioenergy plants (Table 2). LAI_max for PFT14 (eucalypt), PFT15 (poplar and willow), PFT16 (Miscanthus), and PFT17 (switchgrass) are set to 7, 9, 10, and 8, respectively (Ceulemans et al., 1992; Heilman et al., 1996; Whitehead and Beadle, 2004; Heaton et al., 2008; Zub and Brancourt-Hulmel, 2010).

For woody PFTs in ORCHIDEE, the partitioning between aboveground and belowground sapwood biomass is a function of forest age (Krinner et al., 2005):

$$\begin{array}{}\text{(5)}& {\displaystyle }{f}_{\mathrm{ab},t}=f_{\mathrm{ab},\min }+(f_{\mathrm{ab},\max }-f_{\mathrm{ab},\min })\times (\mathrm{1}-e^{-t/\mathit{\tau }}),\end{array}$$

where f_ab, t is the fraction of sapwood allocated to aboveground at age t; f_ab, min and f_ab, max are the minimum and maximum fraction allocated to aboveground (0.2 and 0.8, respectively); and τ is an empirical parameter. This equation implies that more biomass is allocated to belowground sapwood to develop coarse roots in younger forests. The partition between aboveground and belowground biomass is influenced by resource supply like water and nutrient availability (Litton et al., 2007). For example, belowground carbon allocation in eucalypts is observed to be strongly reduced by irrigation (Barton and Montagu, 2006; Stape et al., 2008; Ryan et al., 2010). Fertilized poplars also showed greater shoot growth than control plots (Coleman et al., 2004). We assumed that bioenergy trees should usually be under intensive management (e.g., irrigation and fertilization), especially in the establishment year (Jacobs, 1981; Caslin et al., 2015; Isebrands and Richardson, 2014). A higher water and nutrient availability then implies a lower investment of biomass on roots for bioenergy trees. In the ORCHIDEE version used here, as there is no specific fertilization or irrigation practice included, the idealized approach chosen to partially account for these managements operations was to reduce τ in Eq. (5) from 5 to 2 years (Table 2) to give a maximum allocation of sapwood biomass to aboveground faster than in the standard version. The difference in these two values is illustrated in Fig. S2. Also because the rotation length for bioenergy trees is usually several years only (Karp and Shield, 2008), it is reasonable to assume that these plants allocate more biomass aboveground in the first few years. However, trees like poplar and willow can sprout from the remaining stem or root (Isebrands and Richardson, 2014), which is not accounted for in the model. Last, we also adjusted the factor (β, Table 2) in the exponential function to calculate the soil water stress in ORCHIDEE (Krinner et al., 2005; McMurtrie et al., 1990) to reduce the soil moisture stress on bioenergy trees (Fig. S3).

2.3.3 Phenology parameters

An adjustment of parameters related to phenology was performed for the two herbaceous bioenergy PFTs (PFT16 and PFT17; Table 2) to derive the total biomass production for the whole growing season. Lewandowski et al. (2003) and Zub and Brancourt-Hulmel (2010) reviewed the growth temperature and growing season length of Miscanthus and switchgrass and found that these two crops have higher cold tolerance and a longer growing season than grasses. Compared to maize, Miscanthus has an earlier leaf onset and later leaf fall, and thus its growing season length is 59 % longer (Dohleman and Long, 2009). Some Miscanthus genotypes need fewer cumulative degree days for shoot emergence (60 to 118 degree days) and a high frost tolerance (−9 to −6 ^∘C; Farrell et al., 2006). To account for this frost tolerance and longer growing season, we decreased the growing degree days for leaf onset in the model (GDD_onset) from 700 (standard value for C₄ crop PFT) to 320 degree days (same as the default value for C₄ grass PFT in ORCHIDEE) and the critical temperature for leaf senescence (T_senescence) from 10 to 0 ^∘C for PFT16 and PFT17 (Table 2). Note that we did not set T_senescence as −9 to −6 ^∘C because frost tolerance was only documented for certain Miscanthus genotypes, so we used a conservative value of 0 ^∘C for Miscanthus and switchgrass PFTs. In addition, we increased the critical leaf age beyond which leaves enter senescence (t_leaf) and the minimum leaf age to allow leaf senescence (t_leaf, min) to be the same as the default values for C₄ grass PFT (PFT11 in Table 1) in ORCHIDEE (Table 2).

Figure 2 Harvest age (a) and rotation length (b) in the evaluation dataset. The box plot indicates the interquartile range, and the number of observations is also shown. In this study, the harvest age class is set to 4–6 years (red shading).

2.3.4 Biomass harvest

The harvest index (HI) determines how much aboveground biomass is harvested. Theoretically, all the aboveground biomass of a lignocellulosic crop can be used for energy production. Some IAMs (e.g., GCAM3.0; Kyle et al., 2011) indeed assume an HI of 1 for switchgrass, for instance. In practice, harvesting Miscanthus and switchgrass is usually performed in winter and early spring after drying and nutrient recycling through leaf senescence (Lewandowski et al., 2003; Zub and Brancourt-Hulmel, 2010), which leads to a lower biomass at harvest but enhances nutrient conservation. For example, 18–46 % of the nitrogen in Miscanthus can be recycled through leaf falling to soil and translocation from shoots to rhizomes (Cadoux et al., 2012). Similar seasonal nitrogen dynamics were also observed for switchgrass (Heaton et al., 2009). In fact, Miscanthus is recommended to be harvested between January and March in practice guidelines (Nixon and Bullard, 2001). Otherwise, fertilizers have to be applied to amend the nutrient removal from harvest, which is neither cost effective nor environmentally friendly. For bioenergy trees, current harvesting techniques can hardly harvest 100 % of aboveground biomass (Jacobs, 1981; Isebrands and Richardson, 2014; Caslin et al., 2015). Following Caslin et al. (2010), Richards et al. (2017), and Zhuang et al. (2013), we used an HI of 0.9 (i.e., 90 % of aboveground biomass is harvested) for all the bioenergy PFTs in ORCHIDEE (Table 2). However, for simulations using future land use maps generated from IAMs, we would recommend setting the HI the same as in IAMs to be consistent.

The rotation length for eucalypt, poplar, and willow varies among different tree types, species, locations, and plantation purposes (Ugalde and Pérez, 2001; Keoleian and Volk, 2005; Karp and Shield, 2008; Isebrands and Richardson, 2014; Caslin et al., 2015). For example, eucalypt and poplar for sawlog and veneer utilization are often on rotations of 8–20 years, depending on regions (Ugalde and Pérez, 2001; Isebrands and Richardson, 2014). But short-rotation coppice bioenergy plantations of poplar and willow have shorter cutting cycles of 3–5 years (Keoleian and Volk, 2005; Karp and Shield, 2008; Isebrands and Richardson, 2014; Caslin et al., 2015). A rotation length of 8 years was used in the LPJml model for bioenergy trees (Beringer et al., 2011). In ORCHIDEE, the rotation length for bioenergy tree PFTs is associated with the setting of age classes (see Sect. 2.2). Namely, harvesting starts from the second youngest age class, and thus the age in the second youngest forest age cohort should be set up as the same as the rotation length. For the idealized simulations presented below, we used a rotation length of 4–6 years based on the harvest age and rotation length in the evaluation dataset (Fig. 2; Sect. 3.2). Here, the harvest age (Fig. 2) represents the age at which the biomass of bioenergy trees was harvested or estimated. It is directly reported by the original literature and corresponds to the reported yield. Rotation length (Fig. 2) is the management practice reported in the original literature, and it is the same as the harvest age in most studies. In other studies, however, some trees may be harvested earlier or later than the regular rotation length, e.g., for comparison purposes. In addition, not all literature reported both harvest age and rotation length (see the number of observations in Fig. 2).

3 Model evaluation

3.1 Evaluation dataset

We used the global bioenergy crop yield dataset from Li et al. (2018, see “Data availability”) to evaluate the performance of the modified ORCHIDEE-MICT-BIOENERGY model. This global dataset was compiled from more than 200 field-measurement-based studies with five main bioenergy crop types, i.e., eucalypt, poplar, willow, Miscanthus, and switchgrass (Li et al., 2018). Most of the measurements (> 90 %) are based experimental trials, especially for Miscanthus and switchgrass. About 98 % of the compiled observations are reported as the aboveground biomass, and the rest are reported as the total of aboveground and belowground biomass. We thus did not exclude the observations of the total biomass in the model–observation comparison since their fraction is very low (< 2 %). The biomass yield in this dataset is compiled in a unit of ton DM (dry matter) ha⁻¹ yr⁻¹, corresponding to the mean annual biomass yield. For example, if the original literature reported the total harvested biomass of poplar at a certain age, the total biomass amount is divided by age to get the mean annual biomass yield. If the original literature reported the annual harvested biomass of Miscanthus for several years, each annual yield is taken as one observation. Note that this dataset does not distinguish the utilization of the plantation (for bioenergy use or for timber or pulpwood). In order to evaluate the simulated biomass yields by ORCHIDEE at half-degree resolution, we calculated the median and range of all observations in each half-degree grid cell containing at least one site of the dataset. Each half-degree grid cell may contain observations from different sites or one site with different species, genotypes, and/or treatments (e.g., different irrigation or fertilization levels). Globally, the number of half-degree grid cells containing observations for PFT14 (eucalypt), PFT15 (poplar and willow), PFT16 (Miscanthus), and PFT17 (switchgrass) are 63, 120, 69, and 44, respectively (see maps in Sect. 3.5), giving a total of 296 grid cells (some may have several crops in common).

3.2 Simulation setup

The setup for the site-scale simulations in ORCHIDEE-MICT-BIOENERGY is as follows. The model is forced with 30 min time step climatic forcing data, CRU-NCEP v7 (Viovy, 2017) recycling the period of 1990–2000. The CRU-NCEP forcing data represent a merged product of the CRU TS climate dataset (Harris et al., 2014) and NCEP reanalysis data (Kalnay et al., 1996). Some observation sites have reported mean annual temperature (MAT) and precipitation (MAP), and we verified that these data are consistent with the MAT and MAP from the CRU-NCEP v7 climate forcing data we used (Fig. S4). Thus no bias correction was applied to the CRU-NCEP v7 climate forcing. The soil texture map used in the model is based on the 12 USDA texture classes from Reynolds et al. (2000).

Figure 3 Biomass yields from the observations and simulated by the ORCHIDEE model. The error bars of observations in the left panel represent the range of different observations in this half-degree grid cell caused by different sites, treatments, species, and genotypes. The error bars of modeled yields represent the range of different harvest years caused by the interannual variability of climate. PFT14 is tropical bioenergy tree, eucalypt; PFT15 is temperate bioenergy tree, poplar and willow; PFT16 is C₄ bioenergy grass, Miscanthus; PFT17 is C₄ bioenergy grass, switchgrass. The red line indicates the 1 : 1 ratio line.

We assumed a homogenous coverage (100 %) of one single bioenergy crop PFT in a grid cell covered by the same PFT type as the site observations. We set an annual harvest fraction of 1 % of the grid cell each year. The 1 % annual harvest fraction is just an artificial value to make sure that there is always forest in the second youngest age class available for harvest every year after a stable rotation is established. We compared the annual harvested biomass in the bioenergy harvest pool per area unit, so the harvest area has no influence on our model evaluation. For the bioenergy trees (PFT14 and PFT15), a spin-up of 100 years without harvest was run first to derive biomass evolution in time to define the respective biomass boundaries for age classes in each grid cell (see Yue et al., 2018). The biomass boundaries are grid-cell specific because of the different vegetation growth rates in different grid cells. The six age classes from youngest to oldest are thus set to correspond to 0–4, 4–6, 6–10, 10–30, 30–50, and > 50 years. We set the second youngest age class that is used in priority for bioenergy harvest (Sect. 2.2) to 4–6 years (Fig. 2) based on harvest age and rotation length reported by the original publications in the evaluation dataset (Li et al., 2018). After spin-up, the simulations for PFT14 (eucalypt) and PFT15 (poplar and willow) were run with a bioenergy harvest process for 50 years because we only harvested the second age class (4–6 yr) and 50 years is long enough to establish a stable rotation. The harvested biomass amount for the last 10 years was used to calculate the median and range of the simulated yields. Note that we artificially harvest 1 % of the grid cells each year, and the harvested patches will be planted immediately. After the first 5 years (one rotation length), there is always a fraction reaching a full rotation and ready for harvest. The harvest in the last 10 years thus represents 10 harvest events. We divided the harvested biomass by 5 years (4–6 years in the second youngest age class) to obtain the annual mean yields of PFT14 (eucalypt) and PFT15 (poplar and willow). A carbon-to-dry-matter ratio of 0.5 was used to convert the unit of yields into ton DM ha⁻¹ yr⁻¹.

For the bioenergy grasses (PFT16 and PFT17), simulations were performed directly (without spin-up) with harvest for 50 years, and similarly, the yields of the last 10 years were used for comparison with site-observed values. Note that we aim to assess the performance of simulated biomass yields rather than the state of the carbon pools, including litter and soil organic matter, that depend on site history. As litter and soil carbon pools do not influence vegetation productivity in the model, we did not perform a full long spin-up of carbon pools to their equilibrium values.

3.3 Simulated bioenergy yields at global level

The simulated bioenergy biomass yields in comparison with field observations for the four bioenergy crops are shown in Fig. 3. The model–observation results generally lie around the 1 : 1 ratio line (Fig. 3a, c, e, g). Although the regression between modeled and observed medians is not significant with a low r² value because of the overestimation and underestimation at some sites (Fig. 3a, c, e, g), the difference between the two samples of modeled and observed yields is not significant (t test, p > 0.17) and the percent bias (PBIAS, defined as the sum of biases divided by the sum of observed values; Moriasi et al., 2007) ranges from 2 to 8 % for all PFTs, implying that the global distributions of modeled and observed yields are consistent. ORCHIDEE-MICT-BIOENERGY reproduces the frequency distributions of the observed biomass yields across different grid cells well (Fig. 3b, d, f, h). The median observed and simulated biomass yields in all grid cells are 16.0 and 17.5 ton DM ha⁻¹ yr⁻¹ for PFT14 (eucalypt), 8.4 and 8.3 ton DM ha⁻¹ yr⁻¹ for PFT15 (poplar and willow), 12.7 and 10.8 ton DM ha⁻¹ yr⁻¹ for PFT16 (Miscanthus), and 8.7 and 9.0 ton DM ha⁻¹ yr⁻¹ for PFT17 (switchgrass). PFT14 (eucalypt) shows a large spread in both the observed and simulated biomass yields. Some site observation data with high yield (> 25 ton DM ha⁻¹ yr⁻¹) were not reproduced by the model for eucalypt. By contrast, observed and simulated yields for PFT15 (poplar and willow) and PFT17 (switchgrass) concentrate in a relatively narrow range. In addition, the error bars for most sites (67, 73, 74, and 64 % for PFT14 to PFT17, respectively) reach the 1 : 1 line (Fig. 3a, c, e, g), implying that at least some observations in these grid cells can be represented by the model.

Figure 4 Biomass–age curves at different sites for PFT14 (tropical bioenergy tree, eucalypt). The site number, coordinates, and country for each site are also shown. Biomass at most sites refers to aboveground biomass, except Site nos. 5, 11, and 12 (labeled “total”, i.e., the sum of aboveground and belowground biomass; the same total biomass from the model is used for these sites). The detailed site information is shown in Table S2.

It should be noted that it is impossible to perfectly reproduce observations in all grid cells, i.e., all dots in Fig. 3 on the 1 : 1 line, because of uncertainties in the observation dataset, e.g., treatments, genotypes, and local fertilization or irrigation practices, as well as in soil characteristics and climate forcing variations prescribed in the model. The error bars of modeled yield (y axis) come from the range of different harvest years and represent interannual variability. The error bars of the observations (x axis) represent the range from different sites, crop species, genotypes, and treatments, as well as the observation number in each grid cell. It is difficult to systematically synthesize all these factors to give an optimal observed yield in each grid cell. First, different species and genotypes are impossible to be accounted for in a global vegetation model, and thus a further classification of such information would not help the model evaluation. Second, some management practices are difficult to quantify. For example, some studies report irrigation as an amount per year, while others report only “irrigating when necessary”. The fertilization rates are also difficult to synthesize between different studies because they applied different types of fertilizers, some annually but some in random years. Third, each observation is associated with different managements and treatments, and there is no uniform standard to weight all these different managements. Last, global vegetation models usually run at a half-degree resolution, which may not fully represent the site-level climate variations and soil properties.

Figure 5 Biomass–age curves at different sites for PFT15 (temperate bioenergy tree, poplar and willow). The site number, coordinates, and country for each site are also shown. Biomass at most sites refers to aboveground biomass, except Site no. 17 (labeled “total”, i.e., the sum of aboveground and belowground biomass; the same total biomass from the model is used for this site.). The detailed site information is shown in Table S2.

3.4 Biomass–age relationship at site level

A good representation of biomass–age curves for bioenergy trees in the model is crucial to reproduce the yields, especially in the first several years after planting (≤ rotation length). However, most of the observations in the global evaluation dataset were only mean annual yield (Li et al., 2018). This precludes a detailed analysis of biomass dynamics over time for bioenergy trees. We thus selected 22 studies (Table S2) from the evaluation dataset that reported the biomass amount of multiple ages (at least 2 years) and at the same site for eucalypt, poplar, or willow. We went through the original articles to derive the biomass–age curves and compared them with the same curves from the model simulations (Figs. 4 and 5).

There is a good agreement on the biomass–age relationship of eucalypt between models and observations for some sites in Australia and China (Site nos. 2, 8–12 in Fig. 4). But the model underestimates the biomass evolution of eucalypt at Site no. 13 in New Zealand and overestimates it at Site nos. 5–6 in China (Fig. 4). For poplar and willow, there are two long-term (> 10 years) consecutive observation sites in Wisconsin, USA (Site no. 13 and no. 15 in Fig. 5), where the model captures the biomass–age relationship well. At some other sites (Site nos. 2, 3, 7, 14, and 17 in Fig. 5), however, the model results only agree with observations for the first few years and then deviate from the observations afterwards. The model generally coarsely underestimates the biomass of poplar and willow at all ages at the sites in the eastern (Site no. 1 in Fig. 5) and western (Site nos. 5 and 6 in Fig. 5) coastal regions of the US, in the UK (Site nos. 8 to 11 in Fig. 5), and in Sweden (Site no. 16 in Fig. 5), but overestimates in India (Site no. 12 in Fig. 5) and at one site in China (Site no. 18 in Fig. 5).

Figure 6 The map of relative difference between simulated and observed biomass yields for PFT14 (tropical bioenergy tree, eucalypt). The inset plot shows the frequency of the relative difference between the model and observation.

Possible reasons for the model–observation differences at each site using the information reported in the original studies (see details in Table S2) include the different varieties of species (e.g., genotypes) and management (e.g., fertilization, irrigation, or spacing) in the field, which were not explicitly considered in the model. For example, the model overestimates biomass at Site no. 4 in Fig. 4 because of the large spacing of plantation in the trial at that experimental site (Han et al., 2010), which results in lower biomass yield when converting the unit of ton DM plant⁻¹ yr⁻¹ to ton DM ha⁻¹ yr⁻¹. Site nos. 13, 14, and 15 in Fig. 5 are from the same study (Strong and Hansen, 1993), and the model reproduces it at Site nos. 13 and 15 but underestimates it at Site no. 14. This is because the biomass–age curves at Site nos. 13 and 15 are from the average of several genotypes (some have higher yields and some lower), but only one genotype with relatively high yield was planted at Site no. 14 (Strong and Hansen, 1993), causing a model underestimation at Site no. 14. In addition, our model seems to systematically underestimate the biomass production of willow for the sites in the UK (Site nos. 8–11 in Fig. 5). This observed biomass production in the UK was based on a range of willow varieties in trial experiments, and the authors (Lindegaard et al., 2011) claimed that the trial experiments generate higher yields than large-scale commercial plantations because of the differences in land quality and practice guidelines (e.g., cutting, harvest index). Despite some model–observation differences, we emphasize that the modeled biomass–age curves are consistent with observations for most sites within the rotation length.

Figure 7 The map of relative difference between simulated and observed biomass yields for PFT15 (temperate bioenergy tree, poplar and willow). The inset plot shows the frequency of the relative difference between the model and observation.

Figure 8 The map of relative difference between simulated and observed biomass yields for PFT16 (C₄ bioenergy grass, Miscanthus). The inset plot shows the frequency of the relative difference between the model and observation.

Figure 9 The map of relative difference between simulated and observed biomass yields for PFT17 (C₄ bioenergy grass, switchgrass). The inset plot shows the frequency of the relative difference between the model and observation.

Figure 10 The relative difference between simulated and observed yield in different MAT and MAP intervals. The median values of model–observation differences of all grid cells in each MAT and MAP interval are shown. PFT14 is tropical bioenergy tree, eucalypt; PFT15 is temperate bioenergy tree, poplar and willow; PFT16 is C₄ bioenergy grass, Miscanthus; PFT17 is C₄ bioenergy grass, switchgrass.

3.5 Maps of differences between simulated and observed yields

The spatial distributions of relative differences between simulated and observed biomass yields are shown in Figs. 6 to 9 for each PFT. The observations for eucalypt are mainly distributed in Brazil, tropical Africa, South Asia, and Australia (Fig. 6). ORCHIDEE-MICT-BIOENERGY slightly underestimates biomass yield for PFT14 (eucalypt) in Brazil and overestimates some grid cells in southern China and Australia. Some biomass observations of eucalypts in Australia are obtained from native forests (Li et al., 2018), which may partly explain the overestimation by the model.

Poplar and willow are mainly planted in temperate regions like the United States, Europe, and Central and East Asia (Fig. 7). ORCHIDEE-MICT-BIOENERGY underestimates the biomass yields for PFT15 (poplar and willow) in the western US but overestimates the yields in the eastern US. There is no distinct pattern for the differences between observations and model results in Europe, with both underestimation and overestimation across grid cells. But it seems that the simulated biomass yields are lower than observations in Sweden. In Central and East Asia, biomass yields in the inland grid cells are generally underestimated but those in the coastal areas are overestimated.

Most of the observations for Miscanthus are from Europe although some trial tests are also available in the eastern US and a few in China (Fig. 8). In the US, a very slight underestimation of yield was found in the inland areas, while overestimation occurred more close to the ocean. The model underestimates biomass yields for PFT16 (Miscanthus) in the UK and southern Europe and slightly overestimates it in other areas in Europe. There are only three grid cells with Miscanthus yield observations in China, and they are all largely overestimated in the simulations.

Switchgrass is a native perennial grass in North America (Lewandowski et al., 2003) and thus mainly grows in the US (Fig. 9). There are also very few observations in Europe and East Asia. ORCHIDEE-MICT-BIOENERGY can generally reproduce the biomass yields for PFT17 (switchgrass) in the central US but overestimates in the eastern US, especially in some grid cells around the Great Lakes. The simulated biomass yields are lower than observations in grid cells in Europe and China but fit well with observations in the grid cell in Japan.

3.6 Model–observation difference in different climate bins

The differences between simulated and observed biomass yields for bioenergy crop PFTs in different MAT and MAP intervals are shown in Fig. 10. There is no systematical bias of simulated biomass yields in the climate space except in the climate zones with relatively high MAT and MAP (upper-right grids in Fig. 10b, c) for PFT15 (poplar and willow) to PFT17 (switchgrass). For these PFTs, it seems ORCHIDEE overestimated the yields with MAT > 15 ^∘C and MAP > 1000 mm yr⁻¹. The strong underestimation (darker blue color) seems more aligned to the drier regions, especially for poplar and willow (PFT15, Fig. 10b).

The distribution patterns in Fig. 10 also reflect the different climate conditions of growth for these bioenergy crops. Consistent with their physiological characteristics, eucalypts grow in tropical regions (Fig. 6) with MAT > 10 ^∘C and MAP > 500 mm yr⁻¹ (Fig. 10). By contrast, poplars and willows grow in temperate regions (Fig. 7) and some under low MAT and MAP (Fig. 10). Miscanthus and switchgrass are usually planted in Europe and the US (Figs. 8 and 9) with moderate MAT and MAP.

We further investigated whether other climate forcing variables in the model impact the model–observation differences using the multiple linear regression method (Table S3) and the regression tree method (Breiman et al., 1984; Pedregosa and Varoquaux, 2011; Fig. S5). In these two methods, PFT types and nine climate forcing variables (Table S3) were used as independent variables and the relative model–observation difference as a dependent variable. The multiple linear regression is nonsignificant (p= 0.28) with a very low r² (0.01), suggesting that the variations in the relative model–observation differences are mostly explained by other factors rather than the climate forcing biases used in the model. In the regression tree (Fig. S5), the first discriminator is shortwave radiation but it only split very few samples. Although north wind speed separates a relatively large proportion of samples (Fig. S5), it has little to do with the biomass production in the model. Therefore, results from these two regression methods suggest the model–observation biases are unlikely caused by the model simulation.

Figure 11 Comparison of biomass yields simulated by ORCHIDEE-MICT-BIOENERGY and previous versions. Only median values in half-degree grid cells, some containing multiple sites, are shown for both simulated and observed yields. Red circles represent the simulations using specific bioenergy parameterizations (same as Fig. 3). Gray squares represent using the herbaceous crop PFTs of the previous ORCHIDEE version for bioenergy trees, i.e., PFT13 for both PFT14 and PFT15 (Table 2). Blue triangles represent the simulations using the right PFTs but the parameters of the previous ORCHIDEE version, i.e., parameters of PFT2 (tropical broad-leaved evergreen), PFT6 (temperate broad-leaved summergreen), PFT13 (C₄ crop), and PFT13 (C₄ crop) for PFT14 (eucalypt), PFT15 (poplar and willow), PFT16 (Miscanthus), and PFT17 (switchgrass), respectively (Table 2).

4 Discussion

4.1 Model performance before and after bioenergy crop implementation

In this study, we added four new PFTs to represent the main lignocellulosic bioenergy crops and implemented new parameterizations for each new PFT. As a first step, we evaluated the biomass production from bioenergy crops in ORCHIDEE using a global field measurement dataset. We compared the biomass yields simulated by the new ORCHIDEE-MICT-BIOENERGY with the yields from the previous ORCHIDEE version (Fig. 11). In the previous version, bioenergy crops were all taken as herbaceous C₄ crops (PFT13), and thus severe overestimation (overestimating 60 % on average) occurs for tropical bioenergy trees (i.e., eucalypts; gray squares in Fig. 11a). Although using herbaceous C₄ crops generally reproduces the observed biomass yields of poplars and willows (gray squares in Fig. 11b), different carbon dynamics in litter and soil and water and energy balance can be expected.

Using the right tree PFTs for bioenergy trees and the right herbaceous PFTs for bioenergy grasses but without new parameterizations also results in significant biases in the simulated yields compared to observations (blue triangles in Fig. 11). Specifically, using the default parameters of the previous version is found to largely underestimate the biomass yields of all bioenergy trees (blue triangles in Fig. 11a, b). For bioenergy grasses, a slight underestimation was found for Miscanthus (blue triangles in Fig. 11c), while large overestimation was found for switchgrass (blue triangles in Fig. 11d) with previous default parameters. The large biomass yields of C₄ crops in the previous ORCHIDEE version (blue triangles in Fig. 11c, d) mainly result from the high V_cmax 25 (Table 2), which is not the reason for the high yields of Miscanthus and switchgrass (Fig. 1). We emphasize again that different bioenergy crops achieve high productivities through different pathways based on their plant traits (Sect. 2.3) and it is important to specifically consider these traits by proper parameterizations in the global vegetation models.

4.2 Management impacts on parameters

We adjusted some key parameters (e.g., V_cmax, J_max, and SLA) related to the productivity of bioenergy crops based on a collection of field measurements. We only took the medians and the ranges to validate the parameter values in the model but did not explicitly consider the impacts of management (e.g., fertilization, species) on these parameters, neither in the model nor in the measurements. Here, we summarize some management effects on these parameters for different bioenergy crops based on measurements as follows.

1. Miscanthus. Wang et al. (2012) found that biomass yield of Miscanthus increased under nitrogen addition through elevated SLA, but fertilization did not affect V_cmax, stomatal conductance (g_s), or the extinction coefficient (k). Yan et al. (2015) measured photosynthesis variables of three Miscanthus species in two experimental fields and found significantly higher g_s, J_max, and V_cmax of Miscanthus lutarioriparius than M. sacchariflorus and M. sinensis.

2. Switchgrass. SLA differed significantly among nine cultivars of switchgrass but did not respond significantly to water stress or nitrogen application for individual cultivars (Byrd and May II, 2000). Trócsányi et al. (2009) reported a lower SLA of switchgrass from the early harvest than from the late harvest. Hui et al. (2018) investigated the leaf physiology of switchgrass under five precipitation treatments and found a significantly higher photosynthesis rate and g_s under elevated precipitation but no significant difference under reduced precipitation compared to control plots.

3. Eucalypt. Lin et al. (2013) measured the photosynthesis response of six Eucalyptus species to temperature and found significantly different J_max 25 and V_cmax 25 among species but nonsignificant differences in their ratios (J_max 25 ∕ V_cmax 25) and in the temperature response of J_max and V_cmax. With extra nitrogen supply, J_max and V_cmax of Eucalyptus grandis increased significantly, mainly associated with elevated leaf nitrogen content (Grassi et al., 2002). Sharwood et al. (2017) also found that J_max and V_cmax of Eucalyptus globulus were correlated with leaf nitrogen content and the ratio of J_max ∕ V_cmax was constant under elevated CO₂ or elevated temperature, but SLA is influenced by different CO₂ and temperature treatments.

4. Poplar and willow. In experimental trials of three Populus deltoides clones and two P. deltoides × P. nigra clones, J_max and V_cmax of the former species were significantly higher than the latter hybrid despite some clonal variations (Dowell et al., 2009). Wullschleger (1993) summarized the species-specific estimates of J_max and V_cmax, and the five Populus species displayed large variations. In a poplar free-air CO₂-enrichment (PopFACE) experiment, P. alba, P. nigra, and P. euramericana showed significant differences in g_s but nonsignificant differences in J_max and V_cmax among species, while the elevated CO₂ significantly decreased J_max and V_cmax but had no influence on g_s species (Bernacchi et al., 2003). SLA was also found to differ significantly between the P. deltoides × P. nigra family and the P. deltoides × P. trichocarpa family (Marron et al., 2007). For willows, SLA increased significantly under fertilization and irrigation, but the magnitude of response varied among six varieties of Salix species (Weih and Rönnberg-Wästjung, 2007). Similarly, the response of SLA and g_s to nitrogen fertilization differed among three willow clones, but no significant difference in V_cmax was found between fertilization and control plots for all clones (Merilo et al., 2006).

In general, the values of parameters like V_cmax, J_max, and SLA differ among different species or genotypes within each bioenergy crop type. The parameter responses to management like fertilization and irrigation also show large variations depending on the specific species. Although the effects of management on these parameters seem evident in some cases, a set of quantitative relationships that can be applied in relation to simple management operations in a global vegetation model for large-scale and generalized PFTs is still lacking. Expanding the PFT level to the species level in global vegetation models requires substantial computational resources and, more importantly, there may not be enough measured parameters of each species for all the processes implemented in the models. At this stage, therefore, using the medians and ranges across a great number of observations is a more justified and practical way to tune the parameters in the models. But more field measurements and quantitative reviews of relationships between individual parameters and individual managements as well as interactions between different parameters and managements are highly needed in future research.

4.3 Management impacts on yields

Management like fertilization, irrigation, and species plays an important role in biomass yields. In ORCHIDEE-MICT-BIOENERGY, nutrient limitations, and management by irrigation and fertilization are not explicitly implemented. Instead, we used parameter values in the range that favors a higher productivity (Sect. 2.3, Fig. 1) and compared the simulated yields with the median values of all observations regardless of the management (Fig. 3). We further categorized the observations into three groups (fertilization, non-fertilization, or non-reported) and compared with simulated yields (Fig. S6). There is no systematic bias between simulated yields and yields at fertilized sites for all PFTs (orange dots in Fig. S6). The model seems to overestimate the yields of eucalypt at sites with non-reported information of fertilization (most gray dots above the 1 : 1 line in Fig. S6a, Table S4) and overestimate yields of poplar and willow at sites without fertilization (green dots in Fig. S6b, Table S4). Yields at sites with non-reported fertilization information are underestimated by the model for Miscanthus (gray dots in Fig. S6c, Table S4) but overestimated for switchgrass (gray dots in Fig. S6d, Table S4).

We did not group the observations based on different fertilization rates because there are large variations in the biomass response to fertilization rates. For example, in a quantitative review by Heaton et al. (2004), the relationship between yields of Miscanthus and nitrogen application rates were not significant. Cadoux et al. (2012) reviewed 11 studies that measured Miscanthus yields under fertilization, and the biomass response to nitrogen fertilization was positive in 6 of the studies with no response in the others. Similarly, some studies showed a positive biomass response of poplar to nitrogen fertilization, but others did not (Kauter et al., 2003). Eucalypt also showed a variable response to fertilization, while the general response was positive (De Moraes Gonçalves et al., 2004). In quantitative reviews of fertilization effects on yields of switchgrass (Wang et al., 2010) and willow (Fabio and Smart, 2018), the relationship between biomass yields and nitrogen fertilization rates was significantly positive but the coefficient of determination (r²) was very low. In summary, biomass response to fertilization varied largely, and evidence from field measurements is not conclusive. More importantly, the basic soil characteristics should be taken into account in addition to the fertilization rates, but unfortunately we did not have information on soil nutrient contents or types, nutrient stoichiometry, or the rates and timing of applied fertilizers for each site from observations.

We also separated the observations based on irrigation information (irrigation, nonirrigation, and non-reported) in comparison with modeled yields (Fig. S7). Both underestimation and overestimation were found for sites with different irrigation management for different PFTs. The yields of eucalypt were underestimated at sites with irrigation (blue dots in Fig. S7a, Table S4) but overestimated at sites with non-reported irrigation information (gray dots in Fig. S7a, Table S4). Compared to fertilization, not many sites reported irrigation information and the quantification of irrigation rates is more difficult. For example, some studies reported an irrigation amount per year, while some others only reported descriptive information like “soil moisture maintained to field capacity” or “irregular irrigation when necessary”.

Comparison between simulated yields and observations for the main species of bioenergy crops is shown in Fig. S8. The model overestimated yields of Eucalyptus urophylla × E. grandis, E. globulus, and E. nitens (Fig. S8a, Table S5). For poplar and willow, the model generally overestimated yields of Populus deltoides × P. nigra and P. deltoides but underestimated yields of P. trichocarpa and Salix schwerinii × S. viminalis (Fig. S8b, Table S5). There is an underestimation of yields for Miscanthus × giganteus but overestimation for Miscanthus sinensis. In fact, the observed yields of the former are significantly higher than yields of the latter (t test, p < 0.01). Only four sites reported yields for Panicum pretense, and they were overestimated by the model (Fig. S8d, Table S5).

4.4 Future development

Although the model can generally reproduce the bioenergy crop yields on a global scale, there are still some regional biases of biomass yields for different bioenergy crops. For example, ORCHIDEE-MICT-BIOENERGY underestimates the biomass yields of Miscanthus in the UK by 43 % (Fig. 8) and overestimates the yields of switchgrass in the eastern US by 18 % (Fig. 9). Thus, for regional use of modeled results, slight modifications of the related parameters would be needed.

In addition to the yields from aboveground biomass, the allocation of belowground biomass also needs to be modified, and the resulting soil carbon stocks need to be evaluated. In the current version, the non-harvested parts of biomass go to the litter pool after each harvest. In reality, however, stumps and coarse roots remain alive in coppicing practices of tree species like eucalypt, poplar, and willow, and new shoots grow out of these stumps in the next growing season. Similarly, new shoots grow out of rhizome for perennial grasses like Miscanthus in the next growing season after harvest. Carbon in such live biomass compartments does not transfer to the litter or soil and thus does not contribute to soil carbon stocks. It is necessary to correct the model processes in this respect before applying this model to account for the full carbon cycle involving bioenergy plants. Meanwhile, a global observation dataset of belowground biomass and soil organic carbon for bioenergy crops would be desirable to systematically evaluate the model, but this does not exist to the best of our knowledge. For a long-term perspective, the implementation of explicit managements and interactions between bioenergy yields and nutrient limitations are increasingly important to simulate the carbon reduction potentials of bioenergy crop deployments.

Beside the biogeochemical processes, it is also critical to further parameterize and evaluate biophysical processes, especially in the coupled simulations of global vegetation models with climate models to calculate the biophysical feed-backs. Field measurements on, e.g., leaf traits, heat exchange, and the transpiration of bioenergy crops extend our knowledge of these biophysical processes and need to be integrated adequately in the global vegetation models.

5 Conclusions

Bioenergy crops have been extensively assumed in IAMs and are an important type of future land use. However, most global vegetation models do not have specific representations of these bioenergy crops. It is important to accurately represent the physiology, phenology, and carbon allocation of these crops because it fundamentally impacts the hydrology dynamics, energy balance, and carbon cycle. Especially for woody bioenergy crops like eucalypts, poplars, and willows, not only the biomass yields but also the seasonal variations, biophysical effects, and carbon turnover are impacted by new parameterizations.

In this study, we demonstrated the importance of the proper representative of bioenergy crops in a global vegetation model to reproduce the observation-based biomass yields. We introduced new bioenergy crop PFTs based on their plant characteristics, modified the parameters relevant to productivity based on field measurements and empirical evidence, and added the dedicated harvest process to simulate bioenergy biomass yields. The bioenergy crop simulations in ORCHIDEE-MICT-BIOENERGY generally reproduced the observation-based biomass yields for bioenergy crops at global level. However, it is still difficult to match observations site by site due to the uncertainties in the observation dataset and the lack of explicit managements in the model. Evaluations on soil carbon dynamics and biophysical variables are further needed. Our work improves the performance of ORCHIDEE on bioenergy crop modeling, and the parameters used in ORCHIDEE-MICT-BIOENERGY also provide guidance for other vegetation models on incorporating dedicated bioenergy crops.