Evaluation of Different Crop Models for Simulating Rice Development and Yield in the U.S. Mississippi Delta

The United States is one of the top rice exporters in the world, but warming temperatures and other climate trends may affect grain yield and quality. The use of crop models as decision support tools for a climate impact assessment would be beneficial, but suitability of models for representative growing conditions need to be verified. Therefore, the ability of CERES-Rice and ORYZA crop models to predict rice yield and growing season duration in the Mississippi Delta region was assessed for two widely-grown varieties using a 34-year database. CERES-Rice simulated growth duration more accurately than ORYZA as a result of the latter model’s use of lower cardinal temperatures. An increase in base and optimal temperatures improved ORYZA accuracy and reduced systematic error (e.g., correlation coefficient increased by 0.03–0.27 and root mean square error decreased by 0.3–1.9 days). Both models subsequently showed acceptable skill in reproducing the growing season duration and had similar performance for predicting rice yield for most locations and years. CERES-Rice predictions were more sensitive to years with lower solar radiation, but neither model accurately mimicked negative impacts of very warm or cold temperatures. Both models were shown to reproduce 50% percentile yield trends of more than 100 varieties in the region for the 34-year period when calibrated with two representative cultivars. These results suggest that both models are suitable for exploring the general response of multiple rice cultivars in the Mississippi Delta region for decision support studies involving the current climate. The response of rice growth and development to cold injury and high temperature stress, and variation in cultivar sensitivity, should be further developed and tested for improved decision making tools.


Introduction
Rice (Oryza sativa L.) is a staple food for roughly half the world's population and is the most important cereal food crop [1]. Continued yield increases are needed to meet the needs of a rising global population [2], but production is increasingly constrained in many agricultural areas by land and water availability [3][4][5] and frequency of heat stress episodes, where warmer temperatures occurring during critical stages of rice growth threaten yield stability [6][7][8]. The United States (U.S.) is the 4th largest rice exporter [9] and consistently produces the highest average yield per area in the world, in large part, with experimental data for rice grown in tropical climate conditions [43,44]. As compared with the U.S. Mississippi Delta, this climate has relatively smaller variability in seasonal and diurnal temperatures, which may require re-parameterization of the cardinal temperatures [45]. In contrast, the CERES-Rice model adequately captured the effects of temperature on phenology in the mid-south region of the U.S. [46].
Model ensemble approaches composed of more than one model may be more robust for predicting climate responses [21,22]. In this context, evaluation of both CERES-Rice and ORYZA was relevant for the topic of developing reliable decision support for the Mississippi Delta region. While the majority of rice modeling studies focus on varieties in the indica sub-population, rice growers in the U.S. primarily use temperate or tropical japonica sub-populations [6]. It was not clear whether the models used for Asian country studies, such as ORYZA, are suitable for the U.S. The objectives of this study were, thus, to (1) calibrate CERES-Rice and re-parameterize the cardinal temperatures in ORYZA to simulate growth duration for representative cultivars in the U.S. Mississippi Delta region, (2) evaluate and compare model performance in reproducing observed grain yields, and (3) assess the extent to which each model can duplicate yield trends for a broader set of cultivars grown in the Delta when calibrated with a limited number of varieties. The results of the study will establish the accuracy of utilizing well-known rice models for Southeast U.S., verify the extent to which such models are able to predict the growing season duration, yield trends in comparison with multi-year observed data, and provide insights into model knowledge gaps associated with climate effects.

Observed Data
The Uniform Regional Rice Nursery (URRN) dataset was utilized for this study [24] and included replicated annual field trial results from federal and university rice breeders across multiple states in the U.S. Mississippi Delta region, which covers a rice growing area of about 800,000 ha [4,5]. Trials were conducted in 3.3 m 2 field plots with two to four field replicates each year. Urea was applied in all plots as basal N fertilizer at rates ranging from 364 kg ha −1 to as high as 470 kg ha −1 for most years and locations, except for a lower application of 202 kg ha −1 in 2015 [24]. Trials were conducted that included 101 commercial varieties among others and included data on grain yield, days to heading, days to harvest, and plant height from 1983 to 2016. Not all varieties were grown each year as a result of new cultivar releases and decreased economic importance of older ones. However, cultivars Lemont, Cocodrie, Wells, Cypress, Priscilla, and CL161 were identified as the most well-represented varieties, and, of these, Lemont and Cocodrie were the most consistently grown at all locations and years. Based on availability of the most detailed data records, Lemont for years 1983-1998 and Cocodrie for 1999-2016 were used for this modeling activity at four locations (Stuttgart, Arkansas (AS), Crowley, Louisiana (LA), Stoneville, Mississippi (MS), and Beaumont, Texas (TX)).
Observed weather station data for daily maximum and minimum air temperature, solar radiation and rainfall were obtained from the following website: https://beaumont.tamu.edu/climaticdata/ StationInventory.aspx?index=Inventory. Solar radiation was not available at Crowley, LA and Stoneville, MS, and was obtained from a nearby weather station using the same website. Climate conditions of all four sites were classified as a U.S. warm temperate climate regime with humid and hot summers [47]. However, locational differences with respect to in-season climatology, growing season duration, and day/night temperatures existed (Table 1). April to August is the main growing season of rice at Beaumont, TX and Crowley, LA, and these locations were also more humid and wetter compared to Stuttgart, AR and Stoneville, MS, with average seasonal total rainfall of 744 and 1033 mm during 1983-2016, respectively. In contrast, Stuttgart, AR and Stoneville, MS were relatively dryer in terms of total average seasonal rainfall in the range of 424 to 467 mm from 1983 to 2016 (Table 1). Soil texture and hydrological properties data, such as saturated volumetric water content and hydraulic conductivity, were obtained from the Food Agriculture Organization soil dataset FAO-UNESCO (1974) [48] and Pachepsky and Park (2015) [49].

Development
CERES-Rice includes nine phenological stages in which the duration is calculated using degree days [50,51]. The daily increment in heat units (HU, • Cd (degree-days) d −1 ), or daily thermal time, are as in Equation (1).
where daily minimum (T min ) and maximum (T max ) air temperatures are above a base temperature (T base ) of 9 • C and less than an optimal temperature (T opt ) of 33 • C. When outside this range, HU is calculated by hourly temperature (T d ) at the time of day (h) as in Equations (2) and (3).

Biomass Growth
Potential biomass production for each day is obtained by Equation (4).
where PCARB is daily potential biomass production (g plant −1 d −1 ), RUE is radiation use efficiency (g MJ −1 ), PLTPOP is plant population (plant m −2 ), and IPAR is the fraction of photosynthetically active radiation (PAR, MJ) intercepted by the plants as: where k is the canopy light extinction coefficient and LAI is the green leaf area index. Potential biomass production may be limited by a non-optimal temperature, water stress, or nitrogen stress. Daily biomass is allocated among organs such as leaves, stems, roots, panicle, and grain according to carbon partitioning coefficients, which vary with the stage of development and general growing conditions [50].

Development
There are four main development stages determined by daily mean air temperature, water, nutrient stresses, and photoperiod for photoperiod-sensitive varieties. Cardinal temperatures include base (T base , • C), optimum (T opt , • C), and maximum temperatures (T high , • C), which default to 8, 30, and 42 • C, respectively, for all stages and cultivars [52]. The development rate has a bilinear response to temperature over short time periods (one hour). Hourly temperatures (T d , • C) are estimated by a sinusoidal function of daily minimum (T min , • C) and maximum (T max , • C) air temperatures.
where h is the time of day. Hourly increments in heat units (HUH, • Cd (degree-days) h −1 ) are given by Equation (7).
where T base , T opt , and T high are the cardinal temperatures as previously defined. The daily increment in heat units (HU, • Cd (degree-days) d −1 ) is then calculated as:

Biomass Growth
ORYZA uses a leaf-level photosynthesis algorithm to predict the CO 2 assimilation rate. The total daily rate is affected by daily absorbed radiation, temperature, and LAI. Radiation fluxes, attenuation, and absorption are based on the work of Spitters et al. (1986) [53]. The assimilation rate is calculated for a single shaded and sunlit leaf at each of three canopy depths using Equations (9)- (11).
where A L is the gross CO 2 assimilation rate (kg CO 2 ha −1 leaf h −1 ), A m is the CO 2 assimilation rate at light saturation, ε is the initial light-use efficiency (kg CO 2 ha −1 leaf h −1 (J m −2 leaf s −1 )), I a is the amount of absorbed radiation for either shaded or sunlit leaf (J m −2 leaf s −1 ), and ε340 ppm is ε at a CO 2 concentration of 340 ppm, which follows a function of average daytime temperature as in Goudriaan and van Laar (1994) [54]. Instantaneous rates of leaf CO 2 assimilation over three depths in the canopy are integrated to obtain total canopy CO 2 assimilation using a Gaussian integration at three different times during the day and then integrated into daily total rates of CO 2 assimilation [55].
Assimilate is partitioned between shoots and roots according to partitioning coefficients changing with phenological progress. The gross daily growth rate (G p , kg ha −1 d −1 ) is calculated according to Equation (12).
where A d is the daily rate of gross CO 2 assimilation (kg ha −1 d −1 ), R m is the maintenance respiration (kg ha −1 d −1 ), R t is the amount of available stem reserves for growth (kg ha −1 d −1 ), and Q is the assimilate requirement for dry matter production (kg CH 2 O kg −1 dry matter).

Model Input Data and Calibration
Inputs required to run either model include plant phenotype characteristics, crop management, soil properties, and daily weather data as summarized in Table 2. For calibration, data from Stoneville, MS in 2005 was used for Cocodrie, and 2001 for Lemont [24]. The calibration data included yield, unit grain weight (g/1000 seed), sowing date, flowering, and harvest dates. Table 2. Input requirements and data sources used to calibrate and run CERES-Rice and ORYZA for cultivars Lemont (1983Lemont ( -1998 and Cocodrie (1999-2016) for each location.

Model
Input Type Variables Data Source

Crop cultivar coefficients for Lemont and Cocodrie
Calibrated by Genetic Coefficient Estimator Management planting date, row, and plant spacing [24,57] Calibration inputs yield, grain weight, flowering, and harvest dates [24] ORYZA Weather Tmin, Tmax, SRAD, rainfall, VP, WD iAIMS Climate data (1) Soil soil texture and hydrological properties (2) [48,49] Crop phenology development rate fraction of total dry matter partitioned to organs over growth stages specific leaf area index Calibrated by DRATES Standard crop file from ORYZA v3 Management dates of emergence [24] Calibration inputs dates of emergence, panicle initiation, flowering, and maturity [24,58] Both models Planting dates CERES-Rice was calibrated using the Genetic Coefficient Estimator (Gencalc), which is a software tool contained as part of the Decision Support System for Agrotechnology Transfer (DSSAT) interface [62] with results shown in Table 3. For ORYZA calibration, the DRATES calibration program was used based on procedures defined by Bouman and Van Laar (2006) [34]. Parameters describing biomass partitioning among plant organs and specific leaf area were fixed to that contained within the standard crop file included with ORYZA (v.3), but a difference in cultivar spikelet growth factor (SPGF) was obtained for each location. After this initial calibration, a re-parameterization was done to account for differences in cardinal temperatures between the U.S. and tropical climate conditions in which the model was originally developed for. A sensitivity analysis consisting of factorial combinations of three base (8, 9, and 12 • C), two optimal (30 and 33 • C), and two maximum temperature values (42 and 999 • C) was conducted to obtain an optimal cardinal temperature value set that reduced systematic error with respect to the growing season duration for either cultivar using this same data. Henceforth, this version of ORYZA was referred to as 're-parameterized.' All other calibration values were the same as the original version (Table 3). : Time from seed emergence to the end of juvenile phase (expressed as growing degree days (GDD) in • C). P2O: Critical photoperiod or the longest day length (in hours). P2R: extent to which phasic development leading to panicle initiation is delayed (GDD in • C). P5: time period from grain filling to physical maturity (GDD in • C). G1: Maximum spikelet number coefficient. G2: Maximum possible single grain weight (g) under ideal growing conditions. G3: Scalar vegetative growth coefficient for tillering relative to the IR64 tillering coefficient. G4: Temperature tolerance is the scalar coefficient. DVRJ: Development rate in a juvenile phase ( • C d −1 ). DVRI: Development rate in a photoperiod-sensitive phase ( • C d −1 ). DVRP: Development rate in panicle development ( • C d −1 ). DVRR: Development rate in a reproductive phase ( • C d −1 ).
where N is the number of observations (years), and O and S are observed and simulated values at year i. The index of agreement and mean absolute percentage error (MAPE) were defined as: where d is the agreement index value, with symbols as defined above. After model evaluation, results from the model predictions of yield for these two varieties were compared against the URRN data to explore the distribution of observed versus simulated yields across all cultivars (up to 101) for a given year at each location to evaluate trends and identify off-nominal years with respect to yields.

Comparison of Growth Duration
CERES-Rice generally showed better performance simulating rice growth duration as compared to default-ORYZA in term of the correlation coefficient and RMSE (Figure 1a  ORYZA cardinal temperatures were subsequently re-parameterized to better represent U.S. conditions. Optimal parameter values (with respect to improved correlation coefficients and RMSE) obtained via sensitivity analyses of ORYZA using factorial combinations of cardinal temperatures were 12 • C, 33 • C, and 42 • C (base, optimal, and maximum temperature) at Beaumont and Stoneville for Lemont, and Stoneville for Cocodrie ( Figure 2). Slightly different values of 9, 30, and 42 • C were obtained for Cocodrie at Crowley, LA. However, for the simulations in this study, the same cardinal temperature values across all locations of 12 • C, 33 • C, and 42 • C for base, optimal, and maximum temperature, respectively, were used.  Similar performance metrics were obtained with re-parameterized ORYZA in predicting growth duration as compared to CERES-Rice (Figure 1, Re-ORYZA). These two models were able to capture the annual variation of observed growth duration for most years and locations for either cultivar ( Figure 3). However, both models showed relatively low accuracy for certain years likely as a result of climate extremes. For examples, the years 1983 to 1987, and 1992 in Stoneville, MS for cultivar Lemont (Figure 3a) were most strongly underestimated by ORYZA and CERES-Rice, and were associated with later occurrence of frost in the growing season (minimum temperatures at 0 °C or below) as much as 20 days after the planting date as compared with the 34-year mean of nine days ( Figure S3). Years 2009 and 2010 for Cocodrie (Figure 3d) at Stoneville, MS were associated with higher than average rainfall (>225 mm over the 465 mm annual season mean) and lower solar Similar performance metrics were obtained with re-parameterized ORYZA in predicting growth duration as compared to CERES-Rice (Figure 1, Re-ORYZA). These two models were able to capture the annual variation of observed growth duration for most years and locations for either cultivar (Figure 3). However, both models showed relatively low accuracy for certain years likely as a result of climate extremes. For examples, the years 1983 to 1987, and 1992 in Stoneville, MS for cultivar Lemont ( Figure 3a) were most strongly underestimated by ORYZA and CERES-Rice, and were associated with later occurrence of frost in the growing season (minimum temperatures at 0 • C or below) as much as 20 days after the planting date as compared with the 34-year mean of nine days ( Figure S3). Years 2009 and 2010 for Cocodrie (Figure 3d) at Stoneville, MS were associated with higher than average rainfall (>225 mm over the 465 mm annual season mean) and lower solar radiation, which may have had an indirect response on decreased growth duration not captured by the models. In contrast, growth duration observations at Crowley, LA for Cocodrie from 2003 to 2016 were well reproduced by both models (Figure 3c) possibly due to lower annual variability in mean temperature (Table 1 with a standard deviation of 0.6 • C).

Comparison of Observed and Simulated Yields
Model predicted yields were evaluated over two periods (1983-1998 and 1999-2016) for the two cultivars at four locations. At Stuttgart, AR, there was a large discrepancy between observed and modelled yields for years 1989 for Lemont and 2006Lemont and , 2008Lemont and , and 2009 for Cocodrie for both models (Figure 4a). These outlying years contributed to the low correlation coefficients for either model (between 0.07 and 0.24, Table 4). In other years, simulated yields for Lemont from re-parameterized ORYZA were generally closer to observations than CERES-Rice, especially in 1994 (Figure 4a). Seasonal mean maximum and minimum daily temperatures were 1.1 and 1.3 °C lower in this year compared to the average 34-year value ( Figure S1), and simulated yields in CERES-Rice were more than 15% higher than observed values, as compared with a less than 2% increase for re-parameterized ORYZA. Across all years for Stuttgart, AR, the slope of predicted grain yield versus seasonal

Comparison of Observed and Simulated Yields
Model predicted yields were evaluated over two periods (1983-1998 and 1999-2016) for the two cultivars at four locations. At Stuttgart, AR, there was a large discrepancy between observed and modelled yields for years 1989 for Lemont and 2006, 2008, and 2009 for Cocodrie for both models (Figure 4a). These outlying years contributed to the low correlation coefficients for either model (between 0.07 and 0.24, Table 4). In other years, simulated yields for Lemont from re-parameterized ORYZA were generally closer to observations than CERES-Rice, especially in 1994 (Figure 4a). Seasonal mean maximum and minimum daily temperatures were 1.1 and 1.3 • C lower in this year compared to the average 34-year value ( Figure S1), and simulated yields in CERES-Rice were more than 15% higher than observed values, as compared with a less than 2% increase for re-parameterized ORYZA. Across all years for Stuttgart, AR, the slope of predicted grain yield versus seasonal temperature was negative for both models, but nearly three times lower for CERES-Rice (-863 versus -290 kg • C −1 , Table 5). This implied a higher yield estimation at cooler temperature years as compared to ORYZA. The observed standard deviation for Cocodrie (±1607 kg ha −1 ) was nearly three times greater than simulated from CERES-Rice (±594 kg ha −1 ) and 1.5 times that of ORYZA (±1098 kg ha −1 ) ( Table 4 and Figure 4a). The higher variation in the observed yield can likely be explained by the difference in field management, the impact of biotic factors, and climate variability not duplicated in either model. For example, yields were affected by diseases and pests in 2006 [64] and this may have been responsible for yield fluctuations in other years as well.   Table 4. Correlation coefficient (r), root mean square error (RMSE), standard root mean square error (SRMSE = RMSE/average yield) between observed (Ob) and simulated (S) yields, ratio of simulated average yield to observed average yield (S/ob), the mean absolute percentage error (MAPE), and standard deviation of observed (Ob Std) and simulated (S Std) rice yield for varieties Lemont (1983Lemont ( -1998, Cocodrie (1999-2016), and comparison of 50% percentile yield of 101 varieties for 1983-2016 to simulated yield for cultivar Lemont (1983Lemont ( -1998 and Cocodrie (1999-2016). Models were CERES-Rice and a re-parameterized version of ORYZA.  At Crowley, LA, better agreement (r = 0.52) between simulated yield from re-parameterized-ORYZA and observed yields for Lemont was found for 1983-1998 (Figure 4b), but a somewhat lower correlation (r = 0.41) was observed for CERES-Rice (Table 4). Compared to other locations, there were more years with high rainfall amounts that coincided with reduced solar radiation and lower temperatures that likely constrained the rice yield ( Figure S6). For example, a large yield reduction occurred in 2012 that neither model was able to replicate, and this year was associated with lower radiation, higher rainfall, and lower seasonal temperatures ( Figure S2) compared to other years. This suggested both models may have difficulties simulating yields under a combination of lower radiation and temperatures.

Model
At Stoneville, MS, both models simulated most of the annual variability for Lemont from 1983 to 1998, excluding 1996 when the observed yield was higher (Figure 4c). Similar performance metrics for mean yield and variability of yield were obtained (Table 4). For Cocodrie, simulated yield from re-parameterized-ORYZA was close to observations from 1999 to 2005, but observed yield was slightly underestimated by CERES-Rice during this period (Figure 4c, Table 4). After 2005, there were larger discrepancies in simulated and observed yields for both models, particularly for 2006, 2009, and 2011. As climate conditions did not appear unusual for these years, and models under-predicted observed values, we are not sure why this relatively poor performance occurred.
At Beaumont, LA, correlations from both models were poor for Lemont from 1983 to 1998 (Table 4). In particular, the large yield reductions in 1986, 1988, and 1991 were not captured by models (Figure 4d). After 1992, simulated yields from CERES-Rice were lower than observed yield and re-parameterized-ORYZA, which is a result likely associated with lower solar radiation values measured for those years ( Figure S4). Following this point, in CERES-Rice, yield was more responsive to solar radiation (m = 706 kg MJ −1 m −2 day −1 , Table 5) from 1983 to 1998, as compared with re-parameterized-ORYZA (m = 284 kg MJ −1 m −2 day −1 , Table 5). From 1999 to 2016, observed and simulated yield for Cocodrie was significantly correlated in re-parameterized-ORYZA (r = 0.58) and CERES-Rice (r = 0.56) (Table 4), and the models captured the variability of observed yield for most years except 2010 and 2012 (Figure 4d) in which simulated yields were higher than observations. For 2010, the night temperature (22.7 • C) was higher than average for other years (21.4 • C) ( Figure S4), but no other reasons for the decline in yield were identified, suggesting a non-climatic factor may also be responsible.

Model Skill in Simulating 50% Percentile Yield of 101 Varieties
The ability of the models to capture the observed annual yield trends of the full set of 101 cultivars in the URRN database using only this calibration for these two varieties was assessed. Significant correlations were found between simulations and the 50% percentile yield (median) of observations at three out of four locations ( Figure 5). Both models showed a similar performance in this respect and captured most of the annual variability of observations except for some extreme years (Table 4).
At Stuttgart, AR, there was a relatively small increase trend (63 kg year −1 ) in median yield response versus time (Table 1) (Figure 5a). For these years, solar radiation declined ( Figure S4), and may reveal an over-sensitivity of both models to lower solar radiation. Despite models being calibrated for a single variety in those years, the predictions were within the 25% to 75% distribution of 101 cultivars for most years except 2006 and after 2012.  A trend of increasing median yield for 101 cultivars was also observed for Crowley, LA (Figure 5b, Table 1), which both models replicated. Exceptions were 1996 and 1997, in which yields were under-predicted, and 2007 and 2012, in which yields were over-predicted. These latter two years were associated with more cloudy days with higher rainfall and cooler temperatures, and presumably higher relative humidity, suggesting the models were unable to replicate the combined negative impact on yield ( Figure S2). Observed yields in 2012 for 101 cultivars ranged broadly from 1133 to 9756 kg ha −1 , indicating a wide range of cultivar responses to these wetter, cooler conditions. In this year, rice was also planted 8to 23 days earlier than the average sowing date [60]. Further evaluation of in-season climate data showed an increase of 6.7 days in which mean daily temperature was below 23 • C during reproductive development, as compared with the average of 1.3 days between 1999 and 2016, which can negatively impact yields [60].
Median yields also increased with time (Figure 5c, Table 1) for Stoneville, MS, and variability among observed cultivar yields was relatively larger than simulated ( Table 4). As with results for Lemont and Cocodrie (Figure 4c), the mean yields were simulated well for most years-correlation coefficients were 0.7 or higher for both models (Figure 5c). However, both models under-estimated yields for outlying years 1996, 2006, and 2007-2011 for reasons that were not clear from the available climate data.
From 1985 through 1991, URRN yields were over-estimated by both models at Beaumont, TX (Figure 5d), with similar trends as reported in the prior section. Re-parameterized ORYZA appeared to more closely replicate yields prior to 2007, while CERES-Rice yields were closer to those observed from 2007 through 2016. Correlations were relatively high for both models (0.56 and higher) despite over-predictions. Median yields in this location appeared to have plateaued after 2002, with considerably more variability than the other three locations (Figure 5d, Table 1), which is a trend that both models captured.

Simulated Rice Responses
Simulations of growth duration for CERES-Rice were more accurate as compared with default-ORYZA. This result was expected given that ORYZA was initially developed for indica rice sub-populations in tropical climates with a smaller range of climatic variation than the U.S. [43]. The default base and optimal temperatures for ORYZA were 1 and 3 • C lower than those in CERES-Rice. After re-parameterizing the cardinal temperatures in ORYZA, these values were quite similar between models, even though the model source code and implementation of these values in terms of phenology vary. These changes improved ORYZA accuracy and reduced systematic error (correlation coefficeint increased by 0.03-0.27 and RMSE decreased by 0.3-1.9 days across all locations), and exhibited similar metrics as obtained with CERES-Rice (Figure 1a,b). We applied the same set of cardinal temperatures to the re-parameterized ORYZA across both cultivars and locations, but future studies need to consider the extent to which temperature responses vary within the tropical japonica subpopulation. This is a question requiring more data than is currently available to test for phenotype differences. In years with normal climate conditions, CERES-Rice and re-parameterized-ORYZA showed similar performance in a simulating yield. Standard RMSEs across both models, cultivars, and locations ranged from 0.11 to 0.26 (Table 4), which were lower than those reported in other types of modeling studies [65], and, thus, indicated good model performance. Studies by Espe et al. (2016) [20] and Tsvetsinskaya et al. (2003) [46] observed similar accuracy for both models in the U.S. as well. Espe et al. (2016) [20] reported that ORYZA simulated rice yields for the hybrid variety CXL745 well in the Southern U.S., with 77% of simulated yields falling within ±15% of the observed means. Similarly, in the present study, 74%-75% of simulated yields for Cocodrie and Lemont were within ±15% of observance. Although RMSEs were relatively low and Willmott's index of agreement was above 0.5 for most location-cultivar combinations, low correlation coefficients (Table 4) were observed for both models. This was primarily the result of the limited ability of either model to capture observed standard deviations across all years (Table 4) and was associated with several outlying years skewing the comparisons between observed and simulated values. Common climate conditions for these outlying years were associated with (a) colder early seasonal temperatures where the frequency of night-time temperatures below 14 • C was relatively high, such as at Stuttgart, AR, (b) higher seasonal rainfall amounts associated with cooler temperatures and lower solar radiation as at Crowley, LA, (c) higher night-time temperatures as at Beaumont, TX in 2010 ( Figure S4C), and (d) years with lower solar radiation, with other factors being equal.
Performance between models was similar, except re-parameterized ORYZA that provide better yield responses during lower temperature episodes (e.g., Stuttgart, AR 1994) and years with lower solar radiation (e.g., Beaumont, TX 1993-1998, Figure 4d). When comparing the simulated results for the two cultivars to the range of the 101 varieties in the URRN database, the correlations between the 50% percentile yield of observations and simulations were significant at the 5% level at three out of four sites for both models ( Figure 5, Table 4). Correlations between predicted and observed values were higher for both models for these median yields as compared with results for the two individual cultivars, even though RMSE values were not necessarily lower in comparison (Table 4). This was because an individual cultivar response involves different sensitivities to climate variables and field management that can be challenging to capture over such a broad range of data. However, at a larger scale, such as the case when yield responses are grouped across multiple varieties at each location, such differences are offset by each other, and the median yield is mainly determined by climate variability rather than phenotypic difference. These results suggest that either model, once calibrated with a subset of relevant cultivars, is suitable for exploring the general or relative response of multiple rice cultivars in the Mississippi Delta region for decision support studies relevant to the rice production industry.

Model Knowledge Gaps
A higher frequency of temperatures less than 14 • C was observed at Stuttgart and Beaumont (16 and 17 days, respectively), and rice yields at both locations for some cultivars were negatively associated with cold temperatures ( Figure S5). Minimum daily air temperatures below 12-15 • C can impact yield [66,67], and several years at Stuttgart were associated with more cold events in which the models were unable to replicate yields accurately despite including procedures that simulated the effect of cold on rice yield [25,68,69]. Rice sterility is more sensitive to daily minimum, as opposed to mean temperature between panicle initiation and 50% heading stages [70]. This phenomenon likely needs to be incorporated into the two rice models as well as the impacts of low temperatures on seedling vigor [71]. Espe et al. (2016) [20] also indicated that cold-induced sterility posed a challenge to ORYZA's accuracy in a similar study in the U.S. Alternatively, cool temperatures could have limited stand establishment, development of yield components, or grain filling if they occurred during seedling, vegetative, or post flowering stages respectively.
Exposure to high temperatures near or during flowering (above 33.7 • C or higher for short periods of time) can also reduce yield due to reduced pollination [72][73][74] and spikelet sterility [72]. An improved heat sterility model including the hour at which flowering occurs during the day and transpirational cooling was introduced in the version of ORYZA [75]. However, a recent study of Mohammed et al. (2018) [76] showed that Southern U.S. cultivars had earlier flower opening times (normally open at about 8:00-9:00 a.m. or at about 10:00-11:00 on rainy days) as compared to conditions where this sterility model was developed. High night temperatures (32 • C) were also associated with reduced pollen germination and spikelet fertility in Beaumont, TX [77,78]. Thus, the feasibility of this heat sterility model needs to be further tested in Southern U.S. conditions given warming temperature trends (Table 1).
Other limitations in model performance may be associated with predictions of photosynthesis and its relationship to solar radiation and heat stress. CERES-Rice uses a radiation use efficiency approach (Equation (4)) with simple multipliers to adjust for the impact of high heat on a 24-h basis. This model was observed to be more sensitive to lower solar radiation, and yield was greatly underestimated during these conditions (Figure 4d). ORYZA uses a more mechanistic representation of a leaf photosynthetic rate, which is scaled to canopy at selected hours during the day (Equations (9)-(11)). However, both approaches involve empirical and/or temperature independent relationships applied in calculating CO 2 assimilation, which may underestimate biomass growth in years with excessively high heat. For example, shorter growth duration and lower observed yield were simulated by both models in 1998 at Stuttgart, AR (Figure 4a) in response to a higher temperature (about 2 • C higher than average from 1983-1998, Figure S1a). The rising temperature may increase plant maintenance respiration [79,80], leading to a reduction in the amount of assimilates for growth and yield [81]. Since it is likely that temperature extremes are likely to be more common under future climates, the depiction of the sensitivity of rice CO 2 assimilation rates and yield formation processes to heat need to be re-evaluated. More mechanistic approaches, such as those that integrate coupled biochemical models for photosynthesis and stomatal conductance within an energy balance at the leaf level, are one promising method that can directly simulate the interaction of solar radiation, air temperature, and CO 2 on canopy and leaf temperature and, thus, effects on crop photosynthesis, transpiration, and other factors [23]. The sensitivity of rice growth and development to extreme climate events and warmer night temperature will need to be further developed and tested to improve the skill of ORYZA and CERES-Rice models for application in the Mississippi Delta and evaluation of climate change response. The present results, however, confirm that these two models can provide reasonably accurate responses in comparison with the historical yield record at representative sites across this region.

Conclusions
An increase in base and optimal temperatures improved ORYZA accuracy and reduced systematic error for applications in the U.S. Mississippi Delta region. Both CERES-Rice and Re-parameterized ORYZA rice models showed acceptable skill in reproducing the growing season duration, and exhibited higher accuracy at sites with lower interannual variability in temperature. Both models had similar performance for predicting rice yield for most varieties, locations, and years, but CERES-Rice was more sensitive to years with lower solar radiation. Negative impacts of very warm or cold temperatures were challenging for both models to replicate, which suggests that improving algorithms associated in reproducing the cold injury and high temperature stress are needed for rice model improvement. Similarly, the current representations of photosynthesis need to be improved to better capture interactions of solar radiation, temperature, and CO 2 concentration on the plant growth rate and yield. Integrating a coupled biochemical model of photosynthesis and stomatal conductance within an energy balance approach is one such methodology, which has been successfully implemented for other crops. Broad variation among annual cultivar yields also indicates the need to improve model performance in considering different cultivar sensitivities to temperature stress. However, the present results suggest that these calibrated crop models can be used to simulate the relative response of multiple rice cultivars to historical and current climate conditions in the Mississippi Delta region for decision support studies.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/12/1905/s1, Figures S1-S4: Boxplots of daily distribution of solar radiation, rainfall, minimum temperature (TN) and maximum temperature (TX) per year during the rice growing season from 1983 to 2016 at (S1) Stuttgart, AR, (S2) Crowley, AR, (S3) Stoneville, MS and (S4) Beaumont, TX; Figure S5: Pearson correlation coefficients (p < 0.05) between rice yield for 32 cultivars and seasonal average values of minimum temperature (TN), maximum temperature (TX), rainfall and the total days when minimum temperature is below 14 • C or maximum temperature is above 33 • C. Annual yield records for different varieties range from 11 to 22 years between 1983 and 2016 at each location, so the critical value (dotted line) was different for each cultivar; Figure S6: Seasonal average daily solar radiation (blue) and maximum average daily temperature (red) versus seasonal total rainfall from 1983 to 2016 (p < 0.01) for the Crowley, LA location.

Conflicts of Interest:
The authors declare no conflict of interest.
Disclosure Statement: Mention of a trademark or proprietary product does not constitute a guarantee or warranty of the product by the U.S. Department of Agriculture and does not imply its approval to the exclusion of other products that also can be suitable. USDA is an equal opportunity provider and employer. All experiments complied with the current laws of the United States, the country in which they were performed.