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

Abstract

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 flat, 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 influencing the temperature were classified 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.

1. 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, the air temperature reaches its minimum just before sunrise. In winter, minimum air temperature can decrease due to plant additional radiative cooling, which, in soil, could be colder than air [5]. During the growing season, the presence of a cover crop means that more sunlight is reflected and less energy is stored in the soil. Therefore, the cover crop can increase the frost risk by lowering T_min [6]. This situation was observed in experiments done in managed forests (trees in rows with significant space between), where the average T_min was substantially lower than in unmanaged ones, with a high and dense understory [7,8].

During the day, the bare orchard rows represent an important sink (soil evaporation is a non-negligible component of evapotranspiration) and energy source which contributes to the raising or lowering of the T_max [9]. During the growing season, leaf area density affects heat, energy, and momentum transfer between the underlying soil surface and the above orchard air space, like in sparse forests, lowering the T_max.

The feedback between the energy and water balance of a plant canopy and its development is, in general, the subject of numerous studies, but not many studies deal with orchards [10,11]. These studies are commonly focused on: (a) the impact of the air temperature on plant and yield development [12,13]—not the other way around [10,14]; or (b) specifics of the orchard micro-climatic conditions often on one location and in one season in respect to climatological measurements. Some show very little or non-existant differences in air temperature measured in the orchard and on the climatological station [10].

Key factors affecting a range of the daily temperature variation (DTV) can be classified as static (factors that do not change or change slowly, commonly in one direction) and dynamic (factors that could change during the period of interest) [15]. On synoptic and mesoscales, the static factors are latitude (higher DTV if closer to the Equator), altitude (higher DTV at lower elevations), orography (higher DTV in lowlands), soil type (highest DTV in the case of sandy bare soil or rock), and prevailing weather patterns, while the dynamic ones are cloudiness (higher DTV under a clear sky) and vegetation cover (higher DTV over bare soil) [3]. The T_max, T_min, and DTR variability differs from one season to another and depends on the atmospheric circulation and local physical processes [16,17].

The motivation to deeply explore T_max, T_min, and DTR in orchards came from our awareness about the impact of vegetation on air temperature and humidity and their extreme values. Since in situ micrometeorological measurements are not always present and agrometeorological networks are usually sparse, users frequently ask questions about differences between the data sources. The most popular data sources are climatological measurements and reanalysis data. The reanalysis data are a relatively new product that combines past short-range weather forecasts with observations through data assimilation [18,19]. This data source is getting more popular because of a lack of other data sources. Earlier data are usually not homogenous, are not filled, or in many cases are without proper metadata. Therefore, using reanalysis as the source of weather and climate data in environmental studies is becoming more common [20,21].

The objectives of this evolved study can be summarized through the following research questions.

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

2. Experiments

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

2.2. Data

2.2.1. 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. List of AWSs and climatological stations used in the study. AWS is elevation + short location + cover crop + orchard type to provide information for comparison. The closest city was named as the climatological station.

	Short Name	AWS Location	Latitude/ Longitude	Elev. (m)	Orchard	Interval	Age/ Height	Cover Crop	Density	Soil	Row Orientation	Notes
R1—South Backa	79m_NSC_A	Novi Sad, Cenej	45.36/ 19.81/	79	Apple	9 February 2013 19 November 2018	17 y. 3 m	No	Semi-intensive1.8 m × 4 m	Chernozem calcareous (micelar) on loess terrace (clay loam)	NW-SE
	80m_NSK_A	Novi Sad, Kac	45.29/ 19.90	80	Apple	8 July 2017 31 December 2018	- 1	-	-	Aluvial il (clay loam)	NW-SE
	78m_VRV_Pl	Vrbas, Vrbas	45.56/ 19.64	78	Plum	11 April 2013 8 April 2016	-	-	-	Chernozemlike meadow soil on loess terrace (clay loam)	NW-SE
	135m_NSN_Pr	Novi Sad, Nestin	45.23/ 19.47	135	Pear	2 October 2012 31 December 2018	15 y. 3 m	No	Semi-intensive 2 m × 4 m	Chernozem calsareous brownized (clay loam)	N-S	Danube on 300m Fruska Gora
	151m_NSC_G_Pe	Novi Sad, Cerevic	45.22/ 19.66/	151	Peach	1 January 2016 4 October 2017	12 y. 2.5 m	Yes	-	Chernozem calsareous brownized (clay loam)	NNW-SEE	Danube on 600m Fruska Gora
Climatological station		Novi Sad	19.85/ 45.33	84						Chernozem calsareous brownized (clay loam)
R2—North Backa and Banat	91m_SOR_A	Sombor, Ridjica	45.99/ 19.09/	91	Apple	27 June 2011 31 December 2018	4 y. 2 m 2	No	Semi intensive 2 m × 4 m	Chernozem on sandy loess (sandy loam)	NNW-SEE	Kiđoš river 1.5 km
	91m_SBV_Pe	Subotica, Bac. Vinogradi	46.11/ 19.88	90.5	Peach	1 January 2016 31 December 2018	-	-	-	Soloncak (sandy clay loam)	WSW-ENE
	74m_SEK_G_A	Senta, Kanjiza	46.05/ 20.07	74	Apple	24 February 2011 31 December 2018	20 y. 3 m	Yes	Semi intensive 2 m × 3 m	Black limeless soil (clay)	NNW-SEE	Tisa river 3 km
	119m_SLJ_A	Subotica, Ljutovo	46.08/ 19.52	119	Apple	8 March 2014 31 December 2018	-	-	-	Chernozemlike medow soil on loess plateau (sandy loam)	NNW-SEE
Climatological station		Subotica (Palic)	19.15 45.77	102						Black limeless soil on sand (sandy loam)
R3—South Banat	82m_CRC_A	Crvena Crkva	44.90/ 21.36	82	Plum	1 January 2013 23 December 2016	-	-	-	Chernozem brown (clay loam)	NNE-SSW
	127m_PBK_A	Pancevo, Ban. Karlovac	45.03/ 21.03	127	Apple	16 January 2013 14 December 2016	-	-	-	Chernozem calcareous (micelar) on loess plateau (sandy loam)	NNE-SSW
	80m_BCT_G_A	Bela Crkva, Tirovo	44.89/ 21.42	80	Apple	24 March 2015 31 December 2018	22 y. 2.5 m	Yes	-	Alluvial sandy soil (loamy sand)	W-E
Climatological station		Vrsac	21.32 45.15	83						Black limeless soil (clay)
R4—Central Banat	71m_KIK_Pl	Kikinda, Kikinda	45.86/ 20.53	71	Plum	22 April 2012 31 December 2018	15 y. 3 m	No	-	Black limeless soil on sand (sandy loam)	NE-SW
	74m_ZRS_A	Zrenjanin, Sutjeska	45.39/ 20.69	74	Plum/ Apple (2016)	8 December 2011 31 December 2018	20 y. 25 y.	No	-	Chernozemlike medow soil on loess terrace (sandy clay loam)	N-S	Concrete terrace
	85m_ZNM_A	Zrenjanin, Novo Milosevo	45.71/ 20.37	85	Apple	8 March 2013 1 October 2017	-	-	-	Chernozem calcareous (micelar) on loess terrace (clay loam)	NNE-SSW
Climatological station		Kikinda	20.46 45.83	81						Chernozemlike medow soil on loess plateau (sandy loam)
R5—Srem	94m_SMK_G_Pr	Sr. Mitrovica, Kukujevci	44.92/ 19.75	94	Pear	23 January 2014 31 December 2018	-	Yes	-	Pararendzina soils on loess (sandy loam)	NE-SW
	125m_SMD_A	Sr. Mitrovica, Divos	44.99/ 19.61	125	Apple	6 November 2014 31 December 2018	-	-	-	Chernozem calsareous brownized (clay loam)	WNW-ESE
	145m_RNS_G_Pe	Ruma, N. Slankamen	45.14/ 20.23	145	Peach	1 January 2016 24 March 2017	-	Yes	-	Pararendzina soils on loess (sandy loam)	NNE-SSW	Danube on 100 m, Fruska Gora
Climatological station		Sr. Mitrovica	19.55 45.10	82						Chernozem calsareous brownized (clay loam)

¹. When the sign “-” is present it means that there is no reliable or available information. ². 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.

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

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

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

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

2.2.4. 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. The average seasonal anomaly of T_max and T_min for four selected climatological stations in respect to mean T_max and T_min for the 1981–2010 period (station Vrsac is not included since data were not available for this entire climatological period).

	Spring		Summer		Autumn		Winter
	Tmax	Tmin	Tmax	Tmin	Tmax	Tmin	Tmax	Tmin
2013	0.1	0.4	1.3	0.5	1.0	1.2	1.2	1.9
2014	1.0	1.3	−0.1	0.3	0.9	2.4	3.7	3.3
2015	0.5	0.2	2.6	1.6	0.8	1.3	1.4	2.4
2016	0.9	0.8	0.6	0.8	0.1	−0.1	1.9	1.2
2017	1.5	0.4	3.3	1.2	0.6	0.7	0.4	−0.8
2018	2.4	2.0	1.3	1.9	3.0	1.4	1.1	1.5

). 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 (1951–2015), 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 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].

2.3. 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₀ states that samples were drawn from the same population (medians are equal), and the alternative hypothesis H₁ 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₀ 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.

3. Results

3.1. Factors Driving Variability of Temperature Extremes

Static and dynamic factors selected and analyzed are given in

Table 3. List of static and dynamic variables with the expected casual pathways.

Static/Dynamic	Variable	Selected AWS	Specific to Orchard	Identified Impact	Causal Pathway	Associated Figures
St	Soil texture	74m_SEK_G_A 91m_SOR_A	No	Yes	Daily soil temperature variation decreases in order sand > loam > peat > clay producing air Tmax and Tmin in the same order [36]	Figure 5, Figures S1.1 and S1.2
St	Orography	94m_SMK_G_Pr 135m_NSN_Pr	No	No	South-facing slopes have higher sun exposure	Figure 6, Figures S1.3 and S1.4
St	Row orientation	-	Yes	No	The W-E-oriented rows have higher soil sun exposure
St	Elevation	All stations	No	Yes	Standard temperature decrees with elevation 0.7 °C/100 m	Figure 8 and Figures S2.1–S2.5
Dy	Cover crop phenology	74m_SEK_G_A 91m_SOR_A 94_SMK_G_Pr 135m_NSN_Pr	Yes	Yes	The presence of a short cover crop lowers the Tmax and Tmin. Vineyard and orchards with the cover crop are typically between 0 °C and 0.5 °C colder than one with bare soil [6]	Figure 5, Figures S1.1 and S1.2 Figure 6, Figures S1.3 and S1.4
Dy	Tree phenology	All stations	Yes	Yes	The increase of LAI minimizes Tmax and Tmin variation during the vegetation period	Figure 7, Figures S1.5 and S1.6
Dy	Tree age	74m_SEK_G_A 91m_SOR_A	Yes	Yes	The increase of LAI and closed tree crown minimize Tmax and Tmin variation during the vegetation period	Figure 5, Figures S1.1 and S1.2

. 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 (Figure 5, Figures 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 (Figure 5, Figures 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).

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 after DoY 90 in all years except 2016 for pear and peach) when plants are fast developing, reaching their maximum LAI (Figures S1.5 and S1.6). In this period, T_max and T_min develop constant differences during the summer and until the foliage dies in late autumn. The best examples of this impact are visible in Figures S2.1–S2.5.

The impact of small elevation changes on T_max, T_min, and DTR was inspected based on the intraregional analysis (Table 1 and Table 2, Figure 8 and Figures 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.7 °C per 100 m, this elevation difference should result in a temperature change of 0.1 to 0.4 °C. The small elevation difference impact was identified on the stations over the 110 m altitudes (R1: 135 m, 151 m, R2: 119 m, R5: 125 m, 145 m). This difference was most visible in the values of DTR, where all stations over 110 m (more than 30 m above the climatological station) had negative DTR or slightly positive differences during the seasons (Figures S2.1–S2.5) due to the lower T_max and higher T_min on higher elevation, which corresponds to the standard temperature elevation decrees.

3.2. 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. Average and standard deviation of T_max, T_min, and DTR for all AWS stations under the specific sub-region.

		Spring		Summer		Fall		Winter
		AV	SD	AV	SD	AV	SD	AV	SD
R1	Tmax	18.5	6.9	28.6	4.8	18.5	6.8	6.2	5.3
	Tmin	7.9	4.8	15.9	3	8.2	4.9	−0.1	4.4
	DTR	10.6	4.4	12.7	4.0	10.3	4.5	6.3	3.4
R2	Tmax	18.3	7.1	29.2	4.5	17.6	7.2	5.3	4.9
	Tmin	6.1	5.2	14.6	3.2	6.4	5.1	−1.6	4.6
	DTR	12.1	4.8	14.6	4.3	11.1	5.3	6.9	3.6
R3	Tmax	18.2	6.9	29.3	4.8	17.8	7.4	6.8	5.3
	Tmin	8	5.9	15.4	3.4	8	5.9	0.1	4.8
	DTR	11.0	5.3	13.9	4.8	9.8	5.2	6.6	3.8
R4	Tmax	19.2	6.6	29.1	4.6	18.3	7.1	5.7	5
	Tmin	7.9	4.7	15.8	3.2	8.2	5.1	−0.4	4.6
	DTR	11.3	5.1	13.8	4.3	10.1	4.9	6.1	3.6
R5	Tmax	18.5	7.3	28.8	4.8	17.4	7.3	5.8	5.8
	Tmin	7.5	4.9	16.6	3.4	7.1	5.4	−1	4.9
	DTR	11.0	4.9	12.7	4.3	10.3	5.1	6.8	3.8

, 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 (Table 4 and Table 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.

3.3. 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. Average of standard deviation (SD), centered root-mean-square error (RMSE), and Pearson correlation coefficient (CC) for T_max, T_min and DTR for five regionally representable AWSs, corresponding ERA-Land data, and closest climatological station.

			R1	R2	R3	R4	R5
Tmax	AWS	SD	10.2	10.4	10.3	10.2	9.8
	ERA5-Land	SD	9.1	9.3	9.1	9.2	9.3
		RMSE	2.0	1.9	2.6	2.1	1.9
		R	0.99	0.99	0.99	0.98	0.98
	Clim. Stat.	SD	9.6	10.0	9.5	9.8	9.7
		RMSE	1.4	1.1	1.5	1.5	1.1
		R	0.99	1.00	0.99	0.99	0.99
Tmin	AWS	SD	7.0	7.1	7.1	7.1	7.1
	ERA5-Land	SD	7.6	7.8	7.6	7.7	7.9
		RMSE	2.5	2.9	2.9	2.4	2.1
		R	0.97	0.97	0.95	0.97	0.97
	Clim. Stat.	SD	7.2	7.4	7.6	7.3	7.3
		RMSE	1.3	1.7	2.5	1.2	1.7
		R	0.99	0.98	0.95	0.99	0.97
DTR	AWS	SD	5.3	5.4	5.7	5.3	4.6
	ERA5-Land	SD	3.3	3.1	3.3	3.1	3.3
		RMSE	3.8	4.3	4.9	3.8	2.9
		R	0.84	0.86	0.84	0.85	0.81
	Clim. Stat.	SD	4.4	4.8	4.9	4.4	4.8
		RMSE	2.1	2.0	3.1	2.0	2.1
		R	0.94	0.94	0.88	0.94	0.92

.

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

4. Conclusions

This study emerged from the analysis of the AWS measurements conducted over six years (2013–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.