Cascade prediction of winter wheat developmental stages considering snow cover period in Hokkaido, Japan

Abstract

　We propose a model to predict winter wheat developmental stages in Hokkaido, a snowy region in Japan. The model uses regression equations incorporating day of the year (DOY) data of multiple developmental stages preceding the prediction target stage and air temperature. Using DOY of snow melting dates, the timing of snow disappearance is considered explicitly. We predicted dates of panicle formation, heading, flowering, and maturity of winter wheat ‘Kitahonami’. We also applied the conventional developmental rate (DVR) model, which considers the vernalization requirement of winter wheat, and compared the obtained prediction accuracy against that of our proposed model. For the DVR model, the root mean square errors (RMSEs) were, respectively, 3.6, 1.9, 2.5, and 3.1 days for panicle formation, heading, flowering, and maturity. By contrast, for the proposed model, they were 2.8, 1.9, 1.9, and 2.6 days, respectively. Mean absolute errors (MAEs) for the proposed model were also smaller than those for the DVR model. The proposed model shows improved prediction accuracy, indicating that the use of DOYs of multiple stages and air temperature is beneficial for predicting various developmental stages.

1. Introduction

For obtaining consistent yields with many crops, including wheat, carrying out various cultivation management tasks such as fertilization, chemical spraying, and harvesting according to specific developmental stages is beneficial. Predicting the developmental stages at which these management tasks will be implemented is therefore crucially important. Moreover, weather conditions during periods separated by developmental stages affect yield (Shimoda et al., 2022). Therefore, predicting the developmental stages is also important for yield prediction. Especially in Hokkaido, which produces about 70% of Japan’s wheat, obtaining higher yields by implementing appropriate cultivation management based on developmental stage prediction will be helpful for securing the stable supply of Japanese wheat. In Hokkaido, wheat is covered by several tens of centimeters of snow for about four months, from about two months after it is sown. When air temperatures rise and snow depths lessen, snow melting agents are sometimes applied to hasten the restart of wheat growth after snow has melted completely (Shimoda and Hamasaki, 2021). Therefore, building a prediction model for wheat development in Hokkaido necessitates consideration of the snow cover period, snow depth, and the timing of artificial snow melting when snow melting agents are sprayed.

Many methods have been proposed for predicting wheat developmental stages. These include formulations of the daily developmental speed as a function of temperature and photoperiod. Proposed models include an accumulated temperature model that uses a linear function (Shimazaki et al., 2022) and DVR models able to use both linear and nonlinear functions (Ritchie, 1991; Hunt and Pararajasingham, 1995; Wang and Engel, 1998; Yan and Wallance, 1998; Brooking and Jamieson 2002; Streck et al., 2003; Xue et al., 2004; Maruyama et al., 2010; Nakazono et al., 2014; Zheng et al., 2014; Kawakita et al., 2019; 2020a, b; Zhao et al., 2021). Among these efforts, Maruyama et al. (2010) and Nakazono et al. (2014) proposed a DVR model using air temperature and photoperiod as model inputs to predict developmental stages for Japanese wheat. Kawakita et al. (2020a) reported that adding a vernalization function to the models proposed by Maruyama et al. (2010) and Nakazono et al. (2014) improved their accuracy. By contrast, Shimazaki et al. (2022) achieved the same prediction accuracy as that reported by Maruyama et al. (2010) using an accumulated temperature model that changed the weighting of accumulation of air temperature depending on the season.

Prediction of the developmental stages of wheat in Hokkaido necessitates consideration of environmental factors such as air temperature and photoperiod, as well as the period during which wheat is covered by snow. By contrast, Maruyama et al. (2010), Nakazono et al. (2014), and Kawakita et al. (2020a) do not incorporate consideration of the snow coverage period. Shimazaki et al. (2022) proposed a simple approach in a preliminary study. For that approach, air temperature is not accumulated when the snow depth exceeds a certain threshold. However, snow depth was not incorporated into the final prediction model because the effects of accumulating of air temperature depending on snow depth were not stable. Therefore, the period during which wheat is covered with snow was not considered explicitly.

As described herein, we propose a model that predicts the timings of developmental stages using simple regression equations incorporating air temperature and days of the year (DOYs) of multiple developmental stages preceding the prediction target stage as explanatory variables. This study is intended to make it possible to predict developmental stages of wheat in Hokkaido using the proposed model, with consideration of the air temperature surrounding the wheat under snow and the timing of snow melting, and including artificial snow melting by the spraying of snow melting agents.

2. Models

For comparison, we also applied the DVR model reported by Kawakita et al. (2020a), which is regarded as having the best prediction accuracy of wheat developmental stages among existing prediction models for Japanese wheat. This section explains the equations of our proposed model and the DVR model reported by Kawakita et al. (2020a).

2.1 Proposed model

2.1.1 Model principle

In winter wheat cultivation, panicle formation and heading are the timings for fertilization. These stages, flowering, and the days near these stages are the timings for chemical spraying to defend plants from disease and lodging. Maturity is important for determining harvesting timing. For this study, we developed a model for predicting dates of panicle formation, heading, flowering, and maturity of winter wheat. The model used 1) day of the year (DOY) data of multiple developmental stages preceding the prediction stages and 2) air temperature, which is the most fundamental meteorological factor associated with crop development.

The reason for using DOYs of multiple stages is explained hereinafter. For predicting panicle formation, we used DOYs of the sowing date and snow melting date following the earlier research described below. Nakazono et al. (2014) reported that the sowing date affects the developmental rate up to the jointing stage in winter wheat, and predicted the jointing stage from the sowing date in warmer region at central Japan. By contrast, in Hokkaido, wheat is covered with snow after being sown. Actually, snow has an insulating and shading effect. After snow melts completely, however, wheat is affected directly by the air. Therefore, snow melting dates are useful as starting dates for developmental prediction. For Hokkaido, prediction models based on snow melting dates have been proposed (Hokkaido Research Organization, 2014; NARO, 2021). Following earlier research, we attempted to build a model using both sowing dates and snow melting dates directly. For predicting heading, flowering, and maturity using the proposed model, we also used the DOYs of multiple preceding stages. The stages used for each prediction are shown in Figure 1 and Table 1. Eguchi and Shimada (2000) reported that, after the panicle formation stage, the number of days from one developmental stage (stage 1) to the next (stage 2) was correlated with the DOY of stage 1 for several Japanese wheat varieties in western Japan. We obtained similar results for ‘Kitahonami’ in Hokkaido. Additionally, we found some correlation between the number of days from stage 2 to the next (stage 3) and the DOY of stage 1. These results are presented in Figure 4. Details are described in Subsection 4.1. According to these ideas, our proposed model uses DOYs of multiple preceding stages.

In addition, using the proposed model, we estimated temperatures surrounding the wheat under snow, not only based on air temperature but also based on snow depth, referring to work reported by Ritchie (1991).

Fig. 1.Prediction flow for each developmental stage. Symbols D₁ (days), D₂ (days), and D₃ (days) are defined respectively as the DOYs of stages 1, 2, and 3. Stages 1, 2, and 3 are developmental stages that occur in sequence. Symbol T₁ represents the average effective temperature during the period of D₁ – D₂. Symbol T₂ stands for the average effective temperature during D₂ – D₂+d, where d (days) is a hyperparameter.

Table 1. Target variable D₃–D₂ (days) and explanatory variables D₁ (days), D₂ (days), T₁ (°C), and T₂ (°C) in the proposed model. The value ranges of the respective variable for our data are also shown in “[ ]”

2.1.2 Model flow and equations

Figure 1 presents the predicted flow for each developmental stage in the proposed model. The model uses the following equation, with D₁ (days), D₂ (days), and D₃ (days) defined respectively as the DOYs of stages 1, 2, and 3. Stages 1, 2, and 3 are developmental stages that occur in sequence. The number of days between stages 2 and 3 (D₃-D₂) is predicted as

where T₁ and T₂ respectively represent the average effective temperatures from D₁ to D₂ and from D₂ to D₂+d. In addition, a₀, a₁, a₂, a₃, and a₄ are parameters. Also, d is a hyperparameter. It is not known how long of a period from the previous stage is important to estimate the progress of subsequent development by air temperature. Therefore, the period length is determined in advance using d. The determination methods for a₀, a₁, a₂, a₃, a₄ (a₀, a₁, and a₂ have no units; the units of a₃ and a₄ are days °C^-1) and d (days) are explained in Subsection 3.2.2.

T₁ and T₂ are calculated as

and

where

and

For Equation 4, it was assumed that growth does not proceed at temperatures below 0 °C, similarly to the equation in Kawakita et al. (2020a). In Equation 5, T_i (°C) and SD_i (cm) respectively denote the daily mean air temperature and the daily snow depth on day i. According to earlier studies, temperatures under the snow are influenced not only by air temperature but also by the snow depth (Aase and Siddoway, 1979; Fukuda, 1982; Tsuchiya, 1985; Ritchie, 1991; Hirota et al., 2006; Zheng et al., 2014; Zhao et al., 2021). The latter two terms of Equation 5 represent the snow coverage effect, by which the temperature under snow is higher than the air temperature because of heat stored in the soil referring (Ritchie, 1991). The proposed model uses these two equations for calculating T₁ for panicle formation prediction from the sowing date to the snow melting date.

Using Equations 1-5 for panicle formation prediction, the number of days from snow melting to panicle formation (D₃-D₂) was predicted using the DOY of the sowing date (D₁), the DOY of the snow melting date (D₂), the average effective temperatures from the sowing date to the snow melting date (T₁), and the average effective temperatures for d days after snow melting (T₂). The dates of heading, flowering, and maturity were also predicted. Table 1 presents the target variable D₃-D₂ and explanatory variables D₁, D₂, T₁, and T₂ for prediction for each stage.

For this study, snow melting dates were defined as the day on which more than half the area of the field had no snow. For wheat cultivation in snowy regions, to hasten the restart of wheat growth after snow has melted completely, snow-blackening by a snow melting agent, which facilitates absorption of solar radiation, is sometimes applied to accelerate snow melting (Shimoda and Hamasaki, 2021). Measurement data of snow melting dates used for this study include both natural melting dates by rising air temperatures and artificial melting dates affected by the application of snow melting agents. However, because few actual measurements of snow melting dates were available, the date was determined using the following steps. 1) If the measurements were available, then the date was used considering not only natural snow disappearance by rising air temperatures but also artificial snow disappearance by spraying snow melting agents. 2) If the measurements were not available, then the date was estimated based on the snow depth. We assumed x as the day when the snow depth reaches 0 cm in the National Agriculture and Food Research Organization (NARO) meteorological grid data (https://amu.rd.naro.go.jp) (Ohno et al., 2016; Kominami et al., 2019). If x was earlier than April 1, then the snow melting date was inferred as x. 3) If no measurement was available, and if x was later than April 1, then the snow melting date was inferred as April 1, assuming that the snow melting agent was applied before April 1.

2.2 DVR model considering vernalization requirements of winter wheat

Kawakita et al. (2020a) predicted emergence dates from sowing dates. They then predicted flowering dates from emergence dates. For emergence prediction, the base temperature was set as 0 °C. The emergence dates were determined as dates on which the accumulation of the effective daily mean temperature becomes greater than 110 °C day after the sowing date. For flowering prediction, the DVR model considering the vernalization requirement of winter wheat was used. DVR was represented by the following equations as

where T_i (°C) and P_i (hour) respectively represent the daily mean air temperature and photoperiod on day i. In addition, f (T_i), ɡ (P_i), and f_ver (T_i) respectively stand for the temperature, photoperiod, and vernalization response functions. DVR=0 if P_i < P_b. f_ver (T_i) varies from 0 to 1. T_b (°C) and P_b (hour) respectively signify the base air temperature and photoperiod. Model coefficients C, C_t, C_p, and C_ver are used. The inflection point of the inverse sigmoidal curve, denoted as T_bv (°C), was set as 8 °C for ‘Kitahonami’. Finally, VD_full stands for the number of effective vernalization days for the plant to be fully vernalized; VD_base represents the base effective vernalization days. They were set respectively as 43 and 8.

Kawakita et al. (2020a) predicted emergence dates using the accumulated effective temperature model. They then predicted flowering dates using the DVR model from the predicted emergence. For this study, we also predicted the emergence dates using the accumulated effective temperature model, and then predicted panicle formation, heading, flowering, and maturity dates using the DVR model. We considered the following four developmental stage periods after emergence: 1) from emergence to panicle formation, 2) from emergence to heading, 3) from emergence to flowering, and 4) from emergence to maturity. In the DVR model, the timing of a developmental stage can be described as the day on which the developmental index (DVI), represented by the accumulation of DVR, exceeds a certain value. We set DVI as 0 and 1, respectively, at the beginning and end of each period. The number of days in each period Y (days) was estimated as

and

where DVR_i represents the daily DVR on day i, and where Ŷ is the estimated Y.

3. Materials and Methods

3.1 Data used for this study

The most cultivated winter wheat variety in Hokkaido is ‘Kitahonami’, which accounts for 80% of the winter wheat area in Hokkaido (Hokkaido Prefecture, 2023). After analyzing data of multiple developmental stages of ‘Kitahonami’, we predicted four developmental stages: panicle formation, heading, flowering, and maturity. Table 2 presents the number of data, the minimum, maximum, mean, and standard deviation of the DOYs for each stage. Panicle formation is the day when the panicle length of the main culm reaches an average of 2 mm. The heading stage is the day when 40-50% of all stems have headed. The flowering stage is marked by the day when 40-50% of the whole cultivated area has bloomed. The maturity stage is the day on which the grain moisture content is approximately 40%. These data were collected at the Hokkaido Research Organization (HRO) Chuo (Naganuma Town: 2006-2021), Tokachi (Memuro Town: 2004-2021), and Kitami (Kunneppu Town: 2005-2021) stations, at the National Agriculture and Food Research Organization Hokkaido Agricultural Research Center (NARO/HARC) in Sapporo City and Memuro Town (2011-2020) (Shimoda and Hamasaki, 2021), and at several fields in the Sorachi, Kamikawa, and Tokachi areas (2007-2022) (JA Bibai and JA Minenobu, 2015; Shimoda et al., 2021) (Figure 2). There were 95 data for panicle formation, 103 for heading, 25 for flowering, and 76 for maturity (Table 2).

Table 2. Details of developmental stages data observed in this study. N_data represents the number of data points between 2004 and 2022. D_min (days), D_max (days), D_mean (days), and D_std (days) respectively denote the minimum, maximum, mean, and standard deviation of the DOYs for each stage.

Fig. 2. Map showing locations and municipalities included in this investigation of wheat developmental stage in Hokkaido, Japan. Blue points show the investigation locations. Red, blue, and green areas are investigated municipalities.

We used agrometeorological grid square data from NARO (https://amu.rd.naro.go.jp) (Ohno et al., 2016; Kominami et al., 2019). Data of daily mean air temperature (°C), photoperiod (hour), and daily snow depth (cm) were used for this study. The values of meteorological data in a 1 × 1 km meshed grid, including the latitude and longitude of the respective sites, were used. The average values of meteorological data in the analysis year for HRO Chuo, Tokachi, and Kitami stations, for NARO/HARC Sapporo and Memuro stations, and for Sorachi, Kamikawa, and Tokachi areas are shown in Table 3. The average values of daily mean air temperature and photoperiod in each developmental period, number of days on which the snow depth is greater than 0 cm from sowing to maturity, and daily snow depth when the snow depth is greater than 0 cm are shown in the table.

Table 3. Average values of meteorological data in the analysis year for HRO Chuo, Tokachi, and Kitami stations, for NARO/HARC Sapporo and Memuro stations, and for Sorachi, Kamikawa, and Tokachi areas. Average values of daily mean air temperature (T) (°C) and photoperiod (P) (hour) in each developmental period, number of days on which snow depth is greater than 0 cm from sowing to maturity (NS) (days), and daily snow depth (SD) (cm) when snow depth is greater than 0 cm are shown.

3.2 Training and testing

3.2.1 Training and testing datasets

For testing, we performed a cascade prediction that continuously predicts panicle formation, heading, flowering, and maturity. The prediction uses the following six stages of data: 1) sowing, 2) snow melting, 3) panicle formation, 4) heading, 5) flowering, and 6) maturity data. If 1) to 6) data are available for all survey points and all years, then one can predict the panicle formation, heading, flowering, and maturity in the cascade prediction. However, the stage type of data that is available differs depending on the survey point and year. Furthermore, the difficulty arose of the small number of data for flowering. Therefore, we first set up the testing dataset for flowering. We let $D_{\mathit{\text{stage}}}^t$ be the true (measured) DOY for one stage. Superscript t denotes the true value. First, five data for which $D_{\mathit{\text{stage}}}^t$ of all six stages ( $D_{\mathit{\text{sowing}}}^t$ , $D_{\mathrm{\text{snow melting}}}^t$ , $D_{\mathrm{\text{panicle formation}}}^t$ , $D_{\mathrm{\text{heading}}}^t$ , $D_{\mathrm{\text{flowering}}}^t$ , and $D_{\mathrm{\text{maturity}}}^t$ ) were available were selected for the testing dataset that can evaluate the prediction accuracy of flowering. Next, 15 data for which five stages other than flowering were all available were added to the five previously selected data. Consequently, a total of 20 data were set as the testing dataset for panicle formation, heading, and maturity.

From the remaining data, we chose training data for predicting each stage. We selected data for which $D_{\mathit{\text{sowing}}}^t$ , $D_{\mathrm{\text{snow melting}}}^t$ , and $D_{\mathrm{\text{panicle formation}}}^t$ were available for panicle formation prediction. In addition, $D_{\mathrm{\text{snow melting}}}^t$ , $D_{\mathrm{\text{panicle formation}}}^t$ , and $D_{\mathrm{\text{heading}}}^t$ were available for heading prediction. We used $D_{\mathrm{\text{panicle formation}}}^t$ , $D_{\mathrm{\text{heading}}}^t$ , and $D_{\mathit{\text{flowering}}}^t$ for flowering prediction. In addition, we used $D_{\mathrm{\text{panicle formation}}}^t$ , $D_{\mathrm{\text{heading}}}^t$ , and $D_{\mathit{\text{maturity}}}^t$ for maturity prediction. The numbers of training data for panicle formation, heading, flowering, and maturity were, respectively, 75, 83, 20, and 56.

The testing data prepared using the procedure described above had little yearly bias. However, for the five data for flowering, large fractions of the data were obtained in Sapporo. Furthermore, regarding the 20 data for panicle formation, heading, and maturity, large fractions of the data were acquired in Tokachi area. Therefore, a future challenge is to re-evaluate the model’s generality through data accumulation and retesting.

3.2.2 Method for training, validation, and testing of the proposed model

First, we set hyperparameter d to some values between d₁ and d₂ that respectively signify the minimum and maximum values of the range to search for the estimated d. For set d, we performed five-fold cross-validation. For five-fold cross-validation, data other than the test data were divided into five groups. One group was used as validation data; the other four groups were used as training data (Figure 3(a)). For each of the five combinations of training and validation data, parameters a₀, a₁, a₂, a₃, and a₄ were optimized by linear regression from the training data to build a model. The accuracy was evaluated by calculating the square error with the validation data. Because training and validation were performed on five combinations, five squared errors were obtained. After these errors were averaged, the accuracy at set d was evaluated using the mean squared error. Five-fold cross-validation was performed by setting d to each value between d₁ and d₂. The accuracy is obtained for each value between d₁ and d₂. Then the value of d which showed the best accuracy was adopted. The steps described above can be expressed mathematically as

Fig. 3.Training, validation, and testing datasets for the proposed model (a) and training and testing datasets for the DVR model (b). Five-fold cross-validation was applied to derive the hyperparameters and parameters for the proposed model.

where $\hat{d}$ stands for the estimated d. Regarding other variables, m_v stands for the number of validation data, N denotes the number of days from the previous stage to the target stage, which is equal to the target valuable D₃-D₂ (N=D₃-D₂; D₃-D₂=D_{panicle formation}-D_{snow melting}, D_heading-D_{panicle formation}, D_flowering-D_heading, D_maturity-D_heading, were given, respectively, for the prediction of panicle formation, heading, flowering and maturity), ${N^t}^{(i,j)}$ expresses the true value of jth N at the ith validation, and ${N^p}^{(i,j)}$ signifies the jth predicted value at ith validation. The values of d₁ and d₂ were, respectively, 0 and 50 for panicle formation, heading, and flowering, and were 30 and 80 for maturity. Here $\hat{d}$ is the optimal solution of d that had the highest average validation accuracy in the brute force calculation in the range of d₁-d₂ (Equation 13). Applying $\hat{d}$ to Equation 3, we determined parameters a₀, a₁, a₂, a₃, and a₄ using the entire training and validation dataset. The estimated parameters $\widehat{a_0},\widehat{a_1},\widehat{a_2},\widehat{a_3}$ , and $\widehat{a_4}$ were

where m_tv stands for the total number of training and validation data, ${N^t}^{\left(i\right)}$ represents the true value of ith N, and ${N^p}^{\left(i\right)}$ expresses the ith predicted value. We used R software ver. 4.2.1 and the Fitting Linear Models, “lm” function of the R to solve the optimization.

3.2.3 Method for training and testing of the DVR model

The training and testing datasets were selected as shown in Figure 3(b). The parameters C, C_t, C_p, C_ver, T_b, and P_b were determined using the training dataset. In addition, C_t, C_p, C_ver, T_b, and P_b were estimated as

where m_t stands for the number of training data. This optimization was solved using R software ver. 4.2.1 and the Fitting Linear Models, “optim” function of the R using the method developed by Nelder and Mead (1965). The initial values of C_t, C_p, C_ver, T_b, and P_b were, respectively, 0.50, 0.50, 0.50, 7.00 °C, and 10.00 h. They were determined in the ranges of 0.10 ≤ C_t ≤ 1.00, 0.10 ≤ C_p ≤ 1.00, 0.10 ≤ C_ver ≤ 1.00, 0.00 ≤ T_b ≤ 15.00, and 8.00 ≤ P_b ≤ 11.00. In addition, C was estimated as

3.3 Prediction accuracy evaluation

For the proposed model, we evaluated the prediction accuracy of each stage of DOY when each prediction was made successively. For this analysis, the predicted DOY of the previous stages was used. A cascade prediction of stages was executed as follows. Let D_stage (days) be the DOY of the target stage. It can be represented as

where D_{previous stage} (days) denotes the DOY of the previous stage and N_stage (days) stands for the number of days from the previous stage to the target stage. For example, D_{panicle formation} is predicted using the following equation: D_{panicle formation}=D_{snow melting}+N_{panicle formation}. For heading, flowering, and maturity, the values for D_stage are also calculated respectively with D_{panicle formation} and N_heading, D_heading and N_flowering, and D_heading and N_maturity.

For the DVR model, panicle formation, heading, flowering, and maturity were predicted from emergence dates predicted from sowing dates. Therefore, D_stage can be represented as

where M_emergence (days) denotes the number of days from the sowing date to the emergence date and M_stage (days) stands for the number of days from the emergence date to the target stage.

To evaluate the prediction accuracies of these models, we calculated the root mean square error (RMSE). We also investigated the mean accuracy by the mean absolute error (MAE), which is less sensitive to outliers than RMSE. The values of RMSE and MAE were calculated as

and

where ${D_{\mathit{\text{stage}}}^t}^{\left(i\right)}$ is the true value of ith $D_{\mathit{\text{stage}}}$ , and ${D_{\mathit{\text{stage}}}^p}^{\left(i\right)}$ is its predicted value. For the proposed model, the $D_{\mathit{\text{stage}}}^t$ is represented by

where $N_{\mathit{\text{stage}}}^t$ is the true value of $N_{\mathit{\text{stage}}}$ . For $D_{\mathit{\text{stage}}}^p$ , it is

for panicle formation and

for heading, flowering, and maturity. In Equation 21 and 23, $N_{\mathit{\text{stage}}}^p$ is the predicted value of $N_{\mathit{\text{stage}}}$ . Equation 22 used $D_{\mathit{\text{previous stage}}}^t$ because $D_{\mathit{\text{previous stage}}}$ for panicle formation, $D_{\mathrm{\text{snow melting}}}$ , was not predicted. For the DVR model,

and

where $M_{\mathit{\text{emergence}}}^t$ and $M_{\mathit{\text{emergence}}}^p$ respectively represents the true value and the predicted value of $M_{\mathit{\text{emergence}}}$ . In addition, $M_{\mathit{\text{stage}}}^t$ and $M_{\mathit{\text{stage}}}^p$ are the true value and the predicted value of $M_{\mathit{\text{stage}}}$ , respectively.

4. Results

4.1 Correlation between objective and explanatory variables in the proposed model

Figure 4 presents test results of correlation obtained between the target variable (D₃-D₂) and explanatory variables (D₁, D₂, T₁, and T₂) in the prediction of the respective stages where (D₁, D₂, D₃)=(D_sowing, D_{snow melting}, D_{panicle formation}), (D_{snow melting}, D_{panicle formation}, D_heading), (D_{panicle formation}, D_heading, D_flowering), (D_{panicle formation}, D_heading, D_maturity) were given, respectively, for the prediction of panicle formation, heading, flowering and maturity. Symbol T₁ represents the average effective temperature during the period from D₁ to D₂. Symbol T₂ represents the average effective temperature during the period from D₂ to D₂+d, where d is the hyperparameter. Strongly or moderately negative correlation was found especially between D₃-D₂ and D₂ and between D₃-D₂ and T₂ in all developmental stage predictions. By contrast, correlation between D₃-D₂ and D₁ and between D₃-D₂ and T₁ was not strong, except for D₃-D₂ and D₁ in maturity prediction. Strong correlation found for D₂ and T₂ and weak correlation found for D₁ and T₁ suggest that wheat development can be affected more strongly by the immediate developmental state and by the environment.

4.2 Prediction accuracy

Table 4 presents the RMSEs and MAEs of the cascade prediction for all stages. For the DVR model, the RMSEs were, respectively, 3.6, 1.9, 2.5, and 3.1 days for panicle formation, heading, flowering, and maturity. By contrast, for the proposed model, they were 2.8, 1.9, 1.9, and 2.6 days, respectively. The MAEs for the DVR model were, respectively, 3.1, 1.4, 1.6, and 2.7 days for the stages. For the proposed model, they were 2.3, 1.3, 1.4, and 2.1 days, respectively. Therefore, the proposed model demonstrated improved prediction accuracy for all stages. The histograms of prediction error between the predicted and actual values (Figure 5) also show that the prediction errors of the proposed model were more likely to be zero than those of the DVR model.

Fig. 4.Scatter plots of the target variable (D₃ – D₂ (days)) and explanatory variables (D₁ (days), D₂ (days), T₁ (°C) and T₂ (°C)) for panicle formation (a), heading (b), flowering (c), and maturity (d). Correlation coefficients R and P values are also shown for each chart. Let D_stage be the DOY of the target stage. (D₁, D₂, D₃)=(D_sowing, D_{snow melting}, D_{panicle formation}), (D_{snow melting}, D_{panicle formation}, D_heading), (D_{panicle formation}, D_heading, D_flowering), (D_{panicle formation}, D_heading, D_maturity) were given, respectively, for the prediction of panicle formation, heading, flowering, and maturity. Symbol T₁ represents the average effective temperature during the period of D₁ – D₂. Symbol T₂ represents the average effective temperature during D₂ – D₂+d, where d (days) is a hyperparameter.

Table 4. Root mean square errors (RMSEs) and mean absolute errors (MAEs) for cascade prediction. RMSE and MAE are given in days.

Table 5 presents the correlation coefficients among D₁, D₂, T₁, and T₂ in the proposed model. Some variables were correlated with correlation coefficients greater than 0.4. Therefore, we investigated whether multicollinearity affected the prediction accuracy. To justify the use of D₁, D₂, T₁, and T₂ in the proposed model, we investigated how the prediction accuracy changed when fewer variables were used. The RMSEs and MAEs of the predictions using only DOYs (D₁ and D₂), only temperatures (T₁ and T₂), only stage 1 variables (D₁ and T₁), and only stage 2 variables (D₂ and T₂) were investigated. In addition, the RMSEs and MAEs of the predictions using three variables (D₁, D₂, and T₁; D₁, D₂, and T₂; D₁, T₁, and T₂; D₂, T₁, and T₂) were evaluated. Table 6 presents the prediction accuracies of each stage using the selected parameters and all parameters D₁, D₂, T₁, and T₂ (the results obtained using all the parameters are the same as those in Table 4). Comparison of the results obtained using fewer variables and those obtained using all variables shows that the RMSEs and MAEs using T₁ and T₂, and D₂ and T₂ were smaller for flowering, but almost all developmental stages predicted with D₁, D₂, T₁, and T₂ were more accurate than those with fewer variables.

Fig. 5.Histograms of prediction error between predicted and actual values. Upper panels (a) show results for the proposed model. Lower panels (b) show results for the DVR model.

Table 5. Correlation coefficients between explanatory variables D₁ (days), D₂ (days), T₁ (°C), and T₂ (°C) in the proposed model. Details of D₁, D₂, T₁, and T₂ are presented in the footnotes.

Table 6. RMSEs and MAEs of predictions using explanatory variables of only DOYs (D₁ (days) and D₂ (days)), only temperatures (T₁ (°C) and T₂ (°C)), only stage 1 variables (D₁ and T₁), and only stage 2 variables (D₂ and T₂) in the proposed model. In addition, those using three variables (D₁, D₂, and T₁; D₁, D₂, and T₂; D₁, T₁, and T₂; D₂, T₁, and T₂) were also evaluated. Results obtained using all variables (D₁, D₂, T₁, and T₂) (the same as results in Table 4) are also shown. Details of D₁, D₂, T₁, and T₂ are presented in the footnotes. The units of RMSE and ME are days.

4.3 Model parameters

Table 7 presents the parameters and hyperparameters of the proposed model. It is apparent that the magnitudes of the absolute values of parameters a₁, a₂, a₃, and a₄ differ depending on the stage. However, a₂ and a₄ were negative throughout all stages. Therefore, a later date of the previous stage (D₂) and a higher average effective temperature between the previous stage to the target stage (T₂ from D₂ to D₂+d) are associated with fewer days from the previous stage to the target stage (D₃-D₂).

Table 7. Parameters and hyperparameters of the proposed model. Symbols a₀, a₁, a₂, a₃, and a₄ (a₀, a₁, and a₂ have no units; the units of a₃ and a₄ are days °C^-1) are parameters; d (days) is a hyperparameter.

Table 8 presents the parameters used for the DVR model. Kawakita et al. (2020a) reported that T_b=11.75±1.72 and P_b=8.40±0.28 for flowering prediction of ‘Kitahonami’. In that paper, (T_b, P_b)=(6.48, 11.00), (7.96, 10.99), (7.05, 10.99), (8.77, 11.00) were given, respectively, for the prediction of panicle formation, heading, flowering and maturity. The air temperature in Hokkaido, the target area of this study, is lower than that in either Honshu or Kyusyu, the target areas of Kawakita et al. (2020a). Moreover, T_b for ‘Kitahonami’ in Hokkaido was smaller; P_b was larger.

Table 8. Parameters for the DVR model. T_b (°C) and P_b (hour) respectively represent the base air temperature and photoperiod. C, C_t, C_p, and C_ver are model coefficients.

5. Discussion

5.1 Difference between the proposed model and conventional models

Which explanatory variables are selected is a crucially important consideration for improving the model accuracy. Most conventional models emphasize the use of air temperature and photoperiod as explanatory variables and apply the nonlinear functions to express the relation between development and air temperature, the relation between development and photoperiod, and the relation between development and vernalization. The developmental stages are then predicted by accumulating the developmental degree as calculated using the nonlinear functions (Ritchie, 1991; Hunt and Pararajasingham, 1995; Wang and Engel, 1998; Yan and Wallance, 1998; Brooking and Jamieson 2002; Streck et al., 2003; Xue et al., 2004; Zheng et al., 2014; Kawakita et al., 2019; 2020a, b; Zhao et al., 2021). By contrast, our proposed model uses regression equations incorporating the DOYs of multiple preceding stages (D₁ and D₂) and the air temperature (T₁ and T₂) (Equation 1) to calculate the number of days from the preceding stage to the target stage. DOYs have not been used in conventional models, but all developmental stages predicted using D₁, D₂, T₁, and T₂ were more accurate than those incorporating fewer variables (Table 6). Therefore, the usage of D₁, D₂, T₁, and T₂ in the proposed model was regarded as appropriate for predicting the timing of developmental stages (panicle formation, heading, flowering, and maturity in this study).

We investigated which variables among D₁, D₂, T₁, and T₂ are important for accurate prediction of developmental stages. As one assessment method, we compared the results obtained using two explanatory variables in four combinations. Results show that the prediction accuracy in all stage predictions was best when using D₂ and T₂ (Table 6). The result indicates that the DOY of the stage immediately preceding the prediction target stage and the subsequent air temperature from the preceding stage are important for predicting development stages. Furthermore, our results showed that the prediction accuracy decreased most when T₂ was removed from the used explanatory variables (Table 6), which indicates the air temperature after snow melting as the most important variable for panicle formation prediction. The large negative correlation between the number of days from sowing to panicle formation and mean air temperature during the period was reported by Eguchi and Shimada (2000) in western Japan. Therefore, the importance of air temperature before panicle formation stage shares regardless of whether snow coverage occurs or not.

5.2 Principles for applying the proposed model to Hokkaido

The proposed model handles the air temperature before the snow melting date differently than the conventional developmental prediction model for Japanese wheat. Previous studies by Maruyama et al. (2010), Nakazono et al. (2014), and Kawakita et al. (2020a) do not incorporate consideration of snow coverage into the developmental prediction models. Shimazaki et al. (2022) took a simple approach by not accumulating air temperatures when the snow depth exceeds a certain threshold in a preliminary study. In our proposed model, the temperatures surrounding the wheat under snow are estimated from the air temperature and the snow depth using a very simple method. In this calculation of temperature during the snow coverage period, the actual snow melting dates that include both the natural melting dates and the artificial melting dates because of the application of snow melting agents are also considered.

To examine the importance of using the actual snow melting date, we evaluated the prediction accuracy using the estimated snow melting date based on the simple definition of snow melting shown below, instead of using the actual measured snow melt date. The definition used here is the following: the day on which the snow depth reaches 0 cm after March 1 in the agrometeorological grid square data from NARO is defined as the snow melting date. This definition expresses snow melting naturally because of rising air temperature and does not consider the artificial snow melting which occurs earlier than the natural snow melting because of the application of snow melting agents. Using the snow melting dates estimated using the definition above, the RMSE for panicle formation, heading, flowering, and maturity were 3.1, 2.0, 1.3, and 2.9 days, respectively. The MAEs were 2.7, 1.6, 0.6, and 2.2 days, respectively. This result demonstrates that the prediction accuracy was worse when using the estimated snow melting dates than when using the actual measured dates. Therefore, the actual snow melting date was found to be a useful variable. However, as described in Subsection 2.1.2, we adopted the assumption when no actual measurement of the snow melting date was available. The estimated errors of the snow melting dates cumulatively propagate in the cascade prediction of panicle formation, heading, flowering, and maturity. Therefore, input of the actual snow melting dates is important for predicting these stages more accurately.

For this study, the equations reported by Ritchie (1991) were used to estimate the temperature surrounding the wheat under snow from air temperature and snow depth. The equations were originally proposed for a developmental prediction model based on thermal time (Ritchie, 1991; Zheng et al., 2014; Zhao et al., 2021). Reportedly, the thermal time is calculable with sufficient accuracy based on air temperature, and the equations reported by Ritchie (1991) depending on air temperature and snow depth are used only for verification of frost damage in the major global crop model WOFOST. However, the problem persists of not strictly considering the heat balance between the atmosphere, snow, and soil in the Ritchie’s equations. In the future, it is necessary to verify what kind of temperature estimation equation is suitable for highly accurate prediction of developmental stage of wheat in Hokkaido, where there is a snow coverage period.

In addition to the snow melting dates and the heat balance between the atmosphere, snow, and soil, the most important point of consideration is vernalization (Ritchie, 1991; Hunt and Pararajasingham, 1995; Wang and Engel, 1998; Yan and Wallance, 1998; Brooking and Jamieson 2002; Streck et al., 2003; Xue et al., 2004; Zheng et al., 2014; Kawakita et al., 2020a, b; Zhao et al., 2021). According to the calculation of f_ver (T_i) in Equation 10, vernalization requirement is achieved around DOY = 350 (around December 17) in the year of sowing in Hokkaido. The date of DOY = 350 falls in the period from the sowing date to the snow melting date which is considered in the T₁ calculation for the panicle formation prediction. However, the vernalization function is not applied to the T₁ calculation in the current proposed model. Therefore, it is possible that the prediction accuracy of panicle formation decreased because the change in the developmental rate before and after vernalization was not considered.

The relation between snowmelt and spring vegetation phenology has been investigated in Arctic, Alaska, and North America (Potter, 2020; Kelsey et al., 2021; Zheng et al., 2020; Zheng et al., 2022; Wu et al., 2023). They reported that the timing of snowmelt affects the start of the growing season. In these earlier studies, the timing of snow melting was estimated using remote sensing data. Similarly, in this study, the usage of records along with diagnostic results based on satellite imagery related to natural and artificial snow melting dates is expected to improve prediction accuracy in the future. In addition, inclusion of the vernalization function in the proposed model is expected to improve the prediction accuracy of panicle formation.

6. Conclusion

We constructed a model to predict the developmental stages of winter wheat considering snow cover periods in Hokkaido. In the model, DOY data of multiple developmental stages preceding the prediction target stage and air temperature were used as explanatory variables. To predict panicle formation, heading, flowering, and maturity of winter wheat ‘Kitahonami’ in Hokkaido, the proposed model yielded better accuracy than the conventional DVR model. Therefore, the use of multiple preceding stages and air temperature has been demonstrated as useful for predicting various developmental stages. The DOY of snow melting date was found to be a useful explanatory variable. Comparisons of calculated values revealed that the prediction accuracy using actual measured snow melting dates, including artificial snow melting dates because of the application of snow melting agents, was better than that using snow melting dates that were estimated based on meteorological data. In future studies, further improvement of prediction accuracy can be expected by the comprehensive usage of records and diagnostic results based on satellite images related to natural and artificial snow melting dates. It is also expected that prediction accuracy will be improved with the proposed model by inclusion of a formula expressing the effects of vernalization, although that feature was not incorporated into the current proposed model.

Acknowledgments

Figures 1 and 3, and Tables 1 and 2 were referred from “Cascade prediction of the day of year of winter wheat developmental stages in Hokkaido, Japan” by Iwasaki et al., 2023 with permission of the American Society of Agricultural and Biological Engineers (ASABE). This publication was supported by the Japan Society for the Promotion of Science (23HP2004).

References

• Aase JK, Siddoway FH, 1979: Crown-depth soil temperatures and winter protection for winter wheat survival. Soil Science Society of America Journal 43, 1229-1233.
• Brooking IR, Jamieson PD, 2002: Temperature and photoperiod response of vernalization in near-isogenic lines of wheat. Field Crops Research 79, 21-38.
• Eguchi H, Shimada S, 2000: Phenological development of wheat cultivars in the field and screening for early maturing traits. Japanese Journal of Crop Science 69, 49-53.
• Fukuda M, 1982: The Estimation method of frost penetration depth under snow packs condition. Research reports of grant-in-aid for Special Project of Hokkaido University in Study on Snow Distribution and Characteristics and its Effect on Frozen Soil Depth, Vegetation and Insect Ecology in Hokkaido, 55-74.
• Hirota T, Iwata Y, Hayashi M, et al., 2006: Decreasing soil-frost depth and its relation to climate change in Tokachi, Hokkaido, Japan. Journal of the Meteorological Society of Japan. Ser. II 84, 821-833.
• Hokkaido Prefecture 2023: Wheat, beans, and minor grains handbook.
• Hokkaido Research Organization, 2014: High-quality and stable cultivation method for winter wheat ‘Yumechikara’. Available at: https://www.hro.or.jp/upload/48940/2015104.pdf.
• Hunt LA, Pararajasingham S, 1995: CROPSIM-WHEAT: A model describing the growth and development of wheat. Canadian Journal of Plant Science 75, 619-632.
• Iwasaki C, Shimoda S, Sugiura R, et al., 2023: Cascade Prediction of the Day of Year of Winter Wheat Developmental Stages in Hokkaido, Japan. 2023 ASABE Annual International Meeting, Paper, 2300403.
• JA Bibai, JA Minenobu, 2015: Bibai City, Sorachi Agricultural Extension Center in Hokkaido. Autumn Wheat Project, 1-16.
• Kawakita S, Inaba S, Takahashi T, et al., 2019: Evaluation of non-linear wheat development models and optimization methods for their parameter determination. Journal of Agricultural Meteorology 75, 120-128.
• Kawakita S, Ishikawa N, Takahashi H, et al., 2020a: Winter wheat phenological development model with a vernalization function using sigmoidal and exponential functions. Journal of Agricultural Meteorology 76, 81-88.
• Kawakita S, Takahashi H, Moriya K, 2020b: Prediction and parameter uncertainty for winter wheat phenology models depend on model and parameterization method differences. Agricultural and Forest Meteorology 290, 107998.
• Kelsey KC, Pedersen SH, Leffler AJ, et al., 2021: Winter snow and spring temperature have differential effects on vegetation phenology and productivity across Arctic plant communities. Global Change Biology 27, 1572-1586.
• Kominami Y, Sasaki K, Ohno H, 2019: User’s Manual for the Agro-Meteorological Grid Square Data, In NARO Ver. 4 (in Japanese). NARO, 67.
• Maruyama A, Kurose Y, Ohba K, 2010: Modeling of phenological development in winter wheat to estimate the timing of heading and maturity based on daily mean air temperature and photoperiod. Journal of Agricultural Meteorology 66, 41-50.
• Nakazono K, Ohno H, Yoshida H, et al., 2014: Modeling phenological development in wheat. Japanese Journal of Crop Science 83, 249-259.
• NARO, 2021: A program for predicting wheat growth considering snowmelt. Available at: https://www.naro.go.jp/collab/program/laboratory/harc/141232.html.
• Nelder JA, Mead R, 1965: A simplex method for function minimization. The Computer Journal 7, 308-313.
• Ohno H, Sasaki K, Ohara G, et al., 2016: Development of grid square air temperature and precipitation data compiled from observed, forecasted, and climatic normal data. Climate in Biosphere 16, 71-79.
• Potter C, 2020: Snowmelt timing impacts on growing season phenology in the northern range of Yellowstone National Park estimated from MODIS satellite data. Landscape Ecology 35, 373-388.
• Ritchie JT, 1991: Modeling plant and soil systems, Agronomy Monographs. Madison, 1-545.
• Shimazaki Y, Seki M, Shimoda S, 2022: Examination of phenological developmental prediction based on the accumulated temperature of barley cultivated in snowy regions. The Hokuriku Crop Science 57, 23-28.
• Shimoda S, Hamasaki T, 2021: Potential benefits of promoting snowmelt by artificial snow blacking on the growth of winter wheat and their dependence upon regional climate. International Journal of Biometeorology 65, 223-233.
• Shimoda S, Onodera M, Okumura O, et al., 2021: Effects of snow compaction ‘yuki-fumi’ on soil frost depth and volunteer potato control in potato-wheat rotation system in Hokkaido. Plant Production Science 24, 186-197.
• Shimoda S, Terasawa Y, Nishio Z, 2022: Improving wheat productivity reveals an emerging yield gap associated with short-term change in atmospheric humidity. Agricultural and Forest Meteorology 312, 108710.
• Streck NA, Weiss A, Xue Q, et al., 2003: Improving predictions of developmental stages in winter wheat: A modified Wang and Engel model. Agricultural and Forest Meteorology 115, 139-150.
• Tsuchiya F, 1985: Effect of soil frost and frost penetration on volcanic ash soil in agricultural field in Tokachi, eastern Hokkaido, Japan. Ph.D. thesis, Hokkaido Univ., Japan, pp. 234.
• Wang E, Engel T, 1998: Simulation of phenological development of wheat crops. Agricultural Systems 58, 1-24.
• Wu Y, Xiao P, Zhang X, Liu, et al., 2023: Effects of snow cover on spring vegetation phenology vary with temperature gradient across the Pan-Arctic. Journal of Geophysical Research: Biogeosciences 128, e2022-JG007183.
• Xue Q, Weiss A, Baenziger P, 2004: Predicting phenological development in winter wheat. Climate Research 25, 243-252.
• Yan W, Wallance DH, 1998: Simulation and prediction of plant phenology for five crops based on photoperiod × temperature interaction. Annals of Botany 81, 705-716.
• Zhao HD, Sassenrath GF, Zambreski ZT, et al., 2021: Predicting winter wheat heading date: A simple model and its validation in Kansas. Journal of Applied Meteorology and Climatology 60, 1685-1696.
• Zheng B, Chenu K, Doherty A, Chapman S, 2014: The APSIM-wheat module (7.5 R3008). Agricultural Production Systems Simulator (APSIM) Initiative.
• Zheng J, Xu X, Jia G, 2020: Understanding the spring phenology of Arctic tundra using multiple satellite data products and ground observations. Science China Earth Science 63, 1599-1612.
• Zheng J, Jia G, Xu X, 2022: Earlier snowmelt predominates advanced spring vegetation greenup in Alaska. Agricultural and Forest Meteorology 315, 108828.