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

Ricardo Costa, Helder Fraga, André Fonseca, Inaki Garcia de Cortazar

Atauri, Maria C. Val, Cristina Carlos, Samuel Reis, João A. Santos

► To cite this version:

Ricardo Costa, Helder Fraga, André Fonseca, Inaki Garcia de Cortazar Atauri, Maria C. Val, et al. Grapevine Phenology of cv. Touriga Franca and Touriga Nacional in the Douro Wine Region: Modelling and Climate Change Projections. Agronomy, 2019, 9 (4), pp.210. 10.3390/agronomy9040210. hal-02623900

HAL Id: hal-02623900

https://hal.inrae.fr/hal-02623900

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

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

Ricardo Costa <fnref n="1" />, Helder Fraga <fnref n="1" /> ID, André Fonseca <fnref n="1" />, Iñaki García de Cortázar-Atauri <fnref n="2" /> ID,
Maria C. Val <fnref n="3" />, Cristina Carlos <fnref n="3" />, Samuel Reis <fnref n="1" /> and João A. Santos <fnref n="1" />,* ID

Received: 3 April 2019; Accepted: 23 April 2019; Published: 25 April 2019

<figure><img src="image_2.png" /></figure>

check for
updates

**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^2 < 0.7$) for budburst, while high performances ($R^2 > 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.

**Keywords:** grapevine; phenology; phenology modelling platform; *Touriga Franca*; *Touriga Nacional*; climate change; RCP4.5; EURO-CORDEX; Douro wine region; Portugal

## 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 \times 10^6$ 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

<footnote n="1">¹ Centre for the Research and Technology of Agro-Environmental and Biological Sciences, CITAB, Universidade de Trás-os-Montes e Alto Douro, UTAD, 5000-801 Vila Real, Portugal; rjncosta@utad.pt (R.C.); hfraga@utad.pt (H.F.); andref@utad.pt (A.F.); sreis_1992@hotmail.com (S.R.)</footnote>

<footnote n="2">² INRA, US1116 AgroClim, 84914 Avignon, France; inaki.garciadecortazar@inra.fr</footnote>

<footnote n="3">³ ADVID, Associação para o Desenvolvimento da Viticultura Duriense, Parque de Ciência e Tecnologia de Vila Real-Régia Douro Park, 5000-033 Vila Real, Portugal; carmo.val@advid.pt (M.C.V.); cristina.carlos@advid.pt (C.C.)</footnote>

<footnote>* Correspondence: jsantos@utad.pt; Tel.: +351-259350729</footnote>
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–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–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–33]. Projected increases in temperature are expected to drive earlier development stages and, as a result, a general advancement of grapevine phenology [24,34–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): "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.

<table><caption>**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).</caption><thead><tr><th rowspan="2">Wine Region</th><th rowspan="2">Latitude (°N)</th><th rowspan="2">Longitude (°W)</th><th rowspan="2">Altitude (meters)</th><th rowspan="2">Variety (cv.)</th><th colspan="3">Available Time Period</th></tr><tr><th>Budburst <sup>1</sup></th><th>Flowering <sup>2</sup></th><th>Veraison <sup>3</sup></th></tr></thead><tbody><tr><td>Douro/Porto</td><td>41°14'24" N-41°02'24" N</td><td>7°01'48" W-7°47'24" W</td><td>85-588</td><td>Touriga Franca</td><td>2014-2017</td><td>2014-2017</td><td>2014-2017</td></tr><tr><td>Douro/Porto</td><td>41°14'24" N-41°02'24" N</td><td>7°01'48" W-7°47'24" W</td><td>85-588</td><td>Touriga Nacional</td><td>2014-2017</td><td>2014-2017</td><td>2014-2017</td></tr><tr><td>Lisboa</td><td>39°02'24" N</td><td>9°10'48" W</td><td>85</td><td>Touriga Franca</td><td>1995-2014</td><td>1995-2014</td><td>1995-2014</td></tr><tr><td>Lisboa</td><td>39°02'24" N</td><td>9°10'48" W</td><td>85</td><td>Touriga Nacional</td><td>1990-2000;<br>2006-2014</td><td>1990-2000;<br>2006-2014</td><td>1990-2000;<br>2006-2014</td></tr><tr><td>Vinhos Verdes</td><td>41°48'36" N</td><td>8°24'36" W</td><td>70</td><td>Touriga Franca</td><td>2005-2009</td><td>2005-2009</td><td>2005-2009</td></tr><tr><td>Vinhos Verdes</td><td>41°48'36" N</td><td>8°24'36" W</td><td>70</td><td>Touriga Nacional</td><td>2005-2009</td><td>2005-2009</td><td>2005-2009</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Map of mainland Portugal with wine regions' boundaries. The three wine regions providing phenological data to the present study are grey shaded (Vinhos-Verdes, Lisboa, and Douro/Porto). The location of the vineyard sites EVAG, in Vinhos-Verdes wine region, and DP, in Lisboa wine region, are also outlined. (b) Douro/Porto wine region and its sub-regions (Baixo-Corgo, Cima-Corgo, and Douro Superior), along with the geographical locations of the 12 vineyard sites. The main rivers are also plotted, including the Douro River.</figcaption></figure>

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

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,

<footnote n="1">¹ BBCH 07, ² BBCH 65, ³ BBCH 81.</footnote>
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 describes their corresponding mathematical formulations and tested parameters.

<table><caption>**Table 2.** Phenological models used in the present study, along with their relevant bibliographic references and pre-defined starting dates.</caption><thead><tr><th>Type of Model</th><th>Model Name</th><th>Reference</th><th>Starting Date</th></tr></thead><tbody><tr><td rowspan="3">Chilling model-Dormancy</td><td>Bidabe</td><td>[44]</td><td rowspan="3">August 1st of previous year</td></tr><tr><td>Chuine</td><td>[45]</td></tr><tr><td>Smoothed Utah (SU)</td><td>[46]</td></tr><tr><td rowspan="4">Forcing model, Post-Dormancy, Flowering, Veraison</td><td>Growing degree-days (GDD)</td><td>[47]</td><td rowspan="4">January 1st, previous phenological stage</td></tr><tr><td>Richardson</td><td>[48]</td></tr><tr><td>Sigmoid</td><td>[49]</td></tr><tr><td>Wang</td><td>[50]</td></tr></tbody></table>

<table><caption>**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.</caption><thead><tr><th>Model</th><th>Equation</th><th>Parameters *</th></tr></thead><tbody><tr><td>Bidabe</td><td>$f_{Bidabe} = Q_{10}^{-T_{min}} + Q_{10}^{-T_{max}}$</td><td>$Q_{10} \in [0, 5]$</td></tr><tr><td>Chuine</td><td>$f_{Chuine} = 1/(1 + \exp[a(T_d - c)^2 + b(T_d - c)])$</td><td>$a, b, c$</td></tr><tr><td>GDD</td><td>$f_{GDD} = Max(0; T_d - T_b)$</td><td>$T_b$</td></tr><tr><td>Richardson</td><td>$f_{Richardson} = Max(\text{Min}(T_d - T_{low}, T_{high} - T_{low}), 0)$</td><td>$T_{low}, T_{high}$</td></tr><tr><td>Sigmoid</td><td>$f_{Sigmoid} = 1/(1 + \exp[e(T_d - d)])$</td><td>$e, d$</td></tr><tr><td>Smoothed Utah</td><td>$f_{Smoothed\ Utah} = \frac{1}{1+e^{-\frac{T_d-T_{m1}}{T_{opt}-T_{m1}}}$; If $T_d > T_{m1}$;<br>$\frac{0.5(T_d-T_{opt})^2}{(T_{m1}-T_{opt})^2}$; If $T_{m1} < T_d < T_{opt}$;<br>$1-(1-min)\frac{(T_d-T_{opt})^2}{2(T_{n2}-T_{opt})^2}$; If $T_{opt} < T_d < T_{n2}$;<br>$min + \frac{1-min}{1+e^{-\frac{T_{n2}-T_d}{T_{n2}-T_{opt}}}}$; If $T_{n2} < T_d$</td><td>$T_{m1}, T_{opt}, T_{n2}, min$</td></tr><tr><td>Wang</td><td>$f_{Wang} = Max[(2(T_d - T_m)^\alpha (T_{opt} - T_m)^\alpha - (\frac{T_d-T_m}{T_{opt}-T_m})^{2\alpha}], 0]$;<br>with $\alpha = \ln(2)/\ln(\frac{T_M-T_m}{T_{opt}-T_M})$</td><td>$T_m, T_M, T_{opt}$</td></tr></tbody></table>

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

<footnote>* $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.</footnote>
## 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^2$), 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}(O_i - P_i)^2}{n}}, \qquad (1)$$

$$EF = 1 - \frac{\sum_{i=1}^{n}(O_i - P_i)^2}{\sum_{i=1}^{n}(O_i - \bar{O}_i)^2}, \qquad (2)$$

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

where $O_i$ represents the observed phenological dates, $\bar{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}(O_i - \hat{P}_i)^2}{n}} \qquad (4)$$

where $\hat{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), 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.

<table><caption>**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.</caption><thead><tr><td>GCM</td><td>RCM</td><td>Abbreviation</td></tr></thead><tbody><tr><td>CNRM-CERFACS-CNRM-CM5</td><td>SMHI-RCA4</td><td>CNRMSMHI</td></tr><tr><td>MPI-M-MPI-ESM-LR</td><td>CLMcom-CCLM4-8-17</td><td>MPICLM</td></tr></tbody></table>

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; Table 6 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 (Tables 5 and 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^2$) 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 (Tables 5 and 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^2$ 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 (Tables 5 and 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).
<table><caption>**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.</caption><thead><tr><th rowspan="2">Phenophase</th><th rowspan="2">Model</th><th rowspan="2">Description</th><th colspan="4">Model Performance</th></tr><tr><th>RMSEP</th><th>EF</th><th>R<sup>2</sup></th><th>AIC</th></tr></thead><tbody><tr><td rowspan="8">Dormancy-Budburst</td><td></td><td>With dormancy</td><td></td><td></td><td></td><td></td></tr><tr><td>Bidabe + Wang</td><td>Q<sub>10</sub> = 0.9, T<sub>m</sub> = 5.5, T<sub>M</sub> = 34.3, T<sub>opt</sub> = 14.7</td><td>6.4</td><td>0.36</td><td>0.62</td><td>185.1</td></tr><tr><td>Bidabe + Sigmoid</td><td>Q<sub>10</sub> = 1.0, d = -22.4, e = 7.4</td><td>6.5</td><td>0.35</td><td>0.58</td><td>186.5</td></tr><tr><td>SU + Sigmoid</td><td>T<sub>m1</sub> = -13.4, T<sub>opt</sub> = 29.2, T<sub>n2</sub> = 44.6, min = -0.5, d = -39.9, e = 7.5</td><td>6.2</td><td>0.41</td><td>0.64</td><td>187.4</td></tr><tr><td>SU + GDD</td><td>T<sub>m1</sub> = -2.8, T<sub>opt</sub> = 29.2, T<sub>n2</sub> = 36.6, min = -0.3, T<sub>b</sub> = 0</td><td>6.5</td><td>0.35</td><td>0.59</td><td>190.1</td></tr><tr><td>SU + Wang</td><td>T<sub>m1</sub> = -32.2, T<sub>opt</sub> = 29.3, T<sub>n2</sub> = 35.7, min = -0.3, T<sub>m</sub> = 5.3, T<sub>M</sub> = 34.8, T<sub>opt</sub> = 15.1</td><td>6.4</td><td>0.37</td><td>0.63</td><td>190.4</td></tr><tr><td>Chuine + Sigmoid</td><td>a = 0.3, b = -21.1, c = -2.5, d = -40, e = 7.5</td><td>6.4</td><td>0.38</td><td>0.57</td><td>192.3</td></tr><tr><td>Bidabe + Richardson</td><td>Q<sub>10</sub> = 0.9, T<sub>low</sub> = 0.2, T<sub>high</sub> = 43.2</td><td>6.5</td><td>0.35</td><td>0.54</td><td>201.4</td></tr><tr><td rowspan="10">Budburst-Flowering</td><td></td><td>Without dormancy</td><td></td><td></td><td></td><td></td></tr><tr><td>Sigmoid</td><td>d = -40, e = 7.5</td><td>6.4</td><td>0.38</td><td>0.61</td><td>181.6</td></tr><tr><td>GDD</td><td>T<sub>b</sub> = 6.6</td><td>5.1</td><td>0.69</td><td>0.85</td><td>159.1</td></tr><tr><td>Sigmoid</td><td>d = -0.2, e = 16.9</td><td>5.1</td><td>0.70</td><td>0.86</td><td>159.4</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 6.5, T<sub>high</sub> = 36.9</td><td>5.1</td><td>0.69</td><td>0.85</td><td>161.1</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 40, T<sub>opt</sub> = 29.1</td><td>5.1</td><td>0.69</td><td>0.86</td><td>161.9</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 5, T<sub>high</sub> = 25</td><td>5.2</td><td>0.68</td><td>0.84</td><td>162.6</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 5, T<sub>high</sub> = 20</td><td>5.3</td><td>0.67</td><td>0.83</td><td>164.1</td></tr><tr><td>GDD</td><td>T<sub>b</sub> = 0</td><td>5.9</td><td>0.60</td><td>0.78</td><td>171.1</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 31.6, T<sub>opt</sub> = 25.5</td><td>5.6</td><td>0.62</td><td>0.87</td><td>172.3</td></tr><tr><td rowspan="8">Flowering-Veraison</td><td>Sigmoid</td><td>d = 0.0, e = 13.0</td><td>4.3</td><td>0.83</td><td>0.92</td><td>141.1</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 0.0, T<sub>high</sub> = 21.5</td><td>4.2</td><td>0.83</td><td>0.91</td><td>142.1</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 40, T<sub>opt</sub> = 25.3</td><td>4.2</td><td>0.84</td><td>0.92</td><td>142.6</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 40, T<sub>opt</sub> = 25</td><td>4.2</td><td>0.84</td><td>0.92</td><td>143.1</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 40, T<sub>opt</sub> = 26</td><td>4.2</td><td>0.83</td><td>0.91</td><td>144.2</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 5, T<sub>high</sub> = 20</td><td>4.3</td><td>0.82</td><td>0.91</td><td>144.9</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 36.6, T<sub>opt</sub> = 25.6</td><td>4.4</td><td>0.82</td><td>0.91</td><td>148.9</td></tr><tr><td>GDD</td><td>T<sub>b</sub> = 0</td><td>5.1</td><td>0.76</td><td>0.90</td><td>157.6</td></tr></tbody></table>
<table><caption>**Table 6.** The same as Table 5, but for cv. Touriga Nacional.</caption><thead><tr><th rowspan="2">Phenophase</th><th rowspan="2">Model</th><th rowspan="2">Description</th><th colspan="4">Model Performance</th></tr><tr><th>RMSEP</th><th>EF</th><th>R<sup>2</sup></th><th>AIC</th></tr></thead><tbody><tr><td rowspan="10">Dormancy-Budburst</td><td></td><td>With dormancy</td><td></td><td></td><td></td><td></td></tr><tr><td>Bidabe + GDD</td><td>Q<sub>10</sub> = 0.9, T<sub>b</sub> = 0</td><td>7.8</td><td>0.18</td><td>0.43</td><td>190.5</td></tr><tr><td>Bidabe + Wang</td><td>Q<sub>10</sub> = 0.7, T<sub>m</sub> = -36.7, T<sub>M</sub> = 34, T<sub>opt</sub> = 18</td><td>7.7</td><td>0.20</td><td>0.44</td><td>193.7</td></tr><tr><td>Bidabe + Richardson</td><td>Q<sub>10</sub> = 0.9, T<sub>low</sub> = 0, T<sub>high</sub> = 37.3</td><td>7.8</td><td>0.18</td><td>0.42</td><td>193.9</td></tr><tr><td>SU + Sigmoid</td><td>T<sub>m1</sub> = -34.1, T<sub>opt</sub> = 29.9, T<sub>n2</sub> = 43, min = -0.2, d = -39.7, e = 7.5</td><td>7.6</td><td>0.23</td><td>0.48</td><td>196.0</td></tr><tr><td>SU + GDD</td><td>T<sub>m1</sub> = -31.1, T<sub>opt</sub> = 29.5, T<sub>n2</sub> = 32.1, min = -0.1, T<sub>b</sub> = 0</td><td>8.0</td><td>0.15</td><td>0.41</td><td>197.8</td></tr><tr><td>Bidabe + Sigmoid</td><td>Q<sub>10</sub> = 1, d = -38.3, e = 7.5</td><td>7.8</td><td>0.19</td><td>0.41</td><td>197.9</td></tr><tr><td>SU + Wang</td><td>T<sub>m1</sub> = -31.1, T<sub>opt</sub> = 30, T<sub>n2</sub> = 32.2, min = -0.4, T<sub>m</sub> = -16.1, T<sub>M</sub> = 25.1, T<sub>opt</sub> = 16.4</td><td>7.9</td><td>0.17</td><td>0.42</td><td>201.0</td></tr><tr><td></td><td>Without dormancy</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="10">Budburst-Flowering</td><td>Sigmoid</td><td>d = -40, e = 7.5</td><td>7.2</td><td>0.20</td><td>0.42</td><td>190.9</td></tr><tr><td>GDD</td><td>T<sub>b</sub> = 7.0</td><td>5.4</td><td>0.68</td><td>0.83</td><td>169.6</td></tr><tr><td>Sigmoid</td><td>d = -0.3, e = 14.4</td><td>5.4</td><td>0.69</td><td>0.84</td><td>170.5</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 7.0, T<sub>high</sub> = 20.8</td><td>5.4</td><td>0.69</td><td>0.83</td><td>171.1</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 40, T<sub>opt</sub> = 29.3</td><td>5.3</td><td>0.69</td><td>0.83</td><td>172.5</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 5, T<sub>high</sub> = 25</td><td>5.7</td><td>0.67</td><td>0.82</td><td>173.6</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 5, T<sub>high</sub> = 20</td><td>5.8</td><td>0.67</td><td>0.82</td><td>173.7</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 26.5, T<sub>opt</sub> = 21.8</td><td>5.3</td><td>0.66</td><td>0.84</td><td>177.3</td></tr><tr><td>GDD</td><td>T<sub>b</sub> = 0</td><td>6.8</td><td>0.61</td><td>0.78</td><td>180.1</td></tr><tr><td rowspan="10">Flowering-Veraison</td><td>Sigmoid</td><td>d = -0.4, e = 14.4</td><td>6.0</td><td>0.85</td><td>0.93</td><td>151.3</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 0, T<sub>high</sub> = 23.7</td><td>5.6</td><td>0.86</td><td>0.93</td><td>151.6</td></tr><tr><td>GDD</td><td>T<sub>b</sub> = 0</td><td>5.6</td><td>0.85</td><td>0.92</td><td>154.8</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 40, T<sub>opt</sub> = 26.9</td><td>5.7</td><td>0.86</td><td>0.93</td><td>155.7</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 40, T<sub>opt</sub> = 25</td><td>6.2</td><td>0.85</td><td>0.91</td><td>161.0</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 5, T<sub>high</sub> = 20</td><td>6.1</td><td>0.83</td><td>0.92</td><td>161.7</td></tr><tr><td>Richardson</td><td>T<sub>low</sub> = 5, T<sub>high</sub> = 25</td><td>5.9</td><td>0.83</td><td>0.91</td><td>161.8</td></tr><tr><td>Wang</td><td>T<sub>m</sub> = 0, T<sub>M</sub> = 43.2, T<sub>opt</sub> = 29.9</td><td>5.9</td><td>0.83</td><td>0.91</td><td>164.9</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Chronograms of the departures between simulated and observed phenological timings of each phenophase on cv. *Touriga Franca*: (a) Dormancy–Budburst (D-B); (b) Budburst–Flowering (B-F), and (c) Flowering–Veraison (F-V), obtained from the eight selected models (cf. legends) and assembled for the 14 vineyard sites (Figure 1).</figcaption></figure>

## 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 \,^\circ C, T_M = 40 \,^\circ C, T_{opt} = 29.1 \,^\circ 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 \,^\circ C, T_{high} = 21.5 \,^\circ C$) and the Wang ($T_m = 0 \,^\circ C, T_M = 40 \,^\circ C, T_{opt} = 26 \,^\circ 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 \,^\circ C, T_{high} = 20.8 \,^\circ 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 \,^\circ C, T_{high} = 23.7 \,^\circ C$) and the Wang ($T_m = 0 \,^\circ C, T_M = 40 \,^\circ C, T_{opt} = 26.9 \,^\circ 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 ($T_b = 0 \,^\circ C$) model shows a tendency to have later timings. Conversely, the Wang ($T_m = 0 \,^\circ C, T_M = 31.6 \,^\circ C, T_{opt} = 25.5 \,^\circ 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 \,^\circ 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 ($T_b = 0 \,^\circ C$) model presents later phenological timings, while the Wang ($T_m = 0 \,^\circ C, T_M = 26.5 \,^\circ C, T_{opt} = 21.8 \,^\circ 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 ($T_b = 0$ °C) shows the earlier timings.

<figure><img src="image_5.png" /></figure>

**Figure 3.** Chronograms (2020–2100, RCP4.5) of the phenological timings (days of the year, DOY) of (a) Dormancy–Budburst (D-B); (b) Budburst–Flowering (B-F) and (c) Flowering–Veraison (F-V), and for cv. Touriga Franca (TF), obtained from the eight selected models (cf. legends), averaged for all climate models and the 12 sites in the Douro Demarcated Region (DDR) (Figure 1). The average maximum–minimum ranges for each phenophase (among the different phenology models) are also pointed out (values in days within boxes in the right-upper corner of each panel). The red curves correspond to the (a) January–April, (b) May–June, and (c) July–August mean temperatures (temperature scale inverted).
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Chronograms (2020–2100, RCP4.5) of the phenological timings (days of the year, DOY) of (a) Dormancy–Budburst (D-B); (b) Budburst–Flowering (B-F) and (c) Flowering–Veraison (F-V), and for cv. *Touriga Nacional* (TF), obtained from the eight selected models (cf. legends), averaged for all climate models and the 12 sites in the Douro Demarcated Region (DDR) (Figure 1). The average maximum–minimum ranges for each phenophase (among the different phenology models) are also pointed out (values in days within boxes in the right-upper corner of each panel). The red curves correspond to the (a) January–April, (b) May–June, and (c) July–August mean temperatures (temperature scale inverted).</figcaption></figure>
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–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–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.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/9/4/210/s1, Figure S1: Chronograms of the departures between simulated and observed phenological timings of each phenophase on cv. *Touriga Nacional*: (a) Dormancy–Budburst (D-B); (b) Budburst–Flowering (B-F), and (c) Flowering–Veraison (F-V), obtained from the eight selected models (cf. legends) and assembled for the 14 vineyard sites (Figure 1). Figure S2. Chronograms (2020–2100, RCP4.5) of the phenological timings (days of the year, DOY) of (a) Dormancy–Budburst (D-B); (b) Budburst–Flowering (B-F); and (c) Flowering–Veraison (F-V), and for cv. *Touriga Franca* (TF), obtained from the eight selected models (cf. legends), averaged for the 12 sites in the DDR (Figure 1) and for the climate model CNRM-SHMI. The average maximum-minimum ranges for each phenophase (among the different phenology models) are also pointed out by vertical arrows. Figure S3. Chronograms (2020–2100, RCP4.5) of the phenological timings (days of the year, DOY) of (a) Dormancy–Budburst (D-B); (b) Budburst–Flowering (B-F); and (c) Flowering–Veraison (F-V), and for cv. *Touriga Nacional* (TN), obtained from the eight selected models (cf. legends), averaged for the 12 sites in the DDR (Figure 1) and for the climate model CNRM-SHMI. The average maximum-minimum ranges for each phenophase (among the different phenology models) are also pointed out by vertical arrows. Figure S4. As for Figure S2, but for the climate model MPI-CLM. Figure S5. As for Figure S3, but for the climate model MPI-CLM.

**Author Contributions:** Conceptualization, H.F., I.G.d.C.-A., and J.A.S.; Data curation, R.C., A.F., M.C.V., C.C., and S.R.; Formal analysis, R.C., H.F., and S.R.; Funding acquisition, J.A.S.; Investigation, R.C., H.F., A.F., M.C.V., C.C., S.R., and J.A.S.; Methodology, R.C., A.F., I.G.d.C.-A., and J.A.S.; Project administration, J.A.S.; Resources, A.F., I.G.d.C.-A., and J.A.S.; Software, R.C., H.F., A.F., I.G.d.C.-A., and S.R.; Supervision, H.F., I.G.d.C.-A., and J.A.S.; Validation, I.G.d.C.-A.; Visualization, S.R.; Writing—original draft, R.C.; Writing—review and editing, H.F., A.F., I.G.d.C.-A., M.C.V., C.C., S.R., and J.A.S.

**Funding:** This research and the APC were funded by the European Union's Horizon 2020 Research and Innovation Programme, grant number 810176.

**Acknowledgments:** This study was supported by the Clim4Vitis project—"Climate change impact mitigation for European viticulture: knowledge transfer for an integrated approach", funded by European Union's Horizon 2020 Research and Innovation Programme, under grant agreement n° 810176. This work was also supported by National Funds by FCT—Portuguese Foundation for Science and Technology, under the project UID/AGR/04033/2019. FCT scholarship SFRH/BPD/119461/2016 and contract CEECIND/00447/2017 are acknowledged. The authors also acknowledge the FCT scholarship given to R.C., PD/BD/135326/2017, under the Doctoral Programme 'Agricultural Production Chains—from fork to farm' (PD/00122/2012). Climate data was provided by ADVID, SOGRAPE vinhos SA, Quinta do Noval, and Direcção Regional de Agricultura e Pescas do Norte (DRAPN). Phenological data from DDR was obtained through the support of ADVID members.

**Conflicts of Interest:** The authors declare no conflict of interest.
## References

1. OIV. *State of the Vitiviniculture World Market*; OIV: Paris, France, 2018; p. 14.

2. IVV. *Vinhos e Aguardentes de Portugal, Anuário 2017*; Ministério da Agricultura, do Desenvolvimento Rural e das Pescas, Instituto da Vinha e do Vinho: Lisboa, Portugal, 2017; p. 224.

3. Fraga, H.; Santos, J.A. Vineyard mulching as a climate change adaptation measure: Future simulations for Alentejo, Portugal. *Agric. Syst.* **2018**, *164*, 107–115. [CrossRef]

4. Fraga, H.; Costa, R.; Santos, J.A. Multivariate clustering of viticultural terroirs in the douro winemaking region. *Cienc. Tec. Vitivinic.* **2017**, *32*, 142–153. [CrossRef]

5. Fraga, H.; Santos, J.A. Daily prediction of seasonal grapevine production in the Douro wine region based on favourable meteorological conditions. *Austr. J. Grape Wine Res.* **2017**, *23*, 296–304. [CrossRef]

6. Keller, M. Phenology and growth cycle. In *The Science of Grapevines*, 2nd ed.; Keller, M., Ed.; Academic Press: San Diego, CA, USA, 2015; Chapter 2; pp. 59–99. [CrossRef]

7. Moriondo, M.; Ferrise, R.; Trombi, G.; Brilli, L.; Dibari, C.; Bindi, M. Modelling olive trees and grapevines in a changing climate. *Environ. Model. Softw.* **2015**, *72*, 387–401. [CrossRef]

8. Andreoli, V.; Cassardo, C.; La Iacona, T.; Spanna, F. Description and preliminary simulations with the Italian vineyard integrated numerical model for estimating physiological values (IVINE). *Agronomy* **2019**, *9*, 94. [CrossRef]

9. Winkler, A.J. *General Viticulture*; University of California Press: Berkly, CA, USA, 1974.

10. Huglin, P. *Nouveau Mode D'évaluation des Possibilités Hélio thermiques d'un Milieu Viticole*. COMPTES RENDUS de L'Académie D'AGRICULTURE; Académie d'agriculture de France: Paris, France, 1978.

11. Jones, G.V.; Davis, R.E. Climate influences on grapevine phenology, grape composition, and wine production and quality for Bordeaux, France. *Am. J. Enol. Viticult.* **2000**, *51*, 249–261.

12. Jones, G.V.; Duchêne, E.; Tomasi, D.; Yuste, J.; Braslavska, O.; Schultz, H.R.; Martinez, C.; Boso, S.; Langellier, F.; Perruchot, C.; et al. Changes in European winegrape phenology and relationships with climate. In *Proceedings of the XIV GESCO Symposium, Geisenheim, Germany*, 23–26 August 2005.

13. Van Leeuwen, C.; Garnier, C.; Agut, C.; Baculat, B.; Barbeau, G.; Besnard, E.; Bois, B.; Boursiquot, J.-M.; Chuine, I.; Dessup, T.; et al. Heat requirements for grapevine varieties is essential information to adapt plant material in a changing climate. In *Proceedings of the 7th International Terroir Congress(Agroscope Changins-Wädenswil: Switzerland)*, Changins, Switzerland, 19–23 May 2008; pp. 222–227.

14. Chuine, I.; García de Cortázar-Atauri, I.; Kramer, K.; Hänninen, H. Plant development models. In *Phenology: An Integrative Environmental Science*; Springer: Berlin, Germany, 2013; pp. 275–293.

15. Orlandi, F.; Bonofiglio, T.; Aguilera, F.; Fornaciari, M. Phenological characteristics of different winegrape cultivars in Central Italy. *Vitis* **2015**, *54*, 129–136.

16. Quenol, H.; García de Cortázar-Atauri, I.; Bois, B.; Sturman, A.; Bonnardot, V.; Le Roux, R. Which climatic modeling to assess climate change impacts on vineyards? *Oeno One* **2017**, *51*, 91–97. [CrossRef]

17. Pouget, R. Le débourrement des bourgeons de la vigne: Méthode de prévision et principes d'établissement d'une échelle de précité de débourrement. *Connaiss. Vigne Vin* **1988**, *22*, 105–123. [CrossRef]

18. Riou, C.; Carbonneau, A.; Becker, N.; Caló, A.; Costacurta, A.; Castro, R.; Pinto, P.A.; Carneiro, L.C.; Lopes, C.; Clímaco, P.; et al. *Le Determinisme Climatique de la Maturation du Raisin: Application au Zonage de la Teneur en Sucre Dans la Communauté Européenne*; Office des Publications Officielles des Communautés Européennes: Luxembourg, 1994; p. 319.

19. García de Cortázar-Atauri, I.; Brisson, N.; Gaudillere, J. Performance of several models for predicting budburst date of grapevine (Vitis vinifera L.). *Int. J. Biometeorol.* **2009**, *53*, 317–326. [CrossRef]

20. Nendel, C. Grapevine bud break prediction for cool winter climates. *Int. J. Biometeorol.* **2010**, *54*, 231–241. [CrossRef] [PubMed]

21. Caffarra, A.; Eccel, E. Projecting the impacts of climate change on the phenology of grapevine in a mountain area. *Austr. J. Grape Wine Res.* **2011**, *17*, 5–61. [CrossRef]

22. Fila, G.; Gardiman, M.; Belvini, P.; Meggio, F.; Pitacco, A. A comparison of different modelling solutions for studying grapevine phenology under present and future climate scenarios. *Agr. Forest Meteorol.* **2014**, *195*, 192–205. [CrossRef]

23. Amerine, M.; Winkler, A. Composition and quality of musts and wines of California grapes. *Hilgardia* **1944**, *15*, 493–675. [CrossRef]
24. Duchene, E.; Huard, F.; Dumas, V.; Schneider, C.; Merdinoglu, D. The challenge of adapting grapevine varieties to climate change. *Clim. Res.* **2010**, *41*, 193–204. [CrossRef]

25. García de Cortázar-Atauri, I.; Chuine, I.; Donatelli, M.; Parker, A.; Van Leeuwen, C. A curvilinear process-based phenological model to study impacts of climatic change on grapevine (Vitis vinifera L.). *Proc. Agro* **2010**, *11*, 907–908.

26. Parker, A.K.; García de Cortázar-Atauri, I.; van Leeuwen, C.; Chuine, I. General phenological model to characterise the timing of flowering and veraison of Vitis vinifera L. *Austr. J. Grape Wine Res.* **2011**, *17*, 206–216. [CrossRef]

27. Parker, A.; García de Cortázar-Atauri, I.; Chuine, I.; Barbeau, G.; Bois, B.; Boursiquot, J.-M.; Cahurel, J.-Y.; Claverie, M.; Dufourcq, T.; Gény, L.; et al. Classification of varieties for their timing of flowering and veraison using a modelling approach: A case study for the grapevine species Vitis vinifera L. *Agric. For. Meteorol.* **2013**, *180*, 249–264. [CrossRef]

28. Molitor, D.; Junk, J.; Evers, D.; Hoffmann, L.; Beyer, M. A High-resolution cumulative degree day-based model to simulate phenological development of grapevine. *Am. J. Enol. Viticult.* **2014**, *65*, 72–80. [CrossRef]

29. García de Cortázar-Atauri, I. Adaptation du modèle STICS à la vigne (Vitis vinifera L.): Utilisation dans le cadre d'une étude d'impact du changement climatique à l'échelle de la France. Ph.D. Thesis, École nationale supérieure agronomique, Montpellier, France, 2006.

30. Pieri, P. Climate change and grapevines: Main impacts. In *Climate Change, Agriculture and Forests in France: Simulations of the Impacts on the Main Species. The Green Book of the CLIMATOR Project (2007–2010)*; ADEME: Angers, France, 2010; pp. 213–223.

31. Cuccia, C.; Bois, B.; Richard, Y.; Parker, A.K.; García de Cortázar-Atauri, I.; Van Leeuwen, C.; Castel, T. Phenological Model Performance to Warmer Conditions: Application to Pinot Noir in Burgundy. *J. Int. Sci. Vigne Vin* **2014**, *48*, 169–178. [CrossRef]

32. Molitor, D.; Caffarra, A.; Sinigoj, P.; Pertot, I.; Hoffmann, L.; Junk, J. Late frost damage risk for viticulture under future climate conditions: A case study for the Luxembourgish winegrowing region. *Austr. J. Grape Wine Res.* **2014**, *20*, 160–168. [CrossRef]

33. Moriondo, M.; Bindi, M. Impact of climate change on the phenology of typical mediterranean crops. *Italian J. Agrometeorol.* **2007**, *3*, 5–12.

34. Webb, L.B.; Whetton, P.H.; Barlow, E.W.R. Modelled impact of future climate change on the phenology of winegrapes in Australia. *Austr. J. Grape Wine Res.* **2007**, *13*, 165–175. [CrossRef]

35. Petrie, P.R.; Sadras, V.O. Advancement of grapevine maturity in Australia between 1993 and 2006: Putative causes, magnitude of trends and viticultural consequences. *Austr. J. Grape Wine Res.* **2008**, *14*, 33–45. [CrossRef]

36. Fraga, H.; García de Cortázar-Atauri, I.; Malheiro, A.C.; Santos, J.A. Modelling climate change impacts on viticultural yield, phenology and stress conditions in Europe. *Glob. Chang. Biol.* **2016**, *22*, 3774–3788. [CrossRef]

37. Vivin, P.; Lebon, E.; Dai, Z.; Duchêne, E.; Marguerit, E.; García de Cortázar-Atauri, I.; Zhu, J.; Simonneau, T.; Van Leeuwen, C.; Delrot, S. Combining ecophysiological models and genetic analysis: A promising way to dissect complex adaptive traits in grapevine. *Oeno One* **2017**, *51*, 181. [CrossRef]

38. Van Leeuwen, C.; Seguin, G. The concept of terroir in viticulture. *J. Wine Res.* **2006**, *17*, 1–10. [CrossRef]

39. Wolkovich, E.M.; García de Cortázar-Atauri, I.; Morales-Castilla, I.; Nicholas, K.A.; Lacombe, T. From Pinot to Xinomavro in the world's future wine-growing regions. *Nat. Clim. Chang.* **2018**, *8*, 29–37. [CrossRef]

40. Magalhães, N. *Tratado de Viticultura: A Videira, a Vinha e o Terroir*; Chaves Ferreira: Lisboa, Portugal, 2008; p. 605.

41. Lorenz, D.H.; Eichhorn, K.W.; Bleiholder, H.; Klose, R.; Meier, U.; Weber, E. Growth Stages of the Grapevine: Phenological growth stages of the grapevine (Vitis vinifera L. ssp. vinifera)—Codes and descriptions according to the extended BBCH scale. *Austr. J. Grape Wine Res.* **1995**, *1*, 100–103. [CrossRef]

42. Fonseca, A.; Santos, J. High-resolution temperature datasets in Portugal from a geostatistical approach: Variability and extremes. *J. Appl. Meteorol. Climatol.* **2018**, *57*, 627–644. [CrossRef]

43. Caffarra, A.; Eccel, E. Increasing the robustness of phenological models for Vitis vinifera cv. Chardonnay. *Int. J. Biometeorol.* **2010**, *54*, 255–267. [CrossRef]
44. Bidabe, B. Action of temperature on evolution of buds in apple and comparison of methods used for estimating flowering time. *Ann. Physiol. Veg.* **1967**, *9*, 65.

45. Chuine, I. A unified model for budburst of trees. *J. Theor. Biol.* **2000**, *207*, 337-347. [CrossRef]

46. Bonhomme, M.; Rageau, R.; Lacointe, A. Optimization of endodormancy release models, using series of endodormancy release data collected in france. *Acta Hortic.* **2010**, *872*, 51-60. [CrossRef]

47. De Réaumur, M. Observations du thermometre, faites à Paris pendant l'année 1735 compareés avec celles qui ont été faites sous la ligne à l' Ile de France, à Alger et en quelques-unes de nos îles de l'Amérique. *Memoires de l'Academie Royale des Sciences* **1753**, 545-576.

48. Richardson, E. A model for estimating the completion of rest for 'Redhaven' and 'Elberta' peach trees. *HortScience* **1974**, *9*, 331-332.

49. Hänninen, H. Modelling bud dormancy release in trees from cool and temperate regions. *Acta Forestalia Fennica* **1990**, *213*, 1-47. [CrossRef]

50. Wang, E.L.; Engel, T. Simulation of phenological development of wheat crops. *Agric. Syst.* **1998**, *58*, 1-24. [CrossRef]

51. Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. Equation of state calculations by fast computing machines. *J. Chem. Phys.* **1953**, *21*, 1087-1092. [CrossRef]

52. Chuine, I.; Bonhomme, M.; Legave, J.M.; García de Cortzar-Atauri, I.; Charrier, G.; Lacointe, A.; Améglio, T. Can phenological models predict tree phenology accurately in the future? The unrevealed hurdle of endodormancy break. *Glob. Chang. Biol.* **2016**, *22*, 3444-3460. [CrossRef]

53. Wallach, D.; Makowski, D.; Jones, J.W. *Working with Dynamic Crop Models: Evaluation, Analysis, Parameterization, and Applications*; Elsevier: Amsterdam, The Netherlands, 2006.

54. Jacob, D.; Petersen, J.; Eggert, B.; Alias, A.; Christensen, O.B.; Bouwer, L.M.; Braun, A.; Colette, A.; Déqué, M.; Georgievski, G.; et al. EURO-CORDEX: New high-resolution climate change projections for European impact research. *Reg. Environ. Chang.* **2014**, *14*, 563-578. [CrossRef]

55. Yang, W.; Andreasson, J.; Graham, L.P.; Olsson, J.; Rosberg, J.; Wetterhall, F. Distribution-based scaling to improve usability of regional climate model projections for hydrological climate change impacts studies. *Hydrol. Res.* **2010**, *41*, 211-229. [CrossRef]

56. Ramos, M.C. Projection of phenology response to climate change in rainfed vineyards in north-east Spain. *Agr. Forest Meteorol.* **2017**, *247*, 104-115. [CrossRef]

57. Fraga, H.; Costa, R.; Moutinho-Pereira, J.; Correia, C.M.; Dinis, L.-T.; Gonçalves, I.; Silvestre, J.; Eiras-Dias, J.; Malheiro, A.C.; Santos, J.A. Modeling phenology, water status, and yield components of three portuguese grapevines using the STICS crop model. *Am. J. Enol. Viticult.* **2015**, *66*, 482-491. [CrossRef]

58. Fraga, H.; Amraoui, M.; Malheiro, A.C.; Moutinho-Pereira, J.; Eiras-Dias, J.; Silvestre, J.; Santos, J.A. Examining the relationship between the enhanced vegetation index and grapevine phenology. *Eur. J. Remote Sens.* **2014**, *47*, 753-771. [CrossRef]

59. Cola, G.; Mariani, L.; Salinari, F.; Civardi, S.; Bernizzoni, F.; Gatti, M.; Poni, S. Description and testing of a weather-based model for predicting phenology, canopy development and source-sink balance in Vitis vinifera L. cv. Barbera. *Agric. Forest Meteorol.* **2014**, *184*, 117-136. [CrossRef]

60. Tomasi, D.; Jones, G.V.; Giust, M.; Lovat, L.; Gaiotti, F. Grapevine phenology and climate change: Relationships and trends in the Veneto Region of Italy for 1964-2009. *Am. J. Enol. Viticult.* **2011**, *62*, 329-339. [CrossRef]

61. Fraga, H.; Santos, J.A.; Moutinho-Pereira, J.; Carlos, C.; Silvestre, J.; Eiras-Dias, J.; Mota, T.; Malheiro, A.C. Statistical modelling of grapevine phenology in Portuguese wine regions: Observed trends and climate change projections. *J. Agr. Sci. Camb.* **2016**, *154*, 795-811. [CrossRef]

62. Bindi, M.; Miglietta, F.; Gozzini, B.; Orlandini, S.; Seghi, L. A simple model for simulation of growth and development in grapevine (Vitis vinifera L.). 2. Model validation. *Vitis* **1997**, *36*, 73-76.

63. Ramos, M.C.; Jones, G.V.; Yuste, J. Phenology of Tempranillo and Cabernet-Sauvignon varieties cultivated in the Ribera del Duero DO: Observed variability and predictions under climate change scenarios. *Oeno One* **2018**, *52*, 31-44. [CrossRef]

64. Oliveira, M. Calculation of budbreak and flowering base temperatures for vitis vinifera cv. touriga francesa in the Douro Region of Portugal. *Am. J. Enol. Viticult.* **1998**, *49*, 74-78.
65. Moncur, M.W.; Rattigan, K.; Mackenzie, D.H.; McIntyre, G.N. Base temperatures for budbreak and leaf appearance of grapevines. *Am. J. Enol. Viticult.* **1989**, *40*, 21–26.

66. Jones, G.V. Climate change in the western united states grape growing regions. *Proceedings of the Seventh International Symposium on Grapevine Physiology and Biotechnology*. *Acta Horticult.* **2005**, *689*, 41–59. [CrossRef]

67. Bock, A.; Sparks, T.; Estrella, N.; Menzel, A. Changes in the phenology and composition of wine from Franconia, Germany. *Clim. Res.* **2011**, *50*, 69–81. [CrossRef]

68. Jones, G.V.; White, M.A.; Cooper, O.; Storchmann, K. Climate change and global wine quality. *Clim. Chang.* **2005**, *73*, 319–343. [CrossRef]

69. Malheiro, A.C.; Campos, R.; Fraga, H.; Eiras-Dias, J.; Silvestre, J.; Santos, J.A. Winegrape phenology and temperature relationships in the Lisbon Wine Region, Portugal. *J. Int. Sci. Vigne Vin.* **2013**, *47*, 287–299. [CrossRef]

<figure><img src="image_9.png" /></figure>

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).