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

Abstract

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 CO₂, 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.

1. 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 FAO-AEZ model also only produces yields at maturity, and not the dynamic daily changes in yields, which restricts its applications [12]. The second tef model, the FAO’s AquaCrop model, simulates biomass production using the daily crop transpiration and the normalized crop water productivity [13]. A crop specific harvest index is then used to calculate the grain yield based on the simulated biomass production [14]. The AquaCrop tef model was designed for water limited growing scenarios and therefore does not directly take soil nutrients and fertilizer management practices into account [13]. This limits its application in Ethiopia, where soil degradation [15] and low levels of fertilizer use [16] cause nutrient limited growing conditions. In addition, the AquaCrop tef model was not developed to take differences in daylength into account [1,13,17,18,19], limiting its application to low latitudes.

The Decision Support System for Agrotechnology Transfer (DSSAT) is a globally used crop modeling platform that works on a daily time step and accounts for the effects of nitrogen, CO₂, and daylength on crop growth and development [20,21]. The goal of this paper was to calibrate an existing, well-tested wheat crop module for tef in the DSSAT platform to create a DSSAT-Tef model using published field data for tef and other C4 crops, such as sorghum.

2. Materials and Methods

2.1. Model Choice

The tef crop model was created using the DSSAT-NWheat model as a basis. The NWheat model was selected, instead of a C4 crop model like DSSAT-Sorghum, as a starting point, because of its robust model structure, which has been successfully applied to many different growing conditions around the world [22]. The parameters and structures of DSSAT-NWheat were preserved unless the literature suggested a significant difference between tef and wheat. The major differences between tef and wheat are that tef is a short-day [6] C4 crop [23], while wheat is a long-day [24] C3 crop [25]. In order to adapt the wheat model for tef, changes were made to the phenology, photoperiod response, radiation use efficiency, and transpiration efficiency for both standard and elevated atmospheric CO₂ based on published research for tef and other C4 crops, such as sorghum, maize, and millet. The Agricultural Production Systems Simulator (APSIM)-Sorghum model [26], in particular, was referenced when there were gaps in the tef literature, as sorghum is a C4 crop, like tef, and the parameters were often already in the correct format for the DSSAT-NWheat model. This paper uses model specific nomenclature, but the DSSAT parameter names are explained so that a reader with a general understanding on crop modeling can understand them. For more detailed descriptions of the general DSSAT model structure and of the DSSAT-NWheat model, the reader should refer to Jones et al. [21] and Kassie et al. [22].

2.2. Field Locations

The DSSAT-Tef simulations were compared to observed field data from 9 Ethiopian locations and one location in the United States (

Table 1. Source of experimental data for tef model calibration.

Location Name	Latitude/Longitude	Elevation (m a.s.l.)	Soil Type	Avg. Yearly Rainfall (mm)	Year	Source
Mekelle 1, Ethiopia	13°3′ N/39°6′ E	2130 1	Cambisol	600	2008	[13]
Mekelle 2, Ethiopia	13°3′ N/39°6′ E	2130 1	Cambisol	600	2009	[13]
Ilala, Ethiopia	13°4′ N/39°4′ E	NA	Vertisol	650	2008	[13]
Fallon, NV, USA	39°29′20″ N/118°49′41″ W 2	1215 2	Dia Loam Soil	127 2	2009	[29]
Jamma District, Ethiopia	10°27′ N/39°15′ E	2630	Vertisol	868.2	2011	[16]
Ofla District, Ethiopia	12°31′58″ N/39°30′13″ E	2450	Vertisol	980.5	2011	[31]
Ude 1, Ethiopia	8°48′ N/39°38′ E	1800	Vertisol	775 3	1997	[27]
Kajima, Ethiopia	8°48′ N/39° 38′ E	1920	Andisol	775 3	1997	[27]
Gare Arera 1, Ethiopia	09°03′ N/38° 30′ E	NA	Nitisol	1046 4	2002	[35]
Gare Arera 2, Ethiopia	09°03′ N/38°30′ E	NA	Vertisol	1046 4	2002	[35]

¹ From Araya et al. [1]. ² 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]. ³ Annual rainfall in 1997. ⁴ Annual rainfall in 2002.

). The Ude 2 location from Tulema et al. [27] was discarded due to reported high weed pressure, which the model could not simulate as it only accounts for abiotic stresses. Many tef field experiments from the literature could not be used for model calibration, because they lacked sufficient data on the environmental conditions or management practices. Due to limited field data, there were no separate calibration and validation datasets.

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].

Mekelle 1, Ethiopia						Mekelle 2, Ethiopia
Cambisol						Cambisol
Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5	Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5
0.05	0.14	0.27	0.46	0.72 7	1	0.05	0.14	0.27	0.46	0.72 7	1
0.10	0.14	0.27	0.46	0.72 7	1	0.10	0.14	0.27	0.46	0.72 7	1
0.20	0.14	0.27	0.50	0.60 7	1	0.20	0.14	0.27	0.50	0.60 7	1
0.30	0.22	0.42	0.52	0.60 7	0.61	0.30	0.22	0.42	0.52	0.60 7	0.61
0.40	0.22 6	0.42 6	0.52 6	0.60 6	0	0.40	0.22 6	0.42 6	0.52 6	0.60 6	0
0.50	0.22 6	0.42 6	0.52 6	0.60 6	0	0.50	0.22 6	0.42 6	0.52 6	0.60 6	0
0.60	0.22 6	0.42 6	0.52 6	0.60 6	0	0.60	0.22 6	0.42 6	0.52 6	0.60 6	0
0.70	0.22 6	0.42 6	0.52 6	0.60 6	0	0.70	0.22 6	0.42 6	0.52 6	0.60 6	0
0.80	0.22 6	0.42 6	0.52 6	0.60 6	0	0.80	0.22 6	0.42 6	0.52 6	0.60 6	0
0.90	0.22 6	0.42 6	0.52 6	0.60 6	0	0.90	0.22 6	0.42 6	0.52 6	0.60 6	0
Ilala, Ethiopia						Fallon, NV, USA
Vertisol						Dia Loam
Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5	Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5
0.05	0.19	0.38	0.54	1.02 9	1	0.05	0.16 8	0.17 10	0.47 8	1.16 10	1
0.10	0.19	0.38	0.54	1.02 9	1	0.10	0.16 8	0.17 10	0.47 8	1.16 10	1
0.20	0.19	0.38	0.54	1.02 9	1	0.20	0.14 8	0.17 10	0.45 8	0.66 10	1
0.30	0.19	0.38	0.54	1.02 9	0.61	0.30	0.14 8	0.17 10	0.45 8	0.66 10	0.61
0.40	0.19 6	0.38 6	0.54 6	0.90 9	0	0.40	0.14 8	0.17 10	0.45 8	0.66 10	0.50
0.50	0.19 6	0.38 6	0.54 6	0.90 9	0	0.50	0.14 8	0.17 10	0.44 8	0.44 10	0.41
0.60	0.19 6	0.38 6	0.54 6	0.90 9	0	0.60	0.10 8	0.12 10	0.44 8	0.44 10	0.33
0.70	0.19 6	0.38 6	0.54 6	0.90 9	0	0.70	0.10 8	0.12 10	0.44 8	0.44 10	0.27
0.80	0.19 6	0.38 6	0.54 6	0.87 9	0	0.80	0.10 8	0.12 10	0.44 8	0.44 10	0.22
0.90	0.19 6	0.38 6	0.54 6	0.87 9	0	0.90	0.10 8	0.12 10	0.44 8	0.44 10	0.18
2.10	0.19 6	0.38 6	0.54 6	0.67 9	0	2.10	0.09 8	0.12 6	0.40 8	0.39 6	0.02
Jamma District, Ethiopia						Ofla District, Ethiopia
Vertisol						Vertisol
Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5	Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5
0.05	0.28 11	0.51 11	0.54 11	1.02	1	0.05	0.26 8	0.44 8	0.51 8	2.55	1
0.10	0.28 11	0.51 11	0.54 11	1.02	1	0.10	0.26 8	0.44 8	0.51 8	2.55	1
0.20	0.28 11	0.51 11	0.54 11	1.02	1	0.20	0.26 8	0.44 8	0.51 8	2.55	1
0.30	0.28 11	0.51 11	0.54 11	1.02	0.61	0.30	0.29 8	0.40 8	0.49 8	0.06	0.61
0.40	0.26 11	0.48 11	0.51 11	0.90	0	0.40	0.29 8	0.40 8	0.49 8	0.06	0.50
0.50	0.26 11	0.48 11	0.51 11	0.90	0	0.50	0.29 8	0.40 8	0.49 8	0.06	0.41
0.60	0.26 11	0.48 11	0.51 11	0.90	0	0.60	0.29 8	0.40 8	0.49 8	0.06	0.33
0.70	0.26 11	0.48 11	0.51 11	0.90	0	0.70	0.29 8	0.40 8	0.49 8	0.06	0.28
0.80	0.21 11	0.45 11	0.48 11	0.87	0	0.80	0.29 8	0.40 8	0.49 8	0.06	0.22
0.90	0.21 11	0.45 11	0.48 11	0.87	0	0.90	0.29 8	0.40 8	0.49 8	0.06	0.18
Ude 1, Ethiopia						Kajima, Ethiopia
Vertisol						Andisol
Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5	Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5
0.05	0.37 8,9	0.50 8,9	0.52 8,9	1.25	1	0.05	0.25 7,8	0.41 7	0.49 7	1.91	1
0.10	0.37 8,9	0.50 8,9	0.52 8,9	1.25	1	0.10	0.25 7,8	0.41 7	0.49 7	1.91	1
0.20	0.37 8,9	0.50 8,9	0.52 8,9	1.25	1	0.20	0.25 7,8	0.41 7	0.49 7	1.91	1
0.30	0.37 8,9	0.50 8,9	0.52 8,9	1.256	0.61	0.30	0.25 7,8	0.41 7	0.49 7	1.91	0.61
0.40	0.39 8,9	0.49 8,9	0.54 8,9	1.256	0	0.40	0.25 6	0.41 6	0.49 6	1.91 6	0
0.50	0.39 8,9	0.49 8,9	0.54 8,9	1.256	0	0.50	0.25 6	0.41 6	0.49 6	1.91 6	0
0.60	0.39 8,9	0.49 8,9	0.54 8,9	1.256	0	0.60	0.25 6	0.41 6	0.49 6	1.91 6	0
0.70	0.39 8,9	0.49 8,9	0.54 8,9	1.256	0	0.70	0.25 6	0.41 6	0.49 6	1.91 6	0
0.80	0.39 8,9	0.50 8,9	0.52 8,9	1.256	0	0.80	0.25 6	0.41 6	0.49 6	1.91 6	0
0.90	0.39 8,9	0.50 8,9	0.52 8,9	1.256	0	0.90	0.25 6	0.41 6	0.49 6	1.91 6	0
Gare Arera 1, Ethiopia						Gare Arera 2, Ethiopia
Nitisol						Vertisol
Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5	Depth (m)	LL 1	DUL 2	SAT 3	OC 4	SRGF 5
0.05	0.37 8,12	0.45 8,12	0.56 8,12	3.77 12	1	0.05	0.38 8,13	0.53 8,13	0.55 8,13	2.59 13	1
0.10	0.37 8,12	0.43 8,12	0.56 8,12	3.77 12	1	0.10	0.38 8,13	0.53 8,13	0.55 8,13	2.59 13	1
0.20	0.37 8,12	0.45 8,12	0.56 8,12	3.77 12	1	0.20	0.38 8,13	0.53 8,13	0.55 8,13	2.59 13	1
0.30	0.37 8,12	0.45 8,12	0.54 8,12	2.88 12	0.61	0.30	0.37 8,13	0.52 8,13	0.53 8,13	2.05 13	0.61
0.40	0.37 8,12	0.44 8,12	0.54 8,12	2.88 12	0	0.40	0.37 8,13	0.52 8,13	0.53 8,13	2.05 13	0
0.50	0.34 8,9	0.46 8,12	0.52 8,12	2.20 12	0	0.50	0.37 8,13	0.52 8,13	0.54 8,13	2.09 13	0
0.60	0.34 8,9	0.46 8,12	0.52 8,12	2.20 12	0	0.60	0.37 8,13	0.52 8,13	0.54 8,13	2.09 13	0
0.70	0.34 6	0.45 8,12	0.52 6	2.20 6	0	0.70	0.37 6	0.52 6	0.54 6	2.09 13	0
0.80	0.34 6	0.46 8,12	0.52 6	2.20 6	0	0.80	0.37 6	0.52 6	0.54 6	2.09 13	0
0.90	0.34 6	0.46 8,12	0.52 6	2.20 6	0	0.90	0.37 6	0.52 6	0.54 6	2.09 13	0

¹ Soil lower limit, also known as wilting point (%volume/100). ² Soil drained upper limited, also known as field capacity (%volume/100). ³ Soil saturation point (%volume/100). ⁴ Soil organic carbon (%). ⁵ Soil root growth factor. ⁶ This layer is a copy of the layer above it. ⁷ Based on Habtegebrial et al. [37]. ⁸ This value was calculated using the DSSAT SBuild program [20]. ⁹ Based on the Vertisol description in Getu [16]. ¹⁰ 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]. ¹¹ 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]. ¹² Based on a Gare Arera Nitisol described in Tulema et al. [39]. ¹³ Based on a Gare Arera Vertisol described in Tulema et al. [39].

lists the soil characteristics for the experiments used. In cases where the drainage parameters were not given, they were calculated in DSSAT’s SBuild program [20] using the soil clay, silt, and sand percentages. All soils were standardized to a depth of 90 cm with the same layer depths for each soil. In cases where the properties of the soil layer were unknown, it was assumed that they were the same as the layer above. Root growth for all locations, except for the Fallon, NV, and Ofla District locations, was limited to the top 30 cm of the soil, by setting the root growth factor to 0 below 30 cm. This was based on Araya et al.’s [1] findings that the root depth for early maturing tef varieties was generally less than 30 cm. A 90 cm maximum root depth was assumed for the Fallon, NV, and Ofla District locations due to the later maturity [28,29], which is associated with deeper roots for tef [30]. The high yields at these locations could not have been achieved with shallow roots [29,31].

Table 3. Crop management data.

Location	Planting Date	Plant Density	N Fertilizer Application Dates	N Fertilizer Application (kg N/ha)	Irrigation	Cultivar Name
Mekelle 1, Ethiopia	August 2	900 1	50% at planting and 50% 30 DAP 2	60 3	0,20,47,71, and 95 mm	DZ-974 and Keyh (local)
Mekelle 2, Ethiopia	July 31	900 1	50% at planting and 50% 30 DAP 2	60 3	0, 40, and 90 mm	Keyh (local)
Ilala, Ethiopia	August 5	900 1	50% at planting and 50% 30 DAP2	60 3	No	DZ-974 and Keyh (local)
Fallon, NV, USA	June 4	900 1	No fertilizer applied	0	Yes	Dessie, U.K. Brown, U.K. White, PI 193508, PI 193514, PI 195932, PI 273889, PI 329680, PI 347632, PI 494366, PI 494432, PI 494433, PI 494465, PI 494479, PI 557457
Jamma District, Ethiopia	July 244	900 1	At planting	0, 23 5, 46 5, 69 5	No	Wajera (local)
Ofla District Ethiopia	July 244	900 1	At planting	0, 23 5, 46 5, 69 5	No	DZ-Cr-387 (Kuncho)
Ude, Ethiopia	July 22	900 1	On day of planting 6	64 3	No	DZ-01-354
Kajima, Ethiopia	July 22	900 1	On day of planting 6	64 3	No	DZ-01-354
Gare Arera 1, Ethiopia	July 24	900 1	Before sowing	40 7	No	DZ-01-354
Gare Arera 2, Ethiopia	July 25	900 1	Before sowing	60 7	No	DZ-01-354

¹ Rounded from the 923 plants/m² used by Araya et al. [1]. ² Days after planting. ³ The specific type of fertilizer was not specified. It was assumed to be conventional urea. ⁴ Based on Settie et al. [40]. ⁵ As urea N super granule (USG). ⁶ Assumed, as it was not specified in Tulema et al. [27]. ⁷ As a combination of conventional urea and diammonium phosphate.

specifies the management practices used for the simulations. The DSSAT-Tef model uses post-emergence plant density (plants/m²), 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].

2.3. 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. Initial estimated soil water and mineral nitrogen conditions for each location at the start of the simulation. Simulations started one month before planting.

Location	Initial Soil Water (%FC 1)	Initial Soil N (kg N/ha)
Mekelle 1, Ethiopia	100 2	10
Mekelle 2, Ethiopia	100 2	10
Ilala, Ethiopia	100 2	10
Fallon, NV, USA	100	40
Jamma District, Ethiopia	10	30
Ofla District Ethiopia	10	150 3
Ude, Ethiopia	10	20
Kajima, Ethiopia	10	20
Gare Arera 1, Ethiopia	40	10
Gare Arera 2, Ethiopia	40	10

¹ Field capacity. ² After Araya et al. [13]. ³ 91% of soil N was concentrated in the top 30 cm.

). Simulations started one month before planting to allow surface soil water to reach the equilibrium it would have in the actual field.

2.4. 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. Temperature parameters for phenology for different crop models.

Model	Reference	TBASE 1	TOPT 2	TTOPT 3
DSSAT-NWheat	[20]	0	26	34
APSIM-Sorghum	[26]	11	30	42
APSIM-Millet	[43]	10	33	47
DSSAT-Sorghum	[20]	8	34	NA
DSSAT-Millet	[20]	10	34	NA
DSSAT-Tef		7.8	30	42

¹ Base daily mean temperature below which no crop phenology occurs (°C). ² Optimum daily mean temperature at which the maximum rate of crop phenology occurs (°C). ³ Maximum daily mean temperature at which crop phenology stops (°C).

), 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).

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. Cultivar parameters in the DSSAT model.

Cultivar Name	PPSEN 1	P1 2	P5 3	PHINT 4	GRNO 5	MXFIL 6	STMMX 7
Dessie	16.5	100	500	80	32	1.2	1.0
Keyh (local)	18.0	100	600	110	16	1.0	3.0
U.K. Brown	16.5	100	450	100	16	1.0	1.0
U.K. White	18.0	100	450	85	16	1.0	1.0
Wajera (local)	18.0	100	200	80	36	1.0	3.0
Combined DZ-974 and Keyh (local)	18.0	100	500	105	16	1.0	1.0
DZ-01-354	18.0	100	550	100	16	1.0	3.0
DZ-Cr-387	18.0	100	100	120	32	1.3	3.0
PI 193508	19.0	100	200	100	30	1.1	1.0
PI 193514	19.0	100	200	105	16	1.0	1.0
PI 195932	19.0	100	200	105	16	1.0	1.0
PI 273889	19.0	100	190	105	27	1.0	1.0
PI 329680	19.0	100	200	105	16	1.0	1.0
PI 347632	19.0	100	190	105	32	1.8	1.0
PI 494366	14.0	100	700	70	16	1.0	1.0
PI 494432	15.5	100	500	100	26	1.0	1.0
PI 494433	19.0	100	140	120	16	1.0	1.0
PI 494465	19.0	100	140	120	16	1.0	1.0
PI 494479	19.0	100	200	105	16	1.0	1.0
PI 557457	17.0	100	600	75	16	1.0	1.0

¹ Sensitivity to photoperiod factor. ² Thermal time from seedling emergence to the end of tillering (juvenile phase). ³ Thermal time from beginning of grain fill to maturity. ⁴ Phyllochron interval. ⁵ Coefficient of kernel number per stem weight at the beginning of grain filling (100 kernels/g stem). ⁶ Potential kernel growth rate (mg/kernel/day). ⁷ Potential final dry weight of a single tiller (excluding grain) (g/stem).

). These parameters were individually changed for each cultivar in order to fit observed phenology data.

2.5. Photoperiod Response

DSSAT-NWheat uses the parameter nwheats_ppfac to account for photoperiod response. The NWheat equation was changed to the following for tef:

$${\mathrm{tef}}_{\mathrm{ppfac}}=1-0.002\ast \mathrm{PPSEN}\ast {\left(\mathrm{TWILEN}-11\right)}^2$$

(1)

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 (Table 1 and Table 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.

2.6. 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² 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.

2.7. CO₂ Response

Though there were no published field experiments on the effects of elevated atmospheric CO₂ on tef, it is likely that, as a C4 crop, tef would have a less pronounced response to elevated CO₂ than wheat. It has been found that transpiration efficiency increases with increasing CO₂ levels for C4 crops, while transpiration efficiency and RUE increase for C3 crops [48]. DSSAT considered atmospheric CO₂ concentration as an input [20]. DSSAT-NWheat had two major pathways for dealing with the effects of elevated CO₂, transp_eff_coeff, and rue_factor. The coefficient transp_eff_coeff is the transpiration efficiency coefficient and is calculated using:

transp_eff_coeff = 0.009 * (1 + (0.28/350) * (WEATHER% CO₂ − 350))

(2)

where WEATHER% CO₂ is the actual atmospheric CO₂ 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₂ 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₂ 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₂ 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₂ response was used for DSSAT-Tef.

As RUE for C4 crops are not influenced by increasing atmospheric CO₂ [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.

2.8. 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.

2.9. 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², the root mean squared error (RMSE), and the relative root mean squared error (rRMSE). The r² 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.

3. Results

The performance of the model compared with observed data is summarized in

Table 7. Summary of the DSSAT-Tef model performance at Mekelle, Ilala, Jamma District, Ofla District, Ude, Kajima, and Gare Arera in Ethiopia, and Fallon, NV, USA.

Model Attribute	Number of Paired Data Points	Observed Range	R1,2	M2	RMSE 3	rRMSE 4
Emergence (DAP) 5	9	5–6	1.00	0.876	0.8	14.4
Anthesis (DAP) 5	34	44–98	0.98	1.005	2.5	3.6
Maturity (DAP) 5	34	71–150	0.96	1.008	4.4	3.9
Total Aboveground Biomass (kg/ha)	23	1400–9652	0.80	1.396	2624	49.6
Grain yield (kg/ha)	38	352–3013	0.74	1.159	475	41.0
Crop N Uptake (kg/ha)	10	20–206	0.81	0.623	45	46.2
Grain N (kg/ha)	12	4.3–92	0.79	0.794	15	37.3

¹r² for the 1 to 1 line of simulated versus observed data. ² The slope m of a linear regression (forced through origin) is shown as a quantitative indication of over- or under-estimation of a model’s performance. ³ Root mean squared error. ⁴ Relative root mean squared error. ⁵ Days after planting.

. 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²) 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² of 0.98, an RMSE of 2.5 days, and days to maturity with an r² 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² 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² 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 (Figure 4 and Figure 5). The model consistently captured the increasing trend in grain yield in response to increasing nitrogen and irrigation applications (Figure 4 and Figure 5).

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² 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.

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² 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.

4. 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₂ 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 rRMSE for grain yield for the DSSAT-Tef model was 41.0%. DSSAT-Tef’s grain yield rRMSE was higher than the average 13% rRMSE of the DSSAT-NWheat model [22] and the 20.3% rRMSE of the APSIM-Sorghum model [56]. DSSAT-Tef’s grain yield rRMSE was, however, comparable to the 39.5% rRMSE of the DSSAT-Sorghum model [56].

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].

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₂, 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₂ 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.

5. 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₂ concentration. While the impact of elevated atmospheric CO₂ concentrations on tef growth was based on a general understanding of crop response to CO₂, 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.