Estimation of Blooming Start with the Adaptation of the Unified Model for Three Apricot Cultivars (Prunus armeniaca L.) Based on Long-Term Observations in Hungary (1994–2020)

Abstract

The aim of our research was to adapt Chuine’s unified model to estimate the beginning of blooming of three apricot cultivars (‘Ceglédi bíborkajszi’, ‘Gönci magyar kajszi’, and ‘Rózsakajszi C.1406’) in Hungary in the time period 1994–2020. The unified model is based on the collection of chilling and forcing units. The complexity of the model lies in the high number of parameters necessary to run it. Following the work of other researchers, we reduced the number of relevant model parameters (MP) to six. In order to estimate the six MPs, we used a simulated annealing optimization method (known for being effective in avoiding getting stuck in local minima). From the results, we determined the local optimum of six MPs, and the global optimum parameter vector for three apricot cultivars. With these global optimum parameter vectors, the beginning of blooming could be estimated with a root-mean-square error (RMSE) of less than 2.5 days, using the knowledge of the daily mean temperature in the time period 1994–2020.

1. Introduction

The bud burst and the blooming are the most extensively studied phenophases of temperate trees (for more details see [1]). Such studies found that the dominant environmental factor affecting the growth speed of trees is air temperature [2]. In the last few decades, several temperature-based, mechanistic phenology models were produced for simulating the vegetative (e.g., bud burst) or reproductive (e.g., blooming) phenology of temperate and boreal trees (for more details see [1]). These models are based on the response of trees to the chilling requirement (the accumulation of which breaks endodormancy), and also the response of bud growth to forcing heat requirement (which forces the growth during spring after endodormancy break). Phenological models are important tools for planning agricultural practices [3], particularly for the prevention of frost damage in the short term [4,5] and for projecting the impact of climate change in the long term [6,7,8,9].

Dormancy is an important phase in the life of temperate zone fruit trees. Regular dormancy is essential for surviving unfavorable temperatures in winter [10]. Sufficient chilling units are needed for the breaking of endodormancy of apricot (Prunus armeniaca L.) buds. According to Andreini et al. [11], early cultivars are controlled by minimum temperatures, but intermediate and late cultivars are controlled by mean temperatures.

Blooming time is an important phenological trait of apricot cultivars. It is determined genetically, but modified by environmental factors, mainly by several elements of climate, which is why, for the aspect of blooming time, there are significant differences between cultivars, production sites, and years. Knowledge of blooming time of apricot genotypes is also important for basic research, breeding, and good cultivation practice. It is an important element of the exact description of cultivars [12]. Blooming time is a crucial phenological trait for breeding [13,14]. Apricot is an early blooming species in general, and is endangered by frost damage [15]. Extremely early blooming genotypes frequently suffer under low temperatures [16]. There are self-incompatible cultivars in apricot orchards that need pollinizers cultivars with the same blooming time [17]. Therefore, the knowledge of blooming time of cultivars is essential from a practical aspect as well.

Hungary is situated in the northern zone of apricot production. In Hungary, the first flowers of apricot trees open between the end of February and second half of April, depending on the genotype, the orchard site, and the year [17,18,19,20,21]. In southern countries, and under milder climate conditions, the blooming time of apricot cultivars is earlier, sometimes starting at the end of January [22,23,24].

Based on the long-term observations in different parts of the world, the blooming time of apricot cultivars changed during the last few decades, similar to other temperate zone deciduous tree species. Most studies show advancing spring phenological phases and gradually earlier blooming time because of global warming [21,25,26,27,28,29,30,31,32,33,34]. Only a small number of articles, on a few genotypes and locations, show the opposite, i.e., delayed spring phenological processes and later blooming in response to warming [35,36,37]. The hypothesis of Guo et al. [38] is the following: the spring phenological processes of temperate zone deciduous plants in locations with cold winters should be determined primarily by temperatures during ecodormancy, when the accumulation of forcing units is occurring, whereas at sites with warm winters, temperatures during the endodormancy (chilling requirement) should play a more important role in the regulation of spring phenological processes. Insufficient chilling may cause not just delayed, but irregular, blooming, as well as a lack of crop [14,39,40]. Temperature has a significant impact on the regulation of phenological processes just before blooming and during blooming time [41,42].

In our study, we focused on the beginning of blooming of three apricot cultivars in Hungary, based on daily mean temperatures. Phenology models are very sensitive to local climatic conditions and plant species. Therefore, they must be adapted to local conditions [43]. We used the unified model of Chuine [43]. This model is based on the accumulation of chilling and forcing units. The unified model provides a statistical framework for standard model identification, comparison, and simplification [43,44]. With the appropriate parameter selection, the unified model approximates the sequential model [2,45,46,47], the parallel model [2,46,47,48], and the deepening rest model [49], as well as the four phases model [2,50]. During the dormant period, the biochemical and biophysical processes are not precisely known, therefore, parameter estimation is difficult. There are two methods used to estimate parameters: the response of shoots to temperature is studied under controlled conditions, or statistical models are run on long-term data sets [1]. For the estimation of the unified model parameters, optimization methods such as the simulated annealing method (SAM) are used [1]. The SAM handles the local extrema better than traditional methods, by avoiding getting stuck in them [51,52,53].

2. Materials and Methods

2.1. Meteorological and Phenological Data

We used the blooming dates of three self-fertilizing apricot cultivars (Prunus armeniaca L.; ‘Ceglédi bíborkajszi’ (cb), ‘Gönci magyar kajszi’ (gm), and ‘Rózsakajszi C.1406’ (ro)) recorded by the research team of the Department of Pomology of the predecessor institutes the Hungarian University of Agriculture and Life Sciences in the experimental farm in the time period 1994–2020. The beginning of blooming (blooming start) was the day of the phenological process, when at least 5% of flowers were opened on the trees of the observed cultivars. In addition, we used the time when the string stage occurred in at least 50% of the trees [34]. Our collection was in Szigetcsép (47°15′10.6″ N 18°57′26.6″ E) until 2007. Three years before the excision of this plantation, we set up a new experimental plantation in 2004 in Soroksár (47°23′51.8″ N 19°09′03.2″ E), where we continued our studies. The environmental conditions of the two production sites are very similar, because they are very close to each other. The examination of the phenological process of apricot cultivars was carried out in parallel at both plantations in 2004, 2005, and 2006. No significant difference was found between the two production sites in these processes [34]. The beginning and end of the blooming of cultivars occurred on the same days at both sites. Data evaluated in this paper are from Szigetcsép between September 1994 and April 2006, and from Soroksár from September 2006 to April 2020. The training system at both plantations is compact vase with 5 × 3 m spacing with Prunus cerasifera seedling rootstocks. In the orchard, integrated plant protection was applied. In both plantations, six trees per cultivars were available for testing.

For the model building, we used the daily mean temperatures observed in the synoptic station of Marczell György main observatory of the Hungarian Meteorological Service (47°25′49.2″ N 19°06′43.7″ E) between September 1994 and April 2020. These data were considered to be good estimations of the daily mean temperature in the apricot plantations, taking into account the close distance between the observatory and experimental plantations and the low variance of the temperature layer, as well as the flat terrain of the surrounding area [54].

2.2. The Unified Model

The unified model [43] uses nine parameters ($a_c$, $b_c$, $c_c$, $b_f$, $c_f$, $w$, $k$, $C_{crit}$, $t_c$) to estimate the day of bud burst or the beginning of blooming of woody trees. In what follows, during the description of the unified model (UM), we use the definition of UM in the sense that it estimates the beginning of blooming.

The estimated critical amount of accumulated chilling units ($C_{crit}$) necessary to break the endodormancy and the total accumulated chilling units ($C_{tot}$) are defined as follows:

$${{\displaystyle \sum }}_{t_0}^{t_1}\frac{1}{1+e^{a_c{\left(T-c_c\right)}^2+b_c\left(T-c_c\right)}},$$

(1)

$$C_{tot}={{\displaystyle \sum }}_{t_0}^{t_c}\frac{1}{1+e^{a_c{\left(T-c_c\right)}^2+b_c\left(T-c_c\right)}},$$

(2)

where $a_c$ and $b_c$ are estimated weight parameters, $c_c$ is the estimated parameter of the maximum point, $T$ is the daily mean temperature, $t_0$ is a fix date (1st of September), $t_1$ is the beginning of ecodormancy (beginning of the forcing period), and the estimated parameter $t_c$ is the last day of chilling unit accumulation.

By fixing $t_0$ to 1st of September, we begin our chilling unit accumulation calculations on this date [43], but this does not constrain the model accuracy, because under Hungarian environmental conditions, chilling unit accumulation starts much later, mainly around the end of October [9,55].

The function of the chilling unit accumulation (1) or (2) is a bell-shape curve (Figure 1a) with its maximum value taken at the temperature that is optimal for the plant for chilling unit accumulation ($c_c$). This optimal value is usually a positive temperature close to zero.

The estimated critical amount of forcing units ($F_{crit}$) necessary to stimulate the bud burst is calculated as follows:

$$F_{tot}={{\displaystyle \sum }}_{t_1}^{t_2}\frac{1}{1+e^{b_f\left(T-c_f\right)}},$$

(3)

$$F_{crit}=w{e}^{k{C}_{tot}},$$

(4)

where $b_f$ is an estimated weight parameter, $c_f$ is an estimated parameter of the inflexion point of the function, $T$ is the daily mean temperature, $t_1$ is the beginning of ecodormancy (beginning of the forcing period), $t_2$ is the day of the beginning of blooming, 𝑤 > 0 is an estimated weight parameter, and $k<0$ is an estimated parameter of the relation of F_crit and the total amount of chilling units C_tot. The beginning of blooming is set off when F_tot ≥ F_crit. Considering Equation (4), it is obvious that F_crit and C_tot are in a negative relationship [56]. This relationship expresses that a higher amount of chilling units means a lower amount of forcing units is needed to set the beginning of blooming off, which was proven experimentally [57,58,59,60,61].

The function of the forcing unit accumulation (3) is an S-shape curve (Figure 1b), expressing the higher and higher forcing effect of warming temperatures.

We applied two simplifications to Chuine’s unified model, because the high number of parameters can complicate the optimization, often impairing the accuracy of the estimate [62,63,64]. Moreover, there is a high probability that the underlying biological content disappears during optimization [7,65]. The simplifications are as follows:

(1)

According to Caffarra and Eccel [7], ${\mathit{b}}_{\mathit{c}}$ can be set to 0, so we described the c accumulation in the endodormancy with the following equation:

$${{\mathit{C}}^{\prime }}_{*}\left(\mathit{T}\right)={{\displaystyle \sum }}^{}\frac{\mathbf{2}}{\mathbf{1}+{\mathit{e}}^{{\mathit{a}}_{\mathit{c}}{\left(\mathit{T}-{\mathit{c}}_{\mathit{c}}\right)}^{\mathbf{2}}}}.$$

(5)

(2)

It does not strongly constrain the model accuracy if we assume that the chilling unit accumulation ends at the beginning of ecodormancy (${\mathit{t}}_{\mathit{c}}={\mathit{t}}_{\mathbf{1}}$), [43]. We defined this day as the one when the string stage occurs [66,67]. The observed string stage data were available from the data base of the examined apricot cultivars in the time period 1994–2020. Assuming ${\mathit{t}}_{\mathit{c}}={\mathit{t}}_{\mathbf{1}}$, it follows that ${\mathit{C}}_{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}={\mathit{C}}_{\mathit{t}\mathit{o}\mathit{t}}$. In our study, ${\mathit{t}}_{\mathbf{2}}$ is the day of the beginning of blooming.

We had six parameters left to estimate ($a_c$, $c_c$, $b_f$, $c_f$, $w$, $k$). Based on the literature, we chose a particularly large parameter space for the first runs [7,8,9,43,63,64]. The limits of this parameter space are shown in

Table 1. The boundaries of the parameter space and the step length ¹ used for the first estimation.

Parameter	Minimum Value	Maximum Value	Step Length
$a_c$	0	10	0.01
$c_c$	−50	50	0.10
$c_c$	−10	0	0.01
$b_f$	−30	30	0.10
$c_f$	0	200	0.10
$w$	2	9	0.01

¹ These are not fixed step lengths: the $\sigma$ parameter of the Gaussian distribution was used to select the step length.

. Since $k$ is a negative number very close to zero, we performed the optimization for $k_{exp}$ according to the following equation:

$$k=-{10}^{-k_{exp}}.$$

(6)

2.3. Parameter Estimation with the Simulated Annealing Method

We used the simulated annealing method to estimate the parameters [51,52]. During the simulated annealing calculation process, we searched for the optimal six-dimensional parameter vector, which results in the best-fit to the measured blooming data set. The goodness-of-fit was defined as the root-mean-square error (RMSE) between the estimated and observed data. This measure was minimized within the boundaries of the parameter space. To achieve this, we started a random walk process within the chosen allowed parameter space (Table 1). The starting coordinates were randomized according to a uniform distribution defined between the boundaries of the individual parameters.

At the beginning of the random walk, we set a so-called ‘temperature’ parameter $T^{*}$ to its starting value $T_{start}$, defined as the goodness-of-fit at the starting point. In each random walking step, the temperature parameter $T^{*}$ was lowered (‘cooled’) by a factor of $1-T_{step}$, where $T_{step}=0.001$. The ‘cooling’ process was performed from $T_{start}$ to $T_{crit}=0.01·T_{start}$, so the random walk was finished when $T^{*}\le T_{crit}$.

At the beginning of each walk, we randomly chose a parameter to modify (i.e., the axis along which to step with optimization by ‘cooling’), and a step length from a Gaussian distribution with zero mean and $\sigma$ deviation ($\sigma$ was fixed, Table 1). If a step moved out from the interval of allowed parameters, we stepped towards the opposite direction along the axis of the same parameter. Note that in some exceptional cases when this opposite step also moved out from the boundaries (i.e., the step length was too large), we terminated the whole fitting process, and sent an error message.

A step is accepted with a probability 1 if the measure of goodness-of-fit is smaller than it was previously, and it is accepted by a probability of $p=e^{-\frac{Diff}{T^{*}}}$ if the measure is larger than before. Here, $Diff$ is the absolute difference between the old and the new goodness-of-fit measure. One can see that as $T^{*}$ decreases in every step, the same $Diff$ results in a decreasing probability of the new position being accepted. This step acceptance algorithm ensured that the optimization is not stuck in local minima [53]. As a result of one optimization run (a ‘walk’), we obtain the best-fitting vector of six parameters and the corresponding root-mean-square error (RMSE) value. Altogether, we performed 10,000 simulated annealing walks per cultivar to estimate the parameter vectors of the three apricot cultivars.

We used the MATLAB program system, Microsoft Excel (version: 2015), and Microsoft paint.net 4.2.15 software for our calculations and representations of our results.

3. Results

In reporting our results of 10,000 random walks to estimate the parameters of the three cultivars, we use the following notations for the cultivars: ‘Ceglédi bíborkajszi’ as ‘cb’, ‘Gönci magyar kajszi’ as ‘gm’, and ‘Rózsakajszi’ as ‘ro’.

In what follows, we lead the reader through the optimizing procedure in eight steps.

1.

First, we observe that for each apricot cultivar, the forcing parameters $c_f$ and $w$ are strongly related, although not linearly (see the upper panel in Figure 2). Therefore, these two parameters cannot be optimized independently. Using the empirical relationship between $c_f$ and $w$ obtained from the optimization process, we calculate the optimal values of $c_f$ for all fixed values of $w$;

2.

Based on the resulted parameter vectors of all walks, we plot the histogram of the optimal parameter values of $w$ (see the lower panel in Figure 2). We see that the optimal parameter values of $w$ are dense mainly around three or four small to high values; more exactly, around one small, two medium, and one high value for ‘Rózsakajszi C.1406’, while around one small, one medium, and one high value for the other two cultivars. We immediately exclude the high values, because we obtained $C_f$ values in those cases between −10 °C and −30 °C, which are unlikely during the forcing period in Hungary [9,68];

3.

The results of most walks are dense around the small value for each apricot cultivar, and we obtained the lowest RMSE values here, too. So, we fix the parameter value $w$ at the median of the preferred range of ‘small’ optimal parameters $w$: $w_{cb}$ = 15.6, $w_{gm}$ = 17.9, $w_{ro}$ = 19.4;

4.

As a next step, we narrow the parameter space according to the biologically possible parameter values for Hungary (

Table 2. The boundaries of the parameter space and the step length1 for the second estimation.

Parameter	Minimum Value	Maximum Value	Step Length
$a_c$	0.2	1.0	0.001
$c_c$	1	5	0.005
$b_f$	−0.9	−0.1	0.001
$c_f$	6	14	0.010
$k_{exp}$	2	6	0.005

) [7,9,65,69,70]. In the original parameter space, we find several similarly good parameter vectors, that fit statistically very well to the observed blooming dates, but they are biologically impossible.

With the aim of finding biologically realistic optimal parameter vectors [68], we started another 10,000 walks for each cultivar, applying the fixed parameters $w_{cb}$, $w_{gm}$, and $w_{ro}$. As a result, we obtained 10,000 limit vectors of the six parameters that the optimization converged to. In order to visualize this result, we created fifteen histograms (for five parameters of the three cultivars based on 10,000 limit values, not shown) to find the local optimum parameter values. To detect the significantly more frequent limit values, we divided the histograms into one hundred bins, because it is the square root of the amount of the limit parameter values. Thus, in the random case, an average of one hundred data points falls in each bin ($\mu$ = 100). Unsurprisingly, this happened with natural fluctuation following a Poisson distribution. Having such a large number of limit points, the Poisson distribution could be approximated by the Gaussian distribution with standard deviation $\sigma =\sqrt{\mu }$ = 10. Focusing on a histogram, we define the limit point frequency in a bin as significantly high if it is above the value $\mu +3\sigma$ = 130. From the optimized 15 parameters, in seven cases, namely for $c_f^{cb}$, $k_{exp}^{cb}$, $a_c^{gm}$, $c_c^{gm}$, $c_f^{gm}$,$a_c^{ro}$, and $c_f^{ro}$, we obtain one or more significantly frequented bins next to each other. In another six cases, namely for the parameters $b_f^{cb}$, $b_f^{gm}$, $k_{exp}^{gm}$, $c_c^{ro}$, $b_f^{ro}$, and $k_{exp}^{ro}$, there are more non-adjacent significantly frequented bins. In the remaining two cases ($a_c^{cb}$,$c_c^{cb}$), we do not obtain any bin with at least 130 data points, but there are several bins with more than $\mu +2\sigma$ = 120 optimized parameter values. We define the local optima of the parameters as the medians of the significantly frequented bins. In the case of non-adjacent significant bins, we searched for the local optimum of the parameter in the bin with the largest number of optimized values (

Table 3. Characterization of the local optimum parameter vectors of the unified model: the minimum and maximum values of the bins with significantly higher number of optimized limit values (optimum bins), the number of the limits in them, the median and the standard deviation, together with the maximum, the median, and standard deviation of the RMSE. The values belonging to the cultivars ‘Ceglédi bíborkajszi’, ‘Gönci magyar kajszi’, and ‘Rózsakajszi C.1406’ are denoted by the superscripts ‘cb’, ‘gm’, and ‘ro’ respectively.

	${\mathit{a}}_{\mathit{c}}^{\mathit{c}\mathit{b}}$	${\mathit{c}}_{\mathit{c}}^{\mathit{c}\mathit{b}}$	${\mathit{b}}_{\mathit{f}}^{\mathit{c}\mathit{b}}$	${\mathit{c}}_{\mathit{f}}^{\mathit{c}\mathit{b}}$	${\mathit{k}}_{\mathit{e}\mathit{x}\mathit{p}}^{\mathit{c}\mathit{b}}$	${\mathit{a}}_{\mathit{c}}^{\mathit{g}\mathit{m}}$	${\mathit{c}}_{\mathit{c}}^{\mathit{g}\mathit{m}}$	${\mathit{b}}_{\mathit{f}}^{\mathit{g}\mathit{m}}$	${\mathit{c}}_{\mathit{f}}^{\mathit{g}\mathit{m}}$	${\mathit{k}}_{\mathit{e}\mathit{x}\mathit{p}}^{\mathit{g}\mathit{m}}$	${\mathit{a}}_{\mathit{c}}^{\mathit{r}\mathit{o}}$	${\mathit{c}}_{\mathit{c}}^{\mathit{r}\mathit{o}}$	${\mathit{b}}_{\mathit{f}}^{\mathit{r}\mathit{o}}$	${\mathit{c}}_{\mathit{f}}^{\mathit{r}\mathit{o}}$	${\mathit{k}}_{\mathit{e}\mathit{x}\mathit{p}}^{\mathit{r}\mathit{o}}$
Minimum	0.568	1.400	−0.436	7.600	2.080	0.248	3.160	−0.420	7.280	2.040	0.312	2.640	−0.556	7.360	2.040
Maximum	0.576	1.440	−0.428	9.280	2.480	0.256	3.200	−0.380	9.040	2.400	0.320	2.720	−0.484	8.560	2.600
No. of items	128	129	156	7480	1592	143	138	735	7950	1680	132	270	1318	7180	2210
Median	0.572	1.420	−0.432	8.210	2.260	0.252	3.180	−0.400	7.880	2.220	0.316	2.680	−0.520	7.710	2.310
St. deviation	0.002	0.010	0.002	0.420	0.110	0.002	0.010	0.011	0.430	0.100	0.002	0.020	0.021	0.310	0.160
Maximum RMSE	4.40	4.17	3.40	4.60	4.39	5.04	4.66	4.58	7.18	4.68	4.78	5.63	4.62	5.78	5.05
Median RMSE	2.92	2.94	2.84	2.86	2.66	2.82	2.81	2.79	2.81	2.49	2.47	2.38	2.41	2.41	2.15
St. dev. RMSE	0.33	0.33	0.18	0.21	0.31	0.49	0.43	0.27	0.25	0.33	0.57	0.63	0.31	0.31	0.50

).

5.

Then, we searched for the global optimum of the parameter space for each apricot cultivar. We define the global optimum parameter vector as the parameter vector with the lowest root-mean-square error (RMSE) among the grid of values of

Table 4. The values of the global optimum parameter vectors of the unified model and the corresponding RMSE values for ‘Ceglédi bíborkajszi’, ‘Gönci magyar kajszi’, and ‘Rózsakajszi C.1406’.

Parameter	‘Ceglédi bíborkajszi’	‘Gönci magyar kajszi’	‘Rózsakajszi C.1406’
$a_c$	0.949	0.216	0.608
$c_c$	1.50	2.13	2.42
$b_f$	−0.626	−0.443	−0.365
$c_f$	8.30	9.04	8.84
$w$	15.60	17.90	19.40
$k_{exp} \left(k\right)$	2.14 (−0.0072)	2.08 (−0.0083)	2.07 (−0.0086)
$t_c=t_1$	14th of January	22nd of January	30th of January
$t_2$	27th of March	29th of March	1st of April
$C_{crit}=C_{tot}$	12.73	29.78	19.69
$F_{crit}$	14.24	13.99	16.41
$F_{tot}$	14.57	14.31	16.82
RMSE	2.37	2.10	1.49

. It is seen that, in many cases ($a_c^{cb}$, $c_c^{cb}$, $b_f^{cb}$, $a_c^{gm}$, $c_c^{gm}$, $b_f^{gm}$, $a_c^{ro}$, $c_c^{ro}$, $b_f^{ro}$, and $c_f^{ro}$), the global optimal parameter values do not fall in the local optimum bins. This is most surprising for the parameter $c_f^{ro}$, where more than 70% of the walk limits fall in the local optimal bin, but the global optimum parameter value does not;

6.

Using the global optimum parameter vector, we estimated the blooming date for each apricot cultivar with an average error less than 2.5 days (RMSE < 2.5). For comparison, if we take the mean blooming data calculated over all the years as a constant [64], the average error of the estimation (i.e., the error of the base model) is as high as 9.7–10.6 days, depending on cultivars;

7.

Finally, based on the daily average temperature and the global optimal parameter vectors, we calculated the critical amount of chilling and forcing units, and determined the chilling and forcing process for each cultivar in the period 1994–2020 (Figure 1 and Figure 3). We provide the parameter values that are optimized with the simulated annealing method, applying the unified model and the observed string stage and blooming data of years 1994–2000 (Table 4). The temperature that is optimal for the plant for chilling unit accumulation ($c_c$) is 1.50 °C for ‘Ceglédi bíborkajszi’, 2.13 °C for ‘Gönci magyar kajszi’, and 2.42 °C for ‘Rózsakajszi C.1406’ in the period 1994–2020. According to our calculations, the most chilling units ($C_{crit}=$ 29.8 units) are necessary for ‘Gönci magyar kajszi’, and the least chilling units ($C_{crit}=$ 12.7 units) are required by ‘Ceglédi bíborkajszi’ for breaking the endodormancy (Table 4). The inflection point of the forcing unit accumulation (i.e., $c_f$) is between 8.30 and 9.04 °C, depending on the cultivars. This curve has no maximum point, but the forcing unit accumulation is close to the maximum (1 unit) at 12–15 °C (more than 0.9 units) that could be considered as ‘optimal temperature’ for the plant in their preparation for blooming. The average accumulated forcing units for the blooming are between 14.0 and 16.4 units for each cultivar in the period 1994–2020 (Table 4). Surprisingly, the absolute value of parameter $k$ of our results is larger than is reported in the publications of other researchers (i.e., in between −10⁻⁴ and −10⁻⁸) [8,9,43,64]. This may lead to a conclusion that, in the case of Hungarian apricots, the chilling unit accumulation has a relatively larger effect on forcing unit accumulation.

Note that though the change in the observed and estimated blooming start is not proven as significant in the past 26 years, the slopes of their trends are all negative, below −0.2 (

Table 5. The mean, standard deviation (StDev), range, lower and upper 95% confidence limits (LCI, UCI) of the observed and estimated blooming start (BM), and the parameters of the unified model $F_{crit}$, $C_{crit}=C_{tot}$, and $F_{tot}$, together with their slope in the measured 26 years (1994–2020) with their corresponding significance levels ($p$)of the varieties ‘Ceglédi bíborkajszi’, ‘Gönci magyar kajszi’, and ‘Rózsakajszi C.1406’.

		Observed BM	Estimated BM	${\mathit{F}}_{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}$	${\mathit{C}}_{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}={\mathit{C}}_{\mathit{t}\mathit{o}\mathit{t}}$	${\mathit{F}}_{\mathit{t}\mathit{o}\mathit{t}}$
Ceglédi bíborkajszi	Mean	207.9	208.4	14.2	12.7	14.6
	StDev	10.8	12.1	0.4	3.8	0.5
	Range	48.0	51.0	1.5	14.8	1.7
	LCI	203.8	203.7	14.1	11.3	14.4
	UCI	212.1	213.1	14.4	14.2	14.7
	Slope	−0.224	−0.219	0.004	−0.041	−0.001
	p	0.438	0.500	0.696	0.686	0.933
Gönci magyar kajszi	Mean	210.3	210.2	14.0	29.8	14.3
	StDev	10.6	11.2	0.8	6.5	0.8
	Range	46.0	47.0	2.6	22.3	2.8
	LCI	206.2	205.9	13.7	27.3	14.0
	UCI	214.3	214.5	14.3	32.3	14.6
	Slope	−0.283	−0.213	0.001	−0.005	−0.003
	p	0.316	0.479	0.978	0.977	0.881
Rózsakajszi C.1406	Mean	212.9	212.8	16.4	19.7	16.8
	StDev	10.0	10.2	0.7	4.7	0.7
	Range	41.0	40.0	2.6	18.9	3.2
	LCI	209.0	208.9	16.2	17.9	16.5
	UCI	216.7	216.7	16.7	21.5	17.1
	Slope	−0.305	−0.257	−0.001	0.012	0.003
	p	0.252	0.344	0.933	0.925	0.880

). Since the slopes of the parameters $F_{crit}$, $C_{crit}=C_{tot}$, and $F_{tot}$ are almost equal to zero, i.e., they seem to not be changing with time as the blooming start does, our model can be directly applied in climate change impact studies as well. While the unsatisfactory or slowed down chilling unit accumulation delays the blooming start, the warming temperatures speed up the forcing accumulation. This complex relation with opposite directions can be simulated by our model.

4. Discussion and Conclusions

Phenology models play an important role in planning agricultural practice [3,4,5], and also in research of phenological changes in plants caused by climate change [6,7,8,9]. In our work, we fit Chuine’s unified model [43] to the blooming data sets of three apricot cultivars (‘Ceglédi bíborkajszi’, ‘Gönci magyar kajszi’, and ‘Rózsakajszi C.1406’) in Hungary in time period 1994–2020. Szalay et al. [34] analyzed the same data sets in the period 1994–2018. They observe a shift in string stage of ca. 0.5 days/year later, and in the beginning of blooming of ca. 0.125 days/year earlier. They explained this change by the warming climate in the recent decades in Hungary [71]. The chilling units accumulate more slowly in mild winters, which delays the phenological phases in dormancy [72,73]. On the other hand, the warmer spring causes an earlier beginning of blooming [31].

Due to the higher winter temperatures, and the inhibitory effect on chilling unit accumulation, the phenology models that consider only forcing unit accumulation are no longer suitable for estimating the beginning of bud burst and blooming [74,75,76,77]. That is why we chose the unified model [43], which contains both chilling and forcing unit accumulation. For the estimation of the parameters, we used the simulated annealing method [51,52], which avoids the problem of getting stuck in local minima.

We reduced the nine-parameter model to a six-parameter model with two simplifications [7,43], thus, we could accelerate the parameter estimation [63,64]. We show that there is a strong, non-linear relationship between the parameters $c_f$ and $w$. After obtaining a great number of parameters from our random walk optimization, using the histograms of the limit parameters, we defined the local optimum parameter values as the median of the bins with a significantly high number of limits that the optimization converged to. Finally, we defined the global optimum parameter vector as the one with the lowest root-mean-square error (RMSE). It is an interesting result that two-thirds of the global optimum parameter values do not fall into the local optimum parameter bins.

We estimated the beginning of blooming during the period of 1994–2020, using the global optimum parameter vector with an RMSE of less than 2.5 days, which can be considered as a good estimation compared to the results of several papers (2.6–5.6; [8,9,68,78]).

In our calculations, the maximum of the chilling unit accumulation is between 1.50 and 2.42 °C, depending on the cultivars. The inflection point of the forcing unit accumulation curve is between 8.30 and 9.04 °C, while the forcing unit accumulation is almost maximal (more than 0.9 units) between 12 °C and 15 °C, depending on the cultivars. Other research [68,79,80] determined the maximum of the chilling and forcing unit accumulation for European fruit trees at temperatures that are slightly higher: around 5 °C and around 25 °C, respectively. But there may be differences in response to temperature even within populations of a given species [81]. Thus, it is possible that the results of various studies cover a relatively wide range of chilling and forcing unit requirement values. The maximum of chilling unit accumulation is found between −28 °C and 9 °C for Quercus mongolica [63] and Taxus baccata [43], respectively, and the inflection point of forcing unit accumulation is shown to be between −14.5 °C and 16 °C in the case of Olea europaca [43] and Vitis vinifera L. [7], respectively. In the case of extreme negative values, the limit is probably statistically the best parameter vector, but biologically it may be incorrect or hardly interpretable [68]. In the case of mechanistic phenology models, however, the biological representativeness is essential.

In some cases (when parameter $k$ is around −10⁻⁷ or more), a simplification can be applied with $k=0$ and $w=F_{crit}$. But, in view of our research, this simplification is not justified, because the absolute value of parameter $k$ is larger (between −0.0072 and −0.0086) than it is in publications of other researchers (−10⁻⁸–−10⁻⁴) [8,9,43,64]. It can mean that, in the case of Hungarian apricots, the chilling unit accumulation has a larger effect on forcing unit accumulation.

With our adapted Chuine’s unified model, in the future we can provide a more accurate estimation for the blooming time of the apricot varieties in the studied area. Therefore, the preparation for frost-hazardous conditions can become notably more effective. Since the blooming times of apricot varieties in Hungary are quite close to each other [21], and the climatic conditions of apricot growing areas are also very similar [82], our results can easily be adapted to other Hungarian varieties and regions.

Making estimations for the future using climate model data requires other preliminary studies, because the extension of the phenological models for the future may be questionable [7,9,64,68,83,84], though we also can find very promising new results [85,86].