A new nonlinear method for calculating growing degree days

Precise calculations of growing degree days (GDD) are an important component in crop simulation models and managerial decisions. Traditional methods for calculating GDD assume linear developmental responses to temperature and cannot precisely account for the delay in growth or development at temperatures above the optimal temperature (Topt). A new nonlinear method for calculating GDD was developed. Variations in the prediction of the dates since sowing to various developmental stages and performance measures for describing the accumulation of dry matter by GDD for two widely planted crops (corn and wheat) were used to evaluate the new method in comparison with the traditional methods. The new method predicted the dates of the developmental stages more precisely (date variations reduced by 1 d), and the errors for the predictions of the accumulation of dry matter for winter wheat and corn were smaller. The method was most promising for spring wheat. The new method was more stable and more precise than traditional methods, especially when Topt was lower than the maximum air temperature.

Precise calculations of growing degree days (GDD) are important in models simulating crop growth and for the management of field crops. GDD is also a climatic feature. The use of GDD has vastly improved the description and prediction of phenological events compared with other approaches, such as time of year or number of days, particularly for crop phenology and developmental stage 1,2 .
The relationship between the rate of development and temperature is key for calculating GDD. A linear relationship, which assumes that the rate is proportional to the temperature above a threshold, is used most widely and is often precise for intermediate temperatures [3][4][5][6][7] . However, the assumption of rate-temperature linearity will yield errors when temperatures tend toward extremes under variable conditions 8,9 , i.e., a linear relationship between temperature and plant growth is inappropriate in long-term studies, especially for complete life cycles 10 . Many methods that assume a nonlinear relationship have thus been developed, each with strengths and weaknesses. For example, a bilinear approach has been adopted 11,12 in which the responses to sub-and supra-optimal temperatures are described by different linear equations, and the real response curve is generally smooth. An exponential equation is usually effective for simulating responses at low to intermediate temperatures but not for simulating responses to high temperatures because it does not allow for a rate of development at high temperatures 13 . A quadratic equation is a simple model that can allow a lower rate of development at high temperatures 14,15 . However, the temperature response is rarely a symmetric parabola, and applications of quadratic models may thus be inaccurate. Yin adopted a beta function containing three parameters (the cardinal temperatures) to describe the temperature response and reported successful simulations of the development of several crops (corn, wheat, barley, sorghum, and beans); the method was superior to widely used thermal time approaches in predicting crop developmental stages 16,17 . The approach accounts for the asymmetric temperature response for the developmental rate and the decline in the rate above T opt .
Traditional GDD correlates developmental rate linearly to temperatures above the lower threshold temperature in some applications of the procedure; however, linearization is often criticized for its oversimplification despite being widely used 18 . A more stable and less variable GDD, which should be calculated by a precise method, is needed so that the stages of the crop growth period may be accurately compared and predicted regardless of environmental conditions. However, traditional approaches cannot accurately account for the delay in growth or development at temperatures above the optimal temperature (T opt ). The nonlinear relationships between developmental rate and temperature discussed above are rarely applied to calculate GDD. Therefore, the objective of this study was to develop an improved nonlinear method of calculating GDD and to evaluate the accuracy and applicability of the method by comparison with other methods.

Materials and Methods
GDD, which is cumulative daily thermal time (DTT), is calculated as: Calculating DTT is key to the methods for calculating GDD, and present methods for calculating GDD differ in their method for calculating DTT.
Present widely used methods. McMaster and Wilhelm (1997) proposed two methods for calculating DTT (Methods 1 and 2) needed for calculating GDD, and both methods have been widely used in recent studies [18][19][20][21][22] . Method 1, which is simpler than Method 2, calculates DTT as: where T max is the maximum temperature, T min is the minimum temperature, T avg = (T max + T min )/2, T b is the base temperature, and T u is the upper threshold temperature. Method 2 is an improvement on Method 1. T b is compared with T u before the average temperature (T avg ′) is calculated. T m and T n are adjusted if they are <T b or >T u . In this method, DTT is given by: where T m = min (T max , T u ), T n = max (T m , T b ), and T avg ′ = (T m + T n )/2. Method 3 introduces T opt . Thermal time increases linearly until T opt is reached and decreases rapidly by another linear relationship at supraoptimal temperatures (T > T opt ) 9,11,23 . Method 3 is calculated as: T  T  T  T  T  T  T   T  T   T  T  T  T T  T  T   T where HTT is hourly thermal time and T h is the hourly temperature.
However, method 3 is linear both above and below a sharp T opt , and the value of thermal time may thus fluctuate near T opt . The response curve is also less smooth.
Description of the new method. The beta-distribution function has a density of zero when x ≤ 0 or ≥1 and a maximum density at an optimal x between 0 and 1. Replacing the dependent variable x with temperature T between T b and T u leads to the following expression that can be used to describe the relationship between developmental rate and temperature 24 : where R max is the maximum developmental rate and T h is the hourly temperature. The equation assumes that the developmental rate of the crop is maximal at T opt 11,25 because the contribution of thermal time is T opt −T b °C. The hourly thermal time can then be calculated by a beta-distribution function method (BFM) as: ( ) The equation has three parameters (T b , T opt , and T u ), and the thermal time will be zero if T = T b or if T = T u and will be maximum if T = T opt . Figure 1 shows the comparison of the four methods used to calculate DTT and HTT. The daily GDD of Methods 1 and 2 is the same when (T max − T avg ) = (T avg − T min ), the daily GDD is higher for Method 1 than that for Method 2 when (T max − T avg ) > (T avg − T min ), and the daily GDD is lower for Method 1 than that for Method 2 when (T max − T avg ) < (T avg − T min ). The two methods involve a linearly-increasing segment with temperature T up to T u , beyond which GDD = T u − T b . The curve of the temperature response for BFM is smoother than that for Method 3.
Model evaluation. GDD is frequently used to describe biological processes but no canonical forms are available for calculating GDD. Therefore, we used the coefficient of variation (CV) of predictions of developmental stages and the performance of the accumulation of dry matter described by GDD to evaluate the four methods.
Predicting developmental stages. The observed dates of developmental stages from field data were used to calculate the GDD required to reach a particular stage by the four methods. However, the GDD required to reach a particular stage calculated by a particular method with observed dates in different years was not always the same.
The CV of the dates predicted by one method for calculating GDD was used to test the performance of the method. The lower the CV, the better the prediction. CV was calculated as: where SD GDD is the standard derivation of the annual GDD required for a particular developmental stage since sowing, and GDD m is the mean annual daily GDD during the developmental stage since sowing. Willomtt 26 proposed a refined index of agreement (d r ) for evaluating model performance, defined as: where O i represents the developmental stage of an observation, P i represents the method of prediction of the developmental stage, and O represents the mean value of an observation. A d r of 1 indicates perfect agreement between prediction and observation.   Describing the accumulation of dry matter. The accumulation of dry matter is commonly described by a logistic model as a function of GDD 19,[27][28][29][30] . The normalized logistic model was fitted by GDD to test the performance of the methods for calculating GDD as follows: where Y o is the observed amount of dry matter, Y m is the maximum amount of dry matter, y is the normalized amount of dry matter, and a and b are coefficients. The fitted amount of dry matter for each location was evaluated by the root mean square error (RMSE), defined as: where y oi and y fi are the observed and fitted amounts of dry matter, respectively, and n is the number of observations. Data for calculations. Crop phenological data. Data for wheat (a C3 crop) and corn (a C4 crop), which are widely planted around the world in a large range of cardinal temperatures (T b , T opt , and T u ), were used in this study. T b , T opt , and T u for corn were 8, 33, and 40 °C, respectively 7,11,31,32 . Spring and winter wheat respond differently to temperature, and the three cardinal temperatures for the two wheat crops were thus not identical in the study 33 . T b , T opt , and T u for winter wheat were 0, 24, and 45 °C for the vegetative phase (from emergence to heading) 7,31 and 8, 29, and 40 °C for the reproductive phase (from heading to maturity), respectively [34][35][36] . T b , T opt , and T u for spring wheat were 0, 24, and 42 °C, respectively 7,37-41 .
The phenological and temperature data for the two crops were obtained from agro-meteorological experimental stations maintained by the Chinese Meteorological Administration. The records for both crop phenologies were available from 1991 to 2010, although some records were missing. Data from the stations with records for more than eight years were used for our analysis ( Table 1). The five stations (Turpan, Korla, Hezhang, Xinxiang, and Xianyou) represent a wide range of climatic conditions for the two crops. For example, Turpan station is in the warm-temperate and continental-drought climatic zone in northwestern China, which has extreme aridity, high temperatures, and a large temperature difference during the growth period. In contrast, Xianyou station is in a subtropical-monsoon climatic zone in southeastern China.
Data for total accumulated aboveground dry matter. Data for the total accumulated aboveground dry matter from two locations where the two crops are widely planted and have different climates (northern Xinjiang and central Shaanxi Plain) were reported in previous studies and were used to test the performance of the methods (Table 2) [42][43][44][45][46][47][48][49][50][51] . Northern Xinjiang has a temperate continental climate with a mean annual temperature of 8.1 °C and a mean annual precipitation of 577.8 mm. In contrast, Central Shaanxi Plain has a temperate monsoon climate with a mean annual temperature of 12.9 °C and a mean annual precipitation of 641.3 mm. GDD from sowing to a developmental stage could differ between the calculation methods, even for the same crop at the same station. For example, the GDD requirements from sowing to maturity for spring wheat at the Turpan station were 2167.85 ± 81, 2175.30 ± 77, 957.39 ± 52, and 1833.13 ± 56 °Cd for Method 1, Method 2, Method 3, and BFM, respectively. GDD from sowing to a particular developmental stage calculated by a particular    1 and 2). The standard deviations (SDs) of GDDs calculated by Methods 1-3 for winter wheat and corn were generally lower than those calculated by BFM, and the SDs were the lowest for Method 3. However, a lower SD did not necessarily indicate that one method was better than the others because the mean annual GDD calculated by the method may also be lower (e.g., Method 3). GDD has been used to describe crop development and the duration of a process or the time required to reach a particular stage, and confirming the merit of the four methods based only on GDD is thus difficult.
CV for the prediction of developmental stage. CV for corn. The CVs of dates from sowing to the developmental stages predicted by the four methods were generally similar for corn ( Table 3). The CVs calculated by Methods 1 and 2 were higher than those of Method 3 and BFM for the three locations. The CVs calculated by BFM were the lowest for two locations (Xinxiang and Hezhang). In contrast, the predictions by Method 3 were better from the 7th leaf stage to maturity in Korla than from emergence to the 3rd leaf stage compared with those of the other three methods. The new method (BFM) was better than Methods 1 and 2 for all cases. The   CV for spring wheat. For spring wheat, the CVs calculated by BFM for predicting the developmental stages since sowing were the lowest, followed by Method 3 and Methods 1 and 2 (see Turpan in Table 4). Method 2 was slightly better than Method 1. The advantage of BFM tended to increase; e.g., compared with Method 2, the CV of the jointing date predicted by BFM was 1 d lower and the CV of the milk date was 3.1 d lower. CV from sowing to heading was lower for Method 3 than that for the other methods, but the CVs of the method increased with development, especially by maturity (e.g., CV was 5.7 d for Method 3 and approximately 4 d for the other methods).
CV for winter wheat. For winter wheat, the CVs were the lowest for BFM for predicting the developmental stages since sowing, followed by Method 3 and Methods 1 and 2; however, the advantages of BFM and Method 3 were not larger than those for spring wheat (see Hezhang and Xianyou in Table 4). Methods 1 and 2 did not differ in most cases, as expected. CVs calculated by Method 3 were usually nearly 0.5 d lower than those of Methods 1 and 2. In contrast, CVs calculated by BFM were usually <1 d. The performance of the four methods for predicting developmental stage. We further analyzed the performance of the four methods for predicting developmental stage using a refined index of agreement, d r , which is widely used as a goodness-of-fit indicator ( Table 5). All methods performed well for all locations. For winter wheat, d r for BFM and Method 3 was similar in two locations, and the two methods were generally better than Methods 1 and 2. BFM was better than the other three methods for spring wheat in Turpan, followed by Methods 1 and 2 and Method 3. There was no difference among Methods 1, 2 and 3 in three locations for corn, and BFM had better performance than the three methods.
The performance of describing the accumulation of dry matter. The accumulation of dry matter was described well with the normalized logistic model as a function of GDD calculated by the four methods for the two locations (Fig. 4), with an RMSE between observed and fitted values by the four methods within 0.11 for wheat and 0.16 for corn. The methods agreed the most for northern Xinjiang. RMSE fitted with GDD calculated by BFM was the smallest, especially for spring wheat (RMSE was 0.06 for BFM and >0.08 for the other methods), as expected. Methods 1 and 2 did not generally differ for the two crops and the two locations. Method 3 was not better than Methods 1 or 2 in most cases. The curve for the accumulation of dry matter was steeper for Method 3 than the other methods because GDD was lower for Method 3 in some cases (Fig. 4c,d). and Xianyou (f) for winter wheat. T max , maximum temperature; T avg , average temperature; T u , upper threshold temperature; T opt , optimum temperature; T b , lower threshold temperature. T max − T b is the temperature difference between T max and T b , T max − T opt is the temperature difference between T max and T opt , T avg − T b is the temperature difference between T avg and T b , T avg − T opt is the temperature difference between T avg and T opt , T u − T b is the temperature difference between T u and T b , T opt − T b is the temperature difference between T opt and T b , and T u − T opt is the temperature difference between T u and T opt .

Discussion
Most of the methods for calculating GDD assume that crop development responds linearly to temperature. This calculation of GDD is appropriate for predicting plant development if several conditions are met 9,24 . The use of GDD to describe the duration of growth is reasonable when plant developmental rates for a region are linear over a wide range of temperatures 52 . However, several experiments have indicated that the rates could also respond nonlinearly to temperature 24 . To our knowledge, crop plants do not survive the high temperatures that could stop development, and interest in the production of crops is highest when environmental temperatures are near T opt . However, Eqs (5) and (6) indicate that high temperature remains the largest contributor to crop development and they are thus obviously inappropriate. Temperatures above T opt will retard crop development.
Conditions applicable to the GDD calculation methods. The concept of degree days, with its linear relationship between temperature and developmental rate, is inadequate for simulating field populations under highly variable temperatures 8 . The methods of an optimized developmental response (OR) to temperature (e.g., Method 3 and BFM) are less convenient than the methods of a linear developmental response (LR) (e.g., Methods 1 and 2) for calculating GDD. However, the OR methods generally have higher predictive accuracy, which is also supported by the physiological interpretations of optimal, supraoptimal, and limiting temperatures 11 . OR methods should thus have an advantage over LR methods. The OR methods were more accurate than both LR methods for predicting the dates of the developmental stages for corn, BFM was slightly better than Method 3 in some cases, and Methods 1 and 2 were slightly worse than Method 3. BFM best predicted the dates of the developmental stages for wheat since sowing, especially for spring wheat. T opt is usually well above 30 °C for corn, and T max is lower or slightly higher than usual for most growing stages. T opt is lower for wheat than that for corn, but the growth period of winter wheat is generally from mid-October to the end of May; in addition, T max is usually lower than T opt because temperatures in winter and spring are usually low. The OR methods would have more advantages with more days when temperatures are above T opt throughout the crop growth period. The temperature differences from sowing to the developmental stages for the two crops were thus used for further analysis (Fig. 5).
Three locations (Korla, Xinxiang, and Hezhang) for corn had only a few days with T avg > T opt but T m > T opt for more than one third of the days at Korla and Xinxiang from sowing to the 7th leaf stage and for approximately 5-20% of the days during the other developmental stages. However, a few days at Hezhang had T max > T opt . This analysis demonstrated that the CVs calculated by the four methods were similar for Hezhang but CVs calculated by BFM were lower than those of the other methods for two other locations (Korla and Xinxiang). Figure 5 further illustrates why BFM was more precise at predicting developmental dates than the other methods for wheat, especially for spring wheat in Turpan.
The function of ORs to temperature has a large influence on the performance of the OR methods. Method 3 was superior to both LR methods (Methods 1 and 2) in most cases and opposite in some cases. For example, the CVs from sowing to milk and maturity for spring wheat at Turpan were higher for Method 3 than both LR methods, and RMSE for the accumulation of dry matter predicted by Method 3 was 0.1 for spring wheat for Northern Xinjiang. In contrast, BFM was more stable and precise than Method 3.  53 . Small changes in the cardinal temperatures could affect the accuracy of the methods in predicting developmental stages. A sensitivity analysis was used to determine the uncertainty of the effect of the three cardinal temperatures on the precision of the four methods.
One of the three cardinal temperatures changed by −4, −2, −1, 1, 2, and 4 °C, and the other two cardinal temperatures remained unchanged. The CVs of two predicted developmental stages (heading and maturity for wheat, jointing and maturity for corn) since sowing at three stations (Turpan for spring wheat, Xianyou for winter wheat, and Xinxiang for corn) were used for the analysis (Figs 6-8).
T b sensitivity. A change of T b from −4 to 4 °C increased CVs of the predicted developmental stages for the two crops for Methods 1 and 2; Method 3 was similar to Methods 1 and 2, and the CVs of Method 3 decreased in some cases. In contrast, a change of T b from −4 to 4 °C slightly increased the CVs of the predicted developmental stages for BFM. Sensitivity to T b was highest for Methods 1 and 2, followed by Method 3 and was generally lowest for BFM. T opt sensitivity. Only Method 3 and BFM are discussed in this section because T opt was not used for calculating GDD by Methods 1 or 2. Changes of T opt led to highly variable precision for Method 3 in some cases (e.g., CVs doubled when predicting the maturity stage for spring wheat in Turpan and corn in Xinxiang) (Figs 6 and 8).
In contrast, changes of T opt only slightly increased the CVs of the predicted developmental stages for BFM, which generally remained the most precise of the four methods.
T u sensitivity. Changes of T u had little effect on Methods 1-3 for the two crops, as expected, because the environmental temperatures in most of the growth periods were usually <T u . An increase in T u could improve the precision of predicting the developmental stages for wheat for BFM to some extent, especially for spring wheat in hot environments (e.g., Turpan).

Conclusions
Methods 1 and 2 were generally appropriate if T max < T opt for the crop. Method 2 performed as well as Method 1 in most cases but Method 2 was generally slightly better. For both crops, Method 3 and BFM generally performed better than Methods 1 and 2 at predicting the developmental stages since sowing and predicting the accumulation of dry matter with the normalized logistic model as a function of GDD. However, the stability of Method 3 was unsatisfactory because changes in the cardinal temperatures led to highly variable precision in some cases. BFM was sufficiently stable and more precise than the other methods for T opt < T max . This study used specific examples but the results can be applied to any environment or crop.