Maximum, Minimum, and Daily Air Temperature Range in Orchards: What Do Observations Reveal?

: This study was designed to better understand vegetation’s impact on air maximum (T max ), minimum (T min ), and daily temperature range (DTR), as well as seasonality and variability. We selected a ﬂat, under synoptic-scale, northern Serbian region with an operational network of automated weather stations (AWS) for the study. Data were collected directly from the eighteen AWSs placed in the orchard canopy during 2013–2018. Meteorological data, plant phenological data in the form of the BBCH scale, and orchards’ soil characteristics data were collected. Environmental factors inﬂuencing the temperature were classiﬁed as static (slow or unchangeable) and dynamic (fast-changing). The impact of both factors on maximum, minimum, and daily temperature range and its variability were analyzed. Results show that static factors (like soil texture) affect the annual variation of T max , T min, and DTR rather than its variability over the season. The dynamic factors, mainly coming from the plant’s phenology, substantially affected the seasonal variability of these variables. Studies like this suffer from missing data and sparse spatial coverage by the AWS network. Therefore, the alternatives of orchard micrometeorological data, nearest climatological station, and ERA5-Land reanalysis data are tested. Both data sets showcased limitations in their applicability, while reanalysis data deviated more from the in-situ measurements, both seasonally and regionally.


Introduction
Orchards, from a morphological and aerodynamical point, are usually described as a combination of bare soil (or vegetated if the cover crop is present) and managed short tree forests. Unique orchard structures and permanent changes, particularly during the growing season, create specific atmospheric conditions within the surface layer occupied by trees. When the full leaf stage is reached in a forest and the crown is closed, heat, water vapor, and momentum exchange through the forest top are minimal. This allows for the formation of a specific microclimate within the forest canopy. Forests commonly have a lower maximum (T max ), higher minimum (T min ), and lower daily temperature range (DTR) than nearby crop fields [1,2]. Therefore, if the orchards are a combination of forest and crop, it is expected that DTR is lower than in the crop fields. However, it has been observed that DTR is higher inside orchards than in the crop fields [3]. This difference results from the lower vegetation density and defined tree rows which ensure better energy exchange and the rise of DTR.
The T max and T min can be classified as the stressors (frost, cold spells, heatwave) if the critical values come along with temperature-vulnerable periods of the plant growing season, such as flowering or pollination [4]. During the night, canopy air space experiences radiative cooling by soil and plants. Since the cooling is usually a whole night process,

•
How does air T max , T min, and DTR in orchards relate to micro-environmental conditions and their changes in the field? • How far does T max , T min, and DTR in the orchard deviate from measurements performed at the nearest climate station? • How far does T max , T min, and DTR in the orchard deviate from corresponding reanalysis data?
In order to address the first question, we defined micro-environmental conditions of interest using static and dynamic factors. It is expected that static factors will affect the annual variation of T max , T min, and DTR in orchards, while dynamic ones will contribute to their seasonal variability. The second and third questions are addressed using T max and T min measured at the climatological station and assimilated from the ERA5 reanalysis system for locations of selected orchards, respectively. Orchards were carefully selected from the Forecasting and Reporting Service for Plant Protection of the Republic of Serbia (PIS) phenological and micrometeorological network to secure diversity of defined static factors. More details about this network and measurement methodology can be found in [22]. Strong orographic effects were removed from the analysis by choosing locations in the Vojvodina region (Serbia), the flat northern part of the country.

Study Region
Vojvodina province, the northern part of R. Serbia, is a part of the Pannonian Basin and accounts for 27.9% of the country (21,614 km 2 ). The soil characteristic of the Vojvodina region is suitable for intensive crop production, while small areas are under the orchards and vineyards. The fruit production in this region has had a positive trend but with considerable variability in the last two decades, showing the extent of vulnerability to the weather-related risks [23]. The orchards are primarily located in northern Vojvodina, in and around Fruska Gora mountain, and in the south-eastern part of Vojvodina. Vojvodina is roughly 300 km in diameter, and its different landscape units have diameters below 100 km. Therefore, the variations of atmospheric conditions amongst them are in the mesoscale while sharing a common synoptic situation most of the time.

Micrometeorological Data
The maximum and minimum air temperatures used in this study were measured on the automated weather stations (AWSs), a part of the extensive monitoring network governed by the Forecasting and Reporting Service for Plant Protection of the Republic of Serbia (PIS). PIS monitors agricultural production, the development and occurrence of harmful organisms, and the micrometeorological conditions in which the production occurs [22]. The network of AWSs is created and maintained according to PIS requirements related to the monitoring and forecasting of diseases and pests. Currently, the PIS network contains 166 AWSs located all over the country monitoring air temperature, soil temperature, relative air humidity, precipitation, and leaf wetness.
AWSs presented in this paper were selected based on the data quality control (inspected by the station and data manager), length of the continuous measurement (at least two years of continuous measurements in interval 2013−2018), and location (orchard station being well maintained). We selected 18 AWSs out of 31. They are located in apple (9), plum (4), peach (3), and pear (2) orchards in Vojvodina (Table 1, Figure 1). We used the minimum joint period of 5 years (2014−2018), except for year-to-year intraregional analysis. Orchards with hail nets were excluded from this study since hail nets can impact micrometeorological conditions. The temperature extremes were taken as the maximum (or minimum) temperature in the course of a continuous time interval of 24 h, in which the maximum (or minimum) temperature was defined as the highest (or lowest) temperature reported for a given location during a given period [24,25]. The DTR was calculated as the difference between T max and T min . The hourly data were downloaded directly from AWSs, and they are the results of 12 s measurement averaging. Time series were not additionally corrected or filled. Due to the nature of data collection, missing data and strange values in data are part of the process. These points can be wrongly interpreted in the analysis since it is difficult to eliminate the false data if the disturbance which caused it is not explained in detail. Therefore, the selection was managed through visual inspection and critical control during the analysis of the availability. AWSs were located in the canopy row, close to the fruit tree, to measure the temperature and humidity of the canopy air space, which is the living environment of diseases and pests ( Figure 2).  Table 1).

Figure 2.
The AWS located inside the canopy air space in the fruit tree row.  Table 1).

Figure 2.
The AWS located inside the canopy air space in the fruit tree row.   1 . When the sign "-" is present it means that there is no reliable or available information. 2 . Current orchard on this location is young and still growing. Trees are 4 years old and 2 m in height, placed in 4 m rows 2 m apart. However, measurements from this location are dated from 2011 when the station was in the old apple orchard, which was not more than 500 m away. Smoothing of the row data was done with a simple moving average (SMA) [26]. The best model was obtained for the 17-day window (SMA17). However, since smoothing is done to smooth out the irregular roughness to see a clearer signal, we increased the window to 30 days. While smoothing, the temperature difference outliers larger than +/− 10 • C were eliminated (approximation of +/− 5 standard deviations [27]). Because of the extent of the analysis, data produced during this research are also given in the Supplementary Material.
Locations selected for the study were clustered, according to their position based on minimum-distance criteria (50 km threshold), into five sub-regions denoted according to local district names (Srem, Banat, Backa): R1-Southern Backa, R2-Northern Backa and Banat, R3-Southern Banat, R4-Central Banat, and R5-Srem ( Figure 1, Table 1). One or more specific static parameters defined each sub-region: R1, R2, and R4 havedifferent soil textures and distances from the orographic barriers in the direction of prevailing winds and front passages, R3 and R5 have the presence of orographic elements, river bodies, and different landscapes. Here is an overview of sub-region specifics: • R1-an obstacle for S-W warm front is Fruska Gora mountain (8 km-40 km distance, rise (mountain height-450 m)/run (distance) ranges from 0.06 to 0.01); • R2-wide open for circulations predominantly coming from N-W and N-E; there is not any nearby obstacles (rise/run < 0.0045); • R3-mixed orography (combination of partly hilly (Vrsac mountains) and flat terrain), mixed landscape (agricultural land and semi-deserts), and specific weather patterns (the strong influence of Kosava wind); • R4-obstacle for NW winds is the Carpathian Mountain (100 km-150 km distance, rise/run ranges from 0.010 to 0.014); and • R5-mixed orography (more considerable height and horizontal scale of hill (Fruska Gora), flat terrain, and large river (Danube)).

Phenological Data
Phenological observations made by PIS observers were recorded according to BBCH scale ("Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie") on a weekly basis [22,28]. The BBCH scale is the system for uniformly coding phenologically similar growing stages of all mono-and dicotyledonous plant species. Figure 3 shows the phenology stages in the apple, plum, pear, and peach orchards used in this study. The most critical morphological changes occur in spring and correspond to the following phenological stages: BBCH 10 (green leaf tips 10 mm above the bud scales, first leaves separating), BBCH 60 (first flowers open), and BBCH 73 (second fruit fall, maximum leaf area). The leaf fall in autumn was omitted since most data were not collected (autumn phenology was observed only on 135m_NSN_Pr station orchard).

Soil Data
Soil data were taken from a detailed digitalized soil map [29] (Figure 1). Since these data were more detailed than the World Reference Base (WRB) classification and therefore difficult to compare, we also used Soilgrids, a system for digital soil mapping based on a global compilation of soil profile data (WoSIS) and environmental layers [30]. Soil texture data are common for both classifications (the relative proportion of sand, silt, and clay in soil). Since soil texture strongly affects soil drainage, water holding capacity, heat transfer through the soil, and soil temperature, it was used to depicture soil characteristics in this study.

Soil Data
Soil data were taken from a detailed digitalized soil map [29] ( Figure 1). Since these data were more detailed than the World Reference Base (WRB) classification and therefore difficult to compare, we also used Soilgrids, a system for digital soil mapping based on a global compilation of soil profile data (WoSIS) and environmental layers [30]. Soil texture data are common for both classifications (the relative proportion of sand, silt, and clay in soil). Since soil texture strongly affects soil drainage, water holding capacity, heat transfer through the soil, and soil temperature, it was used to depicture soil characteristics in this study.

Climate and Reanalysis Data
For every region, the reference climatological station was selected from the National Observation Network of the Republic Hydrometeorological Service of Serbia (RHMSS) based on the same criteria of the minimum distance (50 km threshold) for the period of 2013-2018. Stations are listed in Table 1. Measurements on the climatic stations were made over the grass on 2m inside the instrument shelter. The climatological period used in the analysis was 1981-2010 for Tmax and Tmin ( Figure 4). A temperature gradient in the climatological data goes from NW (cold) to SE (warm), which was considered during the analysis.

Climate and Reanalysis Data
For every region, the reference climatological station was selected from the National Observation Network of the Republic Hydrometeorological Service of Serbia (RHMSS) based on the same criteria of the minimum distance (50 km threshold) for the period of 2013-2018. Stations are listed in Table 1. Measurements on the climatic stations were made over the grass on 2m inside the instrument shelter. The climatological period used in the analysis was 1981-2010 for T max and T min (Figure 4). A temperature gradient in the climatological data goes from NW (cold) to SE (warm), which was considered during the analysis.
The seasonal T max and T min anomaly for the 2013-2018 period were used to characterize the seasons as "normal" or "extreme" in respect to the climatological period of 1981-2010. In almost all cases, the seasonal T max and T min anomaly was higher than normal ( Table 2). It is important to emphasize that twelve of the fifteen warmest years in Serbia were registered after the year 2000 (period 1951-2018) (by RHMSS), confirmed by Table 3 [31]. During this period, 2013 was very warm with seven heatwaves; 2014 was the rainiest and the second warmest in the period from 1951-2014; 2015 was the third warmest , with six registered heat waves; 2016 was warm and rainy; 2017 was warm and dry; 2018 was the year of climate records in Serbia, the warmest in the history of meteorological measurements. This analysis shed light on the behavior of the orchard air temperature in the "extreme" years.  The seasonal Tmax and Tmin anomaly for the 2013-2018 period were used to characterize the seasons as "normal" or "extreme" in respect to the climatological period of 1981-2010. In almost all cases, the seasonal Tmax and Tmin anomaly was higher than normal ( Table  2). It is important to emphasize that twelve of the fifteen warmest years in Serbia were registered after the year 2000 (period 1951-2018) (by RHMSS), confirmed by Table 3 [31]. During this period, 2013 was very warm with seven heatwaves; 2014 was the rainiest and the second warmest in the period from 1951-2014; 2015 was the third warmest , with six registered heat waves; 2016 was warm and rainy; 2017 was warm and dry; 2018 was the year of climate records in Serbia, the warmest in the history of meteorological measurements. This analysis shed light on the behavior of the orchard air temperature in the "extreme" years. The reanalysis ERA5-Land dataset is a product of the European Centre for Medium-Range Weather Forecasts (ECMWF). These are gridded data (regular latitude-longitude grid of 9 km, ~0.08° × 0.08°) produced hourly for the period of 1981 to the present. ERA5-Land is designed to characterize trends and anomalies of mass and energy cycles [19]. We must emphasize that: "Reanalysis combines model data with observations from across the world into a globally complete and consistent dataset using the laws of physics" [18,32]. Therefore, those are processed data, completed with a numerical model, not only gridded observations. The ERA5-Land dataset was selected because of its high spatial and tem-  The reanalysis ERA5-Land dataset is a product of the European Centre for Medium-Range Weather Forecasts (ECMWF). These are gridded data (regular latitude-longitude grid of 9 km,~0.08 • × 0.08 • ) produced hourly for the period of 1981 to the present. ERA5-Land is designed to characterize trends and anomalies of mass and energy cycles [19]. We must emphasize that: "Reanalysis combines model data with observations from across the world into a globally complete and consistent dataset using the laws of physics" [18,32]. Therefore, those are processed data, completed with a numerical model, not only gridded observations. The ERA5-Land dataset was selected because of its high spatial and temporal resolution, which is a novelty since previous products have had a much lower spatial resolution of 31 km-ERA5 or 80 km ERA-Interim. More details about the ERA5-Land dataset can be found in [19].

Study Design
In order to reach the addressed objectives, the study design was based on regionally specific static and dynamic factors.
(a) Static factors selected to describe the micro-environment of each location were orography, soil characteristic (soil texture-in the aspect of its thermal properties), orchard row orientation, and elevation.
(b) Dynamic factors affecting selected orchards, energy, and water balance at daily and seasonal scales were orchard tree and cover crop phenology dynamics and tree age.
For this purpose, stations were selected based on location characteristics and data quality. The analysis of the data variability was based on the differences between the selected AWSs and the nearest climatological station. Since air temperature is not normally distributed, we used the non-parametric Kruskal-Wallis (K-W) test [33] to compare the medians of three or more groups with different sizes after the value ranks were placed in ascending order. The null hypothesis H 0 states that samples were drawn from the same population (medians are equal), and the alternative hypothesis H 1 that samples were drawn from different populations (different medians). If the p-value is less than the confidence interval (set to 0.05) the H 0 is rejected, and if a p-value is larger than the confidence interval, it does not contain enough information to conclude that the group medians differ. If the p-values are smaller than 0.05, we can reject the idea that the difference is due to random sampling and conclude that populations have different distributions. This variability is the result of the AWS or climatological station location characteristics [34].
Intraregional analysis of T max , T min, and DTR seasonal variability was carried out on for each station (standard deviation) and concerning reference climatological station (difference). The K-W test and Dunn's test were performed for all stations to examine whether differences between stations result from the local characteristics or lack of data [34]. The intraregional analysis should offer information about key factors affecting extreme temperatures in the region and answer the frequently asked question, "Can observations from one location be used on the locations where no measurements are made?" The second and the third research questions were addressed using T max , T min, and DTR measured in the orchard, on the closest climatological stations and from the closest reanalysis point. The Taylor diagram (TD), commonly used to present differences between modeling and measurements, was adopted for the study by replacing model results with observations from the climatological station and reanalysis data. Namely, the TD provides information about three statistical quantities: (1) the Pearson correlation coefficient (CC), (2) the centered root-mean-square error (RMSE), and (3) the standard deviation (SD) [35]. The CC was used to confirm the similarity in the seasonal pattern of the T max , T min, and DTR and RMSE and SD to account for the deviations between orchard measurements in one hand and climatic station measurements and reanalysis data.

Factors Driving Variability of Temperature Extremes
Static and dynamic factors selected and analyzed are given in Table 3. AWSs used to analyze specific factors were selected from the same sub-region (Table 1) to minimize the influence of the temperature gradient ( Figure 4).
The effect of the soil texture was examined on two AWSs 74m_SEK_G_A-clay soil and 91m_SOR_A-sandy loam (Figures 5, S1.1 and S1.2). Clay soil generally has a higher volumetric heat capacity than sandy soils. Therefore, sandy soil has a higher heating rate, and the air temperature above it during the day is higher than the air above the clay. During the night, sandy soil has a higher cooling rate and loses energy faster than clay soil, leading to the lower air T min (Figures 5, S1.1 and S1.2). On top of that, the apple orchard on clay has a grass cover crop, which is also responsible for lowering temperature through evapotranspiration. During the "normal" summer, the T max and T min anomaly was not more than 1 • C ( Table 2). The difference in the T max median in 2014 resulted from AWS 91m_SOR_A position. In 2014, this AWS was in an old apple orchard. The apple orchard density, high LAI, and soil coverage were the reasons for lower summer T max . In 2016, this AWS was moved 500 m to a young new apple orchard. The position change gave T max , T min, and DTR differences that resulted from the soil texture difference, small elevation change, and tree age (Table 3). Trees were too small to impact the temperature difference. The K-W test confirms this result with no significant difference between locations in 2014 but a significant difference after 2016 when the AWS was moved to the young orchard. During spring 2018, the synoptic-scale disturbance when cold air mass dropped the T max below 0 • C by several degrees during March, and the sudden temperature rise during the hottest April recorded since 1951 (Table 2) did not diminish the slight difference between these two locations (Figures S1.1 and S1.2), but it made it statistically insignificant. Seasonally, the statistical significance and impact of locations were more present in spring and summer values of DTR than in autumn and winter, confirming the impact of the fast-changing dynamical factors (Figures S1.1 and S1.2). The increase of LAI and closed tree crown minimize T max and T min variation during the vegetation period Figures 5, S1.1 and S1.2 but a significant difference after 2016 when the AWS was moved to the young orchard. During spring 2018, the synoptic-scale disturbance when cold air mass dropped the Tmax below 0 °C by several degrees during March, and the sudden temperature rise during the hottest April recorded since 1951 (Table 2) did not diminish the slight difference between these two locations (Figures S1.1 and S1.2), but it made it statistically insignificant. Seasonally, the statistical significance and impact of locations were more present in spring and summer values of DTR than in autumn and winter, confirming the impact of the fastchanging dynamical factors (Figures S1.1 and S1.2).  Orchards selected for the study were commonly oriented in the near N-S direction (NNW-SEE) to extend plant exposure to solar radiation during the daytime. If the rows are W-E oriented, the solar radiation will more effectively heat the soil in the rows, and higher temperature extremes can be expected. We can argue that slight constant differences in T max and T min on W-E-oriented orchards are expected, but unfortunately, we could not extract if from this data.
The orography impact on the formation of cold pools is expected in the measured T min during the winter and springtime due to atmospheric stability at sunrise. We speculate that if orographic impact exists, it should be noticed in the region of Fruska Gora mountain (539 m altitude). Therefore, two AWSs were selected on the northern (135m_NSN_Pr) and the southern (94m_SMK_G_Pr) mountain slopes (Table 1, Figure 6). Both stations are placed in the pear orchards with similar soil texture but on different elevations (135m and 94m), with one of them (94m_SMK_G_Pr) having a grass cover crop. Annual variation of extreme temperatures ( Figure S1.3) indicates almost the same temperatures in both locations. Expected slightly higher temperatures on the southern slope we<re reduced by cover crops evapotranspiration (Table 2), shifting radiation balance partitioning in favor of latent heat. This non-significant difference stayed during the 2018 spring warm episode for all variables except T max between 135m_NSN_Pr and the Sr. Mitrovica climatological station on the other side of the Fruska Gora mountain ( Figure 6). Although expected, the location in this example did not contribute significantly to the T max data variability but significantly impacted the T min seasonally ( Figure S1.4).
Leaf appearance in early spring affects radiation balance partitioning in the surface atmospheric layer occupied by fruit trees. All presented data were used in our attempts to isolate plant phenology impact on the orchard canopy's air temperature. The impact of phenology dynamic was the most evident in the seasonal variation of DTR difference, which corresponds to LAI changes. In all examples, under all five sub-regions, the value of DTR difference between AWS and DTR calculated with data from climatological station goes above 0 in the time of intensive development of the crown and drops under during the time of the leaf fall ( Figures S2.1-S2.5). This change corresponds to the increasing impact of evapotranspiration on both T max and T min, which is then reflected on DTR. In 2016, the warm winter and T max in January and February frequently above 15 • C resulted in the early start of pear and peach growing season, inducing flowering at the end of February and mid-March (Figure 3). In the case presented in Figure 7, the impact of phenology can be noticed after DOY 70 when T max on all locations drops below the values from the climatological station. Consequent humidity increase affects longwave radiation fluxes on the canopy level, reducing daily temperature range and daily tendencies of temperature. A similar effect can be seen in 2017 but one month later. However, in 2016 a significant difference in phenology dynamics among locations was reset by a cold spell (1.9 • C < T max < 13 • C) in the middle of March. After this, the plant development harmonized in all locations and variation of extreme temperatures followed plant phenology, reaching its minimum during late spring and summer (DOY 165-255). The phenology impact is most recognizable between BBCH 60 and BBCH 70 (based on Figure 3  The impact of small elevation changes on Tmax, Tmin, and DTR was inspected based on the intraregional analysis (Tables 1 and 2  The impact of small elevation changes on T max , T min, and DTR was inspected based on the intraregional analysis (Tables 1 and 2, Figures 8 and S2.1-S2.5). The most considerable elevation difference between AWSs and regional climatological stations was 65 m in R1, while others were around 40 m. Compared to standard temperature decrease with elevation of near 0.

Small-Scale Variability of Temperature Extremes
Intraregional analysis of extreme temperature variations was carried out to access differences in Tmax, Tmin, and DTR in nearby orchards. The main regional findings are summarized in the following paragraphs, while figures and tables are given in Supplementary  S2 and S3, respectively. R1. The highest values of Tmax and lowest values of Tmin in this region, followed by the largest difference in DTR (both up to 5 °C), were commonly measured at 79m_NSC_A, particularly during the growing season. This is an example of how a small impact of all static factors can produce a noticeable effect, in this case in favor of higher Tmax and lower Tmin from spring to autumn. Namely, 79m_NSC_A was located in the oldest apple orchard in the region (maximum LAI reached), with the shortest distance among trees, lowest altitude, and without cover crops, which led to the temperature differences presented in Figure 8. This orchard is also closest to the climatological station. The lowest standard deviation of temperatures during summer (Table 4, Supplementary S3) occurs from constant leaf area density canopy in all orchards, contributing to a more uniform transfer of heat and water vapor between canopy air space and atmosphere in comparison to spring and autumn. All stations have higher Tmax, lower Tmin, and higher DTR during the summer with respect to the Novi Sad reference climatological station. R3. Tmax and Tmin differences of 2.5 °C (up to 5 °C in specific cases) in respect to reference Vrsac climatological station were observed on the plum orchard (AWS 90m_CRC_Pl) surrounded by flat crop fields and orchards, witnessing the complexity of the region (Figure S2.2). A similar pattern of difference was observed on the apple orchard station 80m_BCT_G_A. The most considerable differences showcased in DTR should be addressed by the impact of Vrsac mountains on atmospheric circulation, the role of Bela

Small-Scale Variability of Temperature Extremes
Intraregional analysis of extreme temperature variations was carried out to access differences in T max , T min, and DTR in nearby orchards. The main regional findings are summarized in the following paragraphs, while figures and tables are given in Supplementary S2 and S3, respectively.
R1. The highest values of T max and lowest values of T min in this region, followed by the largest difference in DTR (both up to 5 • C), were commonly measured at 79m_NSC_A, particularly during the growing season. This is an example of how a small impact of all static factors can produce a noticeable effect, in this case in favor of higher T max and lower T min from spring to autumn. Namely, 79m_NSC_A was located in the oldest apple orchard in the region (maximum LAI reached), with the shortest distance among trees, lowest altitude, and without cover crops, which led to the temperature differences presented in Figure 8. This orchard is also closest to the climatological station. The lowest standard deviation of temperatures during summer (Table 4, Supplementary S3) occurs from constant leaf area density canopy in all orchards, contributing to a more uniform transfer of heat and water vapor between canopy air space and atmosphere in comparison to spring and autumn. All stations have higher T max , lower T min, and higher DTR during the summer with respect to the Novi Sad reference climatological station. R2. The station 91m_SOR_A was placed in the intensively growing apple orchard surrounded by flat crop fields. Annual changes of reference AWS morphology were indicated by season-to-season changing patterns of T max and T min differences as the young orchard progressed. The difference is more prominent for maximum than for minimum temperature, commonly during the summer when a young orchard cannot produce enough evapotranspiration, increasing latent heat and dumping the maximum temperature. Significant differences among AWSs are not expected since there is no considerable difference in this region's static or dynamic factors and orography. However, Figure S2.1 indicates significant differences among 91m_SBV_Pe and 119m_SLJ_A, which are most probably caused by soil texture differences. All stations have higher T max , lower T min, and therefore higher DTR during the summer with respect to the Subotica reference climatological station.
R3. T max and T min differences of 2.5 • C (up to 5 • C in specific cases) in respect to reference Vrsac climatological station were observed on the plum orchard (AWS 90m_CRC_Pl) surrounded by flat crop fields and orchards, witnessing the complexity of the region ( Figure S2.2). A similar pattern of difference was observed on the apple orchard station 80m_BCT_G_A. The most considerable differences showcased in DTR should be addressed by the impact of Vrsac mountains on atmospheric circulation, the role of Bela Crkva lakes and Danube River as a source of humidity, a mix of agricultural land and semi-desert in orchard surroundings, and high speeds of Kosava wind more affecting locations closer to the Danube River bank (90m_CRC_Pl and 80m_BCT_G_A). Within this reasoning, a lower DTR on 127m_PBK_A was expected due to positions furthest from the Kosava entrance point. All stations have higher T max , lower T min, and therefore higher DTR (except 127m_PBK_A as described above) during the summer with respect to the Vrsac reference climatological station.
R4. The sub-region R4-Central Banat data analysis eliminated all but three AWSs (two plum and one apple) with an elevation from 71 to 85 m. For the referent climatological station, we selected Kikinda as the center of this sub-region. The 74m_ZRS_Pl changed its place in 2015/2016 with a 25-year-old apple orchard and a 20-year-old plum orchard next to each other. This change is visible in Figure S2.3 when the trend of the T max difference changed from negative to positive during the summer. The seasonal standard deviation in the R4 region had the same annual pattern. It is the maximum during the fall for T max and minimum in summer. The maximum average T min SD is in fall, while the minimum is in summer (Tables 4 and S3).
The annual trend of T max and T min over this region was relatively uniform, leaving differences in respect to a reference Kikinda of up to 2.5 • C. The only large deviation in trend can be noticed in AWS 74m_ZRS_Pl in 2016 in T min because of the relocation ( Figure S2.3).
R5. Selected AWSs were located in the region of the Fruska Gora mountain. The vicinity of the Danube RRiver, in some cases, is less than 1.5 km (145m_RNS_G_Pe). Typically, higher temperatures at 145m_RNS_G_Pe and lower at 125_SMD_A are caused mainly by their position in respect to the reference station Sr. Mitrovica. Namely, 125_SMD_A is on the slope of the mountain close to the forest. The terrain slope and vicinity of the forest affects orchard micrometeorological conditions by draining cooler air through the orchard rows and reducing the temperature. On the other hand, 145m_RNS_G_Pe is on the eastern slope oriented towards the Danube River and even at a higher altitude. Here the river provides a permanent supply of humidity for the orchard canopy and slightly higher temperatures, particularly during the summer.
Intraregional analysis (Table 4 and Table S3) singled out region R4 as the warmest one in the spring and the summer (with the lowest standard deviation of summer extremes) and R2 as the coldest. This corresponds to the temperature gradient given in Figure 4. Other regions differ significantly, either in magnitude (up to 0.8 • C) or in standard deviation (SD) (up to 0.9 • C) of T max and T min . The difference of 0.9 • C corresponds to the differences between observed data on climatological stations and climate time series from 1981-2010 ( Table 2). The low standard deviation of extreme temperatures during the summer results from radiation partitioning during the growing season, an increased portion of latent heat flux spent on evapotranspiration, and sensible heat flux spent on plant energy balance.
The K-W test and Dunn's test run for all stations under the separate regions revealed a small number of pairs where we cannot conclude the data difference (p ≥ 0.05). These cases are mostly connected to the lack of data, like in paring of 151m_NSC_G_Pe (569) and 135m_NSN_Pr (1813) from the R1 region (Supplementary S4). All others show that local characteristics are responsible for the significance of the difference in the data distributions.

Observed Temperature Extremes in Field vs. Climate Station and Reanalysis Data
The Taylor diagram, commonly used for model performance evaluation, was adopted to compare three meteorological data sources: AWS located in the orchard, the nearest climatological station, and ERA5-Land reanalysis for grid element corresponding to AWS. Since our focus is on orchard micrometeorology, the AWS measurements are used as observed data in TD. Five stations were selected based on intraregional analysis. Selected parings of AWS-climatological stations are: R1: 79m_NSC_A-Novi Sad, R2: 91m_SOR_A-Subotica, R3: 82m_CRC_A-Vrsac, R4: 71m_KIK_Pl-Kikinda, R5: 94m_SMK_G_Pr-Sr. Mitrovica. Details can be found in Table 1. Here we have focused on one paring R1: 79m_NSC_A-Novi Sad. Let us note that the Novi Sad climate station is located approximately 5 km from the AWS 79m_NSC_A, surrounded by the same landscape (crop fields) and similar elevation (84m). It is a common practice by producers in the area to use data from the nearest climatological station if there are no AWS in the field, while recently researchers have used reanalysis data. In the above analysis, the K-W test showed no significant difference between AWS and climatological data for spring 2016 (Figure 7).
Differences among selected data sources are presented in Supplementary S5 and S6 and on scatter plots for 2017 in Figure 9. Seasonal analysis for this station and the other four can be found in Table 5.    The TD averaged through four seasons for T max states (1) CC is between 0.96 and 0.99 for all data sources, but slightly higher for Novi Sad station than for ERA5-Land data; (2) RMSE is less than 2 • C (app. 1.5 • C), while (3) SD for Novi Sad station and ERA5-Land data is smaller than the SD for the 79m_NSC_A AWS (7.7 • C in average trough the seasons). Listed statistical indices become fully evident on the T max regression plot, indicating that T max measured in the orchard is commonly higher than T max measured at Novi Sad station or associated with ERA5-Land. The difference slightly increases with temperature for the entire 2013−2018 period ( Figure S5.5).
The TD averaged through four seasons for T min states (1) CC for Novi Sad station is near 0.99, while CC for ERA5-Land is close to 0.95 (similarity is confirmed), (2) RMSE is higher for ERA5-Land (2 • C) than for Novi Sad station (1.5 • C), and (3) for both Novi Sad and ERA5-Land, SD is more prominent than SD of 79m_NSC_A AWS data particularly during the growing season (spring, summer) ( Figure S5.2). Throughout the year, T min in the orchard is lower than on the Novi Sad station, and it is lower than the data offered in reanalysis. However, there is a difference between the Novi Sad station and reanalysis data which is not evident in T max . It is fair to say that the Novi Sad station data are closest to the orchard measurements; even the deviation still increases with temperature. ERA5-Land shows deviation from field measurements which can be more than 5 • C for T min above 20 • C.
The seasonal differences between T max and T min are most evident in the DTR: (1) CC is between 0.85 and 0.95 (similarity is confirmed) for all data sources, but slightly higher for Novi Sad station than for ERA5-Land data; (2) RMSE for ERA5-Land is 3.8 • C on average (Table 5) which is 1.7 • C more than RMSE for Novi Sad station, while (3) SD for Novi Sad station and ERA5-Land data is smaller than the SD for the 79m_NSC_A AWS (5.3 • C in average trough the seasons). The scatter plot for DTR showcases an even larger difference between datasets than the T min plot. In some cases, differences in DTR are around 10 • C higher in the observer data (AWS) than in the climatological or ERA5-Land data. These large differences are presented with a smaller R 2 of 0.77 for ERA5-Land in 2017.
Since T min is more important during the winter (plant chilling during the dormancy) and spring (frost), T max is more important during the summer (heatwave and draught), and DTR reflects their changes, the analysis for other locations and years was performed accordingly. The SD of both climatological stations and ERA5-Land data was lower than data from AWS for all years (Table 5, Supplementary S6), while RMSE was in range for T max for ERA5-Land: 2.8-1.7 • C and climatological stations: 2.1-0.7 • C, for T min for ERA5-Land: 3.9-1.7 • C and climatological stations: 3-0.9 • C, and for DTR for ERA5-Land: 5.9-2.7 • C and climatological stations: 3.8-1.4 • C. The differences in data sources are most evident in the extreme seasons ( Table 2). The warm winter in 2014 brought large differences in T min for all locations, where all ERA5-Land data better explained the temperatures in the field than the data from climatological stations. However, during the hot spring 2018, the situation was reversed. Winter and spring temperatures together are important for the successful blooming of the following year. As analysis showed, in that part of the year we cannot rely only on climatological and reanalysis data.
The intraregional comparison of SD and RMSE of ERA5-Land and climatological stations did not show large differences. The larger average SD and RMSE were observed in the R3 region, which is the most diverse in the landscape described in Section 3.2. The RMSE was, in all cases, larger for T min than for other values, on average up to 2.9 • C (Supplementary S6). This result speaks about the necessity of the bias correction of the reanalysis.

Conclusions
This study emerged from the analysis of the AWS measurements conducted over six years (2013)(2014)(2015)(2016)(2017)(2018) in the northern agricultural part of R. Serbia. Through data analysis and collecting knowledge about T max , T min, and DTR in orchards, we came to several findings listed below.
Air T max , T min, and DTR in orchards mirror the micro-environmental conditions and their changes in the field. Static environmental factors (elevation, soil texture, row orientation, and orography) impose low to medium impact on T max , T min, and DTR in the flat area orchards of this region. The largest impact is associated with soil texture (up to 2.5 • C). However, dynamic factors (cover crop phenology, tree phenology, and age) have a high impact on both seasonal variations of T max , T min, and DTR and their magnitude. This impact is pronounced in T max and T min but is more easily spotted in DTR since it includes changes of both.
During the winter, T max , T min, and DTR variation among stations in a particular region are small and driven mainly by static factors. As the fruit trees' phenology progresses, the trend of the T max , T min, and DTR differences increases, reaching its maximum in the time of maximum LAI in the late spring and summer. During the summer, differences are almost constant. The closed tree crown and heating from the soil in the rows produce almost the same temperatures in orchards of similar age. At this time, the difference is a result of orchard structure and static characteristics. It confirms the smallest seasonal SD in the summer on all locations. Closed orchards tend to have lower maximum temperatures in summer. Fall is caring for new phenological changes and an increase in variability as the trees go into dormancy. In selected regions, the annual variability of the T min differences is smaller than in the case of T max . This leads to the conclusion that the most crucial factor affecting the magnitude and variability of T max , T min , and DTR in the orchard located in a flat area is plant phenology dynamics.
Abrupt changes in the daily course of temperature can be the result of precipitation. The impact of rain is assessed by comparing extreme temperature differences in every sub-region associated with days with and without rain. However, the analysis of the available data indicated that rain events do not change the trend of extreme temperatures over the season.
While considering alternative sources of data describing orchard micrometeorological conditions, we conclude that both climatological station data and ERA5-Land reanalysis use in vegetation-related studies can produce errors. Comparative studies and validation of the reanalysis in other cases showed temperature bias in the range of 2.5 • C and R coefficient larger than 0.95 [37]. Our findings of higher T min and lower T max in ERA5-Land, which was also the result of the study done with ERA5 data, can narrow the specter of possible reasons [38,39]. Other studies confirmed the applicability of the ERA5 and ERA5-Land data for the gap-filling of the observed meteorology for regional-scale agricultural modeling [21,40].
Due to its short-term adverse effects, let us shed more light on T min data. Plant vulnerability on low temperatures, as one of the key elements of plant production risk, requests high alert in assessing T min data sources in case air temperature is not measured in the orchard. The results presented indicate that reanalysis data could not be a reliable source of information about low-temperature exposure of plants in the case examined in this study. However, ERA5-Land can be improved with data from climatological stations trough bias correction. This process needs to be performed before its use for local case studies. When the climatological station is close enough, with similar static factors, T min data from the climatological station can be used but with awareness of the possible error.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/atmos12101279/s1. Supplementary S1-Seasonal variation of T max , T min and DTR in specific orchards under the specific sub-region ( Figure S1.1-S1.6). Supplementary S2-Annual and seasonal variation of T max , T min and DTR difference in respect to reference climatological station for different sub-regions ( Figure S2.1-S2.4). Supplementary S3-Seasonal average and standard deviation for the AWSs under the specific sub-region (Table S3.1-S3.5). Supplementary S4-Values of p from K-W test and Dunn's test for pairs where p ≥ 0.05. Supplementary S5-Seasonal Taylor diagram of the T max , T min and DTR for different AWSs and in respect to climatological station and ERA5-Land reanalysis for the closest point ( Figure S5.1-S5.15). Supplementary S6-Average of standard deviation (SD), centered root-mean-square error (RMSE), and Pearson correlation coefficient (CC) for T max , T min and DTR for 5 selected AWSs (observed), corresponding ERA-Land data and closest climatological station (model) ( Table S6).