Spatiotemporal Estimation of the Olive and Vine Cultivations’ Growing Degree Days in the Balkans Region

: Olive and vine cultivations are two of the most important crops in Europe, yielding high quality and value food products. The climate change over the Balkans may elevate the agroecological pressure for the established crops and shift their cultivations areas. One of the widely-used agroclimatic indices is the growing degree days (GDD) which accumulates the necessary thermal units for the selected crops. Despite the advances on the agroclimatic research, there are few available methods for spatiotemporal estimation of this useful index. So, this research is focused on the construction of simple and reliable equations for the calculation and projection of olive and vine cultivations’ GDD over the Balkans. The models’ input parameters are the time, the altitude, the distance from the seashore, and the latitude. Its assembly is made by the extracted spatial data, combined with the Agri4Cast dataset for the period of 1980 to 2018 incorporating the regional climate change trend. The results indicate that the most inﬂuential parameter is the time, followed by the latitude, for both cultivations. According to the projections, as quantiﬁed by GDD, a vast sprawl of olive and vine cultivation areas will have been formed to the northern parts of the studied area. To be more precise, the viticulture could expand spatially by 28.8% (of the Balkans area) by 2040, and by 15.1% to 2060, when the olive cultivations’ area could sprawl 23.9% by 2040 and 20.3% by 2060.


Introduction
Numerous studies have thoroughly investigated the relation of agriculture and climate, since the atmospheric conditions are one of the most critical crops' growth factors [1][2][3]. Climate change implications about crops' productivity are not straightforward. On the other hand, the variation of parameters such as temperature, relative humidity, and precipitation may restrict their growth and productivity [4][5][6][7]. Therefore, the high spatial variability of the climatic parameters drives the necessity of a regionally-focused and cropspecialized research. In order to connect the atmospheric with the agricultural environment in a climatic timescale, a wide variety of agroclimatic indices have been coincided [8][9][10]. Frequently, these indices quantify the thermal environment related to the growing and the phenological stages of the crops, and that is why they are of utmost importance on several scales of the agricultural sector, from yield prediction to policymaking [11][12][13]. One of the most widely-used agroclimatic indexes is the growing degree days (GDD), because it incorporates the air temperature, a common measurement, and an effective atmospheric parameter. The conceptual idea was presented in 1730 when Reaumut introduced the concept of heat units, or thermal time [11]. The accumulated heat was found to be a crucial

Study Area and Data
This research is focused on the Balkans area covering the territories of Albania, Bosnia and Herzegovina, Bulgaria, Croatia, Greece, Montenegro, North Macedonia, Romania, Serbia, and Slovenia which are presented in Figure 1. There are two essential criteria that led to the study of the Balkans region. The first is the fact that these countries have a high potential for agricultural productivity, and the other is that they are located in a transitional zone in terms of climatic change.
For this research, a part of the gridded dataset produced by MARS-AGRI4CAST [46] was utilized. The dataset contains meteorological parameters from weather stations interpolated on a 25 km × 25 km regular meteorological grid over Europe for a time period from 1980 to 2018, including parameters such as air temperature, relative humidity, wind speed, precipitation, and radiation. The AGRI4CAST dataset is well-documented and widely used for agricultural research because of its accuracy and consistency in a spatiotemporal point of view. Its characteristics are presented in Table 1 [47][48][49]. The spatial and the temporal dimensions of the dataset are capable of incorporating the climatic trends of the selected area [50,51]. For this research, a part of the gridded dataset produced by MARS-AGRI4CAST [46] was utilized. The dataset contains meteorological parameters from weather stations interpolated on a 25 × 25 km regular meteorological grid over Europe for a time period from 1980 to 2018, including parameters such as air temperature, relative humidity, wind speed, precipitation, and radiation. The AGRI4CAST dataset is well-documented and widely used for agricultural research because of its accuracy and consistency in a spatiotemporal point of view. Its characteristics are presented in Table 1 [47][48][49]. The spatial and the temporal dimensions of the dataset are capable of incorporating the climatic trends of the selected area [50,51].    Table 2 provides us with some information about the number of grid points inside agricultural areas according to Corine land cover classification nomenclature [52], as a measure of the cultivated land per country. The study area's sampling points used for this research are 1532, covering every area no matter if it is agricultural or not. So, the initial dataset is exceptionally large, with 24,306,653 rows of measurements with ten variables' columns.
The altitude information was extracted by a topographical delineation of the study area which is based on a 30-m digital elevation model (DEM) acquired from the Shuttle Radar Topography Mission (SRTM) dataset, provided by the United States Geological Survey (USGS) EROS archive. The SRTM 1 arc-second global elevation data were processed from raw C-radar signals spaced at intervals of 1 arc-second, approximately 30 m, at NASA's Jet Propulsion Laboratory (JPL). The SRTM mission was formed by using interferometric radar, which compares two radar images or signals taken at slightly different angles. This mission used single-pass interferometry, which acquired two signals at the same time by using two different radar antennas. An antenna located onboard the space shuttle collected one dataset, and the other dataset was collected by an antenna located at the end of a 60-m mast that extended from the shuttle. Differences between the two signals allowed for the calculation of surface elevation [53][54][55]. Moreover, a raster of the latitude and a raster of the Euclidean distance from the seashore were used.

Methods
It is important to point out that due to the large dataset, all the processes and analyses were carried out using the R-Language, which is suitable for such research [48,56]. More specifically, R v.4.0.3 core scripts and essential packages, such as the "dplyr" [57] for data handling and manipulation, "purrr" [58] and "broom" [59] for the multiple regression modeling, along with the "fst" [60] which is very fast for reading and writing large data frames, were utilized. Additionally, part of the spatial analysis and the mapping procedures were performed with the packages "raster" [61], "rgdal" [62], and "tmap" [63], along with the open-source GIS software, QGIS v2.18.3 [64]. This research's general workflow is presented in Figure 2, depicting the major datasets, relationships, and tools.
The whole data handling and analysis process became fully automated via the scripts, except the map drawing, so the method is partially reproducible. The above algorithm will be utilized in future research for estimation of agroclimatic indices focused on different cultivations. The whole data handling and analysis process became fully automated via the scripts, except the map drawing, so the method is partially reproducible. The above algorithm will be utilized in future research for estimation of agroclimatic indices focused on different cultivations.

Growing Degree Days
To calculate the growing degree days, the canonical form of the equation was used, which is: where Tmax is the daily maximum air temperature (°C), Tmin is the daily minimum air temperature (°C), and Τbase is the air temperature threshold below that which the plants' growth does not progress. When ( T max +T min 2 ) < T base , the GDD is set equal to 0. Besides, Τbase depends on the biology of the species and the cultivars [16,17]. Also, the period of the year for the accumulation of the GDD depends on the plant's requirements, and the biological phase was examined. For the examination of the thermal requirements of the olive tree (Olea Europea L.) over the study area, as a threshold, the Τbase = 7 °C was used. Moreover, the necessary heat summation limit, which is relevant for flower structure development, was 700 GDD for the period of January to the end of May [21,65,66]. Accordingly, for the viticulture (Vitis Vinifera L.), we used the Τbase = 10 °C as a threshold for the calculation of GDD units

Growing Degree Days
To calculate the growing degree days, the canonical form of the equation was used, which is: where T max is the daily maximum air temperature ( • C), T min is the daily minimum air temperature ( • C), and T base is the air temperature threshold below that which the plants' growth does not progress.
When T max +T min 2 < T base , the GDD is set equal to 0. Besides, T base depends on the biology of the species and the cultivars [16,17]. Also, the period of the year for the accumulation of the GDD depends on the plant's requirements, and the biological phase was examined. For the examination of the thermal requirements of the olive tree (Olea Europea L.) over the study area, as a threshold, the T base = 7 • C was used. Moreover, the necessary heat summation limit, which is relevant for flower structure development, was 700 GDD for the period of January to the end of May [21,65,66]. Accordingly, for the viticulture (Vitis Vinifera L.), we used the T base = 10 • C as a threshold for the calculation of GDD units for the period from April to the end of October each year. The required GDD to characterize a climatic region as suitable for this cultivation is 2000 units [67,68].

Multiple Linear Regression and Cross-Validation
The primary scope of this research is the investigation of the climatic potential of the Balkans to establish extended viticulture and olive trees cultivation on a broader area, Atmosphere 2021, 12, 148 6 of 14 far beyond the southern coastal regions. To achieve this, the simple and robust method of multiple linear regression was selected. The foremost reason for this choice is the straightforward interpretation of the results and the physical explanation of the input parameters (predictors). Moreover, the equation of this type of analysis can be efficiently utilized to calculate the GDD over a selected point and year. So, the equations can be a useful tool for scientists, farmers, and other professionals. The form of the multiple linear equations is the following (Equation (2)): where GDD represents the value of the GDD units (dependent variable) at a certain point; b 0 is constant, and b 1 to b 4 are the coefficients of each independent variable from every multiple regression. alt is the altitude (m), lat is the latitude (degrees), dist is the distance from the shoreline (km), and time is the number of the year. So, the equation contains time-invariant factors such as altitude, latitude, and the shoreline distance and the year (time-variant) factor, to achieve spatiotemporal, inter-and extrapolation. In order to assess the accuracy of the regression model, a k(10)-fold cross-validation test using the Caret R package was performed [69].
After the multiple linear regression tests and the formation of the equations, the raster calculations were made to visualize the spatiotemporal projections of the GDD over the Balkans area. To be precise, the maps of the regions which fulfil the criterion of GDD for each cultivation were drawn. The multiple regression approach is essential for studies on the agroclimatic conditions with GDD or other related indices [32,33,70,71].

Results and Discussion
The results of this research consist of the output of the multiple linear regression (MLR), which are the equations, their coefficients, and their statistics test. The equations by themselves are useful information to estimate the GDD for each point of the selected area and for every year. On the other hand, the interpretation of their statistics tests and coefficients gives us useful insight into the significance of the predictor's variables.
Moreover, it was presented a combined map for each cultivation with the spatiotemporal pattern of the areas that accumulate the GDD units above the threshold for viticulture and olive tree cultivation over the Balkans region. This cartographic material is a tool to identify the potential climatic sustainability for the expansion of vineyards and olive growths.
In case of implementation of the following models, altitude is in meters (m), the latitude is in decimal degrees (DD), the distance from the shoreline in kilometers (km), and the time is the year (y).

Olive Tree Cultivation Multiple Regression Model Results and Spatiotemporal Distribution of the GDD
In total, the multiple linear regression (MLR) model achieved a high R 2 score of 0.759, which means that it can predict a wide range of GDD variation ( Table 3). The F-statistic of the p-value, which gives the overall significance of the model, was very high. Moreover, p-values indicated a significant correlation between the GDD and each predictor variable. The t value for all the predictors showed that there is a significant association between them and the GDD, since t is always significantly different than zero. According to the coefficients, the time and the distance from the shoreline have a positive influence on GDD when the altitude and latitude have negative. Considering the standardized (scaled) estimates coefficients (in the parentheses), it is evident that the latitude, the distance from the equator, has a high negative impact on GDD increase, followed by the altitude. The predictor with a high positive impact on GDD increment is the time (year) when the positive influence of distance from the shoreline is negligible. Moreover, the time variable represents climate change because it is the only predictor that varies with time. So, the mean added GDD values per year were more than four units. The standard error (std. error) for each predictor variable was low, pinpointing the model's accuracy. Moreover, the 10-fold cross-validation test yielded an R 2 equal to 0.76 and an MAE of 96.89. To illustrate the spatiotemporal distribution of the GDD in the Balkans area according to olive tree requirements for full blooming, the above MLR model for three time periods (2020, 2040, 2060) was implemented. In Figure 3, the light green color for the year 2020 is covering the total area with GDD above the threshold level of 700. The green and dark green colors for the years 2040 and 2060, respectively, are depicting the additional areas with GDD above the level of 700. Moreover, in the map are presented the regions where was established the olive tree cultivation according to Corine land cover in 2018. As it is presented, the GDD in 2020 allows the cultivation of olive tree in the southern part of the Balkans, covering an area from Greece to the northern borders of Bulgaria in the valley of Danube. The only regions which are excluded by this criterion are the higher mountainous areas. At this point, it must be pinpointed that there are several climatic restrictions, such as early or late frosts or extremely high temperatures, that make the cultivation of the olive tree inevitable, which it is not taken into account in this framework. The expansion of the areas which address the criterion of 700 GDD for the year 2040 is negligible in the southern part of the studied area. There is only a small addition in the high-altitude zone. The region where is recorded a broader expansion is in the northern part of the Balkans, especially in Northern Serbia and in the south of Romania. The time projection to the year 2060 indicates a further expansion of areas and fulfils the criterion of 700 GDD to the northern part of the Balkans. So, the zone of olive tree cultivation is covering the northern part of Croatian (and the neighboring part of Slovenia) along with a vast sector in Central Romania.
The findings of this study are corroborated to other research regarding olive cultivation in the same area. Fraga et al. [72] indicate the strong latitudinal gradient on the spatiotemporal projection of related parameters, such as growing season start or growing season length. More precisely, they found that growing season can be increased in lowlatitude areas. Additionally, Orlandi et al. [73] analyzed the apparent influence of latitude on flowering of olive trees in two Mediterranean regions. Moreover, Tanasijevic et al. [74] found similar spatial expansion of the olive tree cultivation zone to the north, but this research does not include Bulgaria and Romania. Finally, Oteros et al. [75] found that the altitude plays a vital role in the phenological related climatic parameters of the olive trees.  The findings of this study are corroborated to other research regarding olive cultivation in the same area. Fraga et al. [72] indicate the strong latitudinal gradient on the spatiotemporal projection of related parameters, such as growing season start or growing season length. More precisely, they found that growing season can be increased in lowlatitude areas. Additionally, Orlandi et al. [73] analyzed the apparent influence of latitude on flowering of olive trees in two Mediterranean regions. Moreover, Tanasijevic et al. [74] found similar spatial expansion of the olive tree cultivation zone to the north, but this research does not include Bulgaria and Romania. Finally, Oteros et al. [75] found that the altitude plays a vital role in the phenological related climatic parameters of the olive trees.

Viticulture Multiple Regression Model Results and Spatiotemporal Distribution of the GDD
The results of the multiple linear regression model analysis (Table 4) describe a statistically significant relationship between the output variable (GDD) and the predictors, giving an R 2 equal to 0.8. Moreover, each predictor has significant statistic correlation with the independent variable. The coefficients of the predictors show a high positive relation of the GDD with the time and a high negative relation with the latitude. Additionally,

Viticulture Multiple Regression Model Results and Spatiotemporal Distribution of the GDD
The results of the multiple linear regression model analysis (Table 4) describe a statistically significant relationship between the output variable (GDD) and the predictors, giving an R 2 equal to 0.8. Moreover, each predictor has significant statistic correlation with the independent variable. The coefficients of the predictors show a high positive relation of the GDD with the time and a high negative relation with the latitude. Additionally, there is a negative trend of the altitude for the GDD. Finally, the distance from the shoreline has a weak negative influence on GDD variation.
Moreover, the t value indicates clearly that all the selected predictor variables are important for the estimation of GDD over the Balkans area. Finally, the relative standard error (RSE) and the standard error pinpoint the high accuracy of the model. Also, the 10-fold cross-validation indicates an R 2 equal to 0.80 with a mean absolute error (MAE) of 176.14.
The following map (Figure 4) illustrates the areas which fulfil the criterion of 2000 GDD for the May to October period. As presented for the year 2020, all the southern part of the Balkans is suitable for viticulture, considering GDD. More specifically, the coastal zone of Albania, Montenegro, Croatia, and Slovenia, along with an extended part of Bulgaria and Southern Romania and the plains of North Macedonia and Serbia are suitable for viticulture.  The time projection of GDD for 2040 indicates an extensive expansion of the climatic suitability over the northern part of the study area. There is also an advancement around the mountainous regions of the south. However, the most extended zonal spreading is in North Serbia and Croatia and the northern part of Donau's basin near the Bulgarian-Ro- The time projection of GDD for 2040 indicates an extensive expansion of the climatic suitability over the northern part of the study area. There is also an advancement around the mountainous regions of the south. However, the most extended zonal spreading is in North Serbia and Croatia and the northern part of Donau's basin near the Bulgarian-Romanian borders. In general, the results corroborate recent research, such as Santos et al. [76] and Malheiro et al. [77], indicating the transition of the viticulture zone of the European and, more specifically, the Mediterranean area.
Moreover, as it is presented in Table 5, there is a high potential for olive cultivations and viticulture in terms of GDD. However, there is a broad expansion for both cultivations in 2040. The olive cultivation, according to the findings, may sprawl the next decades to northern areas and higher altitudes. Herein, there is a change of new olive production areas in case the rest of the climatic and aquatic requirements can be fulfilled (e.g., free of frost growing season, adequate water resources). Probably the most interesting finding is the future expansion of the olive zone to the plain areas of Northern Bulgaria and Southern Romania because, in these areas, the agricultural potential is extremely high. In case of successful cultivation of the olive tree in the Northern Balkans, the worldwide map of the olive and olive oil would change. The major expansion of this cultivation (reaching almost 50% in 2040) could be made over the northern areas with low altitude, such as the above-mentioned regions and the plain in Northern Serbia. In the case of viticulture, it may expand covering almost 60% of the Balkans area because of the physiology of the plant, which is inactive during the cold winter period of the year. So, the thermal threshold of 2000 GDD can be reached in a wide area of the Balkans in the next decades, offering the opportunity of an expansion of the vines far from the Mediterranean seashore.
The presented results (maps and numerical data) can give information about the suitability only from an atmospheric-thermal point of view, which is narrowed to the GDD. So, to form a more realistic image about the ability of vine and olive cultivation, a combination of those results with soil maps and hydrological data, growing seasons estimations could be made. Moreover, there are several restrictions which come from the land's geomorphological condition (e.g., aspect, slope), and the land use policy of the area. So, the total area which is suitable for each cultivation is a small fraction of the above estimated.

Conclusions
According to the results' analysis and their discussion, it could be inferred that for the next six decades, there is a high potential for olive tree cultivation and viticulture in the northern parts of the Balkans (Northern Bulgaria and Southern Romania), in terms of growing degree days. Moreover, it is estimated that a sprawl to the higher altitudes will give opportunities for new terroirs if the other agroecological factors allow their cultivation. The multiple regression models are simple but effective on the spatiotemporal estimation of GDD units, specialized on vine and olive tree growth requirements. Except for the ability to estimate the thermal bioclimate for the cultivations, the models shed light on the significance of the input parameters on the GDD values. So, the most influential (time-invariant) parameter is the latitude, followed by the altitude, both with a negative trend on the GDD variation. Nonetheless, the time-predictor coefficients prove that the climate change trend is positive and dominant on the GDD variation. This means that the geographic factors have a less-significant influence on the regional cultivations' status than the climatic trend. Finally, both multiple linear regression equations are accurate and robust, and can be an effective tool for the spatiotemporal estimation of the GDD for the Balkans area.
Author Contributions: Authors contribution described as: (i) I.C. consolidated paper conceptualization, data analysis, algorithm, scripts, and manuscript preparation; (ii) I.P. provided both general and specific data analysis and manuscript preparation; (iii) E.P. provided spatial data analysis and maps preparation, and (iv) P.N. provided conceptualization and manuscript review and editing. All authors have read and agreed to the published version of the manuscript.