A Crop Simulation Model for Tef (Eragrostis tef (Zucc.) Trotter)

Tef is an Ethiopian staple grain that provides both food security and income for smallholders. As tef is nutritious and gluten free, it is also gaining popularity as a health food. A tef model was calibrated based on the Decision Support System for Agrotechnology Transfer’s (DSSAT) NWheat model and included parameter changes in phenology, photoperiod response, radiation use efficiency, and transpiration efficiency for both standard and elevated atmospheric CO2, based on published literature for tef and other C4 species. The new DSSAT-Tef model was compared with tef field experiments. DSSAT-Tef accurately simulated phenology and responded to changes in N supply and irrigation, but overestimated growth and occasionally yields. Simulation-observation comparisons resulted in an RMSE of 2.5 days for anthesis, 4.4 days for maturity, 2624 kg/ha (49.6%) for biomass, and 475 kg/ha (41.0%) for grain yield. Less data were available for N uptake, and the model simulated crop N uptake with an RMSE of 45 kg N/ha (46.2%) and 15 kg N/ha (37.3%) for grain N. While more data from contrasting environments are needed for further model testing, DSSAT-Tef can be used to assess the performance of crop management strategies, the suitability of tef for cultivation across growing environments, and food security.


Introduction
Tef, or Eragrostis tef (Zucc.) Trotter, is an Ethiopian staple grain [1]. A strong preference for tef grain in Ethiopia means that smallholder farmers can demand a higher price for tef than for other grains, making tef an important source of income [2]. Tef's high calcium and iron content [3,4], as well as its lack of gluten [5] have resulted in a growing demand for tef as a health food in industrialized nations. In both Ethiopia and abroad, tef straw is valued as a high quality, low input, warm season fodder [6,7]. Crop models are computer programs that simulate crop growth and yields, under varying environmental conditions and management practices [8]. Crop models can be used to virtually test management practices and crop traits using a fraction of the time and money needed to run a field trial [9]. Crop models are particularly important for assessing the effects of climate change on agricultural production and for evaluating food security [10,11].
Two published crop models currently exist for tef. The Food and Agriculture Organization (FAO)-AEZ crop growth simulation model estimates yields using a three step process [12]. First, the model calculates the radiation limited yield, then the water limited yield, and finally, the yield affected by soil and management limitations. Each yield is calculated based on the yield of the previous step [12]. The FAO-AEZ model has certain limitations. As growing season length is treated as an input, the model cannot simulate the effects of temperature, or day length, on phenology [12]. The yields at these locations could not have been achieved with shallow roots [29,31]. Table 3 specifies the management practices used for the simulations. The DSSAT-Tef model uses post-emergence plant density (plants/m 2 ), rather than the more commonly reported seeding rate (kg seed planted/ha), so the plant density had to be assumed for most sites (Table 3).
A combination of satellite and station weather data was used. The Climate Hazards Group InfraRed Precipitation with Station (CHIRPS) weather dataset [32] was used for the Jamma District location and for the Ude, Kajima, and both Gare Arera locations. A combination of CHIRPS temperature data and observed rainfall data [13] was used for the Mekelle and Ilala locations. The NASA Prediction of Worldwide Energy Resource (POWER) dataset [33] was used for the Ofla District location. The Fallon Experimental, NV, USA. Weather Station was used for the Fallon, NV, location, though missing data were filled in using data from the Fallon NAAS, NV, USA. Weather Station [34]. 1 From Araya et al. [1]. 2 From a Dia loam, 0 to 1 percent slopes group from the Fallon-Fernley Area, Nevada, parts of Churchill, Lyon, Storey, and Washoe Counties' soil survey [36]. 3 Annual rainfall in 1997. 4 Annual rainfall in 2002.

Model Initialization
Not all locations reported the same types of data, so the observed sample size ranged from 9 to 38. As field data for tef were limited and published field experiments rarely gave full details on the environmental conditions and management practices, assumptions had to be made to fill in some of the gaps. Initial conditions were based on observed measurements. If these were not available, the initial conditions were calibrated by individually changing soil water and N to match the simulated yields with the observed yield of the lowest input (N or irrigation) treatment. The same initial conditions calibrated for the lowest input treatment were then applied to all other treatments within this experiment (Table 4). Simulations started one month before planting to allow surface soil water to reach the equilibrium it would have in the actual field. Table 2. Soil characteristics for a Cambisol soil at Mekelle in 2008 and 2009, a Vertisol soil at Ilala and a Dia Loam soil at Fallon, Vertisol soils in the Jamma and Ofla districts, a Vertisol soil at Ude 1, an Andisol soil at Kajima, and a Nitisol and a Vertisol soil at Gare Arera. Lower limit for plant available soil water (LL), drained upper limit (DUL), saturation (SAT), and soil organic carbon (OC). Data after Araya et al. [13], Davison and Laca [29], Getu [16], Okubay [31], Tulema et al. [27], and Tulema et al. [35].   6 2.09 13 0 1 Soil lower limit, also known as wilting point (%volume/100). 2 Soil drained upper limited, also known as field capacity (%volume/100). 3 Soil saturation point (%volume/100). 4 Soil organic carbon (%). 5 Soil root growth factor. 6 This layer is a copy of the layer above it. 7 Based on Habtegebrial et al. [37]. 8 This value was calculated using the DSSAT SBuild program [20]. 9 Based on the Vertisol description in Getu [16]. 10 From the Dia loam, 0 to 1 percent slopes group from the Fallon-Fernley Area, Nevada, parts of Churchill, Lyon, Storey, and Washoe Counties' soil survey [36]. 11 The saturation soil water content (SAT) of a soil layer was based on the soil porosity, which is estimated as 100-Bulk Density/2.65*100. To consider a general but small pore space that never saturates, 3 Vol% was subtracted from the soil porosity in each layer [38]. 12 Based on a Gare Arera Nitisol described in Tulema et al. [39]. 13 Based on a Gare Arera Vertisol described in Tulema et al. [39].

Phenology
General phenology information for tef was taken from the existing literature. The base temperature used to calculate developmental thermal time accumulation was increased from the 0 • C used in the wheat model [20] to 7.8 • C for tef [6]. The vernalization requirement was reduced to 0 in the ecotype file, as tef does not have a vernalization requirement [6].
There was no clear optimum temperature for tef phenology available in the literature. Reported optimal growing temperature ranges included 10-27 • C [2], 15-21 • C [41], and >26 • C [7]. It is unclear if these ranges referred to minimum, average, or maximum daily temperatures. Tsegay [42] wrote that the maximum temperature for tef production was 35 • C. As the literature did not specify for tef if the above temperatures were related to phenology or growth, the DSSAT-Tef model used the optimum and maximum temperatures from APSIM-Sorghum [26] (Table 5), which is a C4 crop, like tef [23]. APSIM-Sorghum, unlike DSSAT-Sorghum, also used the same three temperature parameters for calculating phenology as DSSAT-NWheat, on which the DSSAT-Tef model was based (Table 5). 1 Base daily mean temperature below which no crop phenology occurs ( • C). 2 Optimum daily mean temperature at which the maximum rate of crop phenology occurs ( • C). 3 Maximum daily mean temperature at which crop phenology stops ( • C).
Phenology differences between the individual cultivars were controlled using the photoperiod sensitivity parameter (PPSEN), the thermal time from seedling emergence to the end of tillering (juvenile phase) (P1) parameter, the phyllochron interval parameter (PHINT), and the thermal time from the beginning of grain fill to maturity (P5) parameter in the cultivar file (Table 6). These parameters were individually changed for each cultivar in order to fit observed phenology data.  2 Thermal time from seedling emergence to the end of tillering (juvenile phase). 3 Thermal time from beginning of grain fill to maturity. 4 Phyllochron interval. 5 Coefficient of kernel number per stem weight at the beginning of grain filling (100 kernels/g stem). 6 Potential kernel growth rate (mg/kernel/day). 7 Potential final dry weight of a single tiller (excluding grain) (g/stem).

Photoperiod Response
DSSAT-NWheat uses the parameter nwheats_ppfac to account for photoperiod response. The NWheat equation was changed to the following for tef: where tef_ppfac is a unitless photoperiod response factor; PPSEN is a unitless photoperiod sensitivity factor that varies between cultivars; and TWILEN is the daylength. For NWheat, the ppfac parameter increases with increasing daylength. Tef is a short-day plant [6], so the photoperiod response needed to be changed so that tef_ppfac decreased with increasing day length, and as a result, the rate of development decreased in response to longer days. Van Delden [6] defined the critical daylength for short day plants as the day length above which day length begins to impact the time to heading. The critical daylength of 20 h used in NWheat was reduced to 11 h, which was the average critical daylength across four tef cultivars [6]. The photoperiod sensitivity of each cultivar (PPSEN) was estimated by using phenology data from field trials (Tables 1 and 6). As tef is extremely photoperiod sensitive, the value for PPSEN used in the tef model was significantly higher than the typical range of 1-5 used by the NWheat model (Table 6). Figure 1 shows the response of tef_ppfac across multiple daylengths and PPSEN values.
development decreased in response to longer days. Van Delden [6] defined the critical daylength for short day plants as the day length above which day length begins to impact the time to heading. The critical daylength of 20 h used in NWheat was reduced to 11 h, which was the average critical daylength across four tef cultivars [6]. The photoperiod sensitivity of each cultivar (PPSEN) was estimated by using phenology data from field trials (Tables 1 and 6). As tef is extremely photoperiod sensitive, the value for PPSEN used in the tef model was significantly higher than the typical range of 1-5 used by the NWheat model (Table 6). Figure 1 shows the response of tef_ppfac across multiple daylengths and PPSEN values.

Growth
Though there were no radiation use efficiency (RUE) data available for tef, it was assumed that, as a C4 crop, tef would have a higher RUE than wheat. Models for the C4 crops maize, millet, and sorghum all have higher aboveground RUE values than wheat models [20,[44][45][46]. In order to calculate the DSSAT-Tef aboveground RUE, the percent difference between the DSSAT-CERES-Wheat and DSSAT-Sorghum aboveground RUE was used to estimate the ratio RUE ratio between wheat and a C4 crop. The percent change in RUE between the DSSAT-CERES-Wheat and DSSAT-Sorghum model was 18.5%. This percent increase was applied to the RUE of DSSAT-NWheat to estimate the RUE for tef. The RUE for DSSAT-NWheat was 3.2 g plant dry matter/MJ PAR, so the RUE value for DSSAT-Tef was set to 4.5 g/MJ. The NWheat model used a weighted mean temperature with emphasis on day time temperature [47]. To account for a higher base temperature for photosynthesis in tef, the threshold minimum temperature for photosynthesis was set to 4 degrees below the base temperature for phenology at 3.8 °C, the same difference, but for a different base temperature used in NWheat. These new parameter values were used in all tef simulations without any additional calibration.
The DSSAT-NWheat model initialized the grain weight at the beginning of grain filling using the initial grain weight parameter (INGWT). Growing conditions throughout the grain fill period and cultivar parameters, such as grain number (GRNO), potential kernel growth rate (MXFIL), photosynthesis, and available water soluble carbohydrates, all affect the change in kernel grain weight from this initial value. Tef grains are much smaller than wheat grains [31], so the initial grain weight needed to be decreased for the tef model. The first step to estimating the initial grain weight for DSSAT-Tef was calculating the ratio of the initial grain weight parameter used by the DSSAT-

Growth
Though there were no radiation use efficiency (RUE) data available for tef, it was assumed that, as a C4 crop, tef would have a higher RUE than wheat. Models for the C4 crops maize, millet, and sorghum all have higher aboveground RUE values than wheat models [20,[44][45][46]. In order to calculate the DSSAT-Tef aboveground RUE, the percent difference between the DSSAT-CERES-Wheat and DSSAT-Sorghum aboveground RUE was used to estimate the ratio RUE ratio between wheat and a C4 crop. The percent change in RUE between the DSSAT-CERES-Wheat and DSSAT-Sorghum model was 18.5%. This percent increase was applied to the RUE of DSSAT-NWheat to estimate the RUE for tef. The RUE for DSSAT-NWheat was 3.2 g plant dry matter/MJ PAR, so the RUE value for DSSAT-Tef was set to 4.5 g/MJ. The NWheat model used a weighted mean temperature with emphasis on day time temperature [47]. To account for a higher base temperature for photosynthesis in tef, the threshold minimum temperature for photosynthesis was set to 4 degrees below the base temperature for phenology at 3.8 • C, the same difference, but for a different base temperature used in NWheat. These new parameter values were used in all tef simulations without any additional calibration.
The DSSAT-NWheat model initialized the grain weight at the beginning of grain filling using the initial grain weight parameter (INGWT). Growing conditions throughout the grain fill period and cultivar parameters, such as grain number (GRNO), potential kernel growth rate (MXFIL), photosynthesis, and available water soluble carbohydrates, all affect the change in kernel grain weight from this initial value. Tef grains are much smaller than wheat grains [31], so the initial grain weight needed to be decreased for the tef model. The first step to estimating the initial grain weight for DSSAT-Tef was calculating the ratio of the initial grain weight parameter used by the DSSAT-NWheat model [20] and the final mature grain weight (mg/grain) at 0% grain moisture of an Ethiopian wheat field experiment [31]. This ratio was then applied to the final grain weight of tef at 0% moisture grown at the same location and under the same management conditions [31], in order to calculate the initial grain weight parameter for tef. Based on the estimated initial grain weights, INGWT was changed from 3.50 mg/grain for NWheat to 2.47 mg/100 grains for DSSAT-Tef. The tef model used mg/100 grains in its input file, but once the INGWT was read into the code, it was divided by 100 to get mg/grain.
The units for the grain number cultivar parameter in NWheat were the number of kernels/g stem weight at the beginning of grain filling and affected the rate of grain filling [20]. As tef has a much higher grain number per m 2 than wheat, the GRNO units in the cultivar file were changed to 100 kernels/g stem weight at the beginning of grain filling for the DSSAT-Tef model.
Lodging is a concern for tef production, especially for high input growing scenarios [37]. Due to a lack of lodging data in the available experiments, except for an indication of no or little occurrence of lodging in these datasets, no attempt was made to include lodging in the DSSAT-Tef model.

CO 2 Response
Though there were no published field experiments on the effects of elevated atmospheric CO 2 on tef, it is likely that, as a C4 crop, tef would have a less pronounced response to elevated CO 2 than wheat. It has been found that transpiration efficiency increases with increasing CO 2 levels for C4 crops, while transpiration efficiency and RUE increase for C3 crops [48]. DSSAT considered atmospheric CO 2 concentration as an input [20]. DSSAT-NWheat had two major pathways for dealing with the effects of elevated CO 2 , transp_eff_coeff, and rue_factor. The coefficient transp_eff_coeff is the transpiration efficiency coefficient and is calculated using: where WEATHER% CO 2 is the actual atmospheric CO 2 concentration in ppm. Hammer and Muchow [49] used a transpiration efficiency coefficient of 0.009 kPa for their sorghum simulation. As Hammer and Muchow [49] based their calculations only on aboveground biomass, their transpiration efficiency coefficient was compatible with the DSSAT-NWheat approach for modeling transpiration [20]. Pembleton et al. [50] used the APSIM-Wheat model as a basis for estimating the functions for modifying transpiration efficiency for various forage crops. Their function for forage sorghum was y = 0.0008x + 1. Based on Pembleton et al.'s [50] equation, a 350 ppm increase in CO 2 would result in a 28% increase in transpiration efficiency in sorghum. Teixera et al. [51] arbitrarily assumed a 37% increase in transpiration efficiency for maize and a 69% increase for wheat when atmospheric CO 2 was increased from 350 to 1032 ppm. Assuming a linear response, this would correspond to a 19% and 35% increase in transpiration efficiency for maize and wheat, respectively, if CO 2 were doubled from 350 to 700 ppm. As the calculations of Pembleton et al. [50] were based on field data, rather than an arbitrary assumption, their CO 2 response was used for DSSAT-Tef. As RUE for C4 crops are not influenced by increasing atmospheric CO 2 [48], the rue_factor was set to 1 for DSSAT-Tef. This was accomplished by setting the parameter RUEFAC in the ecotype file to 0, which set rue_factor to a constant of 1.

Cultivar Parameters
After the changes described above were made to the code, the individual cultivars were parameterized one at a time to account for differences in phenology and productivity. Table 6 gives the parameter values used for each cultivar. Cultivar parameters were calibrated in a multistep process, starting with the time to anthesis and followed by time to maturity, total aboveground biomass, and grain yield. PPSEN, P1, and PHINT affected time to anthesis. P5 affected the time to maturity. STMMX affected biomass yields. GRNO and MXFIL affected grain yields. Only one parameter was changed for each run, in order to avoid unintended interaction effects. Tef has high genetic diversity across cultivars. Dawit and Hirut [52] reported that days to 50% flowering ranged between 43 and 83 days across 506 tef accessions and that days to 50% maturity ranged from 93 to 130 days.
DSSAT did not output the days until heading, which means that observed anthesis dates were needed for calibration. Many of the publications used for the DSSAT-Tef model reported days until heading, rather than days until anthesis, however. For the Ofla District location [31], 10 days were added to the days to 50% heading to approximate the days until anthesis. This estimation was based on the difference between days to heading and days to flowering reported by Getu [16]. No anthesis information was reported for the Ude, Kajima, and both Gare Arera locations [27,35]. The observed days until anthesis were approximated by adding 10 days to the days until heading that were reported for the same cultivar grown under a 12 h daylength by Tadelle Tadesse [53]. The Fallon, NV, location [29] reported days to full seed head emergence, which was assumed to be equivalent to anthesis, so no days were added to the reported dates.

Statistics
All statistical analysis was done using the basic R package [54]. The simulated results were compared to the published observed data using the r 2 , the root mean squared error (RMSE), and the relative root mean squared error (rRMSE). The r 2 was calculated for the 1:1 line to quantify the departure of the model from the ideal reproduction of each observed value [55]. The RMSE and rRMSE were also calculated based on the simulated versus observed 1:1 line.

Results
The performance of the model compared with observed data is summarized in Table 7. Observed days from planting to emergence ranged from 5-6 days, but the observed data were only available for three locations. The model captured this small range very well, with a coefficient of determination (r 2 ) of 1.00, an RMSE of 0.8 days, and an rRMSE of 14.4%. Unlike for emergence, there was a large range in observed days from planting to anthesis and maturity, which ranged from 44-98 days and 71-150 days, respectively. The model accurately simulated days to anthesis, with an r 2 of 0.98, an RMSE of 2.5 days, and days to maturity with an r 2 of 0.96 and an RMSE of 4.4 days. Figure 2 shows the simulated versus observed days from planting to anthesis and days from planting to maturity for each location. There were no reported phenology data for Ude or Kajima, so the observed phenology for the same cultivar (DZ-01-354) from a different study [53] was used for calibration. As both Ude and Kajima shared the same coordinates, and therefore weather data, the sites were lumped together as the Ada Region when graphed.  Total aboveground biomass simulations were less accurate than the phenology simulations with an r 2 of 0.80 and an RMSE of 2625 kg/ha (rRMSE 49.7%). Figure 3a shows simulated versus observed total aboveground biomass for each location. The model tended to overestimate total aboveground biomass. Figure 4 indicates that the model consistently captured the increasing trend in total aboveground biomass in response to increasing applied nitrogen. The model overestimated simulated total aboveground biomass more severely for water stressed treatments, such as the zero irrigation treatments at Mekelle and Iala (leftmost red square and blue diamond in Figure 3a), than for irrigated treatments. Figure 5 shows that the model captured the increase in total aboveground biomass, in response to increasing irrigation amounts, for all treatments, except one.
Observed grain yield, at 0% grain moisture, varied from 352-3013 kg/ha. Grain yield simulations had an r 2 of 0.74 and an RMSE of 475 kg/ha (rRMSE 41.0%). Figure 3b shows simulated versus observed grain yield for each location. Grain yield was, in general, well simulated, except for some overestimations, such as the Ilala and irrigated Mekelle 2009 treatments (Figure 3b). There seemed to be no clear relationship between the amount of water or nitrogen applied and the accuracy of the grain yield simulation (Figure 6), suggesting that there was no model bias towards high or low input scenarios (Figures 4 and 5). The model consistently captured the increasing trend in grain yield in response to increasing nitrogen and irrigation applications (Figures 4 and 5). Total aboveground biomass simulations were less accurate than the phenology simulations with an r 2 of 0.80 and an RMSE of 2625 kg/ha (rRMSE 49.7%). Figure 3a shows simulated versus observed total aboveground biomass for each location. The model tended to overestimate total aboveground biomass. Figure 4 indicates that the model consistently captured the increasing trend in total aboveground biomass in response to increasing applied nitrogen. The model overestimated simulated total aboveground biomass more severely for water stressed treatments, such as the zero irrigation treatments at Mekelle and Iala (leftmost red square and blue diamond in Figure 3a), than for irrigated treatments. Figure 5 shows that the model captured the increase in total aboveground biomass, in response to increasing irrigation amounts, for all treatments, except one.      [16] and Okubay [31]. No error bars are available. Observed data from Getu [16] and Okubay [31]. No error bars are available.
Observed grain yield, at 0% grain moisture, varied from 352-3013 kg/ha. Grain yield simulations had an r 2 of 0.74 and an RMSE of 475 kg/ha (rRMSE 41.0%). Figure 3b shows simulated versus observed grain yield for each location. Grain yield was, in general, well simulated, except for some overestimations, such as the Ilala and irrigated Mekelle 2009 treatments (Figure 3b). There seemed to be no clear relationship between the amount of water or nitrogen applied and the accuracy of the grain yield simulation (Figure 6), suggesting that there was no model bias towards high or low input scenarios (Figures 4 and 5). The model consistently captured the increasing trend in grain yield in response to increasing nitrogen and irrigation applications (Figures 4 and 5).   [13], Getu [16], Okubay [31], and Tulema et al. [27,35]. No error bars are available.
Observed total aboveground nitrogen content ranged from 20-206 kg N/ha. There were no inseason measurements available. The model simulated total aboveground nitrogen with an r 2 of 0.80 and an RMSE of 45 kg N/ha (rRMSE 46.2%). Figure 7a shows simulated and observed total aboveground nitrogen for each location that had observed data. The model underestimated total aboveground nitrogen at higher N application treatments, but treatments with no added nitrogen  [13], Getu [16], Okubay [31], and Tulema et al. [27,35]. No error bars are available.
Observed total aboveground nitrogen content ranged from 20-206 kg N/ha. There were no in-season measurements available. The model simulated total aboveground nitrogen with an r 2 of 0.80 and an RMSE of 45 kg N/ha (rRMSE 46.2%). Figure 7a shows simulated and observed total aboveground nitrogen for each location that had observed data. The model underestimated total aboveground nitrogen at higher N application treatments, but treatments with no added nitrogen were close to the observed values. The model captured the increase in total aboveground nitrogen with increasing nitrogen applications.

Discussion
A new DSSAT-Tef model was created based on DSSAT-NWheat and information from the literature. Phenology, photoperiod response, radiation use efficiency, transpiration efficiency, and atmospheric CO2 response were modified to take the specifics of tef into account. The DSSAT-Tef model could reproduce observed phenology data and the observed yield responses to water and nitrogen, without any structural changes to the original model, using key crop parameters (like RUE and transpiration efficiency) from the literature without additional calibration and some calibration of a limited set of cultivar parameters. The range in parameter values was in accordance with the high variation in agronomic characteristics found in tef [2].
The comparison of the DSSAT-Tef model with a range of published field data showed an RMSE for total aboveground biomass of 2624 kg/ha and for grain yield of 475 kg/ha. The rRMSE for total aboveground biomass for the DSSAT-Tef model was 49.6%. DSSAT-Tef's biomass yield rRMSE was higher than the average 15% rRMSE of the DSSAT-NWheat model [22]. The APSIM-Sorghum and DSSAT-Sorghum models [56] also reported lower rRMSEs for total aboveground biomass (11.5% and 12.8%, respectively) than the DSSAT-Tef model.
The DSSAT-Tef model's grain yield RMSE was comparable to the two existing tef models, but the model overestimated biomass production more than the other two models [12,13,17,18]. The FAO-AEZ tef model was applied at 13 locations [12]. The FAO-AEZ grain yield RMSE was 402 kg/ha (rRMSE 29.1%). The FAO-AquaCrop tef model has been tested at multiple locations in multiple papers [13,17,18]. The AquaCrop final biomass RMSE ranged from 466 kg/ha to 720 kg/ha (rRMSE 23.6%) [13,18]. The AquaCrop grain yield RMSE ranged from 116 kg/ha to 540 kg/ha (rRMSE 19.1%) [13,17]. Observed grain nitrogen ranged from 4-92 kg N/ha. There were no in-season measurements available. The model simulated grain nitrogen with an r 2 of 0.78 and an RMSE of 15 kg N/ha (rRMSE 37.3%). Figure 7b shows simulated and observed grain nitrogen for each location that had observed data. The model did not consistently over-or under-estimate the grain nitrogen. While treatments with no applied nitrogen fertilizer gave more accurate simulations for the Jamma District and Gare Arera locations, the zero nitrogen treatment for the Ofla District location was more off than the fertilized treatments at that location. This was possibly due to the high initial N used for the Ofla District location.

Discussion
A new DSSAT-Tef model was created based on DSSAT-NWheat and information from the literature. Phenology, photoperiod response, radiation use efficiency, transpiration efficiency, and atmospheric CO 2 response were modified to take the specifics of tef into account. The DSSAT-Tef model could reproduce observed phenology data and the observed yield responses to water and nitrogen, without any structural changes to the original model, using key crop parameters (like RUE and transpiration efficiency) from the literature without additional calibration and some calibration of a limited set of cultivar parameters. The range in parameter values was in accordance with the high variation in agronomic characteristics found in tef [2].
The comparison of the DSSAT-Tef model with a range of published field data showed an RMSE for total aboveground biomass of 2624 kg/ha and for grain yield of 475 kg/ha. The rRMSE for total aboveground biomass for the DSSAT-Tef model was 49.6%. DSSAT-Tef's biomass yield rRMSE was higher than the average 15% rRMSE of the DSSAT-NWheat model [22]. The APSIM-Sorghum and DSSAT-Sorghum models [56] also reported lower rRMSEs for total aboveground biomass (11.5% and 12.8%, respectively) than the DSSAT-Tef model.
The biomass RMSE of DSSAT-Tef was higher than that of AquaCrop. The higher RMSE might be the result of DSSAT-Tef being tested with data that came from a greater range of conditions than the data used for AquaCrop. All AquaCrop test data came from fertilized experiments in northern Ethiopia [1,13,17,18] and did not include the variability caused by N stress and photoperiod response that was included in the DSSAT-Tef test data. The DSSAT-Tef grain yield RMSE fell within the RMSE range for AquaCrop, but the DSSAT-Tef grain yield rRMSE was higher than that of the AquaCrop model [13,17]. The DSSAT-Tef grain yield RMSE was slightly higher than that of the FAO-AEZ model, and the DSSAT-Tef's grain yield rRMSE was also higher than that of the FAO-AEZ model. A lack of variability in test data for the FAO-AEZ model could be why the grain yield RMSE of the FAO-AEZ model was comparable to that of DSSAT-Tef, despite the fact that the FAO-AEZ model assumed constant management practices across all locations and used seasonal averages and totals for weather data, rather than daily inputs [12].
Unlike the FAO-AEZ and AquaCrop tef models, DSSAT-Tef could model total aboveground N and grain N, which could be used to evaluate the fodder quality of tef straw and the protein content of tef grain. DSSAT-Tef could also respond to more environmental variation as it accounted for the effects of N, photoperiod, and atmospheric CO 2 , although the later was not tested with field data due to a lack of experiments. Further research is needed to improve the model's ability to simulate N demand and partitioning. Further research is also needed to confirm that tef has a similar response to elevated CO 2 as sorghum. If this research becomes available, the DSSAT-Tef model could be used to assess the impacts of climate change on tef production and food security in Ethiopia. More detailed field data might also confirm the explanations proposed below for some of the DSSAT-Tef model responses.
The DSSAT-Tef model overestimated the days to maturity at the Ilala location [13]. This in turn contributed to the overestimation of total aboveground biomass and grain yield. The reported phenology and yield values for both the Mekelle location in 2008 and the Ilala location were an average of two different cultivars, the improved DZ-974 and the local Keyh variety. It was not reported what the ratio of improved to local variety was for each location, and it is possible that the differences in observed phenology for each location were the result of a difference in the percentage of each cultivar present. The combination of the two cultivars was treated as a single cultivar when calibrating the cultivar parameters. As the same cultivars were used for each location, it was assumed that the cultivar parameters were the same for each location. Therefore, the discrepancy between observed and simulated data at the Ilala location could not be resolved by altering cultivar parameters.
The observed days to anthesis and maturity varied across fertilizer treatments for the Jamma and Ofla District locations [16,31]. The model did not account for this spread, however, as DSSAT used temperature and daylength to determine phenology and not nitrogen availability [20]. Neither publication could confirm the reason for the range in phenology, though it was proposed that tef reduced the time until heading when grown under lower N conditions as an avoidance response. By developing more rapidly, the plant would be able to reproduce before depleting the limited soil nutrients [31].
The simulated total aboveground biomass and grain yield at the Ofla District location followed a strong linear trend for the first three fertilizer treatments, but spiked for the highest N treatment. A possible explanation for this is that the last treatment suffered from lodging (although there was no indication in the report), for which the model did not account. Tef is particularly prone to lodging at fertilizer application rates above 60 kg N/ha [37], and the highest fertilizer rate at the Ofla District location was 69 kg N/ha [31]. Future research should consider including a lodging routine for the DSSAT-Tef model to address this issue.
Other causes for discrepancies between the model and the observed data included differences in initial conditions or crop management practices and biological stresses. Not all studies thoroughly described the initial conditions or crop management practices, so in some cases, assumptions had to be made, like on rooting depth, soil moisture, and plant density. The model only accounted for abiotic stresses, so if there were yield losses due to unreported abiotic stresses, the model would not capture them.

Conclusions
The DSSAT-Tef model offered new modeling capabilities beyond the two existing tef models FAO-AEZ and AquaCrop tef, by including photoperiod sensitivity, crop management, crop N dynamics, and the impact of elevated atmospheric CO 2 concentration. While the impact of elevated atmospheric CO 2 concentrations on tef growth was based on a general understanding of crop response to CO 2 , this was not evaluated with observations due to a lack of field data. Further research is also needed to improve N simulations and to add a lodging routine to the model. Additional field data are needed to do a complete evaluation of the model. The DSSAT-Tef model can be used to assess management practices and the suitability of regions for growing tef both in Ethiopia and other parts of the world.
Funding: K.P. received funding from the University of Florida's Graduate Assistantship Program.