Introduction

Grapevine is one of the most profitable crops in Europe. European winegrowing regions represent about 45 % of the cultivated area worldwide (Fraga et al., 2019), and grapevine cultivation is often the result of site-specific environmental and human interactions that contribute to the quality of the final product. However, the high specificity of the climate niches in which grapevine is cultivated renders the crop particularly susceptible to the impact of climate change, which is reflected in the decrease in suitable premium quality wine regions and/or the shift of varieties from their original environment. Previous studies have already demonstrated that current warming is advancing grapevine phenological stages and shifting the annual growing cycle across Europe (Ramos et al., 2008; Tomasi et al., 2011; van Leeuwen et al., 2019). Future temperature increases, as anticipated by future projections (Jacob et al., 2014), will reinforce this trend, forcing the vine to compete in warmer climates across Europe (Molitor & Junk, 2019; Moriondo et al., 2011; Wolkovich et al., 2018). An advance or delay in phenological phases in response to a changing climate causes the growth cycle of crops to occur in hitherto new environmental conditions, possibly exposing critical phenological phases to extreme weather events (Kartschall et al., 2015; Moriondo & Bindi, 2008). As a matter of fact, many studies have addressed issues related to the occurrence of extreme events during sensitive phenological phases on grapevine (Santos et al., 2020).

Under increasing temperature, budbreak may be determined by physico-chemical cellular changes (Dokoozlian, 1999; Kovaleski et al., 2023). Bud vegetative growth is characterised by an initial phase induced by physiological conditions, in which bud exposure is regulated by cold temperatures (endo-dormancy; Dokoozlian, 1999; Anzanello et al., 2018), followed by a second phase in which bud opening is affected by warm-temperature conditions (eco-dormancy; Sarvas, 1974). In this context, phenological models aim to simulate budbreak using one of two approaches based on temperature accumulation: chilling-forcing (CF models), which is used to estimate both endo- and eco-dormancy phases, and forcing schemes (F models), which are initiated after the expected completion of the endo-dormancy period. Previous studies have shown that the endo- and eco-dormancy phases may not be independent and an overlapping accumulation of chilling and forcing units have been found in many deciduous species (Luedeling et al., 2021; Pope et al., 2014). However, this behaviour has not been clearly observed in grapevine (Prats-Llinàs et al., 2019).

The universally adopted approach when exploring the possible impacts of climate change on grapevine phenology first involves calibrating and validating a phenological model in the present period. Observed phenological datasets are essential for training models to simulate grapevine phenology (e.g., Leolini et al., 2020; Morales-Castilla et al., 2020). The obtained parametrisations in the present period are then used for the spatial model application under future climate scenarios. The use of existing parametrisations, which allow the physiology of the cultivar to be produced, are particularly useful during spatial model application (Sgubin et al., 2023), for which the scarcity of observed data makes it difficult to perform model calibrations. In order to simulate the impact of climate change on phenology, the phenological model is then coupled with the projections of future climates, performed by global climate models/regional climate models (GCMs/RCMs) (Ferrise et al., 2016; Fraga et al., 2020; Grillakis et al., 2022; Leolini et al., 2018; Moriondo et al., 2011). Current literature highlights that the use of different sources of uncertainty limits the reliability of estimations thus far. In many cases, future climate simulations are restricted to a single GCM or RCM (Duchêne et al., 2010; De Cortázar-Atauri et al., 2017; Leolini et al., 2018; Grillakis et al., 2022), while other studies stress the uncertainty related to the use of different RCMs to feed a single phenological model (Molitor et al., 2014; Fraga et al., 2016; Fraga et al., 2019; Morales-Castilla et al., 2020). Only few studies have used different phenological models for the estimation of budbreak under changing climate (Caffarra & Eccel, 2010; Hlaszny et al., 2011; Fila et al., 2014), which highlighted average differences of 5-8 days depending on the climate change scenario. Additionally, these differences vary both temporally and spatially. Fraga et al. (2016) and García de Cortázar-Atauri et al. (2017) found differences of 7-20 days in France, but Grillakis et al. (2022) found lower differences in Crete (< 7 days) for the distant future. These mismatches are anticipated to increase when the different RCMs are considered.

To the best of our knowledge, previous studies on different uncertainty sources have considered single sites and mainly annual crops like cereals and oil seeds (Gao et al., 2021; Jiang et al., 2022; Tao et al., 2018); none have focused on perennial crops (Chapagain et al., 2022). Recently, Sgubin et al. (2023) evaluated the uncertainties in projections of suitable grapevine areas when applying different GCMs, RCM climatic downscaling algorithms, representative concentration pathways (RCPs) and phenological approaches. However, their research did not specifically focus on evaluating uncertainties in phenological stages, and it was limited to a combination of three different GCMs and two downscaling algorithms driving four phenological models.

In this context, the objective of the present study was to evaluate the uncertainty in the simulation of budbreak under changing climate related to four different sources: (i) phenological model approaches, (ii) GCMs, (iii) RCMs, and (iv) greenhouse gas concentration scenarios (anthropogenic radiative forcing) across Europe. We used six phenological models, previously calibrated for eight grapevine varieties as described in Leolini et al. (2020). These models were applied to 35 climate scenarios of the future period that had been produced by different combinations of nine GCMs, eight RCMs and two future RCPs, and they were then compared to simulations from 25 climate scenarios in the historical period.

Materials and methods

1. Climate scenarios

The climate projections used in this study are based on the CORDEX-EUR11 ensemble, consisting of 60 RCM simulations at a very high spatial resolution (~12.5 km spatial resolution, ~57k points around Europe and Mediterranean basin). The 60 climate simulations (see Table S1) were produced by different combinations of nine global circulation models (GCMs) and eight RCM simulations for two RCPs (18 simulations for RCP 2.6 and 17 for RCP 4.5), together with 25 historical simulations. The climate dataset (Menz, 2023), containing mean, minimum and maximum air temperature ($T_{\text{mean}}$, $T_{\text{min}}$, $T_{\text{max}}$; °C) and precipitation (Prec; mm) was bias-adjusted using the ISIMIP3BASD methodology (Lange, 2019) within the Horizon 2020 European project Clim4Vitis (grant agreement no. 810176). The ISIMIP3BASD map is a trend-preserving bias adjustment using parametric quantile mapping. The E-OBS v19.0e dataset was used as a reference in the bias adjustment. For this study we chose two 30-year time-slices (baseline: 1976-2005; near-future: 2026-2055) in order to identify climate change impacts on budbreak in the near-future, and also to be aligned with winegrowers’ needs regarding the adoption of ad-hoc adaptation strategies. RCPs 2.6 and 4.5 were selected according to the expected projections of +0.4 to +1.6 ºC (RCP 2.6) and +0.9 to +2.0 ºC (RCP 4.5) temperature increases in 2046-2065 compared to the pre-industrial decades, as reported in the Global warming +1.5 °C special report of IPCC (2022).

2. Phenological models

Six phenological models (PHEs), three thermal-forcing (GDD, WANG and UNIFORC) and three chilling-forcing (BRIN, UNICHILL and UNIFIED) models were used in this study for estimating the budbreak date of different grapevine varieties in Europe under climate change scenarios.

2.1. GDD

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:

2.2. WANG

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:

2.3. UNIFORC

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.

2.4. BRIN

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.

2.5. UNICHILL

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.

2.6. UNIFIED

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).

3. Grapevine varieties and phenological model calibration

Eight grapevine varieties were selected: Chardonnay (CH), Cabernet-Sauvignon (CS), Gewürztraminer (GE), Grenache (GR), Pinot gris (PG), Riesling (RI), Touriga Franca (TF) and Touriga Nacional (TN). These varieties are dominant locally within the main European wine regions (Leolini et al., 2020) and represent a wide range of precocity levels for budbreak occurrence (CH, GR, TF and TN as warm-adapted varieties, and CS, GE, PG and RI as cold-adapted varieties). Additionally, these grapevine varieties are used for the production of high-quality wines worldwide and can be grouped according to average growing season temperatures and ripening potential (Jones, 2006).

The model calibration used in this study was obtained from Leolini et al. (2020) for each variety and phenological model. Observations from three European countries (France, Luxembourg and Portugal; see Table 1 for the geographical distribution of the data) were used for the calibration and validation of the phenological models, as well as weather data obtained from the closest and most representative weather stations for the vineyard weather conditions.

Additionally, during model calibration, the values of the initial parameters retrieved from the literature were constrained within a range of ± 30%, as determined in the method by Leolini et al. (2020), except for $T_{\text{low}}$ and $T_{\text{high}}$ (WANG), which were set at 0 and 45 ºC, respectively, and $T_b$ (GDD), which was delimited within 0-10 ºC. The starting date was set to 1 January for the F models and 1 September from the previous year for the CF models, respectively.

An additional filtering of the phenological observations was carried out during the calibration of the phenological models following the Tukey method (Tukey, 1977). Thus an initial model calibration was conducted and any differences between observations and simulations that were either lower than $Q_1-1\text{.}5\text{·IQR}$ or higher than $Q_3+1\text{.}5\text{·IQR}$ ($Q_1$ and $Q_3$ being the 25 and 75 % values of the distribution, respectively, and ${\text{IQR=Q}}_3-Q_1$) in at least four phenological models were considered outliers and removed from further calibrations (on average 4.8 % of the budbreak observations). Ten random subsamples were obtained, of which 60 % was assigned to the calibration and 40% to the validation of the models. The calibration was performed using the Phenological Modelling Platform (PMP; Chuine et al., 2013) which is based on the Metropolis-Hastings algorithm (Hastings, 1970; Metropolis et al., 1953). The best model parameterisations of the 10 fitted subsamples were further averaged (AV parametrisation; Leolini et al., 2020) and used to simulate budbreak in the future climate scenarios (Menz, 2023).

4. Sources of uncertainty in budbreak estimation

The analysis of the sources of uncertainty in agriculture, as well as in environmental modelling, has been performed using different methodologies (Sgubin et al., 2023; Wallach & Thorburn, 2017). Analysis of variance (ANOVA) is considered to be one of the most robust of these approaches (Chapagain et al., 2022) when determining the importance of the main sources of uncertainty in simulating crop phenology, yield or biomass (Fila et al., 2014; Jiang et al., 2022; Tao et al., 2018; Wang et al., 2020).

In this study, ANOVA was performed to analyse the total variation of the model (${\text{TSS}}_{\text{TOT}}$; Equation 9), explained by a sum of variations (${\text{TSS}}_n$) of the different sources (i.e., GCMs, RCMs, RCPs and PHEs) and their first order interactions in each 30-year each time-slice; i.e., to estimate the weight of each source in the budbreak simulations. We performed the ANOVA on each grid point separately using the following equation:

where ${\text{TSS}}_{\text{GCM}}$, ${\text{TSS}}_{\text{RCM}}$, ${\text{TSS}}_{\text{RCP}}$, ${\text{TSS}}_{\text{PHE}}$ and ${\text{TSS}}_{\epsilon }$ are the total sum of squares from the PHE, GCM, RCM, RCP, and the residuals of the ANOVA, respectively. Additionally, ${\text{TSS}}_{\text{var}\text{1:}\text{var}2}$ represents the first order interactions between two variables. The grid points included in the ANOVA analysis were those in which budbreak date was simulated in at least 80 % of the cases (i.e., climate scenarios x 30-year window) by a phenological model. Conversely, the grid points in which budbreak date was not simulated in 20 % or more of the cases were excluded from the analysis. The simulation of the budbreak was considered to be successful if it was estimated between the $t_0$ (start of the model simulation) of the current year and the following year.

Finally, the statistical analyses were carried out using the statistical software R v. 4.2 (R Core Team, 2022).

Results

1. Spatial variation in budbreak simulations

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.

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.

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.

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).

2. Sources of uncertainty in budbreak estimation

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.

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).

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).

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.

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).

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.

Discussion

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.

Conclusions

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.

Acknowledgements

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.