The GDD model (Reaumur, 1735) utilises the daily ($i$) mean temperature ($T_{{\text{avg}}_i}$, ºC) accumulation above a plant-specific threshold ($T_{\text{low}}$, ºC). The growing degree unit accumulation (${\text{Forc}}_{\text{GDD}}$) starts at a fixed day-of-year (DOY) ($t_0$) and ends on the budbreak date, when a critical sum of degree-days is achieved ($t_{\text{Bud}}$; Equation 1). Thus, the GDD model describes the eco-dormancy period, assuming that the endo-dormancy has been previously satisfied:

where:

The WANG model (Wang & Engel, 1998) utilises a non-linear and non-symmetric function (${\text{Forc}}_{\text{WANG}}$) based on the daily ($i$) average air temperature ($T_{{\text{avg}}_i}$, °C) accumulation from a fixed starting DOY ($t_0$) to the budbreak date ($t_{\text{Bud}}$). The WANG model is used as a forcing model considering the daily thermal unit rate based on three cardinal temperatures: $T_{\text{opt}}$ (optimum air temperature at which the development rate is maximum, ºC), $T_{\text{high}}$ (maximum air temperature above which the development rate is zero, ºC) and $T_{\text{low}}$ (minimum air temperature below which the development rate is zero, ºC; Equations 2 and 3).

where :

with:

The UNIFORC model (Chuine, 2000) is a sigmoidal-based model for describing the eco-dormancy period by calculating the forcing unit accumulation (${\text{Forc}}_{\text{UNIFORC}}$) from a fixed DOY ($t_0$) until budbreak occurrence ($t_{\text{Bud}}$; Equation 4). The model is based on the daily ($i$) average temperature ($T_{{\text{avg}}_i}$, ºC) and two empirical parameters $d$ and $e$ (being $d$ < 0 and $e$ > 0) related to the sharpness of the curve and the mid-temperature response, respectively.

The BRIN model (de Cortázar-Atauri et al., 2009) is used as a chilling-forcing model due to the sequential combination of two phenological functions: the Bidabe’s Cold Action Model (Bidabe, 1965) and the Richardson model (Richardson et al., 1974) for calculating the endo- and eco-dormancy, respectively. The endo-dormancy phase accumulates chilling (cold) units (${\text{Chill}}_{\text{BRIN}}$; Equation 5) based on the thermal dormancy response ($Q_{\text{10}}$) and accumulates from a fixed DOY ($t_0$) until the chilling requirement ($t_{C_{\text{endo}}}$) is achieved. The eco-dormancy period is based on the forcing unit accumulation (${\text{Forc}}_{\text{BRIN}}$, Equation 6) similar to GDD, but limited by lower ($T_{\text{low}}$, ºC) and upper cardinal temperatures ($T_{\text{high}}$, ºC).

where:

where $T_{{\text{avg}}_i}$, $T_{{\text{min}}_i}$ and $T_{{\text{max}}_i}$ are average, minimum and maximum air temperatures (ºC) at day $i$, respectively.

UNICHILL (Chuine, 2000) is a two-stage sequential model used as a chilling-forcing model and is composed of the Chuine’s function for estimating the chilling accumulation (${\text{Chill}}_{\text{CHUINE}}$, Equation 7) during the endo-dormancy period, and the UNIFORC model (see Equation 4) for the forcing unit accumulation (${\text{Forc}}_{\text{UNIFORC}}$, Equation 4) during the eco-dormancy period. Firstly, the model starts accumulating the chilling units at a fixed DOY ($t_0$) until the chilling requirement is met ($t_{C_{\text{endo}}}$), and then the forcing units are accumulated until the forcing requirement is reached ($t_{\text{Bud}}$).

where $t_0$ is the starting DOY, $T_{{\text{avg}}_i}$ is daily average temperature (ºC, day $i$) and $a$, $b$ and $c$ are empirical parameters of the Chuine’s function.

The UNIFIED model (Chuine, 2000) is a Chilling-Forcing model based on the UNICHILL model, which includes the overlapped effect of chilling units on forcing accumulation (during the eco-dormancy period) by implementing an exponential model (Equation 8). Budbreak occurs when both chilling and forcing requirements are satisfied.

${\text{ForcReq}}_{\text{UNIFIED}}$ is the critical amount of forcing units accumulated and ${\text{Chill}}_{\text{UNIFIED}}$ is the total amount of chilling units at day $i$ (Equation 7), and w and $z$ two empirical parameters (with $w$ > 0 and $z$ < 0).

Table 2 summarises the DOY ranges of estimated budbreak in Europe, shown both spatially and temporally in Figure 1, while Table S2 shows the grid points where budbreak was not reached in at least 80 % of the simulations of each phenological model and period.

Table 2. Average differences in budbreak DOYs (Q95% - Q5%; days) reported in Figure S1 for each grapevine variety and period. The table also shows the percentage of grid points in which budbreak date was not simulated in 20 % or more of the cases by a phenological model. The abbreviations of the varieties are contained in Table 1.

		CH	CS	GE	GR	PG	RI	TF	TN
Range (days)	1976-2005	12.1	18.4	12.6	23	12.4	12.4	15.3	15.7
	2026-2055	16.5	21.6	17.5	24.1	16.7	16.3	18.4	18.2
Excluded (%)	1976-2005	6.9	13.6	10.5	7.6	13.2	35.0	10.0	13.7
	2026-2055	2.9	8.9	10.4	7.3	13.8	33.6	5.3	8.4

Regarding the baseline period, 13.8 ± 9.0 % of the grid points were excluded from further calculation and mainly corresponded to the CF models (${\text{CF}}_{\text{excl}}$ = 13.2 ± 8.6%, $F_{\text{excl}}$ = 5.1 ± 1.0%; Table S2), particularly UNICHILL (9.2 %). Spatially (Figure 1), the non-simulated grid points were predominantly situated in Iceland, the Scandinavian mountains, the Alps, the Pyrenees, the Middle East and North Africa, but their extent varied depending on the variety: the cold-adapted varieties, such as GE, PG and RI, showed a larger excluded area in Southern Europe and North Africa, while the warm-adapted varieties (such as TF and TN) showed a larger excluded area in the Scandinavian mountains and Kola Peninsula. RI stood out for showing the highest proportion of excluded grid points (35.0 %), particularly in Southern Iberia and Eastern Europe.

Figure 1. Differences between the Q5% and Q95% quantiles of the budbreak simulations (range = Q95% - Q5%, i.e., 90 % confidence level of the budbreak simulations; in days) in the baseline (1976-2005) and future (2026-2055) periods for all grapevine varieties. The abbreviations of the varieties are presented in Table 1. The ranges equal to or greater than 100 days are symbolised as 100+. Dark regions indicate the grid points in which budbreak date was not simulated by 20 % or more of the cases by a phenological model. The stars show the location of the sites for the calibration and validation of the phenological models.

The correctly estimated grid points within the baseline period showed a mean range of budbreak estimations of ${\text{∆BB}}_{\text{baseline}}$ = 15.2 ± 3.8 days (average ± standard deviation between varieties; Table 2). The agreement (i.e., lower difference) in the budbreak estimates was higher in Central Europe (Figure 1). By contrast, the widest budbreak range (i.e., where budbreak estimations differed the most), was in the Mediterranean area and, to a lesser extent, in Western Europe, close to the regions where budbreak was not correctly simulated. Among the varieties, GR and CS showed the widest estimated budbreak range, and CH, GE, PG and RI the narrowest (Table S1).

Regarding the near future period, 11.3 % ± 9.6 % of the grid points was excluded (Table S2), which is slightly less than the baseline period. Wrongly simulated grid points in the near-future period mainly concerned CF models (${\text{CF}}_{\text{excl}}$ = 10.8% ± 9.1%, $F_{\text{excl}}$ = 2.7% ± 0.8%), in particular UNICHILL (8.2 %). The spatial distribution of the non-simulated grid points (Figure 1) was similar to the baseline period (i.e., Iceland, the Scandinavian mountains, the Kola Peninsula, the Alps, the Middle East and North Africa). Their extent decreased compared to the baseline period (except for GE and PG), particularly in the Alps and Northern and Eastern Europe. Conversely, the number of non-simulated grid points increased in the southern regions, especially for the cold-adapted varieties.

In addition, the Q5% to Q95% range for the correctly simulated grid points when estimating budbreak increased in the near future scenarios (${\text{∆BB}}_{\text{near}-\text{future}}=$ 18.7 ± 2.8 days; Table 2), with the lowest in Central Europe and the highest in the western and southern regions. In the particular case of GR the agreement in budbreak estimation increased in central and northern Spain, France and central Italy compared to the baseline period.

Figure 2 illustrates the uncertainty (in DOYs) associated with the different subsets of models and varieties in the Q5% to Q95% interval in Europe. In general, the simulated budbreak occurrences ranged between DOYs 20 and 240 with one or two peaks in the distribution. The main peak was around DOY 132.8 ± 5.5 days, while the second, if present, was around 10-20 days later. In general, these plots evidenced the average occurrence of budbreak during the baseline period at DOY 138.7 ± 4.6 days. An earlier occurrence of budbreak was also observed for RCP4.5 than for RCP2.6 relative to the baseline period (${\text{∆BB}}_{\text{RCP}}$ = -7.0 ± 0.6 days), the difference between the RCPs being lower (${\text{∆BB}}_{\text{RCP}4\text{.}5-\text{RCP}2\text{.}6}$ = -0.7 ± 0.6 days). Regarding GCMs and RCMs, the comparison of the peaks of the distribution showed that the Q5% - Q95% range of the GCMs was higher than that of the RCMs: on average ${\text{∆BB}}_{\text{GCM}}$ = 8.9 ± 4.2 days and ${\text{∆BB}}_{\text{RCM}}$ = 5.8 ± 0.5 days between the Q95% and Q5% (both baseline and near-future together), respectively.

Figure 2. Distribution of the frequency rates of the budbreak occurrence (DOY) for each variety in relation to each GCM, RCM, RCP and PHE. The dashed lines show the sources with the lowest and highest peaks of the budbreak DOYs distribution by source. The abbreviations of the varieties and the GCM models are in Tables 1 and S1, respectively.

Differences were also found between the phenological model types (Figure 2). F models tended to simulate a larger DOY-range of budbreak occurrences with lower differences between them (i.e., GDD, UNIFORC and WANG) compared to the CF models. In particular, F models tended to predict the budbreak later than the CF models ($\overline{{\text{BB}}_F}-\overline{{\text{BB}}_{\text{CF}}}$ = 5.3 ± 3.6 days). Specifically, UNIFIED ($\overline{{\text{BB}}_{\text{UNIFIED}}}$ = 127.9 ± 5.0 days) and WANG ($\overline{{\text{BB}}_{\text{WANG}}}$ = 137.6 ± 7.1 days) were the phenological models with the earliest and latest budbreak occurrences, respectively. Estimations of the budbreak of individual grape varieties were found to differ: cold-adapted varieties, such as RI (Figure 2), showed a difference of 95.2 ± 7.6 days depending on model type (F or CF) for the earliest budbreak (left side of the curves) in comparison to warm-adapted varieties (such as CH, with a range of 24.8 ± 13.2 days).

The ANOVAs (Table 3) showed that the average TSS in the baseline period was 5687.9 ± 13188.1. Within the TSS, the PHE contained the largest proportion of uncertainty in the estimation of budbreak occurrence ($\overline{{\text{TSS}}_{\text{PHE}}}$ = 94.7 %, which was the highest for the warm-adapted varieties CS, GR, TF and TN: $\overline{{\text{TSS}}_{{\text{PHE}}_{\text{warm}}}}$ = 97.3 %), followed by GCM ($\overline{{\text{TSS}}_{\text{GCM}}}$= 2.6 %, $\overline{{\text{TSS}}_{{\text{GCM}}_{\text{warm}}}}$= 3.9 %, $\overline{{\text{TSS}}_{{\text{GCM}}_{\text{cold}}}}$ = 1.3 %). The sum of the first-order interactions accounted for 1.5 % of the TSS (slightly higher for the warm-adapted varieties with respect to the cold-adapted varieties: 2.1 % and 0.8 %, respectively), with the highest uncertainty from GCM:RCM ($\overline{{\text{TSS}}_{\text{GCM}:\text{RCM}}}$ = 0.6 %). Lastly, the residuals of the ANOVAs (${\text{TSS}}_{\epsilon }$) accounted for 0.2 % of TSS (Figure 3). TSS varied across Europe according to the distribution of the sources of uncertainty (Figure 3). The highest values of TSS in the baseline period were found in Southern Europe, North Africa and the Middle East, especially for the cold-adapted varieties (GE, PG or RI). CS also showed high TSS values in the north-eastern region. TSS is mainly explained by PHE, which was found to be the highest in the southern (i.e., Iberian Peninsula and Middle East) and northern (Scandinavian Peninsula) regions, while the uncertainties of GCM, RCM and the first-order interactions were higher in central Europe, especially for GE, PG and RI. Similarly, the residuals were slightly higher in the central and eastern regions, especially for the cold-adapted varieties.

Table 3. Average and standard deviation (within brackets) of the uncertainties for $\overline{\text{PHE}}$, $\overline{\text{GCM}}$, $\overline{\text{RCM}}$, $\overline{\text{RCP}}$ and their first order interactions (all points together and calculated as $\text{TSS}{\text{\%}}_n=\frac{{\text{TSS}}_n}{{\text{TSS}}_{\text{TOT}}}\ast \text{100}$) by grape variety and period (baseline = 1976-2005, near-future = 2026-2055) for both phenological model types together (Forcing and Chilling-Forcing). $\overline{\text{RES}}$ are the average residuals of the ANOVA. The abbreviations of the varieties are presented in Table 1.

	Var.	$\overline{\text{TSS}}$	$\overline{\text{PHE}}$	$\overline{\text{GCM}}$	$\overline{\text{RCM}}$	$\overline{\text{RCP}}$	$\overline{\text{PHE}:\text{GCM}}$	$\overline{\text{PHE}:\text{RCM}}$	$\overline{\text{PHE}:\text{RCP}}$	$\overline{\text{GCM}:\text{RCM}}$	$\overline{\text{GCM}:\text{RCP}}$	$\overline{\text{RCM}:\text{RCP}}$	$\overline{\text{RES}}$
1976-2005	CH	2363.6 (5367.3)	93.3 (5.5)	3.4 (3.0)	1.5 (1.8)	-	0.5 (0.4)	0.3 (0.3)	-	0.8 (0.9)	-	-	0.1 (0.2)
	CS	6727.3 (10175.7)	97.9 (2.0)	0.9 (0.9)	0.4 (0.5)	-	0.4 (0.4)	0.2 (0.2)	-	0.2 (0.2)	-	-	0.1 (0.1)
	GE	10730.6 (37928)	90.1 (8.8)	5.1 (4.4)	2.3 (3.5)	-	0.7 (0.9)	0.5 (0.6)	-	1.2 (1.3)	-	-	0.2 (0.4)
	GR	6125.3 (3604.2)	97.9 (3.0)	0.9 (1.4)	0.4 (0.7)	-	0.3 (0.5)	0.2 (0.3)	-	0.2 (0.3)	-	-	0.1 (0.1)
	PG	8407.4 (27588.6)	93.0 (5.9)	3.3 (2.6)	1.5 (2.1)	-	0.7 (0.8)	0.5 (0.6)	-	0.8 (1.0)	-	-	0.2 (0.4)
	RI	5098.1 (13063.6)	91.8 (6.9)	3.8 (3.3)	1.6 (2.1)	-	0.8 (1.0)	0.6 (0.7)	-	1.0 (1.1)	-	-	0.3 (0.5)
	TF	3377.7 (5582.3)	96.8 (2.2)	1.6 (1.2)	0.6 (0.7)	-	0.3 (0.3)	0.2 (0.2)	-	0.3 (0.3)	-	-	0.1 (0.1)
	TN	2673.2 (2194.9)	96.6 (1.7)	1.7 (0.9)	0.7 (0.6)	-	0.3 (0.2)	0.3 (0.2)	-	0.3 (0.3)	-	-	0.1 (0.1)
2026-2055	CH	5723.6 (10636.9)	50.3 (18.7)	29 (12.9)	12.6 (7.3)	1.3 (1.5)	0.8 (0.7)	0.4 (0.3)	<0.1 (<0.1)	2.3 (1.7)	1.5 (2.0)	1.4 (1.5)	0.4 (0.3)
	CS	9864.5 (9862.0)	73.9 (14.4)	14.4 (8.9)	6.0 (4.4)	0.7 (1.0)	1.5 (1.8)	0.6 (0.4)	0.1 (0.1)	1.1 (0.9)	0.7 (0.9)	0.5 (0.6)	0.5 (0.4)
	GE	16105.8 (46568)	46.4 (27.3)	28.5 (15.9)	15.5 (10.2)	1.3 (1.5)	1.6 (1.8)	0.7 (0.6)	0.1 (0.1)	2.8 (2.3)	1.1 (1.1)	1.5 (1.5)	0.6 (0.4)
	GR	9703.5 (5669.9)	76.5 (13.7)	13.2 (8.1)	5.4 (4.0)	0.6 (1.0)	1.1 (1.5)	0.4 (0.3)	0.1 (0.1)	1.0 (0.9)	0.7 (1.1)	0.6 (0.7)	0.4 (0.5)
	PG	10587 (26836.7)	52.1 (22.9)	24.7 (13.1)	13.5 (8.7)	1.1 (1.2)	2.0 (2.3)	0.9 (0.8)	0.1 (0.2)	2.4 (2.0)	1.0 (1.1)	1.3 (1.2)	0.8 (0.5)
	RI	7566.5 (15470.3)	49.5 (24.2)	26.0 (13.7)	14.0 (9.3)	1.1 (1.2)	2.7 (3.0)	0.8 (0.7)	0.1 (0.2)	2.4 (2.0)	1.2 (1.2)	1.3 (1.3)	0.9 (0.5)
	TF	7037.4 (11632.4)	65.1 (14.4)	20.4 (10.2)	8.7 (5.0)	0.9 (1.1)	0.9 (1.1)	0.3 (0.2)	<0.1 (0.1)	1.6 (1.1)	1.0 (1.3)	0.8 (1.0)	0.4 (0.3)
	TN	5482.4 (4212.5)	61.7 (13.7)	22.9 (9.8)	9.5 (5.6)	1.1 (0.9)	0.7 (0.4)	0.3 (0.2)	<0.1 (<0.1)	1.7 (1.2)	1.0 (1.1)	0.7 (0.8)	0.3 (0.2)

Figure 3. Spatial distribution of the Total Sum of Squares (TSS) and the uncertainty (%, calculated as $\text{TSS}{\text{\%}}_n=\frac{{\text{TSS}}_n}{{\text{TSS}}_{\text{TOT}}}\ast \text{100}$) from PHE, GCM, RCM and their first order interactions (var1:var2) in Europe for the baseline period (1976-2005). RCPs are missing in the baseline period. Figure 1 shows the location of the calibration and validation sites of the phenological models. Dark regions indicate the grid points in which budbreak date was not simulated by 20 % or more of the cases by a phenological model.

TSS increased in the near-future period (9008.8 ± 16381.1), with the highest uncertainty in the estimation of budbreak occurrence ($\overline{{\text{TSS}}_{\text{PHE}}}$ = 59.4 %) being associated with PHE, especially for CS, GR, TF and TN varieties ($\overline{{\text{TSS}}_{{\text{PHE}}_{\text{warm}}}}$ = 69.3 %, $\overline{{\text{TSS}}_{{\text{PHE}}_{\text{cold}}}}$ = 49.6 %). The high uncertainty explained by PHE (Table 3) was followed by GCM ($\overline{{\text{TSS}}_{\text{GCM}}}$ = 22.4 %, $\overline{{\text{TSS}}_{{\text{GCM}}_{\text{warm}}}}$ = 27.1 %, $\overline{{\text{TSS}}_{{\text{GCM}}_{\text{cold}}}}$ = 17.7 %), RCM ($\overline{{\text{TSS}}_{\text{RCM}}}$ = 10.4 %, $\overline{{\text{TSS}}_{{\text{RCM}}_{\text{warm}}}}$ = 13.9 %, $\overline{{\text{TSS}}_{{\text{RCM}}_{\text{cold}}}}=$7.4 %). The sum of all first-order interactions accounted for 6.0 % of TSS (slightly higher for the warm-adapted varieties than for the cold adapted-varieties: 7.6 % and 4.4%, respectively), with the largest contributions from GCM and RCM ($\overline{{\text{TSS}}_{\text{GCM}:\text{RCM}}}$ = 1.9 %) and PHE and GCM ($\overline{{\text{TSS}}_{\text{PHE}:\text{GCM}}}$ = 1.4 %), respectively. Finally, the residuals of the ANOVAs (${\text{TSS}}_{\epsilon }$) accounted for 0.5 % of the total sum of squares (Figures 4 and S1).

Figure 4. Spatial distribution of the Total Sum of Squares (TSS) and the uncertainty (%, calculated as $\text{TSS}{\text{\%}}_n=\frac{{\text{TSS}}_n}{{\text{TSS}}_{\text{TOT}}}\ast \text{100}$) from PHE, GCM, RCM and RCP in Europe for the future period (2026-2055). OTHERS summarises the uncertainty from the first order interactions of the previous variables (var1:var2) and the residuals of the ANOVA showed in Figure S1. Dark regions indicate the grid points in which budbreak date was not simulated by 20 % or more of the cases by a phenological model. Figure 1 shows the location of the calibration and validation sites of the phenological models.

Spatially, TSS varied across Europe in the near-future depending on the spatial distribution of the sources of uncertainty. The highest TSS values were located close to the excluded grid points, such as those of the Mediterranean basin, the Alps and the Scandinavian mountains. TSS was also high for CS and GR in north-eastern Europe. Regarding the source, the highest uncertainties associated with PHE when explaining TSS were in southern Europe for all grapevine varieties (i.e., in the Mediterranean basin), as well as in eastern Europe for CS, GR, TF and TN. Conversely, GCM and RCM uncertainties were higher in central Europe, while RCP uncertainties were highest in north-eastern and south-eastern Europe. The first order interactions between the variables were slightly higher in the west due to PHE:GCM, while GCM:RCM increased in the lowland regions from Denmark to the Black Sea. In addition, the spatial distribution of the uncertainties also varied between varieties. For instance, in the case of cold-adapted varieties (e.g., GE, PG or RI), the uncertainty was higher for GCM than for PHE in central Europe, whilst in the case of the warm-adapted varieties (e.g., TF) the uncertainty was higher for PHE than for the other sources of uncertainty throughout Europe.

In addition, budbreak estimates differed depending on phenological model type (Table S3, Figure 5). The percentage of the sum of squares from PHE was slightly higher in model type CF ($\overline{{\text{TSS}}_{{\text{PHE}}_{\text{CF}}}}$ = 83.3 %, $\overline{{\text{TSS}}_{{\text{TOT}}_{\text{CF}}}}$ = 1351.4) than in F ($\overline{{\text{TSS}}_{{\text{PHE}}_F}}$ = 81.2 %, $\overline{{\text{TSS}}_{{\text{TOT}}_F}}$ = 741.6) for the baseline and the future periods together ($\overline{{\text{TSS}}_{{\text{PHE}}_{\text{CF}}}}$ = 49.3 %, $\overline{{\text{TSS}}_{{\text{TOT}}_{\text{CF}}}}$ = 2387.3 and $\overline{{\text{TSS}}_{{\text{PHE}}_F}}$ = 38.9 %, $\overline{{\text{TSS}}_{{\text{TOT}}_F}}$ = 1866.7, for CF and F, respectively). Regarding the varieties, the cold-adapted varieties, especially GE, showed the highest contribution of PHE to explaining TSS in both model types for the baseline and the future periods together ($\overline{{\text{TSS}}_{\text{PHE}}}$ = 78.4 %, $\overline{{\text{TSS}}_{\text{TOT}}}$ = 2347.1 and $\overline{{\text{TSS}}_{\text{PHE}}}$ = 69.5 %, $\overline{{\text{TSS}}_{\text{TOT}}}$ = 1607.1, for CF and F, respectively), compared to the warm-adapted varieties ($\overline{{\text{TSS}}_{\text{PHE}}}$ = 56.3 %, $\overline{{\text{TSS}}_{\text{TOT}}}$ = 1391.6 and $\overline{{\text{TSS}}_{\text{PHE}}}$ = 54.5 %, $\overline{{\text{TSS}}_{\text{TOT}}}$ = 1001.4, for CF and F, respectively). In the case of CS, a warm-adapted variety, PHE contributed the most to explaining TSS in the F models and the least in the CF models for both periods. The opposite was true for TN, with PHE contributing the least to explaining TSS for F in both periods (Table S3). Spatially, TSS did not vary significantly with phenological model. In general, absolute ${\text{TSS}}_{{\text{PHE}}_F}$ was lower in central-western Europe, while ${\text{TSS}}_{{\text{PHE}}_{\text{CF}}}$ was generally lower in central-eastern Europe, except for GR and TF. Therefore, the regions with higher ${\text{TSS}}_{{\text{PHE}}_{\text{FCF}}}$ were characterised by high intra- and/or inter-uncertainty between the phenological model types (such as the south of the Iberian Peninsula and the Middle East), or within the regions with low ${\text{TSS}}_{{\text{PHE}}_{\text{CF}}}$ and ${\text{TSS}}_{{\text{PHE}}_F}$, but distinct budbreak simulations (such as northern Turkey).

Figure 5. Spatial distribution of the sum of squares from each phenological model type: all phenological models (FCF), the Chilling-Forcing models (CF) and the Forcing models (F) in Europe for both the baseline (1976-2005) and future period (2026-2055). Dark regions evidence the grid points in which budbreak date was not simulated by 20 % or more of the cases by a phenological model.

The difference in the uncertainty distribution associated with the phenological models was partly due to the different mechanisms underlying the temperature in the phenological models (Figure 6). More specifically, an increase in mean temperature contributed to a decrease in daily chilling rates in Bidabe’s model (BRIN, chilling). On the other hand, in the UNICHILL and UNIFIED models, only a specific range of mean temperatures (CH: -20.2 to +30.0 ºC, CS: -14.5 to +25.8 ºC, GE: -14.2 to +18.1 ºC, GR: -8.3 to +17.1 ºC, PG: -13.4 to +17.3 ºC, RI: -11.9 to +15.4 ºC, TF: -8.9 to +24.4 ºC, TN: -7.3 to +30.0 ºC) were found to contribute positively to the daily chilling unit rate estimation by excluding the effect of extreme temperatures (low and high temperatures). Regarding the forcing models, temperatures below a specific threshold (GDD: 2.8 ± 3.4 ºC, WANG: 0.1 ºC, UNIFORC: -0.34 ± 2.04 ºC, Richardson – BRIN: 3.4 ± 1.2 ºC) were not found to contribute to the accumulation of forcing units. By contrast, the increase in temperature contributed to either a linear (e.g., GDD or Richardson) or parabolic (e.g., WANG or UNIFORC, whose sigmoid functions are also used in UNICHILL and UNIFIED models) increase in the daily forcing unit rate: GDD did not show an upper threshold, because it simulated a linear increase in forcing units with the increase in mean temperature, whereas the Richardson model showed an upper threshold of 30 ºC, which limited the daily forcing unit rate at high temperature. Moreover, high temperatures decreased the daily forcing rate of the WANG model (above DOY 26.1 ± 4.7), which is described by a sigmoid curve.

Figure 6. Daily chilling/forcing unit rates in relation to mean temperature (from -30 ºC to +30 ºC) for each phenological model. The abbreviations of the varieties are presented in Table 1.

The daily chilling and forcing rates had thus accumulated in order to release of endo- and eco-dormancy phases (Figure 7). Accordingly, the six phenological models were applied to a daily climate database for southern France (43.58º N, 3.96º E; period 1950-2017; mean annual temperature = 14.2 ºC; Leolini et al., 2020) considering different daily delta temperatures from -20 to +20 ºC. The endo-dormancy phase was estimated using the BRIN, UNICHILL and UNIFIED phenological models, which are able to account for chilling unit accumulation. In this context, the BRIN model simulated a delay of the endo-dormancy phase under high temperatures (+6.5 ± 0.9 days/ºC on average), while low temperatures resulted in a delay of the phase. Meanwhile, both high (${\text{∆T}}_{\text{mean}}$ > 3.2 ± 3.4 ºC) and low temperatures (${\text{∆T}}_{\text{mean}}$ < -11.5 ± 3.8 ºC) caused a delay in endo-dormancy in the UNICHILL and UNIFIED models. The plateau of the curve of these models corresponded to a specific range of temperature which showed no effect on the endo-dormancy (same DOY for each temperature increase; -11.5 ± 3.8 ºC < ${\text{∆T}}_{\text{mean}}$ < 3.2 ± 3.4 ºC). Finally, UNICHILL simulated the same endo-dormancy DOY for any delta temperature increase (Figure 7, chilling - CS).

Figure 7. Estimation of endo-dormancy (chilling accumulation) and eco-dormancy (forcing accumulation) DOY at different temperature changes (${\text{∆T}}_{\text{mean}}$ from -20 ºC to +20 ºC) using six phenological models. The ΔTmeans were applied on an average year derived from the daily temperature data of southern France (Latitude = 43.58º N, Longitude = 3.96º E; Leolini et al., 2020) for the period 1950-2017 (annual Tmean = 14.2 ºC). The abbreviations of the varieties are presented in Table 1. The vertical grey dashed line indicates the value of ${\text{∆T}}_{\text{mean}}$ = 0 ºC. Negative DOY values symbolise the end of the phase from the preceding year. $t_0$ in chilling models is DOY -121 (or 244 from the previous year), while the forcing models were initiated on 1 January.

Moreover, low temperatures during the Forcing units’ accumulation (eco-dormancy phase) resulted in an estimated delay in budbreak, with an uneven response in UNIFORC, UNIFIED and UNICHILL. The difference was most significant for CS, GE, GR, PG and RI. Additionally, an increase in temperature was found to anticipate budbreak for F models (i.e., GDD, WANG and UNIFORC; -7.6 ± 0.9 days/${\text{∆T}}_{\text{mean}}$). However, at high temperatures CF models estimated delayed eco-dormancy, even beyond the $t_0$ of the following year, and was therefore deemed unachieved. In such conditions, the uncertainty of the estimation of budbreak increased.

In the present study, we investigated the different sources of uncertainty in the simulations of grapevine budbreak for both the baseline (1976-2005) and near-term future (2026-2055) periods across Europe. The purpose of this analysis was to explore the sources of uncertainty in budbreak simulations, rather than to provide answers about the suitability of the grapevine growing areas (which would need to be verified by considering the entire grapevine seasonal cycle (until maturity; Fraga et al., 2016; Sgubin et al., 2023)). To this end, phenological models were applied across the whole of Europe, including unsuitable grapevine growing areas in which cultivation is limited by temperature and soil and water availability; only the areas unsuitable for the potential occurrence of budbreak were excluded.

Moreover, the lack of comprehensive phenological datasets of specific varieties from diverse climatic conditions meant that existing site-specific model parametrisations had to be used. The accuracy of the spatial application of the phenological models using these parametrisations thus limited their ability to reproduce vine phenology under different climates. Despite this being a critical aspect of spatial model applications and the need for an ad-hoc evaluation of this source of uncertainty, the use of existing parametrisations for spatial models has been adopted in other studies (Sgubin et al., 2023). In contrast to previous studies, which focused exclusively on either specific sample locations (Wang et al., 2020) or regional averages (Ferrise et al., 2016; Ramirez-Villegas et al., 2017), here we utilised gridded climate projections to evaluate spatial uncertainty throughout the whole of Europe. The spatial distribution of the total uncertainty was found to vary depending on the spatial distribution of the fractional uncertainties of the sources (i.e., the GCMs, RCMs, RCPs and PHEs) and their interactions (Figures 3, 4 and S1).

Phenological models did not correctly simulate all the grid points in both the present and near-future periods. In the baseline period, an average of 13.8 % of the grid points, which were mainly located in northern and eastern Europe and the Mediterranean basin and mountain areas, were not correctly simulated by one or more phenological models in more than 20% of the cases (phenological models x climate scenarios x 30-years window). This trend slightly decreased in the near-future period to 11.8 % on average, showing an increase in suitable regions in terms of grapevine budbreak occurrence in the near future period; however, this trend varied spatially. The area where budbreak was not correctly simulated was larger in the near-future scenarios in the Mediterranean basin regions (e.g., southern Iberian Peninsula, the Middle East, Morocco and Italy), being mainly related to the CF models. By contrast, northern and eastern Europe were more correctly simulated than in the baseline period. This finding is in line with previous studies that estimated a future shift of the European wine-producing regions (e.g., Molitor & Junk, 2019; Moriondo et al., 2011; Wolkovich et al., 2018), continuing the trend actually observed (Ramos et al., 2008; Tomasi et al., 2011; van Leeuwen et al., 2019).

Therefore, the phenological model approach was found to be, on average, the main source of uncertainty in budbreak estimation, both in the baseline and near-future periods. This uncertainty is due to the intrinsic structure of the phenological models and the data they require. Firstly, the uncertainty from F models can be partially explained by the identification of the starting date ($t_0$) for the eco-dormancy period estimation. In these models, $t_0$ is generally set (de Cortázar-Atauri et al., 2009; Leolini et al., 2020) or calibrated against the observed database (Fila et al., 2014) when observations about endo-dormancy release are not available. This critically limits the use of these models for future application to different environments (e.g., Mediterranean coast; Camargo-Alvarez et al., 2020) due to the delay in the occurrence of chilling temperatures in warmer conditions, as less chill accumulation during winter delays the accumulation of forcing temperatures (Kovaleski, 2024). However, given that the endo-dormancy release observations are not easily available, the use of a fixed starting date seems to represent a compromise for forcing model application at different spatial scales and climate conditions as reported in other studies (e.g., de Cortázar-Atauri et al., 2009; Costa et al., 2019; Sgubin et al., 2023). By contrast, the use of a fixed starting date for the onset of dormancy in CF models (1 September vs 1 August) showed a lower impact on budbreak date estimation compared to the starting date of the beginning of eco-dormancy (Leolini et al., 2020). In this context, if appropriately calibrated, CF models could potentially provide a better estimation of budbreak than the simplest solution of F models. However, the use of an observed comprehensive dataset, including endo-dormancy and budbreak observations at different temperature conditions, is still crucial for obtaining a robust and plausible parameterisation of forcing and chilling-forcing models (Fila et al., 2012). In our study, the limited spatial behaviour of the observed dataset for each variety was partially overcome by the observations from the long time-series, comprising different temperature conditions and budbreak observations (Figure S2). However, in some cases, this was not enough to accurately define the optimal range for chilling unit accumulation (e.g., the UNICHILL chilling unit rate for the Chardonnay variety is highest in the temperature range of -20 ºC to +30 ºC, Figure 6), influencing the simulation of the endo- and eco-dormancy phases (Figure 7). Furthermore, these phenological models do not consider plant adaptability to different climatic conditions; therefore, grapevines may require fewer chilling units and more forcing units in the warmer regions to reach budbreak (Kovaleski et al., 2023). In warmer climates, the risk of freeze damage to buds increases due to the short bud latency and the more extreme and erratic winter cold events. These conditions were partially determined by the use of field data, which may be lacking in terms of extreme winter temperature effect on budbreak date. This effect has previously been reproduced by other authors via ad-hoc experiments leading to reliable parameterisation (e.g., Fila et al., 2012). Although these parameterisations are currently available in previous reports, the number of referenced varieties and models is limited. In this context, since the aim of the present study was to explore the impact of different parameterisations of phenological models on the uncertainty of budbreak estimation in future scenarios, the lack of plausibility of some parameter sets was an additional source of uncertainty to be taken into account.

Aside from the responses of the phenological models to temperature, budbreak simulations were shown to diverge due to the uncertainties introduced by the climate factors and their interactions (average of 4.3 % and 32.0 % in present and future, respectively). These uncertainties were lower than the those found in previous studies (Fraga et al., 2020; Grillakis et al., 2022; Wang et al., 2020), which estimated phenological phases from biased future scenarios. However, the use of uncorrected climate scenarios in climate impact assessments can often give unrealistic results (Lange, 2019), thus increasing uncertainty due to climate variables and their first-order interactions. According to Lange (2019), the CMIP5-based climate scenarios were initially characterised by future average temperature projections that were biased by around 0.1-1.5 ºC with respect to the observed temperature data from E-OBS within the period 2006-2020. These differences in daily average temperature strongly influence the response of the phenological models on the accumulation of chilling and forcing units, and thereby advance or delay the expected budbreak occurrence.

The spatial pattern was similar over the different climate-related sources and varieties (Figures 4 and S1). Specifically, the region with the highest uncertainty from GCMs was predominantly located in central and northern Europe (e.g., Germany, Finland and Norway), while the RCMs contributed more to total uncertainty in the eastern regions (e.g., Ukraine and Poland). The lowest uncertainty explained by both factors was located in the western and southern regions where PHE dominated (e.g., south Iberian Peninsula, Morocco and Middle East countries). The GCM was capable of explaining more than the 40 % of the total uncertainty of budbreak estimation, especially for CH, GE, PG and RI (most of them being cold-adapted varieties) and in a context of climate change, which is higher than the uncertainty explained by PHE. This corroborates the findings of Sgubin et al. (2023), who found higher fractional uncertainties explained by the GCMs compared to the PHEs, especially in the near-future period. Generally, the uncertainty explained by the fine-grid climate scenarios was shown to be mainly associated with GCMs, which influenced climate change signals in Europe (Christensen & Kjellström, 2020), as they have shown different temperature trends with respect to the baseline period (1976-2005, IPCC 2022). However, the uncertainty from the climatic variables (GCMs and RCMs) during the baseline period (1976-2005) was quite low. The uncertainty associated with RCMs for the baseline period (Figure 3) showed a heterogenous pattern, with mountainous areas, such as the Pyrenees and the Alps, and eastern regions showing higher values. These findings are consistent with those of Christensen and Kjellström (2020), who, for the RCMs, found distinct geographical patterns across Europe close to the eastern and northern boundaries and in mountainous areas due to strong interaction with topography, sea ice, snow or soil moisture. Additionally, Christensen and Kjellström (2020) and Sørland et al. (2018) evidenced that the RCMs differed in the simulated seasonal cycle of temperature; such differences increase the uncertainty in the budbreak simulations and further influence the estimated grapevine growth, development and final yield (Sgubin et al., 2023; Vesely et al., 2019).

The overall lowest contribution to total uncertainty was from the RCP scenarios (only included in the near-future period), which is consistent with van Vuuren et al. (2011). Concerning the greenhouse gas concentration scenarios, the results of our study show low uncertainties in the estimation of budbreak by RCPs (approximately 1.0 %, see Table 3). These findings differ from previous studies (De Cortázar-Atauri et al., 2017; Fraga et al., 2020; Grillakis et al., 2022; Leolini et al., 2018), in which differences of one to three weeks across Europe were attributed to the use of multiple RCPs for explaining budbreak occurrences in the future. These differences can be attributed to three factors: firstly, the greenhouse gas concentration of RCP2.6 was shown to be similar to RCP4.5 especially in the near-future period; secondly the scenarios were initiated in 2006 and were predicted to significantly diverge around 2020-2030 (IPCC, 2022), with the response of the climate system to that change following even later; thirdly, the bias correction applied to the climate projections further reduced the uncertainty associated with the RCPs. Hence, within the timeframe 2026-2055, the simulations did not differ significantly.

Our results show that the selection of adequate GCMs, RCMs, RCPs and PHEs (the most widely used uncertainty factors reported in previous studies, such as Jiang et al. (2022) or Sgubin et al. (2023)) represent the main sources of uncertainty. However, the interactions between the variables accounted for, on average, 6.5 % of the total uncertainty. Specifically, GCM:RCM, GCM:RCP and RCM:RCP, when present, were highest over the central, eastern and northern regions, which also influence the response of the phenological models (i.e., PHE:GCM, PHE:RCM and PHE:RCP). These estimations were lower than those previously found for China and Australia by Wang et al. (2020), who reported that the fractional uncertainty due to interaction GCM:RCP exceeded 30 % in the estimation of wheat yield, especially in the driest regions.

Phenological models are thus the main factor contributing to overall uncertainty, particularly in the context of climate change. Variations in temperature due to climate change increase the estimation range of the phenological model types (F and CF). However, the increase in winter temperatures due to climate change (as occurred in Bordeaux, France in January 2023 with ΔT_JAN = +2 ºC) not only favours the advance of budbreak occurrence but also increases the mismatch between estimations by phenological models (favouring the forcing unit accumulation and delaying the chill unit accumulation). In such conditions, budbreak could occur earlier than the estimations by the selected phenological models. Accordingly, vine growers could delay the application of adequate management practices for winter survival, particularly in areas where such occurrences are expected to be more frequent in the future. This, in turn, could lead to a reduction in number of reproductive buds (Kovaleski et al., 2023) and thus in final grape production (Leolini et al., 2018), as occurred in France in 2021 when national wine production fell by 30 % due to spring frosts.

In the present study, we explored the sources of uncertainty in estimations of grapevine budbreak both baseline (1976-2005) and near-future periods (2026-2055). To this end, we used ANOVA and different combinations of nine GCMs, eight RCMs, two RCPs and six PHEs depending on data availability. These sources of uncertainty were non-linearly linked to temperature, which was shown to be the primary driver for phenological development. Low and high temperatures resulted in high uncertainty in budbreak estimation that depended on model type (forcing or chilling-forcing), and in some cases hindered the simulation of budbreak occurrence. While an increase in temperature during the growing phase increased the length of the endo-dormancy period due to the slow accumulation of chilling units, it reduced the length of the eco-dormancy period due to the rapid accumulation of forcing units. Conversely, low temperatures during the growing phase accelerated the accumulation of chilling units, thus initiating the eco-dormancy phase earlier than 1 January when the forcing models start accumulating ($t_0$) by definition. In addition, the chilling models diverged mainly in colder areas, because BRIN is strongly influenced by low temperatures, while UNICHILL and UNIFIED models are influenced exclusively by a limited range of temperatures. By contrast, the forcing models are mainly limited by low temperatures (below t_b), except for WANG which is also influenced by high temperatures (as seen for CS in the Middle East). This was found to favour an additional uncertainty between the two types of phenological models, mainly in the areas where the vine varieties are not adapted: the lowest uncertainty associated with phenological models F was located within the warmest regions for warm-adapted varieties (e.g., GR, TF, TN), whereas it was located in the coldest regions for cold-adapted varieties (e.g., GE, PG, RI). Thus, the regions with the lowest overall uncertainty (where the phenological models were more consistent) were generally located in central Europe.

Therefore, selecting both the adequate phenological model and adequate climate scenario was found to be crucial for a better estimation of budbreak in many European regions, which in turn influences the estimation of further phenological phases. While in the calibration and validation phase the models do not show any substantial differences in budbreak estimations, the inter- and intra-differences between the phenological modelling approaches (forcing and chilling-forcing) in response to predicted future temperature increases were evident. Therefore, future budbreak simulations in central Europe need few phenological models but many climate models to cover the full uncertainty range: by contrast, the Mediterranean regions need few climate models but all phenological models to cover the full uncertainty range.

The authors would like to acknowledge the Clim4Vitis project (Climate change impact mitigation for European viticulture: knowledge transfer for an integrated approach) and the Montevitis project (Integrating a Comprehensive European Approach for Climate Change Mitigation and Adaptation in Montenegro Viticulture), funded by the European Union’s Horizon 2020 Research and Innovation Programme under the grant agreements 810176 and 101059461, respectively. SC-A and MB acknowledge the EU - Next Generation EU Mission 4 “Education and Research” - Component 2: “From research to business” - Investment 3.1: “Fund for the realisation of an integrated system of research and innovation infrastructures” - Project IR0000032 – ITINERIS - Italian Integrated Environmental Research Infrastructures System - CUP B53C22002150006. HF and JS thank the Portuguese FCT for UIDB/04033/2020.