## Introduction

The Dão winegrowing region is located on a plateau in the central part of northern Portugal. It is surrounded by mountains and is protected from both the Atlantic moist winds and the continental influence from inner Iberia. In this wine region small farms predominate. Approximately 80 % of the vineyards are cultivated with red varieties and 20 % with white varieties. For the period 1950–2000, high-resolution bioclimatic zoning divides the region into three bioclimatic groups (Fraga et al., 2014): i) the peripheral zone comprising the northern, western and south/southeastern areas and characterised by a cold and humid climate and cold nights, ii) the central area with a temperate and humid climate and cold nights, and iii) the south/southwestern area with a temperate and humid climate and warm nights. Climate change projections for 2041–2070 reveal that the climate will become warmer and drier with warm nights in about 85 % of the area, while it will become very warm and dry with warm nights in 14 % of the area. Ninety-seven percent of vineyard soils are derived from granite (Dias, 1995), and are classified as Humic or as Dystric Cambisols according to the FAO Soil Classification (FAO, 2015). These are coarse-textured soils, which have a slightly acid reaction and low organic matter content, but they are reasonably enriched with phosphorus (Almeida, 2005).

Touriga Nacional (TN) and Encruzado (EN) are native varieties (red and white respectively) from Portugal, which are recognised in the Dão region for their aptitude to produce high-quality table wines with high ageing potential.

Grapevine growth, fruit ripening and berry quality are strongly dependent on the environmental conditions of the sites where the vines are cultivated (Martínez-Lüscher et al., 2016). Climate and soil are the two central terroir factors, having the most significant influence in viticulture (Jones and Davis, 2000; Magalhães, 2008). Due to this dependence, many studies have identified different relationships between some climatic parameters and vine growth, yield, fruit ripening and berries and wine quality (Costa et al., 2020; Jones and Davis, 2000; Malheiro et al., 2013). Some of these relationships have been integrated into simulation and forecasting models (Costa et al., 2015).

The grapevine vegetative cycle, morphologically described by its phenological stages, is strongly determined by weather conditions, namely precipitation, radiation and temperature (Malheiro et al., 2013). Grapevine phenological development has been the subject of several studies which aimed to characterize its relationship with the meteorological conditions (Martínez-Lüscher et al., 2016), to develop models to make predictions and to classify the varieties according to the heat requirements of the main phenological stages. Based on phenological data from the Lisbon wine region, Lopes et al. (2008) defined a classification of 19 white and 15 red Portuguese varieties according to their heat requirements for budbreak, flowering and veraison. Using the same methodology, Alves et al. (2013) determined the heat requirements for four red grape varieties in the Douro wine region and evaluated the rootstock effect. Based on previously calibrated phenological models, Parker et al. (2013) carried out a classification of 95 varieties according to flowering precocity, and a classification of 104 varieties in terms of veraison precocity. Following a similar methodology, Reis et al. (2020) classified 36 Portuguese varieties using data from the Lisbon wine region.

Most phenological models commonly used to simulate and forecast phenological events are based on the assumption that air temperature is the preponderant environmental factor which determines vine development. These models are widely applied in viticulture as they are an important tool for planning canopy and irrigation management and harvesting, as well as for defining strategies for disease control (Caffarra et al., 2012; Molitor and Berkelmann-Loehnertz, 2011). Furthermore, they can be used to select the varieties better adapted to specific weather and climate conditions (Parker et al., 2013; Parker et al., 2011; Reis et al., 2020), as well as to assess climate change impacts on grapevine development when coupled with climate change scenarios and temperature projections (Costa et al., 2019; Fraga et al., 2016).

The simplest models, called thermal models (TM), only take into account the forcing effect of heat accumulation after a specific date (onset). More complex models, which describe the entire dormancy period (endo- and eco-dormancy), also take into account the chilling accumulation requirement to break dormancy (Jones, 2013). The chilling models (CM) differ in articulating the thermal effects on dormancy and vegetative phases (Reis et al., 2020). Leolini et al. (2020) evaluated the reliability of these two modelling approaches for budbreak simulation by testing six different phenological models (3 TM and 3 CM) and eight varieties cultivated at different latitudes in Europe. Although CMs have shown a generally higher estimation accuracy, none of them have shown clear supremacy over the simpler models, thus highlighting that the TMs may be a good trade-off between accuracy and simplicity (Garcia de Cortázar-Atauri et al., 2009). As cultivated grapevine, Vitis vinifera, is considered to be a species which has a relatively low chilling requirement compared to other perennial woody horticultural species (Londo and Johnson, 2014), chill-heat models do not usually perform significantly better than heat accumulation models. In the Dão wine region, the chilling requirements are indeed consistently fulfilled, phenological development thus being largely controlled by heat accumulation.

Among the TMs, the degree-day (DD) is one of the most widely used for grapevine phenological predictions. In this model, the temperature effect on phenological dynamics is determined through daily temperature accumulation above a pre-defined threshold (T0) often referred to as base temperature, being the accumulation of temperature between a starting day, or previous phenological state, and the target phenological state.

Using phenological data from twelve vineyards in Chile, Santibanez et al. (2014) presented models for three table grape cultivars based on Mitscherlich´s monomolecular equation, where the dependent and independent variables were the phenological stages and the DD respectively, and T0 = 10 ºC and tx = 1 (1st January). Ortega-Farías et al. (2002) used the same approach for two grapevine varieties. Using phenological observations of 81 varieties collected from 23 different locations (pre-dominantly in France), Parker et al. (2011) obtained better results when T0 = 0 ºC and tx = 60 (1st March). However, using data related to eight varieties collected from fifteen experimental vineyards in three European countries, Leolini et al. (2020) found performance to be higher for three TMs when the starting date was fixed on the 1st January (tx = 1) rather than 1st March (tx = 60).

When considering conventional DD to assumes that above T0 the effect of temperature on phenological development is linear and unlimited, Molitor et al. (2014) compared it with two modifications that incorporated additional thresholds. The first modification comprised a second higher temperature threshold, with which the temperature effect on phenological development was maximum and constant. A second modification comprised a third temperature threshold (Tmax), above which temperature accumulation decreased. In both models, the accumulation of heat began on the budbreak date. For temperature thresholds of 5, 20 and 22 ºC, the latter model showed the best prediction efficiency.

The two-parameter sigmoid model (SM) proposed by Hanninen (1990) overcomes the limitation of imposing an artificial threshold for heat accumulation as it allows for a gradual change in the weights given to different temperatures within a transition range (Reis et al., 2020). A recent study that tested several phenological models applied to two varieties cultivated in Portugal (cv. Touriga Nacional and cv. Touriga Franca) and prevalent in the Douro wine region demonstrated that the DD and SM tended to show higher performance when simulating flowering and veraison than other commonly used models, such as Richardson and Wang models (Costa et al., 2019). Parker et al. (2011) also showed that both these models give the highest performance when simulating veraison dates.

Most models are calibrated, validated or tested to simulate only the main phenological stages; i.e., budbreak, flowering and veraison (Caffarra and Eccel, 2010; Cuccia et al., 2014; Garcia de Cortázar-Atauri et al., 2009; Fraga et al., 2015; Morales-Castilla et al., 2020; Parker et al., 2011). Consistent records of intermediate phenological stages (in between budbreak, flowering and veraison) are generally difficult to obtain, thus explaining why most of the studies are focused on the three aforementioned stages. Nonetheless, other models have been developed and applied to predict the whole sequence of phenological stages in the grapevine growing cycle (Fernández-González et al., 2013; Mariani et al., 2013; Molitor et al., 2014; Ortega-Farías et al., 2002; Rodrigues et al., 2016; Verdugo-Vásquez et al., 2017).

The forecast of the intermediate phenological stages is particularly relevant when models are used as decision support tools for farmers when planning canopy management and vine protection, which involve, for example, shoot thinning, leaf removal and phytosanitary measures against downy mildew disease or Botrytis bunch rot.

An important consideration when developing phenology models is the quantity and quality of phenological observations needed to train and validate the model. Ideally, long time series data (or panel data) collected consistently over time should be used to model phenology, as this will reduce model error due to differences in crop management and will minimise observer errors (Darbyshire et al., 2020). Given the reality of available data in many studies some researchers have compiled data from different sites to construct a larger dataset (Darbyshire et al., 2020; Luedeling et al., 2009; Morales-Castilla et al., 2020; Parker et al., 2013; Parker et al., 2011; Parker et al., 2020; Pope et al., 2014).

The present study aimed to 1) provide a detailed discussion of the phenological observations of two representative winegrape varieties (cv. Touriga Nacional and cv. Encruzado) grown in the Dão wine region (DãoWR), Portugal, and its climatic conditions, 2) assess whether a single set of parameters is sufficient to accurately model different phenological timings in both varieties, and 3) calibrate and validate a Sigmoid model to simulate several intermediate phenological stages between budbreak and veraison. To our knowledge, no previous studies has been carried out on either of these two varieties and their intermediate phenological timings. Therefore, the present study aims at assessing the performance of the SMs when reproducing the observed variability of less studied phenological stages. Some of the results may, however, be extrapolated to other varieties and wine regions worldwide.

## Materials and methods

### 1. Study Area and Datasets

The phenological data were collected from nine commercial vineyards and two varietal vineyards in six locations in the Dão Wine Region, hereafter referred to as DãoWR (Figure 1). All the vineyards are non-irrigated. The characteristics of the experimental plots (variety, plantation density, pruning system and rootstock) and the available periods of data are shown in Table 1.

Figure 1. Map of mainland Portugal with Dão Wine Region and TN vineyard plot (red square), EN vineyard plot (Green circle) and weather station locations (blue triangles).

The phenological stages were recorded through site observations based on the modified Eichhorn-Lorenz scale (Coombe, 1995) until the beginning of veraison (EL35). To evaluate the veraison evolution, after EL35 (“Berries begin to colour and enlarge”), a specific scale was used in which 50V corresponds to 50 % of coloured clusters and represents the mid-veraison. For each varietal plot, the observations were carried out on all buds/shoots (until EL18) or inflorescences/clusters (after EF18) of six plants. On the commercial plots, the observations were made on the two buds/shoots (or on its inflorescences/clusters) of the same fruiting spur. In each repetition (three repetitions per plot), five fruiting spurs with two buds on different positions in the cordon and on different plants were selected. Records were taken twice a week until EL27, and once a week after this stage. All observations were undertaken by the same technicians, thus ensuring a homogeneous dataset. For each plant in the varietal plots, or for each repetition in the commercial plots, a phenological stage was evaluated when at least 50 % of the observed parts of the plant reached the respective stage.

As in the commercial plots the measurements were carried out in three repetitions, the maximum number of observations for each phenological stage is thus 51 for cv. Touriga Nacional and 40 for cv. Encruzado (see Table S1 in Supplementary Material). Although these sample sizes are relatively short for more robust statistical assessments, this is indeed the only available dataset of grapevine phenology in this wine region, thus indicating the need to maintain the current observational framework.

Due to the homogeneity of the observation criteria, the proximity of the weather stations to the plots and to the fact that all the plots were pruned in the second half of January, the implications for the model accuracy of different observations between the commercial and varietal vineyards may be comparable, or even less important, than those that could occur in models calibrated with complied datasets from different locations.

The weather stations are located in the vineyards very close to the selected plots (< 100 meters away), except for CP6 and CP9 in which the weather stations (WS5 and WS6) are located 1.3 km and 3.7 km away from the plots respectively; however, in these two cases, due to the relatively flat terrain, the weather stations are still representative of the atmospheric conditions in the corresponding vineyard plots.

Table 1. Geographical location of the selected plots (latitude, longitude and elevation), grown variety, rootstock, planting density, pruning system, weather station identifier, and available period of phenological data.

Plot

Location

Latitude (N)

Longitude (W)

Elevation (m)

Variety

Planting density
(plant/ha)

Pruning system

Rootstock

Available period

VP1

Viseu

40°38´30´´

7°54´32´´

456

TN
EN

3953

URC

110R

2014–2019
2016–2019

VP2

Nelas

40°31´30´´

7°51´26´´

427

TN

4545

BRC

SO4

2014–2019

CP1

Silgueiros

40°35´19´´

7°55´03´´

332

TN

4000

URC

1103P

2017–2019

CP2

Silgueiros

40°35´20´´

7°55´11´´

334

EN

4000

SG

196/17

2017–2019

CP3

Nelas

40°31´27´´

7°51´25´´

428

TN

4545

BRC

100R

2017–2019

CP4

Santar

40°34´17´´

7°53´05´´

372

EN

4000

URC

a)

2017–2019

CP5

Carregal do Sal

40°25´33´´

7°59´42´´

291

TN

3333

URC

a)

2017–2019

CP6

Carregal do Sal

40°25´45´´

8°00´58´´

293

EN

3333

BRC

a)

2017–2019

CP7

Tábua

40°20´32´´

8°01´25´´

296

TN

4545

BRC

1103P

2017–2019

CP8

Tábua

40°20´27´´

8°01´29´´

278

EN

4545

URC

1103P

2017–2019

CP9

São João de Areias

40°23´19´´

8°04´21´´

238

TN

4545

BRC

a)

2019

VP = Varietal Plot, CP = Commercial Plot, TN = Touriga Nacional, EN = Encruzado, URC = Unilateral Royat cordon, BRC = Bilateral Royat cordon, SG = Single Guyot; a) refers to unknown rootstock.

A preliminary quality check was carried out on the temperature time series at each location. For all weather stations, data gaps corresponded to less than 4 % of the entire time series length. These gaps were filled using linear regression estimations between the daily temperatures recorded at a given weather station and the nearest station. In all cases, very high correlation coefficients were found between both datasets on the daily timescale (> 0.95, statistically significant at a 99 % confidence level).

### 2. Phenological models

The thermal models were adjusted to simulate budbreak (EL4) and several phenological stages between budbreak and veraison (Table 2). These models only consider the effect of forcing temperatures and assume that a given phenological stage occurs when the daily sum of the rate of forcing (R), as from the onset date (tx), reaches a specific critical value (F*) (Parker et al., 2011):

${S}_{t}={\sum }_{{t}_{x}}^{t}R\left({Tm}_{i}\right)\ge {F}^{*}$

with the daily rate of forcing (R) calculated using a sigmoid function (Costa et al., 2019):

where Tmi represents the daily mean temperature, while e and d are the fit parameters, determining the inflection point and sharpness of the curve respectively. Hence, the model contains four parameters: d and e describe the forcing rate function, while tx and F* are related to its integration.

Table 2. Phenological stages simulated and onset date of models (tx).

Phenological stages

EL Number

Morphological description

Onset date (tx)

EL4

4

Budburst; leaf tips visible

1st January

15 January

1st February

15 February

1st March

EL9

9

2 to 3 separated leaves; 2–4 cm long shoots

Budbreak (EL4)

EL12

5 separated leaves; shoots about 10 cm long; clear inflorescence

12

EL17

12 separated leaves; well-developed inflorescence, separated single flowers

17

EL19

19

About 16 separated leaves; beginning of flowering (first flower caps loosening)

EL23

23

17–20 separated leaves; 50 % caps off (= flowering)

EL27

27

Setting; young berries enlarging (> 2 mm diam.), bunches at right angles to stems

Flowering (EL23)

EL29

29

Peppercorn-size berries (4 mm diam.); bunches curling downwards

EL31

31

Pea-size berries (7 mm diam.)

EL32

32

Beginning of bunch closure, berries touching (if bunches are tight)

EL35

35

Berries begin to colour and enlarge

50V

-

50 % of cluster coloured

EL Number – modified Eichhorn-Lorenz scale (Coombe, 1995).

The onset dates (tx), previously fixed for heat accumulation to fulfill each phenological stage, are shown in Table 2. Concerning tx, the budbreak (EL4) date of each year was chosen for stages from EL9 to EL23, and the flowering (EL23) date of each year for stages from EL27 to EL35 and 50V. For budbreak, the model was fitted with several possible onset dates for heat accumulation: 1st January (tx = 1), 15 January (tx = 15), 1st February (tx = 32), 15 February (tx = 47) and 1st March (tx = 60).

As the heat accumulation of each intermediate stage was computed from the previous main phenological state (budbreak or flowering), the respective F* parameters thus are a manifestation of the heat requirements of the corresponding phenophases. Once the onset date for heat accumulation was fixed, the forcing rate function parameters e and d were fitted in the model calibration, also leading to a given F*.

### 3. Model selection

The Phenology Modeling Platform (PMP), version 5.5 (Chuine et al., 2013), was used to calibrate and validate the models. PMP estimates parameters based on the Metropolis annealing algorithm (Chuine et al., 2013). For each specific location, the input data comprise the weather station daily minimum, mean and maximum air temperatures, latitude and the observed phenological stages (in days of the year).

For each variety and phenological stage, taking into account the fixed tx for each phenophase, the model parameterisation was performed in two steps. In the first step, the e coefficient (location parameter) was fixed in successive integer values, between a lower and an upper threshold, the remaining parameters (d and F*) being fitted for each e value separately. In the second step, for the e coefficient selected from the analysis of the first step, the responsiveness of the models (variation of the performance) as a function of the d coefficient was analysed. A similar model parameterisation was used by Parker et al., (2020). In present study for budbreak (EL4), the integer values of the e coefficient set in each fit ranged between 6 and 14; for the other phenological stages, the e coefficient ranged between 8 and 15. These ranges were defined according to the results obtained by Reis et al., (2020) for Portuguese varieties.

The root-mean-squared-error (RMSE) and the Nash-Sutcliffe efficiency coefficient (EFF) were used to assess model performance:

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

$EFF=1-\frac{{\sum }_{i=1}^{n}{\left({O}_{i}-{P}_{i}\right)}^{2}}{{\sum }_{i=1}^{n}{\left({O}_{i}-{O}_{m}\right)}^{2}}$

where Oi represents the observed phenological date on the i-th year (vintage), Pi the corresponding simulated values, Om the average of all observed values and n the sample size.

In order to take into account potential model overfitting when assessing model performances, model calibration should be followed by model validation. Preferably the validation should be carried out using independent subsets; e.g., randomly selected sub-samples from the original dataset. Nonetheless, owing to the short sample sizes available for each phenological stage and variety, this methodology cannot be robustly applied. Therefore, the leave-one-out cross-validation method was carried out for all of the selected models. This cross-validation method is applied once for each data point, using all the other points as a training set and the selected one as a single-item test set. The RMSE metric was accordingly adapted to the Root Mean Square Error of Prediction (RMSEP), defined as (Morales-Castilla et al., 2020; Reis et al., 2020):

$RMSEP=\sqrt{\frac{{\sum }_{i=1}^{n}{\left({O}_{i}-{p}_{i}^{v}\right)}^{2}}{n}}$

where Pvi is the predicted value obtained from a leave-one-out cross-validation approach.

To test whether the models predicted better than the null hypothesis of average dates, RMSEa (Pope et al., 2014) was calculated from each observed data point (Oi) with the average of observed dates of all the other points (Omi) being used as a predicted value:

$RMSEa=\sqrt{\frac{{\sum }_{i=1}^{n}{\left({O}_{i}-{O}_{i}^{m}\right)}^{2}}{n}}$

## Results

### 1. Climatic characterisation

By analysing the ombrothermic diagram from WS1 (Viseu) for 2014–2019, two different periods of the year can be identified (Figure 2a). The months between October and March are relatively cool and rainy. As expected, December and January are the coldest months, whereas November is the rainiest. In the growing season, which lasts from April to September, precipitation is significantly lower, particularly from June to September (dry season), July and August being the warmest and driest months. According to the Köppen-Geiger climate classification, this climate is the Csb type; i.e., Mediterranean climate with cool summers.

In order to assess the representativeness of the selected years (2014–2019), their corresponding meteorological conditions are compared with a 30-year time period (1990–2019), namely for the temperatures averaged in the periods of January–March and April–September (Figure S1, Supplementary Material). The results hint at a fairly good coverage of the regional inter-annual climate variability, and can thus be considered to be robust and representative.

Figure 2. (a) Ombrothermic diagram comparing monthly minimum (Tmin), mean (Tmean) and maximum (Tmax) temperatures, and total precipitation from the weather station of Viseu (WS1), 2014–2019; (b) Mean growing season temperature (TGS) per year (averaged over all stations); (c) Mean growing season temperature (TGS) recorded at each weather station, 2017–2019; (d) Mean temperature before budburst (January to March) per year (averaged over all stations); (e) Mean temperature before budburst (January to March) recorded at each weather station, 2017–2019. The variation across sites (b and d) and the inter-annual variation (a, c and e) are also pointed out (only their lower halves in panels b-e).

During the study period, the mean growing season temperature (April to September; TGS) ranged from 17.7 °C at WS1 in 2019 to 20.2 °C at WS2 in 2017. In terms of the average of all the weather stations, the lowest TGS was recorded in 2014 and 2019, while the highest was recorded in 2017 (Figure 2b). For the period 2017 to 2019, the lowest TGS value was recorded at WS1 and the highest at WS6 (Figure 2c). Before budbreak (January to March), the average air temperature was highest in 2019 (Figure 2d). The lowest mean temperature before budbreak was recorded at WS1, and the highest at WS6 (Figure 2e). Overall, mild temperatures (daily mean values around 10 ºC) are typical during the dormancy period, whist growing season mean temperatures of approximately 18–19 ºC highlight the prevailing temperate warm conditions.

### 2. Characterisation of the phenological stages

In Figure 3, the box plots of the observed date of each phenological stage and each studied variety used for the calibration models are depicted. The corresponding descriptive statistics are shown in Table S1 (Supplementary Material). The cv. Encruzado (white variety) tends to reach budburst (EL4) earlier than cv. Touriga Nacional (red variety). In terms of the mean and median (2nd quartile), the difference between the budburst dates was 1 and 3 days respectively. However, Encruzado shows higher variability than Touriga Nacional, with interquartile-ranges (IQR) of 21 and 19 days respectively.

Figure 3. Box plots (mean, minimum, maximum and 1st, 2nd and 3rd quartiles) of the observed distributions of the outlined phenological stage dates (EL4 to EL35 and 50 % V), for cv. Encruzado (EN) and cv. Touriga Nacional (TN).

The inter-annual variability of this phenological event is particularly evident when comparing 2018 with 2019, when the average temperature between January to March was the lowest and the highest respectively (Figure 2d). On average, in 2019 budburst was 19 (EN) and 20 (TN) days earlier. The regional variability of the budburst timings was much more evident for TN in 2018; the January-March mean air temperature was at its lowest, suggesting that cooler conditions in the dormancy period may strengthen spatial variability of budburst dates.

Although the flowering (EL23) of the EN variety also occurred on average two days earlier than TN, the difference in the corresponding inter-annual variability is noteworthy (TN with IQR = 15.4 days, and EN with IQR = 20.5 days). Furthermore, the inter-annual variability was higher in this phenological stage. In agreement with the mean air temperature differences recorded between EL4 and EL19 in 2017 (15.9 °C) and 2018 (15.3 °C), EL23 occurred for both varieties about one month later in 2018. Albeit with significant inter-annual variation, the duration of the budbreak-flowering phase (EL4-EL23) was similar for both varieties: on average, 53 days for EN and 54 days for TN.

As for the previous phase, the duration of the period from EL23 to EL35 was similar for the two varieties, with significant inter-annual variability: in 2018, with a mean air temperature of 21.6 °C it lasted 55–56 days, while in 2019, with a mean air temperature of 19.5 °C it lasted 60–64 days.

This preliminary analysis indicates the high sensitivity of the phenological development of both varieties to air temperature, thus suggesting some predictability potential based on temperature accumulation models.

### 3. Phenological model optimisation

As previously mentioned, the model calibration was performed in two steps. This procedure allowed us to evaluate the feasibility of adjusting a single model to all phenological states of each phenophase (budbreak, budbreak-flowering and flowering-veraison), as well as to investigate whether it can be used to simulate the phenology of both varieties.

Table S2 (Supplementary Material) shows the EFF and RMSE of the models calibrated for EL4 with different e and tx. For both varieties, the highest EFF values correspond to the models with tx = 1st March and > 10. Thus, we can conclude that for both varieties the onset data (tx) for budbreak should be 1st March.

Figure 4 shows the changes in EFF and RMSE of the models fitted in the first step of our procedure; i.e., in response to the variation in the e coefficient.

For budbreak (Figures 4a and 4b), the efficiency of the model increased significantly for increasing e values between 6 and 12, but the increments are not significant for > 10 (< 0.008 for cv. Encruzado and < 0.01 for cv. Touriga Nacional). For the models with e = 12 (which for both varieties corresponds to the lowest RMSE), the d coefficient is -0.35 for EN (EFF = 0.63; RMSE = 6.52) and -0.33 for TN (EFF = 0.84; RMSE = 3.71) (Table S2). For the whole range of tested e values, the EFFs are systematically lower for EN than TN.

Since the aim was to develop a single model for both varieties, the effect of setting the same d coefficient in the second step on the model's performance was investigated. When the same model fitted for EN was used for TN (e = 12, d = -0.35), 0.83 and 3.73 days were obtained for EFF and RMSE respectively. Therefore, the results show that the use of the same model SM (12; -0.35) with tx = 60) to simulate the budbreak of both varieties is feasible without significantly lowering model performance.

In terms of the budbreak-flowering and flowering-veraison phenophases, the efficiency of the models calibrated for the different e coefficient changed in different ways, depending on the phenological stage and variety. Between budbreak and flowering, EFF and RMSE for EL9 (EN) and EL19 (TN) stays more or less constant within the whole range of the e coefficient. For the other phenological stages, the efficiency tends to increase (RMSE decreases) until the e values are between 10 and 12, stabilising for higher values (Figure 4c and Figure 4d). Between flowering and veraison, and for some phenological stages, EFF and RMSE also remain more or less constant within the whole range of the e coefficient. However, for certain phenological stages of EN (i.e., EL31, EL32, EL35 and 50V), the EFF of the models tends to decrease (RMSE tends to increase) for higher e values (Figure 4e and Figure 4f).

### Figure 4. Variation in the model’s efficiency (EFF) and root-mean-squared-error (RMSE) for optimal solution as a function of the e coefficient. a) and b) Budbreak, c) and d) phenological stages between budbreak and flowering, and e) and f) phenological stages between flowering and veraison, for Touriga Nacional (TN) and Encruzado (EN). For Budbreak (a and b) only the EFF and RMSE of models fitted with tx = 1st March are shown.

The threshold of the daily forcing rate sum, which corresponds to phenological stage occurrence (F*), depends on the forcing function parameter values and the period duration after the onset date for heat accumulation. In the model parameterisation procedure, the fitted F* values were higher for lower e and d parameters and for longer heat accumulation periods starting from tx (Table S3).

Based on the range of values for the d coefficient fitted in the first step of the procedure for models with e = 12, the effect on model performance was tested by using the two values for this parameter (d = -0,5 and d = -1). For this analysis, the performance metrics (EFF and RMSE) were compared and the differences with the best-fit model (dEFF and dRMSE) were determined. The results are presented as Supplementary Material (Table S3).

Compared to the best fit model (bfSM), the use of the model with e = 12 and d = -0.5 (SM (12; -0.5)) for all the phenological stages between budbreak and flowering resulted in a maximum EFF reduction (dEFF) of 0.03 (EL12, EL17 and EL19) and a maximum RMSE increase (dRMSE) of 0.45 (EL17) for EN. On the other hand, the use of the model with e = 12 and d = -1 (SM (12; -1)) resulted in, overall, a similar change in model performance when compared to the best-fit model (mean dRMSE = 0.24 for SM (12; -1) vs. mean dRMSE = 0.26 for SM (12; -0.5). Overall, for Touriga Nacional, SM (12; -0.5) showed better performance for the majority of the phenological states (EL9, EL17, EL19 and EL23) than SM (12; -1), and lower performance changes with respect to the bfSM (mean dRMSE = 0.18 vs. mean dRMSE = 0.33).

As was verified in the previous development phase for Touriga Nacional between flowering and veraison, SM (12; -0.5) also showed, globally and in each of the phenological stages, better performance (mean RMSE = 4.89 for SM (12; -0.5) vs. mean RMSE = 5.07 for SM (12; -1)). For Encruzado, the performance of both SM (12; -0.5) and SM (12; -1) was similar (mean RMSE = 5.02 vs. RMSE = 4.93 and mean dRMSE = 0.33 vs. mean dRMSE = 0.25). Therefore, these findings show that SM (12; -0.5) can be used to simulate all the phenological stages between budbreak and veraison of both target varieties, as the prediction performance was very close to those of the best-fit models; there is therefore a good trade-off between accuracy and simplicity.

### 4. Phenological model validation

By applying the previously described methodology, it was possible to establish the same set of parameters to simulate the phenological development (PDM) of both selected varieties (EN and TN). Although the models do not correspond to the best-fit models, the performance metrics are still globally high. RMSE ranged from 3.2 to 6.2 for TN, and from 3.9 to 6.8 for EN. For both varieties and in all phenological stages, RMSE was significantly lower than RMSEa (Figure 5) and the standard deviation of the phenological observations (Table S2). In terms of TN, the model efficiency was greater than 0.71 for all phenological stages (Figure 5). In order to take into account any model overfitting in the assessment of their performances, a leave-one-out cross-validation was applied. As expected, a slight increase in RMSEP was verified when compared to RMSE and for each phenological stage, although most of the values are still below seven days (Figure 5).

Figure 5. (a) Sigmoidal model response functions for phenological development simulation of cv. Encruzado (EN) and Touriga Nacional (TN) in the Dão wine Region. Efficiency coefficient (EFF), root-mean-squared-error (RMSE), Root Mean Square Error of Prediction (RMSEP) and Root Mean Square Error of Prediction with average (RMSEa) for each phenological stage and variety. (b) Critical value of the daily thermal forcing (F*) of each phenological stage for TN and EN.

The F* values for each phenological stage revealed slight differences in the phenological development of the two varieties. Overall, the heat requirements of Touriga Nacional was shown to by slightly lower between budbreak and flowering, but slightly higher between flowering and veraison.

Lastly, Figure 6 depicts the observed and simulated dates of all the stages together, using the selected phenological development models SM (12; -0.35) for EL4 and SM (12; -0.5) for EL9 to 50V. By adjusting the linear regression equation and applying the RMSE equation to the observed and simulated timings in all phenological stages (RMSEg), it can be concluded that the models have a slightly lower performance when predicting the phenological development of EN (RMSEg = 5.4 days) than TN (RMSEg = 5.0 days). These findings are nonetheless quite satisfactory, given the very low complexity of the model and its straightforward application using temperature records from local weather stations.

Figure 6. Scatterplots of the observed vs. simulated dates of all stages using the selected unified phenological development model for (a) Encruzado, (b) Touriga Nacional.

## Discussion

In this study, a phenological development model (PDM) for simulating several intermediate stages between budbreak and veraison for two representative winegrape varieties (cv. Touriga Nacional and cv. Encruzado) grown in the DãoWR was calibrated and validated. This was a thermal model, with the daily sum of the rate of forcing (R) calculated using a sigmoid function. For this purpose, a high-quality and comprehensive dataset which combines phenology data with weather station data in several vineyard sites spread over the target region (DãoWR) was used. To our knowledge, no previous study has been carried out for these two varieties and their intermediate phenological timings.

The performance metrics applied to the observed and simulated timings of all phenological stages (EFFg and RMSEg) show that PDM has good prediction accuracy (Figure 6), which is similar to that of other models developed to predict the whole sequence of phenological stages in the grapevine growing cycle (Molitor et al., 2014; Verdugo-Vásquez et al., 2017). Although the model showed an overall high performance, its simulation accuracy also depended on the phenological stage and variety (ranging from 0.71 (EL35) to 0.91 (EL27) for Touriga Nacional, and 0.58 (EL35) to 0.86 (EL9) for Encruzado). For flowering (EL23) and veraison (EL35), the performance of the model was similar to those obtained with SM for cv. Touriga Nacional by Costa et al. (2019) and for cv. Touriga Nacional and cv. Encruzado by Reis et al. (2020). For budbreak (EL4), the PDM performed better than those obtained by these authors.

The F* values for each phenological stage revealed slight differences between the phenological development of the two varieties: Touriga Nacional develops slightly earlier than Encruzado between budbreak and flowering, but slightly later between flowering and veraison. As expected, these results reflect the differences in observed dates of phenological stages between varieties. The differences in the duration observed between varieties for the phenophases EL4-EL23 and EL23-50V were indeed close to those observed in the Lisboa Wine Region (Lopes et al., 2008; Reis et al., 2020). Based on 29 years (1990 to 2018) of phenological observations in that region, no differences between varieties were found for the duration of the EL4-EL23 period and more or less 5 days for the EL23-50V period. However, both phenophases were longer in the Lisboa Wine Region.

Being a thermal model – different from a Chilling Model (CM), which describes the entire dormancy period (endo- and eco-dormancy) – our model only resolves the eco-dormancy period via the accumulation of forcing units under the assumption that the chilling requirements have already been satisfied (Chuine, 2000). The higher accuracy of CMs for budbreak timing prediction has been linked to the ability of these models to estimate the endo-dormancy release and to better account for the effect of temperature in this phase, thus avoiding the use of a fixed onset date for the beginning of the eco-dormancy period, as is the case in the TMs (Garcia de Cortázar-Atauri et al., 2009; Nendel, 2010). This is indeed a limitation of the model approach described in this study, as the beginning of heat accumulation for budbreak was fixed on 1st March, which does not take into consideration the fact that, depending on the year and site, different grapevine varieties might be in different phases of their development cycle when thermal summation begins. However, this approach was also used by Gutierrez et al. (1985). Leolini et al. (2020) demonstrated that the beginning of the eco-dormancy period was dependent on the fulfilment of specific physiological requirements, which can change with the climatic conditions of the different environments. In the present study, the SM model was found to have higher accuracy when 1st of March was used as a starting date to predict the budbreak timing of the grapevines varieties. These results suggest that although cultivated grapevines are a species with low chilling requirements (Londo and Johnson, 2014), the eco-dormancy period only starts at the end of winter in the study region. Another possibility is that eco-dormancy starts very early with chilling requirements already fulfilled, but without significant heat accumulation. The results of the present study reinforce the conclusion that the starting date of the TM depends on the physiological features of the grapevine variety and their interaction with the environment (Leolini et al., 2020). To overcome this limitation, models that include the endo-dormancy period would need to be tested in future studies.

Another practical limitation is the need for budbreak (EL4) or flowering (EL23) dates (input variables) before simulations can be implemented of intermediate states within the budbreak-flowering or flowering-veraison phenophases respectively. Due to their prediction error, when the simulated dates of EL4 and EL23 as tx are used, the uncertainty in the forecasts of the subsequent phenological stages can increase significantly, especially for the stages at the end of each phenophase; it is therefore recommended that they be determined by direct in situ observations. If the models are to be used as a vineyard management support tool, it will also be necessary to determine whether the accuracy of the models can be significantly improved via one-off readjustments in the F* values when an intermediate phenological dates are available. This also limits their direct application to the assessment of climate change effects on the phenology. In forthcoming research, it will be necessary to access their performance using simulated or fixed onset dates for heat accumulation.

The model developed in this study can be used to predict the dates of several phenological stages, and could therefore have an important role in promoting the sustainability of the wine industry in the Dão wine region. It can be used to better planning canopy management and supporting disease control. For example, the implementation of new strategies for crop protection requires advanced knowledge of phenological timings (Valdés-Gómez et al., 2017). On the other hand, having been calibrated and validated based on a combined phenology-weather dataset from different and representative sites in the Dão wine region, the model's application on a regional scale seems feasible.

Given that the present model does not include the berry ripening period, in forthcoming research it will be combined with specific models that simulate the evolution of the quality indicators used for planning the harvest of wine grapes (e.g., total soluble solids, titratable acidity, pH and phenolic compounds), as suggested by Ortega-Farías et al. (2002) and Fernández-González et al. (2013).

## Conclusions

A phenological development model (PDM) for predicting several intermediate stages between budbreak and veraison of two grapevine Portuguese varieties (Touriga Nacional and Encruzado) growing in the Dão wine region was developed in this study. A thermal model with a daily rate of forcing calculated with the sigmoid equation was used. The model revealed significant predictability, but it was dependent on the phenological stage and variety. In future studies, the model will be combined with specific models that simulate the evolution of winegrape berry quality indicators commonly used for harvest decision support. The relatively low complexity of the model makes it easy to use as a crop management and decision support tool. Similar approaches could be adopted in other wine regions worldwide; the present study is therefore an illustration of conceivable model developments under diverse environmental conditions.

## Acknowledgements

This research was funded by the operation nº CENTRO-04-3928-FEDER-000001 – Strategic Project to Support Wine Sector in the Central Region (Portugal), research activity: Comparative Study of grapevines varieties in the Dão wine Region, sponsored by CENTRO2020 and by National Funds through the FCT (Foundation for Science and Technology, I.P), within the scope of the project Refª UIDB/00681/2020. We also would like to thank the CERNAS Research Centre and the Polytechnic Institute of Viseu for their support. The study was also supported by FCT – Portuguese Foundation for Science and Technology, under the project UIDB/04033/2020, and 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 no. 810176.