Temperature-Based Grapevine Ripeness Modeling for cv. Touriga Nacional and Encruzado in the Dão Wine Region, Portugal

In the present study, we calibrated and validated thermal models to predict the DOY date at which the grape maturity index, potential alcohol/total acidity (PA/TA), reaches 0.75 (MS0.75), 1.0 (MS1), 1.5 (MS1.5), and 2.0 (MS2) for two grapevine Portuguese varieties, Touriga Nacional (TN) and Encruzado (EN), growing in the Dão wine region, Portugal. Daily rates of forcing calculated with the Sigmoid (SM) function and the Degree Day (DD) function were used. The outcomes show that the best performance of the models was obtained for the heat accumulation starting at flowering (tx = EL23). The analysis of model sensitivity to changes in forcing rate coefficients (T0, e, and d) enabled the selection of the same models for all maturity stage of each variety. The selected models revealed significant predictability, though dependent on the grape maturity stage and variety (EFF > 0.81 for TN and EFF > 0.75 for EN). The non-linear regression analyses of sugar concentration (SC) and total acidity (TA) with heat accumulation, calculated using the select models, demonstrated that a high fraction of SC and TA variance was explained by the variation of these temperature-based indices. Comparatively to SC and TA, the results highlight that the thermal conditions accumulated from flowering had a lower influence on pH juice variance.


Introduction
The grape berry quality at harvest mainly depends on the content of water, sugars, organic acids, amino acids, phenolic compounds, and aroma precursors [1]. Following fruit set, the increases in berry weight, volume, or diameter during berry development are typically characterized by a double sigmoid curve, resulting from two consecutive stages of rapid growth separated by a lag phase with slow or no growth [2][3][4]. During the first phase of growth, the sugar content remains low, while several organic acids are accumulated. The principal organic acids of grape berry are tartaric and malic acids, making up approximately 90% of total fruit acidity [3]. Tartaric acid is accumulated during the initial stages of berry development and the malic acid at the end of the first growth phase [2]. At the lag phase, the berries remain firm, while the chlorophyll content begins to decrease, and organic acid concentration reaches its highest level. The most dramatic changes in grape berry composition occur during the second growth phase, or the ripening phase [2,5]. During this phase, the berries soften, the sugar concentration increases reaching maximum levels, and the acid concentration decreases, whereas the compounds responsible for berry color, aroma and flavor are accumulated in the berry skin [4]. One of the main features of the grape ripening process is the accumulation of sugars in the form of glucose and fructose within the berries. The sugar accumulation also follows a sigmoidal pattern, with slow accumulation at the onset of veraison, followed by a rapid accumulation phase from just after veraison to several weeks onwards, then reaching a plateau phase [6]. At the end of the ripening period, the sugar content stabilizes, but sugar concentration may increase or decrease due to berry dehydration [2] or dilution [6]. The decrease in the concentration of tartaric acid verified during the ripening period can be mostly attributed to dilution effects. In contrast, reductions in malic acid concentration result from respiration and enzyme degradation besides dilution [3].
The stage of grape berry maturity has a marked effect on wine quality [7,8]. Consequently, the berry ripeness level at harvest is one of the most important parameters for obtaining high-quality wines [9]. Traditionally, the decision on harvesting is taken based on the sugar concentration (SC), titratable acidity (TA), and pH of the grape juice (technological maturity). However, these parameters only provide information concerning the pulp ripeness [9]. The use of these single parameters to define optimum wine grape maturity has not been very successful in this regard [8,10]. Several studies have been performed to isolate grape berry maturity indices associated with optimum quality. Du Plessis and Van Rooyen [8] found good results for the indices calculated by the product SC × pH or by the ratio SC/TA. Van Rooyen et al. [10] concluded that the range wherein maximum wine quality occurred was too wide in the case of the SC/TA index, while the SC × pH index gave a narrower optimum quality range and similar results for both cultivars.
Grape berry ripening and, consequently, its chemical composition at harvest, are influenced by the grapevine cultivar, agricultural management practices, pedoclimatic conditions, and their interactions. Cameron et al. [11] analyzed the trends in the rate of ripening of 24 cultivars over 20 years in four Victorian vineyard regions and verified significant differences among cultivars. Suter et al. [6], using data collected from 36 cultivars and over 7 years, verified that the variation in sugar accumulation traits was well explained by the cultivar, year, and their interaction. They verified that the sugar accumulation traits were affected by climate factors and grapevine water status, before or after midveraison, by mid-veraison date, and by berry weight at mid-veraison, though the relative importance of these factors varied significantly with cultivar. Using berry composition data at technological maturity for varieties Touriga Nacional and Tempranillo, from three Portuguese wine regions, Costa et al. [12] showed that summertime temperature (June to August) was negatively correlated to berry weight, titratable acidity, anthocyanins, and total phenols index, but was positively correlated to pH and potential alcohol. The influence of precipitation (April-June and July-September) was dependent on the location and variety. Fernández-González et al. [13] demonstrated that sugar accumulation in berries of several varieties was principally correlated with growing-degree days (GDD) calculated using a threshold temperature of 10 • C and accumulated from 1 January. Lucchetta et al. [14] investigated the effects of different pre-harvest canopy management techniques on berry ripening in cv. Sauvignon blanc and verified that the antitranspirant applications delayed ripening by 10 to 15 days, whereas the defoliation at veraison allowed postponement of the ripening process by about 5 days, without altering the sugar/acid ratio. Zheng et al. [15] showed that severe shoot trimming after fruit set could delay berry ripening, giving rise to a better organic acid composition by increasing the tartaric acid while reducing malic acid. Feller et al. [16] concluded that, although deficit irrigation has a strong influence on berry ripening, its positive effect on berry composition only occurred when high temperatures were not a limiting factor.
Although stage 38 on the modified Eichhorn and Lorenz scale of phenology [26] represents "Harvest, berries ripe", there is no universal definition for "berries ripe" [27]. The harvest is not a true reflection of growth or phenological stage, because the decision to undertake the vintage is also influenced by the aimed wine style and quality, the logistic requirements at the winery, labor force availability, as well as weather conditions or pest/disease pressures. Therefore, the models able to forecast the evolution of berry quality attributes, as a complement of other tools used in ripening monitoring, could be of utmost relevance to the vintage date decision and to promote more efficient winery management. These models may also be an important tools to the cultivars classification [6,27] used to support cultivar choice in response to the climate change [28].
The berry ripening dynamics and its relationship with climate conditions have been the aim of several studies. Sadras et al. [29] presented a quantitative model of accumulation of soluble solids that provide a baseline for comparison among varieties. Cortázar-Atauri et al. [30] proposed a model to simulate the dry matter growth based on thermal time and final potential dry weight. In this model, the water content dynamics (or • Brix, using the relationship obtained between them) was modeled defining two components: One related to the berry phenological stage and the other depending on the plant water status. In the Italian Vineyard Integrated Numerical (IVINE) model, the berry sugar content is estimated by its maximum value at harvest and a normalized coefficient, less than 1, calculated using a double sigmoid curve as a function of thermal time [31]. Fernández-González et al. [13] show that during the ripening period, the sugar accumulation of the cv. Treixadura and cv. Godello was mostly correlated with GDD (10 • C) and that the heat required to reach 20 • Brix ranged between 1132 and 1316 for Treixadura, and between 1059 and 1233 for Godello. Ortega-Farías et al. [23] found high R 2 for the non-linear regression between SC, TA, pH, and GDD (10 • C) calculated from budbreak. Based on a time series of sugar concentrations of 65 grapevine varieties, Parker et al. [27] demonstrated that temperaturebased models could be used to predict the time to target sugar concentrations from 170 to 220 g/L. Michelini et al. [28] demonstrated that the GDD-based model may be used to predict malic acid dynamics in cv. Pinot blanc and improved its accuracy with the addition of two other temperature-based parameters.
Along the previous lines, the present study aims to (i) calibrate and validate a linear (degree days-DD) and a non-linear (sigmoid model SM) temperature-based model, using time series of SC and TA data for the varieties cv. Touriga Nacional and cv. Encruzado, collected in several vineyards spread over the Dão wine region (DãoWR), central Portugal, to predict the date at which the grape maturity index (Potential alcohol/Total acidity-PA/TA) reaches 0.75, 1.0, 1.5, and 2.0; and (ii) investigate the relationship between the SC, TA, pH dynamics, and heat accumulation during the growing cycle. To our knowledge, no previous study was carried out for these two varieties and their grape maturity dynamics. Some of the results may, however, be extrapolated to other varieties and wine regions worldwide.

Study Area
The phenological and grape ripening data were collected from nine commercial vineyards and two varietal vineyards from six locations in the Dão Wine Region, Portugal, hereafter DãoWR ( Figure 1). This study was carried out with the cultivars Touriga Nacional (TN) and Encruzado (EN), red and white grapes, respectively. The DãoWR is located on a plateau in central-northern Portugal, surrounded by mountains, and protected from both the Atlantic moist winds and the continental influence from inner Iberia. Fraga et al. [32] categorized the region into three bioclimatic groups. The peripheral zone is characterized by a cold and humid climate, with cold nights, whereas the central area features a temperate and humid climate, with cold nights, and the south/southwestern area experiences a temperate and humid climate, with warm nights. The vineyard locations represent these three bioclimatic groups of the DãoWR. Overall, it is a wine region where small farms predominate, and vineyards are typically rainfed. The characteristics of the experimental plots (variety, plantation density, pruning system, and rootstock), as well as the available data periods, are shown in Table S1.

Phenological Data
Phenological stages were recorded through site observations, based on the modified Eichhorn-Lorenz scale [26] and until the beginning of veraison (EL35). For each varietal plot, the observations were carried out on all buds/shoots (until EL18), or inflorescences/clusters (after EF18), of six plants. At commercial plots (with three replicates per plot), the observations were made on two buds/shoots (or on its inflorescences/clusters) of five plants per replicate. Records were taken twice a week, until EL27, and once a week after this stage. All observations were undertaken by the same technicians, thus warranting a more homogeneous dataset. For each plant, in varietal plots, or for each replicate, in commercial plots, a phenological stage was evaluated when at least 50% of the observed organs reached this stage. From these phenological observations, the budbreak (EL4), flowering (EL23), and beginning of veraison (EL35) dates were determined for each plot/replicate and year.

Berry Quality Monitoring and Maturity Stages
The berry-ripening process was monitored weekly from the veraison to the harvest. Three samples (one from each replicate) of 200 berries were sampled from 30 plants at each commercial plot. At varietal plots, one sample of 100 berries was sampled from 6 plants. The berries' sampling was made alternately collecting three berries at the top, the middle, and the bottom of clusters located at both sides of the plants. On the sampling day, each berry sample was weighted, crushed at constant pressure (hydraulic grape pressing machine 1040 cm 3 , Agro-Moderna, Vila Nova de Famalicão, Portugal), and the grape must volume was recorded. The must sample was immediately filtered, being the sugar concentration (SC), the titratable acidity (TA), and the pH subsequently analyzed. The SC was measured with a refractometer (HANNA Inst. HI996813, Woonsocket, Rhode Island, USA) and subsequently converted to potential alcohol (PA). Immediately before SC measurements, the must temperature was recorded. The TA was measured by titration with 0.1 M NaOH of a dilute solution of 10 mL of must with 25 mL of distilled water, using bromothymol blue as a color indicator. Results were expressed in equivalent of tartaric acid content (g/L). The pH was measured with a potentiometer (Metrohm, 691 pH meter, Herisau, Switzerland). All laboratory determinations were made with three replications. For each sample, the grape maturation index (MI) was calculated as the ratio between the Potential Alcohol (PA) and the total acidity (TA), i.e., MI = PA/TA. The number of berry sampling days of each year and plot, and the total number of berry samples per plot are shown in Table S2.
Each time series (berry attribute measurements of each replicate/plot per year) of both cultivars was assessed using the following criteria: (1) After veraison, the sugar content in berries increases until reaching a plateau at approximately the maximum content, and (2) the sugar content stabilizes at the end of the ripening period, though its concentration may increase due to berry dehydration [2]. To avoid using data where sugar concentrations may result from dehydration rather than physiological ripening, data values were excluded when SC suddenly increased after a plateau of several days. Furthermore, if on a given sampling date, the SC decreased compared to the preceding sample date, the time series period from this value onwards was eliminated [27].
Taking into account the range of MI for both cultivars (Figure 2), four grape maturity stages (MS i ) corresponding to MI equal to 0.75 (MS 0.75 ), 1.0 (MS 1 ), 1.5 (MS 1.5 ) and 2.0 (MS 2 ) were pre-defined. For each time series (plot or replicate × year × cultivar), the day of year (DOY) on which the MI reached each MS i was determined as follows: (1) If on any sampling date, the MI was within the range [i − 0.025; i + 0.025], this date was considered as the DOY that this MS i was reached; (2) alternatively, the DOY was determined by a two-point linear interpolation of the data points on either side of the MS i . The number of data points, the mean, and standard deviation of PA, SC, and TA corresponding to each MS i and cultivar is also presented in Figure 2.

Temperature Dataset
Weather stations are located in the vineyards, very close to the selected plots (<100 m), except for CP6 and CP9, for which the weather stations (WS5 and WS6) are located at 1.3 km and 3.7 km apart from the plots, respectively. In these two cases, however, owing to the relatively flat terrain, the weather stations remain representative of the atmospheric conditions in the corresponding vineyard plots. A preliminary quality checking was carried out in the temperature time series at each location. For all weather stations, data gaps correspond to less than 4% of the entire daily temperature time series length. These gaps were filled in 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). Therefore, the phenological and berry ripening data are complemented by local meteorological data, which is a very important feature to warrant the representativeness of the analysis.

Models Selection and Evaluation
Thermal models were herein adjusted to predict the DOY to reach the grape maturity stages MS 0.75 , MS 1 , MS 1.5 , and MS 2 . These models have been successfully applied as phenological models [17,19,20,22] and were also used by Parker et al. [27] to predict the time to target sugar concentrations between 170 and 220 g/L for a wide range of varieties. Thermal models only consider the effect of forcing temperatures and assume that a given grape maturity stage (MSi) occurs when the daily accumulation of the daily rate of forcing (RF), from the onset date (t x ), reach a specific critical value (F): Two daily rates of forcing types were tested: (a) The Degree Day function (DD) that determine the RF through the daily temperature (T i ) and a fitted temperature threshold (T 0 ), often referred to as base temperature: The Sigmoid function (SM) that calculates the RF from the T i and two fit parameters, e and d, which correspond to the location and sharpness of the curve, respectively: In order to find the timing within the grapevine development cycle from which the temperature has a more significant effect on the ripening process, the models were fitted considering the three main phenological stages, namely budbreak (EL4), flowering (EL23), and beginning of veraison (EL35) as the start date for heat accumulation (t x ), separately. In contrast to the options of setting a date for t x or assuming t x as an adjustment parameter, the use of a phenological stage as t x has the advantage of allowing a better adaptation of the ripening dynamics modeling to the previous stages of the development cycle, also determined by the weather conditions of each year and site.
For each MSi of each cultivar, the Phenology Modeling Platform (PMP), version 5.5 [33], which estimates the model parameters by the Metropolis annealing algorithm, was used to fit the most accurate model. In the model-fitting process, the PMP estimated, for the DD model, the T 0 coefficient (base temperature) and the critical value (F), whereas it estimated the F, e, and d coefficients for the SM model. For the DD model, a sensitivity analysis was performed to assess the model responsiveness to a changing T 0 to find a single parameter that could be used in all MSi. For the SM model, this sensitivity analysis was performed in a two-step approach: In the first step, the coefficient e (location parameter) was fixed in successive integer values, between a lower and an upper threshold. The remaining parameters (d and F) for each fixed e value were subsequently automatically adjusted. In the second step, for the selected e coefficient from the analysis of the first step results, the sensitivity of the models (variation of the performance) to changes in the d coefficient was analyzed.
The root-mean-squared-error (RMSE) and the Nash-Sutcliffe efficiency coefficient (EFF) were used to assess the model performance: where O j represents the observed MSi date, P j is the corresponding simulated values, O m is the average of all observed values, and n is the sample size. The model performance was also indirectly assessed through non-linear and linear regression analyses (R, R 2 , adjusted-R 2 , Std. Error of the Estimate, and p-value) between SC, TA, and pH of each plot and variety, as well as heat accumulation calculated using the select models.
To take into account potential model overfitting, the calibration was 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 pair of grape maturity 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 follows: where P v j is the predicted value obtained from a leave-one-out cross-validation approach.

Phenological, Ripening and Climate Data
For the three main phenological dates, budburst (EL4), flowering (EL23), and beginning of veraison (EL35), the cv. Encruzado tends to reach them slightly earlier than the cv. Touriga Nacional. The average differences between dates were 1 day at EL4 and 2 days at EL23 and EL35. Based on the median, these differences were slightly higher at budbreak (3 days) and veraison (6 days). The interquartile ranges (IQR) show higher variability of observed dates for Encruzado at EL4 and EL23 but lower at EL35 (Table S3).  As expected, the phenological variability was the result of the thermal conditions in each year and plot (Table 1). Before budbreak, in 2017 and 2019, the daily mean temperature was about 10 °C, while in 2016 and 2018, it was approximately 9 °C. During this period, at WS5 and WS6 the mean daily temperature was barely 1 °C higher than at WS1. Between EL4 and EL19, 2017 and 2019 were the years with the highest and the lowest mean temperature, respectively (16 °C vs. 14.3 °C). During the flowering period (EL19-EL27), the mean temperature ranged from 15.6 °C in 2014 to 19.8 °C in 2015. For these two phenophases, the highest temperature was recorded at WS6. For the fruit development period (EL27-EL35) and the ripening period, the highest temperature was recorded in 2016 and 2018, respectively. During the veraison and ripening period (EL35-harvest) the temperature recorded at W2 was, on average, 1.8 °C higher than those As expected, the phenological variability was the result of the thermal conditions in each year and plot (Table 1). Before budbreak, in 2017 and 2019, the daily mean temperature was about 10 • C, while in 2016 and 2018, it was approximately 9 • C. During this period, at WS5 and WS6 the mean daily temperature was barely 1 • C higher than at WS1. Between EL4 and EL19, 2017 and 2019 were the years with the highest and the lowest mean temperature, respectively (16 • C vs. 14.3 • C). During the flowering period (EL19-EL27), the mean temperature ranged from 15.6 • C in 2014 to 19.8 • C in 2015. For these two phenophases, the highest temperature was recorded at WS6. For the fruit development period (EL27-EL35) and the ripening period, the highest temperature was recorded in 2016 and 2018, respectively. During the veraison and ripening period (EL35-harvest) the temperature recorded at W2 was, on average, 1.8 • C higher than those recorded at WS1. Considering the growing season, the warmest years were 2017 and 2018 (19.8 • C and 19.9 • C, respectively), being the highest temperature recorded at WS2 and WS6.
At harvest, the SC (  (Table 1). Furthermore, the average TA reached the lowest value for EN in 2018, and the highest for EN in 2017 and TN in 2019. Table 1. Average and standard deviation of daily mean air temperature before budbreak (1 January to EL4), for the shoot and inflorescences development (E4 to EL19), flowering (EL19 to EL27), berry development (EL27 to EL35), and ripening (EL35 to harvest) phenophases, as well as for the growing season length, GSL (EL4 to harvest), recorded in each (a) year and (b) weather station.   Table 2 shows the model coefficients and the performance statistic of the best-fit models (with all forcing rate parameters fitted) to simulate the DOY of the four grape maturity stages (MSi). Except for MS 0.75 of TN variety, in which the better performance was obtained with the phenological stage EL35 as onset date for heat accumulation (t x ), all of the others MSi reveal better performance with t x = EL23. The models with highest EFF are shown in the grey boxes.

Model Selection
To identify a single parameter that could be used for all MSi, a sensitivity analysis was performed to assess the model responsiveness to changes in the T 0 and e coefficients. The variation in the global efficiency (EFFg, mean of the model EFF for the four maturity stages) as a function of those coefficients ( Figure 4) shows that for TN the best global performance was obtained with e = 15 using SM, and T 0 = 2 using DD, in both cases with t x = EL23. For EN, the best global performance was obtained with e = 13 for SM, and with T 0 = 0 for DD, also with t x = EL23.
Based on the range values of the d coefficient fitted to SM with e = 15 (TN variety) and e = 13 (EN variety) (Table S4)

Model Validation
The methodology described above allowed us to establish the same simulation models (DD and SM models) for the technical grape ripeness (TGR) of the two selected varieties (EN and TN). Although it does not correspond to the previously obtained best-fit models, the performance metrics are still globally high ( Figure 5). Using the SM model, the RMSE ranged between 3.7 (MS 0.75 ) and 5.4 (MS 2 ) for TN, and between 3.5 (MS 0.75 ) and 5.4 (MS 2 ) for EN. Overall, for both varieties, the DD showed slightly better performance than SM, with the RMSE ranging from 3.5 (MS 0.75 ) to 5.0 (MS 2 ) for TN and from 3.4 (MS 0.75 ) to 5.3 (MS 2 ) for EN. For both varieties and in all grape maturity stages, the RMSE of both models was significantly lower than the standard deviation of the observed DOY of each MS i . The EFF of SM is greater than 0.81 for TN and greater than 0.75 for EN. On the other hand, the EFF of DD is greater than 0.84 for TN and greater than 0.76 for EN. The results also show that the accuracy of both models decreases as grape maturity advances. Compared with the results obtained by Parker et al. [27], the selected models of our study presented a significantly better performance.
In order to take into account the model overfitting in the assessment of their performances, a leave-one-out cross-validation was applied. As expected, a slight increase of the RMSEP was verified when compared to the RMSE and for each MS i , though most of the values are still low and below six days ( Figure 5).
The model performance was also indirectly assessed through non-linear and linear regression analyses (R, R 2 , adj-R 2 , SEE, and p-value) between SC, TA, and pH of each plot and variety, and heat accumulation (HA) calculated using the select models ( Figure 6 and Table S6). In general, the results show that a high fraction of SC and TA variance was explained by the variation of HA calculated by the two models. Considering all plots and SC, the adj-R 2 ranged between 0.79 for EN, with HA calculated by DD Several studies also highlighted the key role played by thermal effects on these berry quality attributes during or at the end of the ripening period [12,23,28,[34][35][36][37]. On the other hand, by comparing the results of each plot, the relationships between these grape quality parameters and the thermal conditions are site dependent. These findings demonstrate that, depending on the location, the other forcing factors have different relative importance on the berry ripening dynamics.

Discussion
In this study, the DD and SM models were adjusted, calibrated, and validated to predict the timing of four grape maturation stages (MS i ) of two representative winegrape varieties (cv. Touriga Nacional and cv. Encruzado) grown in the DãoWR. For this purpose, a high-quality and comprehensive dataset, which combines phenology data, time-series of SC and TA data with weather station data in several vineyard sites spread over the target region (DãoWR), was used. Quality checking of the raw data revealed their high consistency and homogeneity, thus warranting the representativeness and robustness of the outcomes of the present study. The four grape maturity stages (MS i ) were defined as being the correspondent to the PA and TA ratio (MI) equal to 0.75 (MS 0.75 ), 1.0 (MS 1 ), 1.5 (MS 1.5 ), and 2.0 (MS 2 ). The models were fitted considering the three main phenological stages, namely budbreak (EL4), flowering (EL23), and beginning of veraison (EL35) as the start date for heat accumulation (t x ), separately. To our knowledge, no previous study was carried out for these two varieties and their grape ripening dynamics.
The results show that the best performance of the models was obtained for the heat accumulation starting at flowering, i.e., t x = EL23 (Table 2 and Figure 4), demonstrating that compared to the whole vegetative cycle (after budbreak) or the period after veraison, the period after flowering was the one for which the thermal conditions have a more significant effect on the technical ripening process. This finding is not totally in agreement with Parker et al. [27], who found better performances of DD and SM models with t x = 91 and t x = 87, respectively. However, Barnuud et al. [34] verified the strength of the associations between TA at grape maturity of 22 • Brix for cv. Cabernet Sauvignon and Shiraz, and climate variables increase steadily from the early growing season towards the berry ripening period, indicating that the thermal conditions during the ripening period are more influential in determining berry composition at maturity. On the other hand, Jones and Davis [36] found a higher correlation of SC and TA with flowering and veraison dates than with budbreak dates. Assessing the relationship between interannual variability in atmospheric conditions on grape berry quality attributes (cv. Touriga Nacional and cv. Aragonez/Tempranillo) in three Portuguese wine regions (Douro, Dão and Alentejo), Costa et al. [12] verified that temperatures from June to August (therefore after flowering) were more influential in determining grape berry composition at harvest. The analysis of model sensitivity to changes in forcing rate coefficients (T 0 , e, and d) allowed us to establish some simulation models for all MS i of each variety. Regarding the best-fit model, the use of selected models (DD [2;EL23] (Figures 2 and 5). As expected, the cross-validation performance metric (RMSEP) slightly increases compared to the RMSE (differences of up to 0.11 for TN and up to 0.95 for EN). Compared with the results obtained by Parker et al. [27], the selected models of our study presented a significantly better performance. For all MS i of TN, the DD performed only marginally better than the SM (differences of up to 0.03 in EFF). For EN, the DD performed marginally better than SM for MS 0.75 and MS 2 , but worst for MS 1 and MS 1.5 (differences of up to 0.01 in EF). Parker et al. [27]) also found slight differences in the performance of these models (up to 0.03 in EF). Comparing the main RMSE or RMSEP of all MS i , both DD and SM performed better for TN than EN.
Other environmental factors, as well as cultivar and management practices, are also known to influence berry ripening dynamics [6,[12][13][14][15][16]38]. The reduction of performance verified in both models and varieties from MS 0.75 to MS 2 suggests that other environmental factors may gain more relevance as ripening progresses. This finding suggests that integrating other environmental factors into the models, such as grapevine water status and precipitation during the ripening period, may be important to improve its accuracy, mainly at the end of ripening. Costa et al. [12] demonstrated that precipitation from April to June (before flowering) had a significant negative effect on potential alcohol at harvest of TN in the DãoW region. Jones and Davis [36] found a similar effect of precipitation during flowering and the ripening period on the sugar concentration of Cabernet Sauvignon and Merlot varieties. Michelini et al. [28] concluded that the DD-based model to predict the malic acid dynamics improved its performance adding other temperature-based parameters.
One of the main features of the grape-ripening process is the accumulation of sugars in the form of glucose and fructose within the cellular medium [39], which are converted into alcohol by the action of the yeast during the fermentation process, releasing carbon dioxide and heat [2]. Besides sugar content, the acid content of grape juice is also a critical quality parameter since it affects wine quality [2]. Berry sugar accumulation and acid degradation are influenced by climatic variables, such as air temperature [27,28]. The non-linear regression analyses of SC and TA with heat accumulation (HA) calculated using the select models demonstrated that a high fraction of SC and TA variance was explained by the variation of these temperature-based indices (Table S6). Considering each plot, the adj-R 2 of the SC-HA relationship ranged from 0.8 to 0.92 for TN and from 0.67 to 0.92 for EN. For the TA-HA relationship, adj-R 2 ranged from 0.81 to 0.96 for TN and from 0.77 to 0.90 for EN. This finding suggests that the temperature effect on sugar accumulation and acid degradation dynamics was site dependent. Comparing the Grapevine Sugar Ripeness (GSR with T 0 = 0 • C), proposed by Parker et al. [27], with the Growing Degree Day (GDD with T 0 = 10 • C), Michelini et al. [28] obtained an R 2 of 0.77 and 0.57, respectively, for the regression with sugar concentration. For the regression of the malic acid content, they obtained R 2 = 0.44 and R 2 = 0.66, respectively, and concluded that the malic acid degradation was influenced by site elevation. With data of one vintage, Ortega-Farías et al. [23] found R 2 of 0.96 (Chardonnay) and 0.99 (Cabernet Sauvignon) for the nonlinear regression between SC and GDD, and R 2 of 0.88 (Chardonnay) and 0.92 (Cabernet Sauvignon) for the regression between TA and GDD.
The influence of climate on grape ripening has also been investigated through the regression analyses between some berry composition attributes at harvest and weather conditions in the different periods of the vegetative cycle. Analyzing climate influences on berry composition at maturity (TTS = 22 • Brix) of Cabernet Sauvignon, Chardonnay, and Shiraz, Barnuud et al. [34] verified that only 59% to 63% of the variations in berry TA were described using generic models based on two temperature-based variables. However, the cultivar-specific models explained significantly higher proportions of the variation in TA (adj-R 2 from 0.7 to 0.82). In the regression models proposed by Jones and Davis [36], the TA variability of Cabernet Sauvignon was explained by rainfall during flowering (positive effect) and by the relative number of days with temperatures greater than 30 • C (negative effect) during flowering and veraison (adj-R 2 = 0.66). In turn, the potential evapotranspiration during the flowering and the relative number of days with temperatures greater than 25 • C (negative effect) during veraison explained the TA variability of Merlot (adj-R 2 = 0.77). Concerning the sugar concentration at harvest, the regression models showed that climatic variables, including the number of days with temperatures greater than 25 • C for cv. Merlot or 30 • C for cv. Cabernet Sauvignon (positive effect), collectively explained 68% of its variability for Cabernet Sauvignon and 79% for Merlot.
Comparatively to the SC and TA, the linear regression analyses between pH and HA highlighted that the thermal conditions from flowering had a lower influence on its variance (adj-R 2 = 0.56, for TN, and adj-R 2 = 0.63, for EN). Barnuud et al. [34] verified that a generic model containing growing season degree days and air vapor pressure deficit (VPD) described only 50% of the pH variations. They concluded that juice pH at maturity was generally weakly associated with the climate variables rather than anthocyanins concentration and TA. The juice pH measures the hydrogen ion concentration in the berry and is related to total acidity [3]. Although lower acidity levels are usually correlated with higher grape pH, this relationship is affected by potassium accumulation, which depends on the temperature itself [40]. With the increasing of juice pH verified during the ripening being the result of the combined effect of the potassium accumulation and the total acidity reduction, factors other than temperature may play a key role on its dynamics.
A practical limitation of the models developed in this study is the need for flowering (EL23) dates (input variables) before implementing simulations of maturity states. Due to the prediction error of phenological models, using flowering simulated dates as t x may increase the uncertainly in the maturity stages forecasting, mainly at the end of ripening. Thus, it is recommendable that they are determined by direct in situ observations. In forthcoming research, it will be necessary to assess their performance using simulated onset dates for heat accumulation. The use of the flowering date has the advantage of allowing a better adaptation of the ripening dynamics modeling to the previous stages of the development cycle, also determined by the local weather conditions.
The models developed in this study that allow predicting the date at which four grape maturity stages (MS 0.75 , MS 1 , MS 1.5 , and MS 2 ) are reached can be a helpful tool to support the vintage date decision and to promote more efficient winery management, as a complement of ripening monitoring. On the other hand, having been calibrated and validated based on a combined phenology-ripening-weather dataset from different and representative sites in the Dão wine region, its application at a regional scale seems feasible.

Conclusions
A thermal model to predict the DOY to reach the grape maturity stage MS 0.75 , MS 1 , MS 1.5 , and MS 2 of two grapevine Portuguese varieties (Touriga Nacional and Encruzado), growing in the Dão wine region, were herein developed. Daily rates of forcing calculated with the Sigmoid function (SM) and the Degree Day function (DD) were used. The outcomes show that the best performance of the models was obtained for the heat accumulation starting at flowering (t x = EL23). The analysis of model sensitivity to changes in forcing rate coefficients (T 0 , e, and d) enabled the selection of the same models for all MS i of each variety. The selected models revealed significant predictability, though dependent on the grape maturity stage and variety. The non-linear regression analyses of SC and TA with heat accumulation, calculated using the select models, demonstrated that a high fraction of SC and TA variance was explained by the variation of these temperature-based indices. Comparatively to SC and TA, the results highlight that the thermal conditions accumulated from flowering had a lower influence on pH juice variance. The relatively low complexity of the models facilitates their implementation as a tool to support the vintage date decision and to promote more efficient winery management, as a complement of ripening monitoring. Similar approaches can be adopted in other wine regions worldwide, thus being the present study an illustration of conceivable model developments under diverse environmental conditions. Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy11091777/s1, Table S1. Geographical location of the selected plots (latitude and longitude), grown variety, rootstock, planting density, pruning system, weather station identifier, and available period of phonological and ripening data, Table S2. The number of berry sampling days of each year and plot, and total number of berry samples per plot, Table S3. Descriptive statistics of observed dates for main phenological stage of Touriga Nacional and Encruzado varieties, Table S4. Efficiency coefficient (EFF) and the root-mean-squared-error (RMSE) of the best-fit SM, SM with t x = EL23, e = 15 and adjusted d, SM with t x = EL23, e = 15 and d = −0.4 for all maturity stages of cv. Touriga Nacional (TN) and of best-fit SM, SM with t x = EL23, e = 13 and adjusted d, SM with t x = EL23, e = 13 and d = −0.35 of cv. Encruzado (EN). dEEF and dRMSE are the differences of EFF and RMSE with respect to the best-fit model, Table S5. Efficiency coefficient (EFF) and the root-mean-squared-error (RMSE) of the best-fit DD, DD with t x = EL23, and adjusted T 0 , GDD with t x = EL23, and T 0 = 2, for all maturity stages of cv. Touriga Nacional (TN) and of the best-fit DD, DD with t x = EL23, and adjusted T 0 , GDD with t x = EL23, and T 0 = 0, for all maturity stages of cv. Encruzado (EN). dEEF and dRMSE are the differences of EFF and RMSE with respect to the best-fit model, Table S6. Statistical quality indicators of the non-linear and linear regression between Sugar content, Total acidity, pH and heat accumulation using SM and GDD models for each plot and variety.
Funding: 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 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. The study was also supported by National Funds through the FCT (Foundation for Science and Technology, I.P), under the projects Refª UIDB/00681/2020 and Refª UIDB/04033/2020.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.