Introduction

Temperature affects several physiological processes in plants. Generally, this influence is due to a thermal action on enzymatic activity, with temperature inducing changes in the conformation of enzymes and consequently on their functionality (^{Sharpe and DeMichele, 1977}; ^{Johnson and Thornley, 1985}; ^{Higley et al., 1986}; ^{Bonhomme, 2000}). The phenological response of plants to temperature is observed as changes in developmental rates (the occurrence of certain phenological events per time unit) (^{Raworth, 1994}). ^{Jaramillo and Guzmán (1984)} found that a significant correlation exists among the number of days elapsed from planting to first harvest, the average air temperature (°C) and the thermal units.

Heat accumulation is referred to as "heat units", "thermal time", "cumulative growing degree days" or "physiological time" and reflects the concept of a quantifiable relationship between temperature and the rate of crop development (^{Gordon and Bootsma, 1993}; ^{McMaster and Wilhelm, 1997}; ^{Bonhomme, 2000}; ^{Snyder et al., 2013}). This relationship has been successfully used in agricultural science to predict and quantify the time between the plant phenological stages (^{Gilmore and Rogers, 1958}; ^{Cross and Zuber, 1972}; McMaster, 1993; ^{Cardina et al, 2007}).

In the study of phenology and crop development, the concept "physiological time" can be expressed in quantifiable terms, with this time measured in growing degree-days (GDD, °C-days) (^{McMaster and Wilhelm, 1997}; ^{Rodríguez et al., 2012}) and representing the thermal units recorded between the minimum (Tb) and maximum (TB) threshold temperatures at the hours of a day (^{Snyder et al., 2013}). Mathematically, the daily GDD are obtained by the integration of growing degree-hours (GDH, °C-hours) for each hour of the day divided by the number of hours in a day (24 h) (^{Roltsch et al., 1999}; ^{Cesaraccio et al., 2001}; ^{Souza et al., 2011}; ^{Rodríguez et al., 2012}; ^{Snyder et al, 2013}).

However, hourly data are not always available, for example, when historical records are used or automatic weather stations are scarce. In such cases, the GDD are estimated using mathematical methods that use only daily maximum (TM) and minimum (Tm) temperatures (^{Snyder et al., 1999}). Among the methods most frequently used for calculating degree-days are the rectangle (^{Arnold, 1959}), triangle (^{Lindsey and Newman, 1956}) and sine wave (^{Baskerville and Emin, 1969}) methods, including their variants (^{Villa-Nova et al., 1972}; ^{Ometto, 1981}; ^{Snyder, 1985}).

According to ^{Roltsch (1999)}, when the TM exceeds the TB, these methods can use two approximations of accumulated thermal units: (1) the vertical cutoff, which assumes that there has been no accumulation of thermal units and (2) the horizontal cutoff, which assumes that the thermal units continue to accumulate until the Tb is reached. The threshold temperatures are the values under which (Tb) or above which (TB) the development rate is zero.

However, the threshold temperature is only a statistical value, which may be distant from the "physiological temperature" for which the development rate is close to zero, essentially because its value may vary depending on the method used for its calculation, on the growth stage or on the physiological process analyzed (^{Wang, 1960}; ^{Durand et al., 1982}; ^{Bonhomme, 2000}; ^{Litschmann et al., 2008}).

In crops such as Brussels sprouts, cabbage, parsley, legumes, fodder, corn, soybeans and tomatoes, the Tb may oscillate between 0 and 10°C (^{Gordon and Bootsma, 1993}). For Coffea arabica L. in Brazil, ^{Lima and Silva (2008)} estimated the Tb and TB at 12.9 and 32.4°C, respectively, from transplanting to the first flowering event, whereas ^{Pezzopane et al. (2008)} obtained a Tb value of 10.2°C from flowering to harvesting based on several harvest cycles using the "Mundo Novo" variety. For most tropical plants, including coffee, the Tb and TB values have been defined, respectively, as 10°C (^{Pedro-Junior et al., 1977}; ^{Jaramillo and Guzmán, 1984}) and 32°C (Jaramillo and Guzmán, 1984; ^{Hatfield and Prueger, 2015}).

^{Roltsch et al. (1999)}, ^{Souza et al. (2011)}, ^{Rodríguez et al. (2012)} and ^{Kean (2013)} have found variation in performance when comparing different methods to estimate GDD based on maximum and minimum daily temperatures and when applying those methods to different growing regions. These disparities become relevant, for example, should the estimated accumulation of degree-days be greater than the "real" accumulation; this example would imply that a shorter chronological time than estimated would be necessary to reach a certain phenological stage, or vice versa (^{Bryant et al., 1998}), which can affect the predictions of phenological models based on "physiological time" (^{Kean, 2013}).

Given the importance of the degree-day calculation in phenological models, the objective of the present research was to evaluate the performance of four horizontal cutoff methods in estimating degree-days within eight Colombian coffee-producing areas, using thermal thresholds associated with tropical plants.

Materials and methods

Study area and meteorological data

Eight localities with contrasting environmental conditions were selected, espatial distribution and economic importance for Colombian coffee crop was considered. Each locality was equipped with an automatic, RAWS-F, Fire Weather, Campbell Scientific^® Remote Automated Weather Station from the Coffee Meteorological Network, a part of the Colombian National Federation of Coffee Producers (Federación Nacional de Cafeteros de Colombia) (Tab. 1, Fig. 2). The air temperature was recorded 24 h a day between February 1, 2014, and January 31, 2015, with five-minute intervals between measurements. The differences in the number of daily logs for the meteorological stations, due to the days with missing data are shown in Table 1. The temperature sensors were previously calibrated, by means of the homologation in parallel between conventional and authomatic stations, with amaximum error of ± 0.5 °C.

TABLE 1 Geographic location of the weather stations used in the study; average annual historical values for precipitation (P), maximum temperature (T _MAX), temperature (T _MED), minimum temperature (T _MIN) and mean relative humidity (RH); and the number of days used in calculating the degree-days for each locality (ND). 

FIGURE 1 Location in Colombia of the weather stations of the Coffee Meteorological Network used in the study. 

FIGURE 2 Description ofthermal units (TU) integrated in the AUC for GDD _TT (a), GDD _Tb (b) and GDD _TB (c)  

TABLE 2 Indirect methods of degree-day estimation used in the analysis. 

Where: TM = maximum daily temperature (ºC), Tm = minimum daily temperature (ºC), Tb = minimum threshold temperature (ºC), TB = maximum threshold temperature (ºC), GDD = growing degree-days according to the indirect method (°C-days), GDD1 = degree-days above Tb (°C-days), and GDD2 = degree-days above TB (°C-days).

Degree-day calculation

For all the developed methods, the Tb and TB values were considered as 10 and 32°C, respectively, following the proposal by ^{Pedro-Junior et al. (1977)}, ^{Jaramillo and Guzmán, (1984)}, ^{Camargo and Pereira (1994)} and ^{Hatfield and Prueger, (2015)} for coffee cultivation.

Numerical integration by trapezoidal rule (reference method)

The integration method was used as the reference method against which the indirect methods were evaluated due to its greater estimation precision and its comparatively frequent data collection at relatively short intervals (^{Cesaraccio et al., 2001}; ^{Souza et al., 2011}; ^{Rodríguez et al., 2012}).

The estimate of GDD using the reference method (GDD_r) was calculated as follows:

Where GDD_r = degree-days of the reference method (°C-days), GDD_TT = total degree-days (°C-days), GDD_TB = degree-days above TB (°C-days), and GDD_Tb = degree-days below Tb (°C-days).

To obtain the values that compose expression (1), we used expression (2), in which the areas of individual trapezoids constructed from the temperature records were calculated to apply the "trapezoid rule". These areas were integrated into an area under the curve (AUC), which was divided by the total number of seconds per day (86400), thus obtaining the daily thermal units, representing the total degree-days per day (GDD_TT) in this case. To calculate the GDD_TB and GDD_Tb, the AUC was integrated, assuming that the temperatures (the m_i values in expression 2) were, respectively, higher than or equal to the TB and below or equal to the Tb at corresponding times (Fig. 2).

Where TU = thermal units, m¡ = temperature in the i^th measurement, and t_¡ = i^th time in the daily logs.

Indirect methods evaluated

The daily "cumulative growing degree days" or "physiological time" in GDD was estimated using the daily maximum and minimum temperature and four indirect methodologies formerly described in the literature (^{Arnold, 1959}; ^{Villa-Nova et al, 1972}; ^{Ometto, 1981}; ^{Snyder, 1985}) (Tab. 2), some of them previously applied to coffee production (^{Lima and Silva, 2008}).

Where: TM = maximum daily temperature (°C), Tm = minimum daily temperature (°C), Tb = minimum threshold temperature (°C), TB = maximum threshold temperature (°C), GDD = growing degree-days according to the indirect method (°C-days), GDD₁ = degree-days above Tb (°C-days), and GDD₂ = degree-days aboveTB (° C-days).

Statistical analysis

A linear regression model of the observed values was adjusted to each loca lity and indirect method, fitting the reference method data (dependent variable, Y) as a function of the data obtained with the indirect method (independent variable, X). The significance of the regression coefficients was evaluated using a t-test with α=5%. To identifying the best adjustment, we also evaluated whether the regression slope coefficient (β₁) differed statistically from one to determine whether the indirect method overestimated or underestimated the GDD. Additionally, the coefficient of determination (R²) was calculated.

The bias of the degree-days accumulated by each conventional method was es timated as the percent bias (PB) of the residual (expressed as percentage), as shown by expression (24).

Where: GDD_r = degree-days of the reference method (°C-days); GDD = growing degree-days according to the indirect method (°C-days).

The statistical analyse s were performed using SAS software version 9.4 (SAS Institute, 2012).

Results and discussion

The results of adjusting the regression model with the data from the reference method expressed as a function of the data obtained from each of the indirect methods are presented in Table 3 for each weather station. The regression coefficients (β₁) were significantly different from zero for all the methods by the t-test at the 5% probability level, indicating that the values obtained with each of the indirect methods help explain the results observed with the reference method. The proportion of the data variability observed with the reference method that was explained by the indirect method varied from 32 to 82%, according to the coefficients of determination (Tab. 3).

Table 3 Regression coefficient (β1), standard error (STE), coefficient of determination (R2) and percent bias (PB) obtained in the relation between accumulated degree-days in each indirect method and the accumulated degree-days estimated by numeric integration during the evaluated period in eight coffee-growing areas of Colombia. 

SE = standard error; a = t-test with alpha 5% for β₁ values significantly different from zero; b = t-test with alpha 5% for β₁ values significantly different from one.

The regression co efficient was equal to one (it means, a 1:1 ratio existed between the indirect method and the reference method) for the i ndirect method of ^{Villa-Nova et al. (1972)} and for the localities of Naranjal (β₁ = 0.960°C-days) and El Rosario (β₁ = 1.053°C-days) (t-test, α = 5%). The coefficient was significantly different from 1 for all the other localities.

Regarding goodness-of-fit, R² for the methods proposed by ^{Arnold (1959)}, ^{Ometto (1981)} and ^{Snyder (1985)} exceeded 0.73 for the El Rosario, El Sauce, Jorge Villamíl, Naranjal and Paraguaicito localities and varied between 0.59 and 0.68 for Julio Fernández and Pueblo Bello. For every locality assessed, the R² obtained for the method of ^{Villa-Nova et al. (1972)} was lower compared to all the other methods. The lowest adjustment values for R² were recorded for a locality in Bertha (between 0.32 and 0.44).

The PB showed that when the degree-day errors accumulated, for the methods of ^{Arnold (1959)}, ^{Ometto (1981)} and ^{Snyder (1985)} overestimated the GDD relative to the reference method by between 8.3 and 17.2% at all eight localities. Under the thermal thresholds evaluated, these three methods were also consistently observed to show no significant differences. The method described by ^{Villa-Nova et al. (1972)}, with negative PB values, was the only one to underestimate the reference method results, ranging between 2.0 and 20.6%.

In the present study, the methods of ^{Arnold (1959)}, ^{Ometto (1981)} and ^{Snyder (1985)} showed no descriptive differences by locality regarding the regression coefficient, R ² or PB values (Tab. 3). ^{Souza et al. (2011)} also reported that these three methods produce similar results, which is due to the similar geometric forms usedby these methods for estimating the degree-days based on Tb and TB.

According to the above, these methods did not present differences, because most of the minimum temperatures were higher than the Tb in the daily records. However, in a very few cases, the maximum temperatures exceeded the TB. This occurrence causes the methods of ^{Ometto (1981)} and ^{Snyder (1985)} to consider expressions based only on the Tb, a situation bearing similarity to the arising results with the method of ^{Arnold (1959)}, which only includes Tb in its mathematical expression. According to ^{Higley et al. (1986)}, even though omitting a maximum development threshold could decrease the precision of the degree-day estimation, the introduced error is low as long as the daily maximum temperatures are generally lower than the maximum threshold temperature for development (TB).

The evaluation criteria did indicate different performance among the localities for themethods of ^{Arnold (1959)}, ^{Ometto (1981)} and ^{Snyder (1985)}. For example, higher coefficients of determination (R ² ) were observed for the localities of El Rosario, El Sauce, Jorge Villamíl, Naranjal and Paraguaicito compared to those on Bertha, Julio Fernández and Pueblo Bello (Tab. 3), this implies a better fit for the first. ^{Carlson and Hancock (1991)} and ^{Roltsch et al. (1999)} highlight the importance of R² as an evaluation criterion to determine the degree of fit between methods.

The PB also varied among the localities, although the variation was inconsistent with that described for the adjustment. Comparatively high PB values were obtained for Bertha and El Sauce, where as those for El Rosario, Jorge Villamíland Pueblo Bello were low (Tab. 3). The greater PB for Bertha could be attributed to the lower temperatures recorded in this area (Tab. 1). By contrast, the greater PB for El Sauce could be mostly due to variation in cloudiness and wind speed, which affect the daily temperature fluctuation. ^{Roltsch et al. (1999)} and ^{Worner (1988)} found that when degree-days are estimated by indirect methods, a greater error is found for cold areas or climates in which fast thermal changes occur. Both error and bias are important parameters to consider because the prediction of phenological events is linked to the accumulation of physiological time (^{Rodríguez et al., 2012}).

The accumulated error evaluated by the PB did not exceed 17.1°% compared to ^{Arnold (1959)}, ^{Ometto (1981)} or ^{Snyder (1985)} methods in any locality (Tab. 3). This error could be considered acceptable for the prediction as regards to the PB found by ^{Rodríguez et al. (2012)} for Colombia (between 7.13 and 30.57%). The different errors involved may cancel each other to some extent, but the greatest source of error, which is thermal accumulation, results from inaccurate temperature estimation (^{Pruess, 1983}). Those methods with positive PB values overestimated the degree-days relative to the ones calculated by the reference method, which was also reported by ^{Souza et al. (2011)}.

An overestimation can be convenient, considering that an underestimation of accumulated degree-days can lead to erroneous conclusions regarding a phenological event in poikilothermic organisms by inaccurately predicting a later event than actually occurs (^{Bryant et al., 1998}). Plant development may be further affected by not only nutritional status and water stress but also photoperiod in long-day plants and, to a lesser extent, in short-day plants (^{Bonhomme, 2000}). Therefore, any deleterious effects of these factors should be reduced for cash crops, with all required substrates present in amounts adequate to ensure growth (^{Higley et al, 1986}).

The method proposed by ^{Villa-Nova et al. (1972)} was the only one to present regression coefficients statistically equal to 1 and to underestimate the reference method, contrasting with the results reported by ^{Souza et al. (2011)}.

Conclusion

Considering the PB, the methods proposed by Arnold, Ometto and Snyder tend to overestimate thermal time, whereas the method proposed by Villa-Nova underestimates thermal time, but with a lower performance regarding to the previous ones. However, the performance of each method can vary between zones due to the agri-environmental conditions typical of each locality. Arnold's method can be taken into account when daily temperatures do not exceed the maximum or minimum threshold considered, as in the present study. The use of the different methods depends on the available information and the objective of the thermal time estimation.