Grapevine Phenology of cv. Touriga Franca and Touriga Nacional in the Douro Wine Region: Modelling and Climate Change Projections

Abstract

Projections of grapevine phenophases under future climate change scenarios are strategic decision support tools for viticulturists and wine producers. Several phenological models are tested for budburst, flowering, and veraison and for two main grapevine varieties (cv. Touriga Franca and Touriga Nacional) growing in the Douro Demarcated Region. Four forcing models (Growing degree-days, Richardson, Sigmoid, and Wang) and three dormancy models (Bidabe, Smoothed Utah and Chuine), with different parameterizations and combinations, are used. New datasets, combing phenology with weather station data, widespread over the Douro wine region, were used for this purpose. The eight best performing models and parameterizations were selected for each phenophase and variety, based on performance metrics. For both cultivars, results revealed moderate performances (0.4 < R² < 0.7) for budburst, while high performances (R² > 0.7) were found for flowering and veraison, particularly when Growing degree-days or Sigmoid models are used, respectively. Climate change projections were based on a two-member climate model ensemble from the EURO-CORDEX project under RCP4.5. Projections depicted an anticipation of phenophase timings by 6, 8 or 10–12 days until the end of the century for budburst, flowering, and veraison, respectively. The inter-model variability is of approximately 2–4 days for flowering and veraison and 4–6 days for budburst. These results establish grounds for the implementation of a decision support system for monitoring and short-term prediction of grapevine phenology, thus promoting a more efficient viticulture.

1. Introduction

Viticulture is one of the most important socioeconomic sectors in Portugal. From a global perspective, Portugal is currently the 11th-largest wine producer and the 9th wine exporter in the world [1]. Portugal accounts for approximately 194,000 ha of vineyard area and with a wine production around 6.7 × 10⁶ hL [1]. Almost half of this production is exported, accounting for nearly 2% of national export income [2]. The country is divided into 14 viticultural regions (12 in the mainland and 2 in the islands), with 25 Denominations of Origin (DOs). Douro/Porto, Vinhos-Verdes, Alentejo, and Lisboa are the main regions in terms of wine production. Although the Alentejo wine region is the most productive [3], the Douro/Porto wine region has the largest vineyard area and is the main region in terms of production, accounting for ca. 25% of all wine produced in Portugal [2].

The Douro/Porto wine region, also known as Douro Demarcated Region (DDR), has unique features, where vineyards are traditionally grown on terraces over steep slopes along the mountainous Douro valley and in an almost monocultural system, resulting in a UNESCO (United Nations Educational, Scientific and Cultural Organization) human heritage vineyard landscape. A large number of grapevine varieties (Vitis vinifera L.) are grown in the region, being many of them native. The complex physiographic and pedoclimatic characteristics of the DDR not only provide a wide range of terroirs for grapevine growth [4], but also make wine production very costly in terms of labour force and vineyard management [5]. Given the heterogeneous conditions in which grapevines are grown within the DDR, it is imperative to better understand the climatic influences on grapevine development in order to improve vineyard management in such a diverse region. In effect, monitoring the phenological phases and timings are some of the most useful information for planning the necessary viticultural practices in each developmental stage. This knowledge may promote more adjusted cultural practices, eventually optimizing production and quality at harvest.

It is known that grapevine development and growth are mainly driven by meteorological and climatic factors [6]. Therefore, factors such as temperature, solar radiation and precipitation may have a great influence on phenology, yields, and grape berry quality and attributes [7,8]. In particular, temperature is a fundamental factor in controlling grapevine phenological development and ripening [9,10,11,12,13]. Phenology models have been developed for a wide range of species, including grapevines, and using observational data from many different countries and regions [14]. Several previous studies have been performed in order to model the main phenological stages of grapevines, e.g., enabling a classification of varieties for technical proposes [15], to predict phenophases or to assess the impact of climate change on phenology [16].

Many different models have been developed during the last 60 years to simulate grapevine phenology. The modelling of budbreak was explored by several authors assuming different hypothesis [17,18,19,20,21,22]: taking or not into account the dormancy period, interactions between dormancy and post-dormancy periods, starting at a fixed date (e.g., January 1st, February 15th, March 15th). These models were tested for many different varieties under different climatic conditions. For flowering and veraison, the conventional linear growing degree-days model (GDD) was initially proposed by Amerine and Winkler [23]. This simple temperature accumulation model uses daily mean temperatures, with a fixed base temperature of 10 °C, and onsets temperature accumulation at April 1st. Others models have been calibrated and tested to simulate flowering and veraison, such as the linear model proposed by Duchene, et al. [24], with a formulation similar to the Winkler model, and the curvilinear model of Wang and Engel, later adapted by García de Cortázar-Atauri, et al. [25]. The latter model was improved by incorporating an optimal temperature (25 to 30 °C) and an upper critical threshold temperature (40 °C). The GFV model (Grapevine Flowering and Veraison model) was proposed by Parker, et al. [26] and Parker, et al. [27]. This model is similar to GDD, but with a 0 °C base temperature and starting temperature accumulation on March 1st. This last model was developed to classify different grapevine varieties within them. Finally, Molitor, et al. [28] developed a high-resolution model to simulate all the phenological stages described in the BBCH (Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie) scale using a capped model. As such, the performance of several phenology models should be tested at a regional level, as there is no best model for all varieties and conditions.

Phenological models may also be coupled with climate change scenarios to project the potential impacts of future climates on phenological timings and grapevine development [21,24,25,29,30,31,32,33]. Projected increases in temperature are expected to drive earlier development stages and, as a result, a general advancement of grapevine phenology [24,34,35,36,37]. Those changes will affect the ripening period, commonly occurring in the warmest part of summer, and may lead to aromas and acidity losses [38]. Consequently, currently planted varieties, which are grown under quite specific conditions today, may no longer thrive in the same place under modified environmental conditions in the future [39]. Therefore, it is critical to understand to what extent temperature may influence the timings of the reproductive and vegetative cycles, as well as identifying varietal differences in phenology and maturity [26]. While several studies have assessed the phenological timings in many regions worldwide, a comprehensive study of phenology in the DDR was not undertaken so far, due to the overall lack of phenological and meteorological data. Hence, the present study attempts to overcome this drawback, by presenting the first grapevine phenological study for the DDR.

Along the previous lines, the aims of the present study are twofold: 1) to identify suitable phenological models for three main phenological timings (budburst, flowering, and veraison) of two leading grapevine varieties (cv. Touriga Franca and cv. Touriga Nacional) currently grown in the DDR; 2) to assess the climate change impacts on these phenological timings and varieties. The datasets and methodologies are described in Section 2. The results are presented in Section 3. Section 4 will provide a discussion and a summary of the main conclusions.

2. Materials and Methods

2.1. Study Area and Datasets

The present study targets the Douro/Porto Wine Region (DDR), located in Northern Portugal (Figure 1a). The DDR is characterized by Mediterranean-type climates, with typically dry and warm summers, though a wide range of mesoclimates and microclimates do exist throughout the region, mostly driven by the multifaceted orography and distance to the North Atlantic. Important east–west gradients in both temperatures and precipitation can be found within the region, and are typified by its 3 sub-regions: from Baixo-Corgo (westernmost sector) to Cima-Corgo (middle sector) and to Douro Superior (easternmost sector) (Figure 1b). This spatial variability is further enhanced by important differences in soil properties, though soils are mostly originated from schist bedrocks [40].

In the present study, phenological data recorded at 12 vineyards (“Quintas”) in the DDR, with data available over the period 2014–2017 (4 years), were used (Figure 1b). The phenological timings were recorded through observations based on the BBCH scale adapted to grapevines and when at least 50% of a pre-defined set of homogenous plants reach the corresponding stage. Budburst (BBCH 07), flowering (BBCH 65) and veraison (BBCH 81) were considered herein [41]. Data for two main grapevines varieties grown in Portugal were selected: cv. Touriga Franca (TF) and cv. Touriga Nacional (TN). As information for other varieties (cv. Tempranillo, cv. Tinta Amarela and cv. Moscatel-Galego) was indeed very limited and inconsistent among vineyards, they were not taken into consideration. Furthermore, in order to enhance the statistical robustness of the selected models for a given grapevine variety, the phenological timings for the same variety for all sites (and terroirs) were considered for model calibration.

Owing to the relatively small sample size of the DDR dataset for model calibration (12 sites × 4 years = 48 records), the same phenological timings for TF and TN from two other Portuguese wine regions were also used: Lisboa (central-west) and Vinhos Verdes (north-west) (Figure 1a). Phenological records were available for the following locations (

Table 1. Available time periods of the phenological data used in the present study (budburst, flowering and veraison dates) by wine region (Douro, Lisboa and Vinhos-Verdes) and for each grapevine variety (cv. Touriga Franca and cv. Touriga Nacional).

Wine Region	Latitude (°N)	Longitude (°W)	Altitude (meters)	Variety (cv.)	Available Time Period
					Budburst 1	Flowering 2	Veraison 3
Douro/Porto	41°14′24″ N–41°02′24″ N	7°01′48″ W–°47′24″ W	85–588	Touriga Franca	2014–2017	2014–2017	2014–2017
Douro/Porto	41°14′24″ N–41°02′24″ N	7°01′48″ W–7°47′24″ W	85–588	Touriga Nacional	2014–2017	2014–2017	2014–2017
Lisboa	39°02′24″ N	9°10′48″ W	85	Touriga Franca	1995–2014	1995–2014	1995–2014
Lisboa	39°02′24″ N	9°10′48″ W	85	Touriga Nacional	1990–2000; 2006–2014	1990–2000; 2006–2014	1990–2000; 2006–2014
Vinhos Verdes	41°48′36″ N	8°24′36″ W	70	Touriga Franca	2005–2009	2005–2009	2005–2009
Vinhos Verdes	41°48′36″ N	8°24′36″ W	70	Touriga Nacional	2005–2009	2005–2009	2005–2009

¹ BBCH 07, ² BBCH 65, ³ BBCH 81.

): “Dois Portos—DP” (Lisboa Wine Region), with 20 years of data within the full period of 1990–2014, and “Estação Vitivinícola Amândio Galhano—EVAG” (Vinhos Verdes Wine Region), with records over the period of 2005–2009 (5 years). Therefore, the total sample size is extended to 73 records for each phenological timing and variety, thus enabling a higher statistical robustness in model calibration. A summary of the entire phenology dataset is provided in Table 1.

For the DDR, daily mean temperature records from the meteorological stations installed at the different vineyard sites and over the period of 2014–2017 (4 years) were used (Figure 1b). For the other two wine regions (Vinhos Verdes-EVAG and Lisboa-DP), the local temperature time series were estimated from the nearest-point of a very-high resolution (~1 km) gridded daily temperature dataset for Portugal [42]. Phenology and temperature data for the Douro Wine Region were provided by the “Associação para o Desenvolvimento da Viticultura Duriense—ADVID”, after collecting information among their members, while phenology data for Lisboa-DP and Vinhos Verdes-EVAG were provided by the “Instituto Nacional de Investigação Agrária e Veterinária”, INIAV, and by EVAG, respectively.

2.2. Phenological Models

Phenological models were selected based on previous studies, which have already shown moderate-to-high performance in simulating grapevine phenological timings in other regions [19,26,43]. Two main groups of phenological models were used: chilling models (Bidabe, Chuine, and Smoothed Utah), used to simulate the dormancy phase; and the forcing models (GDD, Richardson, Sigmoid, and Wang) applied to simulate thermal accumulation during the three phenophases: budburst, flowering, and veraison (

Table 2. Phenological models used in the present study, along with their relevant bibliographic references and pre-defined starting dates.

Type of Model	Model Name	Reference	Starting Date
Chilling model–Dormancy	Bidabe	[44]	August 1st of previous year
	Chuine	[45]
	Smoothed Utah (SU)	[46]
Forcing model, Post-Dormancy, Flowering, Veraison	Growing degree-days (GDD)	[47]	January 1st, previous phenological stage
	Richardson	[48]
	Sigmoid	[49]
	Wang	[50]

).

In the case of budburst, two approaches were tested: sequential models that allow integrating a dormancy phase (described using chilling models) and models without dormancy period (starting at a fixed date—January 1st—and only using forcing models). For sequential models, dormancy was herein calculated starting from August 1st of the previous year, as proposed by García de Cortázar-Atauri, Brisson and Gaudillere [19]. For flowering and veraison only the forcing models were applied. Table 2 describes the tested models and their relevant bibliographic references, while

Table 3. Grapevine phenology models selected in the present study, accompanied by their transfer function equations as a function of daily mean ($T_d$), maximum ($T_{max}$) and minimum ($T_{min}$) temperatures. The corresponding parameters are also outlined.

Model	Equation	Parameters *
Bidabe	$f_{Bidabe} = Q_{10}{}^{-T_{min}}+ Q_{10}{}^{-T_{max}}$	$Q_{10}\in \left[0, 5\right]$
Chuine	$f_{Chuine} = 1/(1 +\exp [a{\left(T_d -c\right)}^2 + b \left(T_d -c\right)])$	$a$$, b$$, c$
GDD	$f_{GDD} = Max \left(0;T_d - T_b\right)$	$T_b$
Richardson	$f_{Richardson} = Max\left(Min\left(T_d - T_{low}, T_{high}- T_{low}\right), 0\right)$	$T_{low}$, $T_{high}$
Sigmoid	$f_{Sigmoid} = 1/\left(1 +\exp \left[e\left(T_d -d\right)\right]\right)$	$e$$, d$
Smoothed Utah	$f_{Smoothed Utah} = \frac{1}{1+e^{-4\frac{T_d-T_{m1}}{T_{opt}-T_{m1}}}}$$; \mathrm{If} T_d$$> T_{m1}$; $\frac{0.5 {\left(T_d-T_{opt}\right)}^2}{{\left(T_{m1}-T_{opt}\right)}^2}$$; \mathrm{If} T_{m1}$$< T_d$$< T_{opt}$; $1-\left(1-min\right) \frac{{\left(T_d-T_{opt}\right)}^2}{2{\left(T_{n2}-T_{opt}\right)}^2}$$; \mathrm{If} T_{opt}$$< T_d$$< T_{n2}$; $min+ \frac{1-min}{1+e^{-4\frac{T_{n2}-T_d}{T_{n2}-T_{opt}}}}$$; \mathrm{If} T_{n2}$$< T_d$	$T_{m1}$,$T_{opt}$,$T_{n2}$, min
Wang	$f_{Wang} = Max[(2{(T_d- T_m)}^{\alpha } {(T_{opt}-T_m)}^{\alpha }-{\left( \frac{T_d-T_m}{T_{opt}-T_m}\right)}^{2\alpha }, 0]$; $\mathrm{with} \alpha = \ln \left(2\right)/\ln \left(\frac{T_M-T_m}{T_{opt}-T_m}\right)$	$T_m$$, T_M$$, T_{opt}$

* $T_{opt}$—optimum temperature; $T_b$—base temperature; $Q_{10}$—Bidabe parameter; a, b, c—Chuine constant parameters; $T_{low}$—Richardson lower plateau temperature; $T_{high}$—Richardson upper plateau temperature; e, d—Sigmoid constant parameters; $T_{m1}$, $T_{n2}$, min—Smoothed Utah parameters; $T_m$—Wang lower threshold temperature; $T_M$—Wang upper threshold temperature.

describes their corresponding mathematical formulations and tested parameters.

For all these models, the parameters were fixed at specific thresholds, some of them defined in the literature taking into account previous research: T_m = 0 °C, T_M = 40 °C and T_opt = 25–26 °C for Wang, T_low = 5 °C and T_high = 20–25 °C for Richardson, and T_b = 0 °C for GDD [19,20,26,43]. Besides these default models, the freely-adjusted best-fit parameters for each of these four models were also selected in the present study, as their performances can be significantly better than using the default parameters. The best eight performing models for each phenological phase were selected based on the performance metrics explained below. The models were considered for each of the three phenophases: dormancy to budburst (D-B), budburst to flowering (B-F), and flowering to veraison (F-V). The calibration of different models was performed independently for each phenophase and for the two selected grapevines varieties (TF and TN).

2.3. Modelling Tools and Performance Verification

In this study, the Phenological Modelling Platform, PMP [14], was used to test and calibrate different phenological models applied to grapevines. PMP is an intuitive platform that allows users to apply and test different phenological models for several purposes. In this platform, climatic and phenological data were used in order to calibrate the models, i.e., estimate the best-fit model parameters, for a specific location and a pre-defined time period. PMP estimates best-fit parameters following an iterative optimization procedure, based on the simulated annealing algorithm of Metropolis, et al. [51].

For each specific location, the input data used were daily mean (T_d), maximum (T_max) and minimum (T_min) temperatures, latitude and observed phenological timings for budburst, flowering, and veraison (in days of the year, DOY). The model performance is assessed herein by four commonly used metrics: the determination coefficient (R²), root-mean-square error (RMSE), Nash–Sutcliffe coefficient of efficiency (EF) and the Akaike Information Criterion (AIC). Their definitions are as follows:

$$RMSE = \sqrt{\frac{{\sum }_{i = 1}^n{\left(O_i- P_i\right)}^2}{n}} ,$$

(1)

$$EF = 1- \frac{{\sum }_{i = 1}^n{\left(O_i- P_i\right)}^2}{{\sum }_{i = 1}^n{\left(O_i- {\overline{O}}_i\right)}^2} ,$$

(2)

$$AIC = n\times \ln \left(\frac{{\sum }_{i = 1}^n{\left(O_i- P_i\right)}^2}{n}\right) +2k+\left(\frac{2k\left(k+1\right)}{n-k-1}\right)$$

(3)

where $O_i$ represents the observed phenological dates, ${\overline{O}}_i$ represents the average observed value, $P_i$ represents simulated phenological dates, $n$ is the sample size, and $k$ is the number of parameters.

EF varies between –∞ and +1, +1 corresponds to a perfect fit, 0 corresponds to the performance of the null model (average), while a negative value corresponds to a worse prediction than the null model. The $AIC$ takes into account the number of model parameters, with the lowest value corresponding to a model that represents the highest variance ratio in the observed dataset with the fewest parameters (parsimony principle). Further, in order to take into consideration model overfitting, a leave-one-out cross-validation scheme was applied to every simulation, following the same methodology as described by Chuine, et al. [52], and the RMSE metric was accordingly adapted to Root Mean Square Error of Prediction (RMSEP), defined as:

$$RMSEP = \sqrt{\frac{{\sum }_{i = 1}^n{\left(O_i-\widehat{P_i}\right)}^2}{n}}$$

(4)

where $\widehat{P_i}$ is the predicted value obtained from a leave-one-out cross-validation approach [53].

2.4. Future Climate Projections

PMP also allows running simulations using climatic data generated from different anthropogenic forcing scenarios and climate model experiments, applying the previously calibrated phenology models. For the simulations under future climates, the daily temperature time series generated by a two-member ensemble of climate model chains (

Table 4. List of the two EURO-CORDEX climate model chains selected for the present study. For each model chain the global climate model (GCM) and the regional climate model (RCM) are outlined, along with their abbreviation and relevant bibliographic reference.

GCM	RCM	Abbreviation
CNRM-CERFACS-CNRM-CM5	SMHI-RCA4	CNRMSMHI
MPI-M-MPI-ESM-LR	CLMcom-CCLM4-8-17	MPICLM

), produced within the framework of the EURO-CORDEX (Coordinated Downscaling Experiment—European Domain) project, were used [54]. A distribution-based scaling method was applied as a bias correction method by EURO-CORDEX [55]. The nearest-point time series for each vineyard site was extracted from the original grid (~12 km spacing). Data for the future period of 2020–2100, under the Representative Concentration Pathway 4.5 (RCP4.5), was selected. This scenario corresponds to an intermediate/moderate anthropogenic (greenhouse gas) forcing and corresponds to an increase of the global mean temperature slightly above 2 °C, thus marginally fulfilling the Paris Agreements.

Although the model calibration was undertaken using data from three wine regions, which enabled an improvement of the statistical robustness of the model fits, the future projections were only developed for the Douro wine region, as this study was never carried out previously. The outputs for the eight selected phenology models and for each phenological phase were applied to each climate model output, thus obtaining 12 sites × 3 phases × 8 phenology models × 2 climate models (576 simulations). In the future simulations, the D-B phase started at August 1st (with dormancy) or January 1st (without dormancy, Sigmoid model), as previously described. As in the observational period, for the other two phenophases (B-F and F-V) the starting date was set at the previous phenological stage (budburst or flowering), but now using the average time series of simulated budburst or flowering dates by the corresponding eight phenology models.

3. Results

3.1. Model Performance Verification

As previously mentioned, the sample size of the observed phenology dataset used for model calibration is of 73 records for each phenological timing and variety. For TF, the minimum/mean/maximum DOY is 61/77/98 for budburst (from March 2nd to April 8th), 112/138/158 for flowering (from April 22nd to June 7th) and 177/201/222 for veraison (from June 26th to August 10th). For TN, the observed minimum/mean/maximum DOY is 64/78/100 (from March 5th to April 10th) for budburst, 112/140/158 for flowering (from April 22nd to June 7th) and 170/213/235 for veraison (from June 19th to August 23rd). Despite the large spread in the dates of each phenophase, these results show that phenophase timings are commonly earlier for TF than for TN, particularly for veraison. In addition, the correlation coefficients between the time series for DP and EVAG, in their overlapping period (2005–2009), are very high (>0.80), thus showing a strong coherency in the inter-annual variability of the phenophase timings, despite the regional differences. No correlations can be computed with the DDR, as only 2014 overlaps with DP.

The best eight models for each corresponding phenophase are shown in

Table 5. List of the selected phenological models for cv. Touriga Franca and for each phenophase (D-B, B-F, and F-V), accompanied by the corresponding model performance parameters (Root Mean Square Error of Prediction (RMSEP), Nash–Sutcliffe coefficient of efficiency (EF), determination coefficient (R²), and Akaike Information Criterion (AIC)). Models are sorted according to their model performance based on AIC.

Phenophase	Model	Description	Model Performance
Dormancy-Budburst	With dormancy		RMSEP	EF	R2	AIC
	Bidabe + Wang	Q10 = 0.9, Tm = 5.5, TM = 34.3, Topt = 14.7	6.4	0.36	0.62	185.1
	Bidabe + Sigmoid	Q10 = 1.0, d = −22.4, e = 7.4	6.5	0.35	0.58	186.5
	SU + Sigmoid	Tm1 = −13.4, Topt = 29.2, Tn2 = 44.6, min = −0.5, d = −39.9, e = 7.5	6.2	0.41	0.64	187.4
	SU + GDD	Tm1 = −2.8, Topt = 29.2, Tn2 = 36.6, min = −0.3, Tb = 0	6.5	0.35	0.59	190.1
	SU + Wang	Tm1 = −32.2, Topt = 29.3, Tn2 = 35.7, min = −0.3, Tm = 5.3, TM = 34.8, Topt = 15.1	6.4	0.37	0.63	190.4
	Chuine + Sigmoid	a = 0.3, b = −21.1, c = −2.5, d = −40, e = 7.5	6.4	0.38	0.57	192.3
	Bidabe + Richardson	Q10 = 0.9, Tlow = 0.2, Thigh = 43.2	6.5	0.35	0.54	201.4
	Without dormancy
	Sigmoid	d = −40, e = 7.5	6.4	0.38	0.61	181.6
Budburst-Flowering	GDD	Tb = 6.6	5.1	0.69	0.85	159.1
	Sigmoid	d = −0.2, e = 16.9	5.1	0.70	0.86	159.4
	Richardson	Tlow = 6.5, Thigh = 36.9	5.1	0.69	0.85	161.1
	Wang	Tm = 0, TM = 40, Topt = 29.1	5.1	0.69	0.86	161.9
	Richardson	Tlow = 5, Thigh = 25	5.2	0.68	0.84	162.6
	Richardson	Tlow = 5, Thigh = 20	5.3	0.67	0.83	164.1
	GDD	Tb = 0	5.9	0.60	0.78	171.1
	Wang	Tm = 0, TM = 31.6, Topt = 25.5	5.6	0.62	0.87	172.3
Flowering-Veraison	Sigmoid	d = 0.0, e = 13.0	4.3	0.83	0.92	141.1
	Richardson	Tlow = 0.0, Thigh = 21.5	4.2	0.83	0.91	142.1
	Wang	Tm = 0, TM = 40, Topt = 25.3	4.2	0.84	0.92	142.6
	Wang	Tm = 0, TM = 40, Topt = 25	4.2	0.84	0.92	143.1
	Wang	Tm = 0, TM = 40, Topt = 26	4.2	0.83	0.91	144.2
	Richardson	Tlow = 5, Thigh = 20	4.3	0.82	0.91	144.9
	Wang	Tm = 0, TM = 36.6, Topt = 25.6	4.4	0.82	0.91	148.9
	GDD	Tb = 0	5.1	0.76	0.90	157.6

;

Table 6. The same as Table 5, but for cv. Touriga Nacional.

Phenophase	Model	Description	Model Performance
Dormancy-Budburst	With dormancy		RMSEP	EF	R2	AIC
	Bidabe + GDD	Q10 = 0.9, Tb = 0	7.8	0.18	0.43	190.5
	Bidabe + Wang	Q10 = 0.7, Tm = −36.7, TM = 34, Topt = 18	7.7	0.20	0.44	193.7
	Bidabe + Richardson	Q10 = 0.9, Tlow = 0, Thigh = 37.3	7.8	0.18	0.42	193.9
	SU + Sigmoid	Tm1 = −34.1, Topt = 29.9, Tn2 = 43, min = −0.2, d = −39.7, e = 7.5	7.6	0.23	0.48	196.0
	SU + GDD	Tm1 = −31.1, Topt = 29.5, Tn2 = 32.1, min = −0.1, Tb = 0	8.0	0.15	0.41	197.8
	Bidabe + Sigmoid	Q10 = 1, d = −38.3, e = 7.5	7.8	0.19	0.41	197.9
	SU + Wang	Tm1 = −31.1, Topt = 30, Tn2 = 32.2, min = −0.4, Tm = −16.1, TM = 25.1, Topt = 16.4	7.9	0.17	0.42	201.0
	Without dormancy
	Sigmoid	d = −40, e = 7.5	7.2	0.20	0.42	190.9
Budburst-Flowering	GDD	Tb = 7.0	5.4	0.68	0.83	169.6
	Sigmoid	d = −0.3, e = 14.4	5.4	0.69	0.84	170.5
	Richardson	Tlow = 7.0, Thigh = 20.8	5.4	0.69	0.83	171.1
	Wang	Tm = 0, TM = 40, Topt = 29.3	5.3	0.69	0.83	172.5
	Richardson	Tlow = 5, Thigh = 25	5.7	0.67	0.82	173.6
	Richardson	Tlow = 5, Thigh = 20	5.8	0.67	0.82	173.7
	Wang	Tm = 0, TM = 26.5, Topt = 21.8	5.3	0.66	0.84	177.3
	GDD	Tb = 0	6.8	0.61	0.78	180.1
Flowering-Veraison	Sigmoid	d = −0.4, e = 14.4	6.0	0.85	0.93	151.3
	Richardson	Tlow = 0, Thigh = 23.7	5.6	0.86	0.93	151.6
	GDD	Tb = 0	5.6	0.85	0.92	154.8
	Wang	Tm = 0, TM = 40, Topt = 26.9	5.7	0.86	0.93	155.7
	Wang	Tm = 0, TM = 40, Topt = 25	6.2	0.85	0.91	161.0
	Richardson	Tlow = 5, Thigh = 20	6.1	0.83	0.92	161.7
	Richardson	Tlow = 5, Thigh = 25	5.9	0.83	0.91	161.8
	Wang	Tm = 0, TM = 43.2, Topt = 29.9	5.9	0.83	0.91	164.9

for TF and TN, respectively. For the D-B phenophase, seven sequential models and one forcing model were chosen. The results show that for D-B, the RMSEP is of 6–7 days for TF, while it reaches values of almost 8 days for TN (Table 5 and Table 6). The EF varies from 0.35 to 0.41 (TF) or from 0.15 to 0.23 (TN). The corresponding coefficients of determination (R²) show values ranging from 0.54 to 0.64 (TF) or from 0.41 to 0.48 (TN). It is important to highlight that most of the models chosen present different sets of parameters, showing very different temperature responses. As the sample size of observed data does not allow choosing a specific model with high statistical confidence, it is important to maintain this diversity in order to cover a wide range of possibilities. Further, there is no clear difference between the experiments with or without dormancy simulation. In fact, the use of a sequential model does not seem to improve over the single-model even if the number of model parameters is higher in the sequential-model [19]. The AIC values hint at this conclusion, revealing higher values in the sequential models compared to the Sigmoid model. Overall, the model performances for D-B are thus moderate and the best performances are obtained from the Sigmoid model, for both TF and TN, or, alternatively, from the Bidabe + Wang (TF) or Bidabe + GDD (TN) sequential models.

In the case of B-F (Table 5 and Table 6), the RMSEP is of 5–6 days (TF and TN), i.e., less one day than for the D-B phase. The EF also presents much higher values than in the previous phenophase, between 0.60 to 0.70 in TF and from 0.61 to 0.69 in TN. The R² are also significantly higher, varying from 0.78 to 0.87 (0.78 to 0.84) in TF (TN). The AIC values are within the range 159–172 (170–180) in TF (TN). Among the three phenophases, the best results are achieved for F-V (Table 5 and Table 6), with RMSEP of 4–5 days in TF and 5–6 days in TN. The EF varies from 0.76 to 0.84 in TF and from 0.83 to 0.86 in TN. The AIC shows relatively low values, from 141 to 158, in TF, and from 151 to 165, in TN.

The errors between simulated and observed dates for TF, over the period from 1995 to 2017, show relatively high values for D-B (up to ca. 20 days), while they are mostly lower than 10 days for B-F and F-V (Figure 2). As expected, the errors are much more scattered for D-B than for the other two phenophases, with an error compensation between the periods of 1998–2003 and 2005–2009. There is also some non-stationarity in errors, which may be attributed to non-homogeneities in phenophase and/or climate data. Similar considerations can be made for TN (Figure S1).

3.2. Future Projections

The simulated temperatures from the ensemble of two climate models under RCP4.5 were performed for 8 phenological models and for 12 vineyard sites in the DDR. In order to avoid an excessively large number of figures, only the mean curves over the 12 sites in the DDR are presented in Figure 3; Figure 4. The corresponding panels for each climate model chain are presented separately in Figures S2–S5. Apart from minor differences, site-specific figures are very similar (not shown). For all models, a gradual anticipation of the three main phenological phases (downward long-term trends) is clear for both TF (Figure 3) and TN (Figure 4), despite not being linear and revealing noteworthy fluctuations. Analysing the chronograms for both TF and TN, it is evident that the eight selected models for each phenophase render coherent results (low inter-model ranges), mostly for flowering and veraison. In general, the inter-model range is of nearly 4–6 days for budburst and 2–4 days for flowering and veraison. The correlations between D-B phase and January–April mean temperature, between B-F phase and May–June mean temperature and between F-V phase and July–August mean temperature are clear, thus highlighting that the projected changes are largely driven by changes in mean temperatures over specific time periods. However, some differences between the curves of temperature and mean phenophase timings are found. This can be explained by the fact that the temperature curves correspond to multi-month averages, while the phenological timings are determined by 1) the accumulation of daily temperatures over longer time periods and 2) by the preceding phenophase timings. Furthermore, apart from GDD, model responses are non-linear functions of temperature.

In the case of D-B for TF, a mean anticipation of 6 days is expected from 2020 to 2100, with important inter-decadal variability, without a clear trend until 2070, followed by a more pronounced long-term decrease from 2070 to 2100 (Figure 3a). The Smoothed-Wang and the Smoothed-Sigmoid are the closest models to the mean. For B-F in TF, a clear anticipation is expected after 2040, despite the delay in the period of 2060–2070 (Figure 3b), which can be attributed to a cooler period in the future climate (Figure 3b). Overall, a mean anticipation of 8 days in flowering is projected until 2100. The Sigmoid (d = −0.2, e = 16.9) and the Wang (T_m = 0 °C, T_M = 40 °C, T_opt = 29.1 °C) models show the closest correspondence with the mean. For F-V, an anticipation of 8–10 days over the full period is projected (Figure 3c), with analogous temporal variability to B-F (Figure 3b). The Richardson (T_low = 0 °C, T_high = 21.5 °C) and the Wang (T_m = 0 °C, T_M = 40 °C, T_opt = 26 °C) models are very close to the mean.

Similar considerations can be made for TN (Figure 4). For D-B of TN, the Smoothed-GDD and Smoothed-Wang are the closest to the mean. For B-F of TN, the Richardson (T_low = 7 °C, T_high = 20.8 °C) and the Sigmoid (d = −0.3, e = 14.4) models show the most central behaviour. For F-V of TN, the Richardson (T_low = 0 °C, T_high = 23.7 °C) and the Wang (T_m = 0 °C, T_M = 40 °C, T_opt = 26.9 °C) models are the closest to the mean. In addition, the anticipation of veraison for TN is more pronounced than for TF (8 to 12 days). These anticipations in phenological timings are in clear agreement with other studies worldwide [21,32,34]. These shifts towards earlier timings may result in changes to the currently established cultural practices and potentially affect the wine characteristics of the Douro.

3.3. Inter-Model Spread

Regarding the inter-model variability, some differences can be highlighted. For D-B of TF (Figure 3a), all the selected models show very similar results, with the exception of the sequential models Bidabe-Richardson and Smoothed-GDD, which reveal some significant departures from the model average. These two models simulate systematically earlier budburst timings than the other models (Figure 3a), despite their strong agreement, as expected taking into account that their transfer functions are very similar (Table 5). Further, the Sigmoid, Bidabe-Sigmoid, Smoothed-Sigmoid, and Chuine-Sigmoid models are in close agreement, as the second model is the same, also suggesting that dormancy does not play a key role in the final outcome. Similarly, the Bidabe-Wang and Smoothed-Wang also reveal similar behaviour throughout the future period. For B-F of TF (Figure 3b), the GDD (Tb = 0 °C) model shows a tendency to have later timings. Conversely, the Wang (T_m = 0 °C, T_M = 31.6 °C, T_opt = 25.5 °C) model reveal earlier timings throughout the analysed period. For F-V of TF (Figure 3c), the selected models highlight a very strong coherency, apart from some light discrepancies found for GDD (T_b = 0 °C), with earlier timings.

For D-B in TN (Figure 4a), the Bidabe-Wang presents the earliest phenological timings throughout the study time period. For B-F in TN (Figure 4b), the GDD (Tb = 0 °C) model presents later phenological timings, while the Wang (T_m = 0 °C, T_M = 26.5 °C, T_opt = 21.8 °C) model shows earlier phenological timings along all the study period. The other models show very strong agreement between them and with values closer to the mean. However, it should be stated that models showing higher agreement are not necessarily more realistic in reproducing timings under future conditions. Lastly, regarding F-V in TN (Figure 4c), the Wang (T_m = 0 °C, T_M = 40 °C, T_opt = 25 °C) shows the later timings, while the GDD (Tb = 0 °C) shows the earlier timings.

4. Discussion and Conclusions

Phenological models are a key tool for viticulturists, as they provide valuable information for management and planning in vineyards [7]. The accurate prediction phenological timings promotes best practices and the timely implementation of suitable measures to optimize grapevine yields and grape berry quality [56,57,58,59]. Moreover, future projections of potential changes in the phenological timings may also deliver important information for medium-to-long term planning [60]. Despite the aforementioned high economic value of the DDR, no previous studies were conducted for this region on phenological modelling and simulation. Unique datasets, combining phenological timings with weather station records, in the Douro wine region were used herein. The weather stations used were installed at the vineyard locations, which is of utmost importance for a proper phenological modelling [61]. Furthermore, two of the most important grapevine varieties cultivated in the Douro wine region (TF and TN) were selected (cv. Touriga Franca—TF and cv. Touriga Nacional—TN). Despite the relatively large number of sites in the Douro valley (12 Quintas), data is only available over a maximum period of 4 years. This dataset was thereby complemented with data from two other Portuguese wine regions (Lisboa and Vinhos Verdes) and for the same grapevine varieties (TF and TN). This approach led to more robust results, as the sample sizes for model calibration were significantly enlarged.

It was herein demonstrated that some specific phenological models are particularly suitable to the simulation of the two selected grapevine varieties. For all phenophases, specific parameter values were chosen within each model (e.g., Tb in GDD), either by selecting reference values, recommended from previous studies, or best-fit values. Although the recommended values did not always lead to the best performances in our study, no definite conclusions can be withdrawn, as only two varieties were selected herein and the corresponding sample sizes are relatively low. As a whole, the present study allowed the identification of the best phenological models (a set of eight models per phenophase), thus warranting their application both in an operational mode (real-time monitoring and short-term prediction) and in future climate change projections (long-term prediction). Nevertheless, more data will allow improving the selection of the best model and the accuracy of the predictions.

In general terms, the results highlight that the phenological models are better in simulating flowering or veraison than budburst. This may be due to the fact that dormancy and other relevant physiological processes are not sufficiently understood, but might be important drivers of budburst occurrence. As also described by García de Cortázar-Atauri, Brisson and Gaudillere [19], no clear improvements were found concerning the incorporation of a dormancy period in the model structure. Moreover, the model performances are generally higher for TF than for TN, which may suggest that the latter is slightly less sensitive to thermal forcing, though this needs to be properly addressed by further research, including field experiments. Overall, the Bidabe-Wang and Sigmoid models showed the best results for budburst in TF, while the Bidabe-GDD and Sigmoid models present the highest performances in TN. For flowering, in both grapevine varieties (TF and TN), the very simple GDD model, with a base temperature of ca. 7 °C, revealed the best performances. Lastly, the Sigmoid model provided the best performances for veraison in both varieties. Hence, the GDD linear model is not always the best approach to model and simulate grapevine phenology, as the phenology–temperature relationship may be of non-linear nature. The present study results are comparable, in terms of model performance, to other studies in other winemaking regions worldwide [19,26,27,43,62,63,64,65]. Although the aforementioned models provide the best approaches, it is important to bear in mind that they are bound to this particular phenology dataset.

Climate change projections of the phenophase timings, based on a two-member ensemble of state-of-the-art climate model chains, hinted at anticipations of the phenophase timings up to 10–12 days until the end of the twentieth first century, though with some dissimilarities amongst phenophases and varieties. In addition, there is a strong coherency between phenological models, as the results do not change significantly with the model choice, apart from a few exceptions, thus highlighting the robustness of the projections. Nonetheless, different grapevines varieties exhibit greater phenological sensitivity to temperature changes than others, which limits the extrapolation of the results found herein to other locations and varieties. These projections for the Douro are in agreement with observed trends worldwide, such as for Alsace (France), [Jones [66]] for the California (USA) and Bock, et al. [67] for Franconia (Germany). These shifts towards earlier phenophase onsets can potentially result in changes to the currently established wine characteristics and typicity. In effect, the expected warming may result in unbalanced wines, with high alcoholic content and excessively low acidity and altered colour and aroma [68,69]. These impacts can be partially overcome through the adoption of suitable adaptation measures, such as selection of grapevine varieties that tend to produce berries with lower sugar contents and higher acidity levels in warm climates. Suitable selections of variety-clone-rootstock combinations are among the most effective medium-to-long term measures to cope with climate change. The cultivation of late maturation cultivars is also key to mitigate the anticipation of phenological stages in the vineyards. Short term measures, like irrigation, application of sun screens, changes in training systems, soil management and mulching application, may also be envisioned by producers and winemakers.

The present study aims at providing some initial insights that can be of use not only to future research, but also to regional stakeholders. In the near future, it is expected that the outcomes of this study will be implemented on a web platform, maintained by ADVID, which will deliver information on phenological timings to regional wine companies and grounded on a network of operational weather stations. Future research should also envision the improvement of the current models, such as by extending the phenological datasets with new data from other regions in Portugal or in other countries where the selected varieties also grow. Other grapevine varieties may also be considered in future studies, but there is a critical lack of information for them, particularly on the Douro native/autochthonous varieties, which may considerably hamper their simulation.