Monitoring of emerging water stress situations by thermal and vegetation indices in different almond cultivars

agronomy Article Monitoring of Emerging Water Stress Situations by Thermal and Vegetation Indices in Different Almond Cultivars Saray Gutiérrez-Gordillo 1,*,†, Javier de la Gala González-Santiago 2,3, Emiliano Trigo-Córdoba 1, Alfredo Emilio Rubio-Casal 4, Iván Francisco García-Tejero 1and Gregorio Egea 2,†   Citation: Gutiérrez-Gordillo, S.; de la Gala González-Santiago, J.; Trigo-Córdoba, E.; Rubio-Casal, A.E.; García-Tejero, I.F.; Egea, G. Monitoring of Emerging Water Stress Situations by Thermal and Vegetation Indices in Different Almond Cultivars. Agronomy 2021,11, 1419. https://doi.org/10.3390/ agronomy11071419 Academic Editor: Pasquale Losciale Received: 27 May 2021 Accepted: 12 July 2021 Published: 15 July 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Center “Las Torres”, Andalusian Institute of Agricultural and Fisheries Research and Training (IFAPA), Carretera Sevilla-Cazalla km 12.2, 41200 Sevilla, Spain; [email protected] (E.T.-C.); [email protected] (I.F.G.-T.) 2Area of Agroforestry Engineering, Technical School of Agricultural Engineering (ETSIA), Universidad de Sevilla, Ctra. Utrera km 1, 41013 Sevilla, Spain; [email protected] (J.d.l.G.G.-S.); [email protected] (G.E.) 3Agrocampo Ingenieros, Avda. Eduardo y Felipe Osborne, 2-11, 11500 Puerto de Santa María, Spain 4 Department of Plant Biology and Ecology, University of Seville, Ap 1095, 41080 Sevilla, Spain; [email protected] *Correspondence: saray.gutierrezgor[email protected]; Tel.: +34-671-532-852 † These authors had an equal contribution in this work. Abstract: In recent years, the area dedicated to modern irrigated almond plantations has increased significantly in Spain. However, the legal irrigation allocations are lower than the maximum water requirements of the crop in most cases. Therefore, almond growers are forced to implement regulated deficit irrigation strategies on their farms, applying water stress in certain resistant phenological periods and avoiding it in sensitive periods. Given the need to monitor the water status of the crop, especially in the most sensitive periods to water stress, the objective of this work was to evaluate the sensitivity of two UAV-based crop water status indicators to detect early water stress conditions in four almond cultivars. The field trial was conducted during 2020 in an experimental almond orchard, where two irrigation strategies were established: full irrigation (FI), which received 100% of irrigation requirements (IR), and regulated deficit irrigation (RDI), which received 70% of IR during the whole irrigation period except during the kernel-filling stage when received 40% IR. The UAV flights were performed on four selected dates of the irrigation season. The Crop Water Status Index (CWSI) and the Normalized Difference Vegetation Index (NDVI) were derived from thermal and multispectral images, respectively, and compared to classical water status indicators, i.e., stem water potential ( Ψstem ), stomatal conductance (g s ), and photosynthetic rate (A N ). Of the four flights performed, three corresponded to mild water stress conditions and a single flight was performed under moderate water stress conditions. Under mild water stress, CWSI was not able to capture the differences between FI and RDI trees that were observed with Ψstem . Under moderate stress conditions, CWSI was sensitive to the water deficit reached in the trees and showed significant differences among both irrigation treatments. No differences were observed in the CWSI and NVDI response to water stress among cultivars. Although NDVI and CWSI were sensitive to water stress, the low signal intensity observed in NDVI makes this index less robust than CWSI to monitor crop water stress. It can be concluded that UAV-based CWSI measurements are reliable to monitor almond water status, although for early (mild) levels of water stress, Ψstem seems to be the preferred option. Keywords: CWSI; NDVI; deficit irrigation; crop water status; unmanned aerial systems; emerging water stress 1. Introduction In recent years, the area devoted to almond crop (Prunus dulcis Mill. (D.A. Web)) has rapidly increased in Spain. While the area dedicated to this crop in 2019 was 822,878 ha, in 2020, there were 844,244 ha, which represents an annual increase of 3% [ 1 ]. In Andalusia, Agronomy 2021,11, 1419. https://doi.org/10.3390/agronomy11071419 https://www.mdpi.com/journal/agronomy Agronomy 2021,11, 1419 2 of 16 the most southern region of Spain, this crop occupies almost 190,000 ha, of which 31,800 ha are under irrigation [ 1 ]. Due to the excellent productive response of almond trees to irrigation, most of the new plantations are irrigated, although the scarcity of water in the region does not allow them to be irrigated to meet their full water requirements. Legal water allocations for this crop in the region are around 2000–2500 m 3 /ha [ 2 ], which makes it necessary to irrigate almond trees with deficit irrigation (DI) strategies. In this line, studies carried out by Lopez-Lopez et al. [ 3 ] established the water needs of the almond tree in the Guadalquivir River basin around 8000 m 3 / ha, confirming the need to define appropriate deficit irrigation strategies for this species. These strategies are designed to minimize the impact of the reduction in irrigation supplies on crop yield, maximizing water use efficiency and crop productivity for the available water resources, although at the expense of limiting crop development and its potential productions [ 4 ]. Although DI strategies may reduce yield, recent studies have shown that they can improve crop quality in several woody species [ 5 , 6 ], including almond [ 7 ], adding value to the crop. However, to ensure that the expected outcomes are obtained when applying DI strategies, proper monitoring of the crop’s water status is normally required to maintain water stress levels within the desired values [8]. Physiological variables such as leaf or stem water potential ( Ψleaf/stem ), stomatal conductance (g s ), or photosynthesis (A N ) have been traditionally used to monitor the crop physiological response to water deficit and to develop appropriate irrigation scheduling criteria. These leaf-level variables require a wide number of replications to be measured to avoid the scarce representativeness of one or two monitoring leaves. Thus, a proper assessment of crop physiological status would require sampling a representative number of leaves per plant and a representative number of plants per orchard to extract robust information [ 9 ]. Moreover, the use of this kind of measurements for irrigation scheduling in commercial orchards is even more difficult than in research trials because of the high soil and crop spatial heterogeneity frequently observed in commercial orchards [10]. Remote sensing techniques based on the use of unmanned aerial vehicles (UAVs) are an interesting alternative to manually determined physiological measurements, providing much more representative information on the existing spatial variability of crop performance and with a temporal resolution adaptable to management needs [ 11 ]. The physiological information obtained remotely by UAVs usually comes from thermal infrared sensors and optical sensors that measure the radiation reflected by the crop in either a reduced number of discrete wavebands (multispectral sensors) or multiple (>100) wavebands from a broader portion of the electromagnetic spectrum (hyperspectral sensors). Although this information can also be obtained by satellite-based remote sensing, UAVs have the advantage of being able to take images at the most appropriate time of the day (particularly relevant for thermal remote sensing) and to discriminate between vegetation and other types of surfaces (e.g., bare soil, cover crops grown in orchard alleys, etc.) in the captured image. As its main disadvantages, it can be pointed out the high cost of the equipment and the time required to obtain images of large farms. Another disadvantage is related to the legal regulations on the use of airspace for this type of vehicles (State Aviation Safety Agency (AESA) in Spain) [12], which prevents their use in certain areas. Thermal infrared sensors mounted on UAVs can provide information on the spatial– temporal variation of canopy temperature (T c ). This variable is related to crop water status when transformed in thermal indices that counteract the effect of certain environmental factors on T c [ 13 ]. The most widely used thermal index in agriculture is the Crop Water Stress Index (CWSI), which was introduced by Jackson et al. [ 14 ] and Idso et al. [ 15 ]. This index varies between 0 and 1 (denoting no water stress and maximum water stress conditions, respectively), and it is calculated from the difference between T c and the air temperature (T air ) on the day of measurement and two reference values, i.e., T wet and T dry , corresponding to the T c at maximum and minimum crop transpiration rates, respectively [ 16 ]. Some authors, such as Bellvert et al. [ 17 ], have already demonstrated the potential of CWSI derived from UAV-based thermal imaging as an appropriate indicator Agronomy 2021,11, 1419 3 of 16 for water stress monitoring in peach and citrus trees. González-Dugo et al. [ 18 ] also demonstrated that CWSI is a valuable method to assess the water status in citrus species but pointed out that short-term fluctuations in T c observed in these species must be carefully accounted for to avoid misinterpretation of CWSI values. Multispectral cameras mounted on UAV measure radiation reflectance in certain wavebands within the visible near-infrared (VIS-NIR) region with the aim of obtaining vegetation indices that can be related to crops performance, including crop water status [ 19 ]. One of the vegetation indices most commonly derived from remote sensing applications is the Normalized Difference Vegetation Index (NDVI), which varies between − 1 and 1 being the values closest to 0 those corresponding to extremely stressed/dead vegetation, those closest to 1 representing very healthy vegetation, and negative values denoting other nonvegetated surfaces (water, snow, etc.) [ 20 ]. In addition, to evaluate crop vegetative vigor, this index has also been used to estimate the crop water status in different woody crops such as peanuts [21], almonds [22], vineyards [23], or apple and citrus [24], among others. The development of tools and methodologies that allow early detection of water stress symptoms is essential to establish irrigation strategies that meet, at the same time, both productive and sustainable objectives [ 25 ]. This is especially true when applying regulated deficit irrigation (RDI) strategies [ 4 ], since there are phenological periods in which uncontrolled water stress can have negative consequences on fruit yield and quality, and over-irrigation can compromise the farm’s sustainability objectives. In the case of almond cultivation, there are numerous scientific evidences on the suitability of using RDI strategies when the water allocated for irrigation is lower than almond water needs [ 26 ]. However, for those phenological periods that are particularly sensitive to water stress and where water stress is sometimes recommended to be avoided, the development of tools for detecting incipient levels of water stress in almond could reduce the production losses associated with the application of RDI strategies [ 27 ]. The usefulness of thermal cameras mounted on UAVs to determine the water status of almond trees has been successfully evaluated [ 28 ]. However, neither their sensitivity for detecting incipient water stress episodes in almond trees, to be avoided in certain phenological periods, nor their comparison with indices obtained from multispectral images (e.g., NDVI) has been evaluated in field trials. Based on the above, the objective of this work was to evaluate the potential of infrared thermography and multispectral images taken from UAVs to detect situations of incipient water stress in almond trees that may compromise crop production when applying RDI strategies in almond orchards grown in drought-prone regions. 2. Materials and Methods 2.1. Experimental Site The experiment was carried out during 2020 in the Andalusian Institute of Agricultural and Fisheries Research and Training (IFAPA), Center “Las Torres” in an experimental almond orchard (Prunus dulcis Mill. (D.A. Webb)) cvs. Guara, Marta, Lauranne, and Marcona (Figure 1) grafted onto GN15 rootstock and located in the Guadalquivir River basin (SW Spain, 37 ◦ 30 0 38.55 00 N; 05 ◦ 57 0 44.98 00 W). Trees were planted in 2016, spaced 8 m ×6 m , and a sub-superficial drip irrigated using two pipelines at a depth of 15–20 cm with emitters of 2.3 L·h−1at 0.75 m intervals. The soil was a 2 m deep silty loam typical of Fluvisol, very homogeneous with an organic matter content of around 1.5%. Roots are situated predominately in the first 0.5 m layer of the soil profile, corresponding to the intended wetting depth. Soil water content values at field capacity ( − 0.033 MPa) and permanent wilting point ( − 1.5 MPa) are close to 0.40 and 0.15 m 3· m −3 , respectively. The climatic classification of the area is attenuated meso-Mediterranean, with an annual ET 0 rate of 1400 mm and an accumulated rainfall of 540 mm (average data corresponding to the last 15 years; obtained from the Andalusian Weather information Network [29]). Agronomy 2021,11, 1419 4 of 16 Agronomy 2021, 11, x FOR PEER REVIEW 4 of 17 layer of the soil profile, corresponding to the intended wetting depth. Soil water content values at field capacity (−0.033 MPa) and permanent wilting point (−1.5 MPa) are close to 0.40 and 0.15 m3·m−3, respectively. The climatic classification of the area is attenuated meso-Mediterranean, with an annual ET0 rate of 1400 mm and an accumulated rainfall of 540 mm (average data corresponding to the last 15 years; obtained from the Andalusian Weather information Network [29]). Figure 1. Experimental plot. Vegetative (A) phenological period and winter dormancy (B). 2.2. Irrigation Treatments Two irrigation treatments were designed: (i) a full-irrigated treatment (FI), which received 100% of irrigation requirements (IR), and (ii) a regulated deficit irrigation (RDI), which received 70% of IR during the whole irrigation period except during the kernelfilling stage (June to August) when it was irrigated with 40% IR. The irrigation doses were calculated using the methodology proposed by Allen et al. [30], for which the values of reference crop evapotranspiration (ETo) were multiplied with the crop coefficient (Kc) and a reduction coefficient (Kr) that depends on the horizontal projection of the tree shade (ground cover). ETo values were obtained according to the Penman–Monteith equation and the climatic variables recorded by a weather station (Davis Advance Pro2, Davis Instruments, Valencia, Spain) installed in the same experimental orchard. The crop coefficients obtained by García-Tejero et al. [31], and the reduction coefficients calculated with the expression proposed by Steduto et al. [32] were used (Table 1). The irrigation frequency was daily during the period of study. Table 1. Local crop reduction (Kr) and crop coefficient (Kc) values used in the experiment. Coefficients March April May Jun July August September October Kc 0.4 0.6 0.9 1.1 1.2 1.1 0.8 0.7 Kr 0.4 0.6 0.7 0.7 0.7 0.7 0.7 0.4 According to the climatic conditions, the irrigation doses applied were 476 and 184 mm for FI and RDI treatments, respectively. At the end of the season, these irrigation doses meant a water saving of 61% for the deficit treatment. 2.3. Plant Measurements and Indices Calculations During the period of maximum evaporative demand of the atmosphere (June to August), the monitoring of the water status of the trees was carried out through punctual measurements of stem water potential (Ψstem), net photosynthetic rate (AN), and stomatal conductance (gs), these measurements being taken between 13:00 and 14:30 (local time). Additionally, coinciding with these punctual measurements, four UAV flights were conducted on days of year (DOY) 176, 192, 199, and 216, which allowed us to obtain mean Figure 1. Experimental plot. Vegetative (A) phenological period and winter dormancy (B). 2.2. Irrigation Treatments Two irrigation treatments were designed: (i) a full-irrigated treatment (FI), which received 100% of irrigation requirements (IR), and (ii) a regulated deficit irrigation (RDI), which received 70% of IR during the whole irrigation period except during the kernel-filling stage (June to August) when it was irrigated with 40% IR. The irrigation doses were calculated using the methodology proposed by Allen et al. [ 30 ], for which the values of reference crop evapotranspiration (ETo) were multiplied with the crop coefficient (Kc) and a reduction coefficient (Kr) that depends on the horizontal projection of the tree shade (ground cover). ETo values were obtained according to the Penman–Monteith equation and the climatic variables recorded by a weather station (Davis Advance Pro2, Davis Instruments, Valencia, Spain) installed in the same experimental orchard. The crop coefficients obtained by García-Tejero et al. [ 31 ], and the reduction coefficients calculated with the expression proposed by Steduto et al. [ 32 ] were used (Table 1). The irrigation frequency was daily during the period of study. Table 1. Local crop reduction (Kr) and crop coefficient (Kc) values used in the experiment. Coefficients March April May June July August September October Kc0.4 0.6 0.9 1.1 1.2 1.1 0.8 0.7 Kr0.4 0.6 0.7 0.7 0.7 0.7 0.7 0.4 According to the climatic conditions, the irrigation doses applied were 476 and 184 mm for FI and RDI treatments, respectively. At the end of the season, these irrigation doses meant a water saving of 61% for the deficit treatment. 2.3. Plant Measurements and Indices Calculations During the period of maximum evaporative demand of the atmosphere (June to August), the monitoring of the water status of the trees was carried out through punctual measurements of stem water potential ( Ψstem ), net photosynthetic rate (A N ), and stomatal conductance (g s ), these measurements being taken between 13:00 and 14:30 (local time). Additionally, coinciding with these punctual measurements, four UAV flights were conducted on days of year (DOY) 176, 192, 199, and 216, which allowed us to obtain mean canopy temperature (Tc) and Normalized Difference Vegetation Index (NVDI) readings of the trees. The Ψstem was measured using a pressure chamber (Soil Moisture Equipment Corp., Sta. Barbara, CA, USA) in four trees per irrigation treatment and cultivar and one leaf per tree (n= 4). The selected leaves for Ψstem determination were mature, located at 1.5 m of height approximately, and NW exposed. About half an hour before readings, selected leaves were covered with aluminum–plastic bags to equalize the leaf-to-stem water potential. Agronomy 2021,11, 1419 5 of 16 In the same trees where Ψstem was measured, leaf gas-exchange measurements were conducted using a portable IRGA (ADC BioScientific Ltd., LCPro, Herts, England). Net photosynthetic rate (A N ) and stomatal conductance (g s ) were measured with an air CO 2 concentration of 400 ppm and a photosynthetic photon flux density between 1300 and 1500 µmol·m−2·s−1 . Similarly, to Ψstem readings, gas-exchange measurements were obtained in one leaf per monitored tree (n= 4), completely exposed to the sun, at 1.5 m of height and preferably south-east facing. A multispectral camera (model P4 multispectral, SZ DJI Technology Co., Ltd., Shenzhen , China) mounted on a multirotor UAV model Phantom P4 (SZ DJI Technology Co., Ltd., Shenzhen, China) equipped with a GNSS receptor was used to acquire multispectral images on the same dates and time as the crop physiological variables were measured. The flight missions were planned with the software DJI GS Pro (SZ DJI Technology Co., Ltd., Shenzhen, China), an iPAD (Apple Inc., California, USA) application for drone operations. The P4 multispectral camera includes one RGB camera and five cameras with a resolution of 2 megapixels that form a multispectral array covering the blue (450 nm ± 16 nm), green (560 nm ± 16 nm), red (650 nm ± 16 nm), red edge (730 nm ± 16 nm), and near-infrared, NIR (840 nm ± 26 nm) bands. The Normalized Difference Vegetation Index (NDVI) was calculated as Equation (1): NDVI =NIR −Red NIR +Red (1) The NDVI images (Figure 2) were then segmented with QGIS software (QGIS Development Team, 2009) to extract the pixels corresponding to pure vegetation. As no cover crop was grown in the orchard inter-rows, pixels with NDVI values lower than 0.45 were considered to belong to either bare soil or a combination of vegetation and bare soil (mixed pixels). The segmented images were transformed into vector layers in which each tree crown was represented with a polygon. Using the QGIS zonal statistics tool, the mean NDVI value of each tree in the experimental plot was calculated. Agronomy 2021, 11, x FOR PEER REVIEW 5 of 17 canopy temperature (Tc) and Normalized Difference Vegetation Index (NVDI) readings of the trees. The Ψstem was measured using a pressure chamber (Soil Moisture Equipment Corp., Sta. Barbara, CA, USA) in four trees per irrigation treatment and cultivar and one leaf per tree (n = 4). The selected leaves for Ψstem determination were mature, located at 1.5 m of height approximately, and NW exposed. About half an hour before readings, selected leaves were covered with aluminum–plastic bags to equalize the leaf-to-stem water potential. In the same trees where Ψstem was measured, leaf gas-exchange measurements were conducted using a portable IRGA (ADC BioScientific Ltd., LCPro, Herts, England). Net photosynthetic rate (AN) and stomatal conductance (gs) were measured with an air CO2 concentration of 400 ppm and a photosynthetic photon flux density between 1300 and 1500 µmol·m−2·s−1. Similarly, to Ψstem readings, gas-exchange measurements were obtained in one leaf per monitored tree (n = 4), completely exposed to the sun, at 1.5 m of height and preferably south-east facing. A multispectral camera (model P4 multispectral, SZ DJI Technology Co., Ltd., Shenzhen, China) mounted on a multirotor UAV model Phantom P4 (SZ DJI Technology Co., Ltd., Shenzhen, China) equipped with a GNSS receptor was used to acquire multispectral images on the same dates and time as the crop physiological variables were measured. The flight missions were planned with the software DJI GS Pro (SZ DJI Technology Co., Ltd., Shenzhen, China), an iPAD (Apple Inc., California, USA) application for drone operations. The P4 multispectral camera includes one RGB camera and five cameras with a resolution of 2 megapixels that form a multispectral array covering the blue (450 nm ± 16 nm), green (560 nm ± 16 nm), red (650 nm ± 16 nm), red edge (730 nm ± 16 nm), and nearinfrared, NIR (840 nm ± 26 nm) bands. The Normalized Difference Vegetation Index (NDVI) was calculated as Equation (1): NDVI = NIR − Red NIR + Red (1) The NDVI images (Figure 2) were then segmented with QGIS software (QGIS Development Team, 2009) to extract the pixels corresponding to pure vegetation. As no cover crop was grown in the orchard inter-rows, pixels with NDVI values lower than 0.45 were considered to belong to either bare soil or a combination of vegetation and bare soil (mixed pixels). The segmented images were transformed into vector layers in which each tree crown was represented with a polygon. Using the QGIS zonal statistics tool, the mean NDVI value of each tree in the experimental plot was calculated. Figure 2. NDVI orthomosaic of the experimental orchard printed on PNOA RGB orthomosaic (https://pnoa.ign.es, accessed on 6 April 2021). Day of year 216. Figure 2. NDVI orthomosaic of the experimental orchard printed on PNOA RGB orthomosaic (https://pnoa.ign.es, accessed on 6 April 2021). Day of year 216. Regarding the thermal readings, a thermal infrared (TIR) camera (Tau 2 324, FLIR Systems, Inc., Wilsonville, OR, USA) was mounted on a multirotor UAV model Phantom 4 Pro V.2.0 (SZ DJI Technology Co., Ltd., Shenzhen, China) equipped with a GNSS receptor. As with the multispectral flights, the flight missions were planned with the software DJI GS Pro (SZ DJI Technology Co., Ltd., Shenzhen, China). The TIR camera was installed to aim vertically downward (nadir view) at the bottom of the UAV. The camera spectral Agronomy 2021,11, 1419 6 of 16 range is 7.5–13.5 µ m with a resolution of 324 × 256 pixels, a focal length of 9 mm, and a field of view of 49 ◦ (H) × 39 ◦ (V). The UAV equipped with the thermal camera was also flown on the same dates and at the same times as the crop physiological variables were measured. The flight height was 30 m above the ground level, delivering thermal images with a ground spatial resolution of about 8 cm. The thermal images were stored on board in a raw format with 14-bit radiometric resolution. At the time of each flight, ground targets varying over a wide range of surface temperature were used for indirect calibration of the thermal imagery [ 16 , 33 ]. In particular, multiple ground-level temperatures were collected with a hand-held infrared thermometer (MI-210, Apogee Instruments Inc., Logan, UT, USA) from bare soil, selected almond canopies, and hot references (40 cm × 50 cm black plastic panels) located in the center of several inter-rows. The thermal images were segmented with QGIS software to extract the pixels corresponding to pure vegetation, for which the polygons representing the tree crowns previously created from the NDVI images were used. Mean canopy temperature (Tc) of each individual was computed using the QGIS zonal statistics tool. The temperature-based crop water stress index (CWSI) was then calculated for each tree crown as follows Equation (2): CWSI =(Tc−Ta)−(Tc−Ta)LL (Tc−Ta)UL −(Tc−Ta)LL (2) where T c− T a denotes the measured canopy-air temperature difference; (T c− T a ) LL is the lower limit of (T c− T a ) for a given vapor pressure deficit (VPD), which is equivalent to a canopy transpiring at its potential rate; (T c− T a ) UL is the maximum (T c− T a ), which corresponds to a non-transpiring canopy. The term (T c− T a ) LL is obtained from a linear equation that must be obtained empirically, the so-called non-water-stressed baseline (NWSB). In this study, the NWSB derived by González-Dugo et al. [ 27 ] for almond trees grown under similar environmental conditions to those prevailing in our experimental site has been used. The value of (T c− T a ) UL was set at 5 ◦ C, in accordance with the (Tc−Ta) observations of González-Dugo et al. [ 28 ] for severely stressed almond trees and as suggested for a great deal of species [16]. 2.4. Experimental Design and Statistical Analysis The experimental design consisted of two blocks, one of them being fully irrigated (FI) and the other one subjected to RDI strategies. In each of these blocks, the four cultivars (Guara, Marta, Lauranne, and Marcona) were randomly arranged, with three replications and four trees per replication and cultivar. Statistical analysis was done by using the Sigma Plot statistical software (version 14.0, Systat Software, Inc., San Jose, CA, USA). For each measurement day, an exploratory and descriptive analysis of physiological measurements ( Ψstem , A N and g s ) and indices (CWSI and NDVI) was made after applying a Levene’s test to check the variance homogeneity of the studied variables. For each cultivar and day of measurement, one-way ANOVA was performed to analyze the effect of irrigation treatment on the studied crop water status indicators. Two-way ANOVA was performed to analyze the interaction of cultivar and irrigation treatment on the crop water status indicators. The SigmaPlot statistical software (version 14.0, Systat Software, Inc., San Jose, CA, USA) was used in the analyses. 3. Results and Discussion 3.1. Weather Conditions Regarding the climatic conditions prevailing during the entire irrigation period (March to October) in which the trial was conducted (Table 2), the average temperature was 20.8 ◦ C with a cumulative rainfall and reference evapotranspiration (ETo) of 243 mm and 1032 mm , respectively. The cumulative crop evapotranspiration (ETc) was 758 mm. Agronomy 2021,11, 1419 7 of 16 Table 2. Monthly average values of the weather variables registered during the experimental irrigation season. Month Air Temperature Relative Humidity Radiation Rainfall ETo ETc March 14.0 71.7 15.9 66.2 80.6 13.0 April 15.9 75.5 16.1 71.6 84.9 34.4 May 21.5 61.5 24.5 45.0 153.1 130.8 June 23.4 52.2 27.9 0.0 167.0 148.7 July 28.9 46.2 27.3 2.0 190.5 182.6 August 25.0 32.9 25.3 0.8 167.1 143.7 September 20.1 58.6 19.8 19.8 118.5 80.2 October 17.3 65.8 15.1 37.4 70.2 24.1 Air temperature ( ◦ C); relative humidity (%); radiation (MJ · m −2 day −1 ); rainfall (mm); ETo, reference evapotranspiration (mm); ETc, crop evapotranspiration (mm). Table 3shows the weather conditions prevailing at the time of measurement for the 4 days when the punctual measurements and flights were carried out. The average air temperature was 34.78 ◦ C, 34.80 ◦ C, 38.50 ◦ C, and 34.56 ◦ C for DOY 176, 192, 199, and 216, respectively. The corresponding relative humidity values were 25.40, 31.47, 17.68, and 24.97, respectively. Finally, the average VPD values were 4.13, 3.82, 5.59, and 4.11, respectively, which may classify the measurement days into low (DOY 192), medium (DOY 176 and 216), and high (DOY 199) evaporative demand of the atmosphere. Table 3. Relative humidity (RH), air temperature (T), and vapor pressure deficit (VPD) on the measurement’s days. DOY Hour RH T VPD 176 11:00 27.6 34.47 3.94 11:30 25.03 34.82 4.16 12:00 24.43 34.84 4.20 12:30 24.52 34.98 4.23 192 11:00 36.25 33.24 3.24 11:30 32.73 34.34 3.64 12:00 30.96 35.1 3.89 12:30 25.93 36.52 4.51 199 11:00 18.72 38.06 5.39 11:30 18.45 38.22 5.45 12:00 17.02 38.62 5.67 12:30 16.52 39.08 5.84 216 11:00 27.1 33.58 3.78 11:30 25.55 34.35 4.03 12:00 23.51 34.9 4.26 12:30 23.71 35.39 4.37 DOY, day of the year; hour (GTM); relative humidity (%); temperature (◦C); VPD (kPa). 3.2. Charactetization of Cultivars Water Status Through Leaf-Level Physiological Measurements Figure 3shows the Ψstem values measured on the four sampling dates for the two irrigation treatments and the four cultivars evaluated. Ψstem values ranged from −0.99 MPa to − 1.31 MPa on DOY 176, − 0.71 MPa to − 1.13 MPa on DOY 192, − 1.25 MPa to −1.56 MPa on DOY 199, and − 0.99 MPa to − 1.38 MPa on DOY 216. The degree of water stress severity reached in the RDI treatment was determined by comparing the mean Ψstem values measured in the RDI trees ( Ψstem_RDI ) with the Ψstem values corresponding to wellwatered almond trees ( Ψstem_REF ) obtained from the reference baseline proposed by Shackel Agronomy 2021,11, 1419 8 of 16 et al. [ 34 ] ( Ψstem_REF = − 0.41–0.12 VPD). The ratio Ψstem_RDI / Ψstem_REF , i.e., the so-called water stress signal as defined by Goldhamer and Fereres [ 35 ], was used to determine the level of water stress severity reached in the RDI trees over the experimental period ( Table 4 ). Goldhamer and Fereres [ 35 ] found a mean stress signal of 1.37 in adult almond trees under mild water deficit. Using this value as a reference to distinguish between mild and moderate water stress conditions, it can be seen in Table 4that the water stress reached in the first three sampling dates can be classified as mild, while the stress level observed in the fourth sampling date can be classified as moderate. Agronomy 2021, 11, x FOR PEER REVIEW 9 of 17 Figure 3. Stem water potential (Ψstem) values for each cultivar and treatment. FI: control treatment, RDI: regulated deficit irrigation treatment. * Denotes significant differences (p < 0.05) between treatments for a particular day. Table 4. Mean vapor pressure deficit (VPD) during the measurement time, stem water potential of well-watered almond trees (Ystem_ref) calculated with the Ystem vs. VPD baseline proposed by Shackel et al. (1998), mean stem water potential of RDI trees (Ystem_RDI), and water stress signal calculated as Ystem_RDI/ Ystem_REF (Goldhamer and Fereres, 2001). DOY VPD Ψstem_REF Ψstem_RDI Water Stress Signal 176 4.13 −0.91 −1.24 1.37 192 3.82 −0.87 −1.05 1.20 199 5.59 −1.08 −1.45 1.34 216 4.11 −0.90 −1.35 1.49 DOY, day of the year; VPD (KPa); Ψstem_REF (MPa); Ψstem_RDI (MPa). The decrease in Ψstem with mild water stress levels has already been described for this species. For example, Barzegar et al. [36], in a trial with six almond cultivars (Azar, Marcona, Mission, Nonpariel, Sahand, and Supernovoa) reported that Ψstem decreased rapidly due to the osmotic adjustment performed by the crop to avoid a reduction of its photosynthetic capacity in the early stages of stress. Larsen et al. [37] also verified that almond leaves reduced the values of Ψstem during progressive water stress to avoid turgor losing. Despite the observed differences in Ψstem between FI and RDI treatments (Figure 3), gs values did not differ among both irrigation treatments in any of the four cultivars analyzed (Figure 4). The gs values measured in both treatments, with minimum values of around 0.25 mmol H2O·m−2·s−1 in both treatments, are typical of well-watered almond trees [38,39] indicating that the level of stress reached in the trees was not sufficient to generate stomatal closure in any of the cultivars. Figure 3. Stem water potential ( Ψstem ) values for each cultivar and treatment. FI: control treatment, RDI: regulated deficit irrigation treatment. * Denotes significant differences (p< 0.05) between treatments for a particular day. Table 4. Mean vapor pressure deficit (VPD) during the measurement time, stem water potential of well-watered almond trees (Y stem_REF ) calculated with the Y stem vs. VPD baseline proposed by Shackel et al. (1998), mean stem water potential of RDI trees (Y stem_RDI ), and water stress signal calculated as Ystem_RDI/ Ystem_REF (Goldhamer and Fereres, 2001). DOY VPD Ψstem_REF Ψstem_RDI Water Stress Signal 176 4.13 −0.91 −1.24 1.37 192 3.82 −0.87 −1.05 1.20 199 5.59 −1.08 −1.45 1.34 216 4.11 −0.90 −1.35 1.49 DOY, day of the year; VPD (KPa); Ψstem_REF (MPa); Ψstem_RDI (MPa). Although the RDI irrigation strategy applied a severe water deficit (40% ETc in the kernel-filling stage and 70% over the rest of the irrigation season), the level of water stress reached in the deficit treatment was mostly mild, with the exception of some days with Agronomy 2021,11, 1419 9 of 16 high evaporative demand in which the level of stress could be moderate (e.g., DOY 216). The reason why the water stress severity in RDI was not higher, as would have been expected with the volumes of water applied, is most likely due to the orchard proximity to the Guadalquivir River (<300 m) and the subsequent existence of relatively shallow groundwater resources. In any case, these conditions have been favorable to achieve the desired experimental conditions, i.e., trees submitted to mild to moderate water stress levels to evaluate the sensitivity of UAV-based remote sensing indices to capture incipient levels of water stress in almond trees. With the exception of Marta cultivar, which showed significant differences in Ψstem between FI and RDI only on the second day of measurement (DOY 192), the rest of the cultivars showed significant differences in Ψstem on DOY 176, 192, and 216 (Figure 3). None of the cultivars showed significant differences in Ψstem between FI and RDI on the third day of measurement (DOY 199), the one with the highest evaporative demand of the four sampling dates (Table 4). According to the Shackel et al. [ 33 ] reference baseline, the Ψstem_REF for DOY 199 is − 1.08 MPa (Table 4). In all cultivars, FI showed Ψstem values lower than Ψstem_REF , with water stress signals of 1.31 (Lauranne), 1.16 (Marcona), 1.19 (Guara), and 1.34 (Marta), which indicates that FI trees also presented mild water stress conditions on DOY 199. The reason why FI presented a mild water stress level on DOY 199 could be due to differences in the climatic conditions of the days used to estimate the irrigation needs and those actually existing on DOY 199. The decrease in Ψstem with mild water stress levels has already been described for this species. For example, Barzegar et al. [ 36 ], in a trial with six almond cultivars (Azar, Marcona, Mission, Nonpariel, Sahand, and Supernovoa) reported that Ψstem decreased rapidly due to the osmotic adjustment performed by the crop to avoid a reduction of its photosynthetic capacity in the early stages of stress. Larsen et al. [ 37 ] also verified that almond leaves reduced the values of Ψstem during progressive water stress to avoid turgor losing. Despite the observed differences in Ψstem between FI and RDI treatments ( Figure 3 ), g s values did not differ among both irrigation treatments in any of the four cultivars analyzed (Figure 4). The g s values measured in both treatments, with minimum values of around 0.25 mmol H 2 O · m −2· s −1 in both treatments, are typical of well-watered almond trees [ 38 , 39 ] indicating that the level of stress reached in the trees was not sufficient to generate stomatal closure in any of the cultivars. Regarding the first option, it is known that almond trees [ 38 , 39 ] and other species of the Prunus genus [ 40 ] show anisohydric behavior in their response to water deficit, characteristic of species with significant drops in Ψstem in response to water deficit as a consequence of poor stomatal control of plant water status [ 41 , 42 ]. However, in the case of almond trees, some authors also reported an isohydric behavior for this species [ 43 ], although our results suggest that the behavior of almond trees has followed an anisohydric strategy. Anisohydric plants need to reach a certain leaf water potential threshold so that stomata begin to regulate water loss [ 44 ], probably due to a low response of the stomata to the chemical signals involved in stomatal regulation (e.g., ABA) at the onset of water stress, or the existence of other mechanisms involved in stomatal regulation (e.g., hydraulic pathways) independent of chemical signaling [ 45 , 46 ]. In this line, Garcia-Tejero et al. [ 13 ] defined for almond trees two leaf water potential thresholds (shaded leaves) related to the response of g s to water stress, the first one established in about − 1.4 MPa denoting the onset of decreases in gs, and a second one set at about −2.0 MPa where severe reductions in g s are observed. The Ψstem values measured in the RDI treatments in our experiment were around the values of the first threshold defined by García-Tejero et al. [ 13 ], so it is suspected that the water stress levels reached by RDI trees are at the threshold at which stomatal regulation begins in this species. However, we should note that the water potential values measured by García-Tejero et al. [ 13 ] may differ from those measured in this work, since [ 13 ] measured water potential on shaded leaves while stem water potential was measured in this study (i.e., on covered leaves to allow leaf and stem water potential to Agronomy 2021,11, 1419 16 of 16 40. Steinberg, S.L.; Miller, J.C.; McFarland, M.J. Dry Matter Partitioning and Vegetative Growth of Young Peach Trees under Water Stress. Funct. Plant Biol. 1990,17, 23–36. [CrossRef] 41. Fernández, J.E. Understanding Olive Adaptation to Abiotic Stresses as a TOOL to increase Crop Performance. Environ. Exp. Bot. 2014,103, 158–179. [CrossRef] 42. Fu, X.; Meinzer, F.C. Metrics and Proxies for Stringency of Regulation of Plant Water Status (Iso/Anisohydry): A Global Data Set Reveals Coordination and Trade-Offs among Water Transport Traits. Tree Physiol. 2018,39, 122–134. [CrossRef] [PubMed] 43. Tardieu, F.; Simonneau, T. Variability among Species of Stomatal Control under Fluctuating Soil Water Status and Evaporative Demand: Modelling Isohydric and Anisohydric Behaviours. J. Exp. Bot. 1998,49, 419–432. [CrossRef] 44. Hamlyn, G.J. Irrigation Scheduling-Comparison of Soil, Plant and Atmosphere Monitoring Approaches. Act. Hortic. 2008 ,792, 16. 45. Hamlyn, G.J.; Tardieu, F. Modeling Water Relations of horticultural CROPS: A Review. Sci. Hortic. 1998,74, 21–46. 46. Limpus, S. Isohydric and Anisohydric Characterization of Vegetable Crops. The Classification of Vegetables by Their Physiological Responses to Water Stress. Project Report, 1st ed.; The Department of Primary Industries and Fisheries: Brisbane, QLD, Australia, 2009. 47. Espadafor, M.; Orgaz, F.; Testi, L.; Lorite, I.J.; González-Dugo, V.; Fereres, E. Responses of Transpiration and Transpiration Efficiency of Almond Trees to Moderate Water Deficits. Sci. Hortic. 2017,225, 6–14. [CrossRef] 48. Egea, G.; Verhoef, A.; Vidale, P.L. Towards an Improved and More Flexible Representation of Water Stress in Coupled Photosynthesis—Stomatal Conductance Models. Agric For. Meteorol. 2011,151, 1370–1384. [CrossRef] 49. Suarez, L.; González-Dugo, V.; Camino, C.; Hornero, A.; Zarco-Tejada, P.-J. Physical Model Inversion of the Green Spectral Region to Track Assimilation Rate in Almond Trees with an Airborne Nano-Hyperspectral Imager. Remote Sens. Environ. 2021 ,252, 112147. [CrossRef] 50. Caturegli, L.; Corniglia, M.; Gaetani, M.; Grossi, N.; Magni, S.; Migliazzi, M.; Angelini, L.; Mazzoncini, M.; Silvestri, N.; Fontanelli, M.; et al. Unmanned Aerial Vehicle to Estimate Nitrogen Status of Turfgrasses. PLoS ONE 2016 ,11, e0158268. [CrossRef] [PubMed] 51. Julie, C.A.; Zinnert, C.; Polo, M.-J.; Young, D.-R. NDVI as An Indicator for Changes in Water Availability to Woody Vegetation. Ecol. Indic. 2012,23, 290–300.