Modeling Phenological Phases of Winter Wheat Based on Temperature and the Start of the Growing Season

: The phenological phases of ﬁeld crops have shifted to earlier times in the Czech Republic in recent decades; additionally, they have shown correlations with temperatures from previous spring months. Using a thermal time model called PhenoClim, the correlations between temperatures and phenophases allow us to evaluate the strongest predictors (i.e., maximum temperature) and indicators of base temperatures and growing degree days for the selected phenophases of winter wheat ( Triticum aestivum L.). With the help of this model, it is possible to explain 0.6–0.82% of the phase variability and to estimate the onset of phenophases for the selected time period and stations (with the RMSE values of 9.4 days for jointing, 4.3 days for heading, and 5.3 days for full ripeness). To further reﬁne the modeled onsets of phenophases, we used satellite data, speciﬁcally the normalized difference vegetation index and the enhanced vegetation index 2 from MODIS; based on these vegetation indices, the start of the growing season (SOS) was determined. After including SOS to model PhenoClim, we modeled the onsets of phenophases, with average accuracies ranging from 6.2 to 15.2. By combining the thermal time model and remote sensing data, speciﬁcally the data concerning the determination of SOS, we can reﬁne the modeling of the onset of full ripeness in some locations.


Introduction
Phenology is an important indicator for assessment of the impacts of climate change [1], owing to its relationship with air temperature [2,3]. Other factors influencing phenology include the length of the photoperiod and the availability of water [4]. Each phenological phase (phenophase) has optimal temperature limits that are required to reach specific phenophase [5]. Temperature models are usually based on the base temperature (T base ) and growing degree day (GDD) values for each phenophase [6]. A challenge in studying phenophases and their onsets estimation of the onsets of phenophases with the highest possible accuracies [6][7][8] for further use in growth models [9] under current and future climatic conditions [10,11]. Based on models of the onset of phenophases, it is possible to evaluate the length of the vegetation period [12] or the length of the period between the onsets of individual phenophases [8].
In recent decades, the use of satellite imagery has come to the forefront in phenology assessments. Both optical [13] and radar satellite data [14] are used for this purpose. Some studies have combined both sets of satellite data to identify phenophases and achieve more accurate results than the results obtained using each set separately [15]. Satellite images enable evaluation of the growth and development of vegetation on a much larger spatial scale than can be achieved with ground observations. However, ground observations (by humans or cameras) play an important role and serve to calibrate satellite data [7].
Phenological metrics (phenometrics) obtained by the functional analysis of vegetation indices (VIs) are used As indicators of the timing of phenophases. Basic phenometrics include the start of the growing season (SOS), also called the green-up date, onset of greenness, or spring phenology; the length of the growing season (LOS); and the end of the growing season (EOS), also called the end of senescence, end of greenness, dormancy, or autumn phenology [16].
Phenological studies most often describe shifts in the onsets of phenophases to an earlier time and their relationships with climatic variables (most often temperature); however, recent studies have dealt with chilling units and the length of the day [17]. Many studies have evaluated wild plant species [18][19][20], and some have reviewed field crops, e.g., rape [21,22] and rice [23]. Fewer publications have examined the relationships between phenophases and satellite data in field crops [13,24,25] than those that have examined wild plants. Therefore, in this study, we focused on the phenology of a field crop (winter wheat), in particular, the phenophases of jointing, heading, and full ripeness. The main objectives of this study are (1) to evaluate the trends in the onsets of selected phenophases for field crops (winter wheat), (2) to establish the most accurate temperature model for determination of selected phenophases, (3) to determine the dates of phenophases (jointing, heading, and full ripeness) based on temperature models for selected areas, and (4) to determine the start of the growing season (SOS) using VIs and its relationship with the temperature model.

Stations and Phenological and Meteorological Data
Phenological data from 28 stations in the Czech Republic at varying altitudes (from 171 to 647 m a.s.l.) were used in the present study, covering the period from 1961 to 2021. The experimental stations represent the soil and climatic conditions of the region; the same set of experimental fields was used in almost all cases. These data were visual ground-based phenological observations of winter wheat (Triticum aestivum L.) from official state trials of the registered and considered cultivars at small-plot experimental stations of the Central Institute for Supervising and Testing in Agriculture (Figure 1). Phenological dates were assessed by highly trained staff. Because the same set of 30-40 winter wheat cultivars (depending on the year) was grown at each site, the so-called standard cultivar phenological dates were used. The standard for the variety trials was always selected with respect to the previous standard cultivar. The dates of onset of three phenophases were evaluated: jointing (BBCH31), heading (BBCH51), and full ripeness (BBCH89). The completeness of data in individual years depended on each station and phenophase. For some stations, data were available for the entire period of 1961-2021, whereas for some stations, data were unavailable in some years. In most cases, data were available for the last thirty years or more, but at several stations, only a short-term series of observations (approx. 11-17 years) was available. Long-term series of observations are generally used as part of phenological studies (approx. 20-30 years), but in this study, a short-term series of observations was also used in order to include as much accurate information as possible about the onsets of phenophases in the model in an effort to increase the spatial resolution of the stations. The geographical coordinates of all stations and the availability of phenological data at each of them are shown in Table 1. A detailed overview of the available years is shown in Table A1.
Correlation coefficients (r, Pearson coefficient) were used as the primary indicators of the strength of the relationships between given variables. The trends represent the slope of the linear regression between the phenological date and the year. Any significance in the observed trends was assessed using a t-test. All tests were performed with the statistical/programming tool R 3.6.1. [26] and AnClim software [27].

PhenoClim Model
PhenoClim is a model that allows for modeling of the onsets of phenophases in locations or time periods with only meteorological parameters [10,32]. It is a temperature model or a so-called thermal time model that works with one climate variable (T max , T mean , and T min ) and uses statistical variables to determine the most suitable values of GDD and T base . The beginning of temperature summation during model calibration was determined using temperature conditions (T min = 0 • C, T max = 5 • C, and T mean = 2.5 • C).
For each phenophase, the best meteorological predictor was searched for, which was determined using the root mean square error (RMSE). The values of T base and GDD are needed to reach a certain phenophase and were simultaneously calculated with the RMSE. Different calibration and validation datasets were used for the analysis, and a total of 168 model runs were calculated for each phenophase. Using the lowest RMSE value, the 10 best model runs were determined; these values were averaged, and the most suitable model was determined for each phenophase. PhenoClim made it possible to include the length of the photoperiod in the calculations and simultaneously consider the snow cover. According to the basic settings of the model, the photoperiod for winter wheat was defined as a short day of 7 h and a long day of 13 h following the approach used by Trnka et al. (2014) [33]. The beginning of the temperature summation was adjusted based on the day length coefficient. When the day length was less than 7 h, the value of the coefficient was 0, and there was no temperature summation. This phenomenon occurred only after crossing the 7 h mark. The snow cover was calculated using seven parameters, which were key for snow accumulation and melting [34]. Both of these functionalities reduced the resulting RMSE values and thus improved the accuracy of the model. PhenoClim was further used (based on the calculated models with the smallest error and the established SOS) to model the dates of the phenophases for sites for which there was no ground-based phenological monitoring.

Remote Sensing Data
The satellite data came from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor that was carried on the polar-orbiting Terra satellite, i.e., the MOD09GQ product, with a spatial resolution of 250 m, which was available through the LP DAAC Data Pool [35]. For the Czech Republic territory, four images were downloaded, and three layers were extracted from each image: two bands (red and infrared) and a quality layer. Each band was covered with a layer of quality, and the grids with the highest qualities were selected. After selection, the values of the normalized difference vegetation index (NDVI) and enhanced vegetation index 2 (EVI2) were calculated as follows: where ρRed and ρNIR are the reflectance values in the red and near-infrared (IR) bands, respectively. The factor G is determined by the c value, and c is derived by linearly fitting ρRed = c × ρNIR [36]. Based on the filter, which works with a ± seven-day time window, the values of the vegetation indices were interpolated in a daily step. Pixels containing only winter wheat vegetation in the period from 2015 to 2020 were selected near stations where visual ground-based phenological observations were available. Based on information from the Land Parcel Identification System (LPIS) database, the plots and pixels containing only with winter wheat were identified. A minimum of five and a maximum of ten available pixels within a radius of five kilometers from the station at approximately the same altitude were required for further analysis of a station. NDVI and EVI2 values were calculated for these pixels. The index values were averaged, and the average values for each station were smoothed using the Gaussian ordinate method with a period = 10. AnClim software [27] was used to smooth the curves, facilitating the used of this method. The adjusted daily index values were converted into Z-score values. A Z-score threshold between 0.00 and 0.60 (in the 0.05 step) was determined at eight stations because data were available for these stations each year between 2015 and 2020. Each of these threshold values corresponded to a threshold for determination of a possible SOS. The conversion of VIs values to Zscore values made it possible to take into account the interannual variability of VI values when determining the SOS. The SOS was determined separately for each station and year. This phenomenon indicated that for each year at one station, 13 days were determined as possible days for the SOS. The established SOS for the individual thresholds was used in the PhenoClim model as the day of the start of temperature summation to calculate the onset of the phenophase. The onsets of phenophases for each threshold in the time period of 2015-2020 were modeled based on the best predictor, GDD, T base , and SOS. The onsets of the phenophases modeled in this six-year time period were compared to ground observations. For modeled onsets of phenophases, the RMSE value was calculated for each threshold, and the best threshold of the Z-score coefficient was determine based on the lowest value.

Results
The first outputs revealed trends in the phenophase dates of winter wheat. At almost all studied stations (28 stations overall, with three phenophases each), phenophases were shifted to earlier times, and at some stations, these shifts were statistically significant. The lowest change and significance levels in trends were detected for the phenophase jointing, which was mostly caused by short observation periods (e.g., Kromeriz: 12 years or Sedlec: 13 years). However, the heading and full ripeness phenophases were observed for longer periods (reaching 61 years), and trends were frequently significant. The dates of observed phenophases (jointing, heading, and full ripeness) advanced by an average of 2.1-3.4 days per decade ( Table 2). The shifting of phenophases to an earlier date differed according to the given phenophase but was observed at almost all stations. Shifts to later dates were observed at only three stations (nonsignificant). The largest significant trends (by 5.9-6.9 days per decade) were calculated for the jointing phenophase (at three stations), as well as the full ripeness phenophase (at Stupice station). The most significant trends were calculated for the full ripeness phenophase (trends were detected for 18 stations among a total of 28 stations); however, for the jointing phenophase, only 11 stations presented with significant trends (Table 2).
According to the lowest RMSE value, the best predictor for the onsets of phenophases was T max calculated by the model PhenoClim for all three phenophases. The value of RMSE (based on the average of the 10 best models) moved in a range of 4.3-9.4 days. The lowest RMSE (4.3 days) was detected by the set of models for the heading phenophase, and the largest RMSE (9.4 days) was detected for the jointing phenophase. Two main parameters subsequently used for the modeling of the terms of phenophases were T base and GDD (Table 3).
Using the model setting (the best predictor, GDD, and T base ) for each phenophase, the onsets of phenophases were modeled for the 28 experimental stations for the period from 1961 to 2021. The lowest correlation (between observed and modeled terms of phenophases) was calculated for the jointing phenophase (r = 0.45-0.87), with only one station with no significance. The correlations for the remaining two phenophases (heading and full ripeness) between observed and modeled onsets were high and significant at all stations (r = 0.68-0.96) ( Table 2). A comparison of observed and modeled phases is displayed for station Chrlice for illustration in Figure 2. The line graphs on the left side of Figure 2 show the time course of observed and modeled onsets of individual phenophases in the period from 1961 to 2021. The scatterplots on the right side of the figure show the relationship between observed and modeled onsets of individual phenophases. Table 2. Linear trends for observed jointing, heading, and ripening phenophases among the 28 observation stations. Positive values of trends indicate a shift in the onsets of phenophases to a later date, and negative values indicate a shift to an earlier date. The values of the correlation coefficients describe the relationships among observed and modeled terms of the phenophases and, which were statistically significant in all cases. * Significant trend at α = 0.05; ** significant trend at α = 0.01; *** significant trend at α = 0.001.  Table 3. RMSE values for given phenophases with the corresponding growing degree day (GDD) and base temperature (T base ) values for the best predictor (T max ) as calibration results of the PhenoClim model. R 2 is the value of the coefficient of determination, and r is the value of the correlation coefficient. Based on the NDVI and EVI2 values, the SOS was determined for eight experimental stations in individual years during the period from 2015 to 2020 for each phenophase. After including SOS in the model, the RMSE value was calculated for the modeled onsets of phenophases, based on which the best threshold of Z-score was determined. Among the average values of eight stations, the best threshold for both vegetation indices was found to have a Z-score of 0. The average RMSE value for this threshold moved within the range of 6.2 to 15.2 ( Table 4). The lowest average RMSE value was found for the full ripeness phenophase by using the SOS determined by EVI2. For the jointing and heading phenophases, improved results were obtained by using SOS from NDVI; nevertheless, the average RMSE values were higher than for the full ripeness phenophase. The best threshold Z-score and corresponding RMSE value differed for individual stations; an overview is shown in Tables A2 and A3, highlighting that for the jointing phenophase and full ripeness phenophase at some stations, the accuracy of the model was improved when SOS was included relative to when it was excluded. A comparison of the observed and modeled values with and without the start of the season is displayed for the Chrlice station for illustration in Figure 3.

Discussion
Ground observations are among the most accurate methods for determining phenophases and are necessary for the calibration and evaluation of the model. In order to ensure that the model to considered as much accurate information as possible about the onsets of phenophases and that this was high-spatial-resolution information, we decided to use phenological data from stations with available short-term observations (e.g., Sedlec (13 years)), in addition to phenological data from stations with available long-term observations. Previous studies dealing with the phenology of winter wheat showed changes in trends comparable to those reported here, e.g., app. by 1.1-2.7 days/decade during the period from 1981 to 2009 [37] or by 0.8-1.8 days/decade during the period from 1935 to 2004 [38]. We observed an average shift of 2.1-3.4 days during the period from 1961 to 2021. Which may be the reason for a more pronounced shift compared to the mentioned studies given the global increase in temperature in recent decades? One reason could be a different time period that was evaluated. Compared to the aforementioned studies, in this study, we included data from the last 30 years in our analysis at most stations, with a mor pronounced increase in the temperature than in previous years; temperature is considered one of the main factors influencing the onset of phenophases. Earlier onsets of phenophases have also been reported in other field crops. In Germany [22], in the period from 1979 to 2013 shifts in winter rape of 4.8 days/decade and in winter rye by 1.3

Discussion
Ground observations are among the most accurate methods for determining phenophases and are necessary for the calibration and evaluation of the model. In order to ensure that the model to considered as much accurate information as possible about the onsets of phenophases and that this was high-spatial-resolution information, we decided to use phenological data from stations with available short-term observations (e.g., Sedlec (13 years)), in addition to phenological data from stations with available long-term observations. Previous studies dealing with the phenology of winter wheat showed changes in trends comparable to those reported here, e.g., app. by 1.1-2.7 days/decade during the period from 1981 to 2009 [37] or by 0.8-1.8 days/decade during the period from 1935 to 2004 [38]. We observed an average shift of 2.1-3.4 days during the period from 1961 to 2021. Which may be the reason for a more pronounced shift compared to the mentioned studies given the global increase in temperature in recent decades? One reason could be a different time period that was evaluated. Compared to the aforementioned studies, in this study, we included data from the last 30 years in our analysis at most stations, with a mor pronounced increase in the temperature than in previous years; temperature is considered one of the main factors influencing the onset of phenophases. Earlier onsets of phenophases have also been reported in other field crops. In Germany [22], in the period from 1979 to 2013 shifts in winter rape of 4.8 days/decade and in winter rye by 1.3 days/decade were reported; in a study of oats in the period from 1959 to 2009, researchers reported a shift in an average of 0.2-3.3 days/decade [39]. The results of these studies showed that the shifts in the onsets of phenophases are species-specific. The abovementioned changes in the onsets of phenophases and in the lengths of time between individual phenophases led to reductions in the growth period and negatively impacted the yields of field crops [40]. Given that temperature affects the onsets of the phenophases of field crops [37], knowledge about the effects of climate change on their phenological development is needed to optimize the farming system and increase the productivity of field crops [1]. Additionally, owing to climate change, increasing temperatures, and potentially more frequent occurrences of extreme temperature conditions, it is necessary to adapt to these conditions and mitigate their impact. One possible solution is the inclusion of varieties with higher temperature requirements [40,41] and increased resistance to drought and heat [42] in the sowing procedure. The results of studies dealing with the temperature requirements of field crops or changes in the onsets of phenophases can thus be crucial for the decision making of breeding companies. Another possible measure to mitigate the negative impacts of climate change is shifting the date of sowing [41,42]. The mentioned measures are mainly based on the relationship between phenology and temperature; temperature is generally considered to be one of the main factors influencing the onsets of phenophases. In our case, T max was identified as the best predictor of the onsets of phenological terms. In addition to air temperature, photoperiod plays an important role in plant phenology, which is important in spring phenology to prevent damage from spring frosts and in autumn phenology to stop plant growth [2]. Considering these two main factors influencing the onsets of phenophases and in an effort to obtain more robust results, photoperiod was also included in the calculations during the calibration of the PhenoClim model. In addition to climatic conditions, the phenologies of field crops are influenced by changes in management, although few studies have been conducted on this topic. Such influences include the type of cultivar (early vs. late), which is related to the density of planting and the date of sowing [22]. For example, a later sowing date can have distinct effects on winter (positive) and spring (negative) wheat [41]. In this case, we did not include other factors in the calculations, but they could be the subject of follow-up studies to provide more realistic information about plant behavior.
The results of the PhenoClim model based only on temperature thresholds and photoperiod showed relatively high accuracies in terms of predicting heading and full ripeness, similar to the results of a study conducted in southeastern USA [6]. Using the same temperature thresholds, researchers modeled three phenological stages for eight winter wheat cultivars during the period from 1999 to 2010 and achieved a similar accuracy (average RMSE in the range of 5.2-7.1 days), as in our study. The results of studies focused on predicting the onset of key phenophases based on GDD and T base could be helpful for crop management practices, such as fertilization, irrigation, pesticide applications, and harvest scheduling. Possibilities of further use of predicting the prediction of the onsets of phenophases of field crops were presented in a study conducted in Denmark [43] based on field experiments of spring barley and winter wheat during the period from 1991 to 2018 to assess the viability of growing cover crops on the basis of harvest date prediction in the current season. In the first step, the ability of the model to predict the maturity date was determined with an accuracy RMSE of 5.5 days for both evaluated crops. Then, based on a calibrated sum of temperatures between maturity and harvest (400 • C for barley and 450 • C for wheat), the model predicted the harvest date and, with an RMSE error of 6.6 days for spring barley and 9.3 days for winter wheat. The model was capable of predicting the date and harvests, which could be helpful for planning autumn work in the field and the possible establishment of cover crops [43].
Finally, a combination of the PhenoClim thermal model and the SOS derived from remote sensing showed that it was possible to model the onset of the full ripeness phenophase with a relatively high accuracy. At half of the assessed stations, even higher accuracy was achieved when SOS was considered than when only the thermal model was used. A similar accuracy in estimating the maturity date of winter wheat was achieved by combining the sum of temperatures and EVI2, achieving an RMSE of 7.3 with a coefficient of determination value of 0.48 (i.e., r = 0.69) [25]. Furthermore, a method for monitoring and predicting heading and flowering dates of winter wheat was proposed based on a combination of temperature sums and NDVI [24]. Researchers used the dynamic threshold method to determine phenological metrics. The RMSE value for the prediction of the onset of heading was in the range of 4.76-6.13 days, and that for the estimation of flowering was in the range of 5.30-6.41 days. In this case, for the heading phenophase, by combining both methods, the accuracy of the model was one day less than that for the full ripeness phenophase, which can be considered a relatively satisfactory accuracy. The lowest accuracy of the model was for the jointing phenophase, probably as a result of the limited amount of phenological input data for this phenophase, given the overall poorer results in the other analyses. With information derived from remote sensing data, it is possible to include another variable in the temperature models reflecting the actual state of the vegetation on a large spatial scale. Additionally, this information can be derived practically in real time, enabling estimation of the onsets of phenophases in a given season. The results of this study could be of practical use for planning field work (e.g., fertilization and pest control).

Conclusions
Analysis of long-term phenological data of winter wheat in the period from 1961 to 2021 showed shifts in the onsets of phenophases to earlier dates for all shifts with statistical significance. The maximum air temperature was shown to be the best predictor of phenophase onsets for model tools. Using models based on temperature thresholds to determine phenophases, it was possible fundamentally approach to methods based on traditional ground observations. Combining remote sensing data with temperature models resulted in decreased accuracy of phenophase onset determination in most cases; however, in the case of the heading and full ripeness phenophase, this deterioration did not exceed an average of four days (for full ripeness, even better results were obtained at four stations), indicating a still high ability of the model to predicting the onsets of phenophases. Acknowledgments: The contributions of Petra Dížková and Jakub Bohuslav were financially supported by grant no. AF-IGA2021-IP036. The authors are grateful to Reinhart Ceulemans, (University of Antwerp, Belgium) for a thorough prereview of the paper and critical comments that considerably improved the manuscript. The Central Institute for Supervising and Testing in Agriculture provided the ground-based phenological data set.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. An overview of experimental stations of the Central Institute for Supervising and Testing in Agriculture with phenological observations and their geographic coordinates (altitude (m.a.s.l.), latitude, and longitude). In the right half of the table, the time periods in which ground observations were available for each phenophase and each station are listed. The period column indicates a specific time period or year. The year column is the total number of years of ground observations for each station.