Diurnal and Seasonal Mapping of Water Deﬁcit Index and Evapotranspiration by an Unmanned Aerial System: A Case Study for Winter Wheat in Denmark

: Precision irrigation is a promising method to mitigate the impacts of drought stress on crop production with the optimal use of water resources. However, the reliable assessment of plant water status has not been adequately demonstrated, and unmanned aerial systems (UAS) offer great potential for spatiotemporal improvements. This study utilized UAS equipped with multispectral and thermal sensors to detect and quantify drought stress in winter wheat ( Triticum aestivum L.) using the Water Deﬁcit Index ( WDI ). Biennial ﬁeld experiments were conducted on coarse sand soil in Denmark and analyses were performed at both diurnal and seasonal timescales. The WDI was signiﬁcantly correlated with leaf stomatal conductance (R 2 = 0.61–0.73), and the correlation was weaker with leaf water potential (R 2 = 0.39–0.56) and topsoil water status (the highest R 2 of 0.68). A semi-physical model depicting the relationship between WDI and fraction of transpirable soil water ( FTSW ) in the root zone was derived with R 2 = 0.74. Moreover, WDI estimates were improved using an energy balance model with an iterative scheme to estimate the net radiation and land surface temperature, as well as the dual crop coefﬁcient. The diurnal variation in WDI revealed a pattern of the ratio of actual to potential evapotranspiration, being higher in the morning, decreasing at noon hours and ‘recovering’ in the afternoon. Future work should investigate the temporal upscaling of evapotranspiration, which may be used to develop methods for site-speciﬁc


Introduction
The aim of precision irrigation is to balance the plant water needs with the water supply [1] in space and time by sustaining transpiration at a desired level near the potential rate and ensuring all irrigation water is utilized by the plant [2]. Despite improvements in water use efficiency by improved varieties and agronomic managements, it remains challenging to assess crop water demand with regard to within-field variation of plant and soil properties, as well as microenvironmental conditions [3][4][5]. The initial crop response to lack of water, i.e., drought, is to save water by decreasing leaf expansion rate while activating a complex stomatal control [6]. Both reactions limit transpiration, but also reduce The aim of this study was to quantify temporally and spatially resolved crop water status on sandy soil under various levels of irrigation using a method with potential for practical on-farm application. The specific objectives were to: (i) Collect UAS thermal data and derive WDI maps for winter wheat on seasonal and diurnal scales, (ii) model the relationship of WDI with crop hydric physiology and root zone soil water status, (iii) validate the ability of WDI to quantify crop drought stress, and (iv) estimate the energy balance and improve the calculations based on ET partitioning between the crop and soil. The test crop was winter wheat, which is particularly interesting due to its long growth season spanning in the northern hemisphere from autumn until the following summer, and large consumption of irrigation water worldwide [31].

Experimental Design and Field Measurements
Field experiments with winter wheat (Triticum aestivum L.) were performed in 2018 and 2019 in central Denmark (56 • 31 59.6 N 9 • 24 39.4 E) on a sand soil (according to the USDA soil textural classification) containing clay, silt, fine and coarse sand of 3.8%, 6%, 32% and 56%, respectively, and organic matter of 2.2% in the topsoil (0-30 cm). The climate is temperate and humid, with cool to cold winters and relatively cool summers. Meteorological variables were measured in 10-min intervals at the nearby station located less than 1 km from the experimental site. Atmospheric conditions at the time of the UAS flight are shown in Table A1 in the Appendix A. The year of 2018 was unusually dry (the 2018 European drought) compared to 2019, which was meteorologically more normal. The winter wheat cv. Hereford was sown on 19 October 2017 and 18 September 2018, and phosphorus (P) and potassium (K) were supplied in early spring according to standard practice. The crop was grown in 30 × 30 m plots, which were replicated four times ( Figure A1 in the Appendix A) and arranged in a randomized split-plot design with irrigation as the main plot and nitrogen (N) fertilization as the subplot factor. Irrigation was either minimum (I0), optimum (I1) or full (I2) depending on the measured soil water deficit. However, in 2018, the difference in crop yield between I1 and I2 irrigations was not significant according to the mixed-effects statistical modelling. Further details about the irrigation amounts and dates are shown in Table A2 in the Appendix A. The N fertilization treatment consisted of basic (N0, 90 kg N ha −1 ) or full (N1, 203 kg N ha −1 ) irrigation for a target yield of 7 Mg ha −1 , split in equal doses in April and May. Crop yields were determined in early August by harvesting four areas of 1.5 × 10 m (15 m 2 ) inside each plot with a Haldrup combiner (Haldrup, Ilshofen, BW, Germany), oven-drying of straw and grain for 48 h at 60 • C and weighing for dry matter (DM) calculation.
Volumetric soil water content during the seasons was measured, covering all treatments by time domain reflectometry (TDR) probes installed vertically to 20 cm, 50 cm and 75 cm depth ( Figure A1) [32]. Field capacity (FC) was measured in March after a period with rain and subsequent dry weather for 2-3 days. For the 0-50 cm soil depth, the FC ranged between 66 mm and 100 mm. The plant total available water (TAW) was calculated as the plot-specific value of FC minus a fixed value of unavailable water for the soil type based on Hansen (1976) [33]. The depletion of available water (water deficit) during the season in each plot equipped with TDR was calculated by subtracting the measured soil water content from the TAW. The fraction of transpirable soil water (FTSW) was calculated by dividing the measured soil water content by TAW.
Stomatal conductance (g s ) was measured on the youngest fully developed and surfacedry leaf on three plants per replicate with a leaf porometer (Model SC-1, Decagon, Pullman, WA, USA) on both the adaxial and the abaxial surfaces at full sunlight. The LWP was measured after each g s measurement on the same leaf using the pressure chamber method (Soil Moisture Equipment, Santa Barbara, CA, USA) [34].

Unmanned Aerial System and Acquisition of Multispectral and Thermal Images
The UAS consisted of a Quadrotor DJI Matrice 100 equipped with DJI Zenmuse X3 (12 Mpix true colour RGB), RedEdge multispectral (1.2 Mpix with RGB, red edge and nearinfrared bands centered at, 475 nm, 560 nm, 668 nm, 717 nm and 840 nm, respectively) and DRONExpert thermal sensors (640 × 480 pixels). These were used to collect images over the field from May until harvest in August at flight altitude of 50-70 m. The multispectral and the thermal flights were performed separately.

Data Acquisition Times
The seasonal response of winter wheat to treatments was studied in 2018 by engaging the UAS on 5 days, i.e., 15, 22, and 31 May; 7 June; and 2 July. The flights on 31 May and 7 June were conducted during irrigation to assess the reliability of the "wet" threshold of the WDI. The flight on 2 July was conducted closer to the maturity period. In order to study diurnal responses, single-day flights were conducted in June in both years over four plots. In June 2018, six flights were conducted between 10:00 and 14:00 (10:15, 11:00, 11:30, 12:00, 12:30, 13:00, 13:30, 14:00), and in 2019, four flights were conducted between 11:00 and 15:00 (11:49, 12:13, 13:20, 14:49). In both years, the winter wheat was at a comparable growth stage of late flowering/early grain filling (65-75 on the BBCH scale) [35]. The meteorological conditions during the season in 2018, as well as the flight days in 2018 and 2019, are shown in the Appendix A.

Data Processing
The multispectral and the thermal images were processed in Pix4Dmapper software (Pix4D SA, Lausanne, Switzerland) using "Ag Multispectral" and "Thermal camera" mode, respectively. The outputs were high-resolution GeoTIFF images. For all images, radiometric calibration was performed with reference panel for the multispectral data and four reference plates (white, gray, dark gray and black to simulate different temperatures) for the thermal data recorded prior to the flights. For the thermal flights without recorded calibration plates, an equation derived from the previous flight dates was used (y = x*1.428 − 9.193, where x is calibration panel digital number and y is corresponding surface temperature). All images were georeferenced with ground control points fixed in the field, resampled to 10 cm resolution and projected in the same grid and coordinate system (ETRS89/UTM zone 32N, code EPSG:25832).

Calculation of the Water Deficit Index
The VIT trapezoid was used as a reference to distinguish between fully saturated and dry surfaces regarding the fraction of vegetation cover [21]. The method combines the multispectral vegetation index (VI) with the temperature anomaly (surface-air difference) to locate pixels within a trapezoid limited by four baselines points, i.e., (1) well-watered and fully transpiring crop, (2) totally dry and non-transpiring crop, (3) saturated bare soil and (4) dry bare soil. The upper and the lower VI boundaries should be assessed from the maximum and the minimum values of the VI maps. Assuming that surface-air temperature difference (T S − T a ) is a linear function of vegetation cover, the line connecting the well-watered crop to saturated soil forms the "wet" baseline, and the line connecting the non-transpiring crop with dry bare soil forms the "dry" baseline (e.g., Figure 1 with the corresponding values of baseline points in Table 1).  Table A2 in the Appendix A for irrigation time and amount).  According to Moran et al. (1994), the WDI relates with ET a and ET p according to Equation (Equation (1)): where ET a and ET P are actual and potential evapotranspiration rate (mm h −1 ), respectively; T S and T a are surface and air temperature ( • C_, respectively; and subscripts "min", "max" and "mes" refer to minimum, maximum and measured, respectively. For all pixels not fully covered by canopy, the WDI equals to the horizontal distance ratio of a given pixel coordinate (T S − T a , f c ) to the "wet" and "dry" baselines. The equation for the "wet" baseline is x = a + b*y and for "dry" is x = c + d*y, where x is (T S − T a ), y is fc. Thus, it follows for the WDI: where f c is fraction of vegetation cover, a and b are slope and intercept of "wet" baseline, and c and d are slope and intercept of the "dry" baseline. The calculations for points 1 and 2 are identical to those of the CWSI "wet" (Equation (4)) and "dry" (Equation (5) where r a is aerodynamic resistance (sm −1 ), R n is net radiation (Wm −2 ), G is soil ground heat flux (Wm −2 _, C v is volumetric heat capacity of air (1200 J • C −1 ), ρ is air density (kg m −3 ), VPD is vapor pressure deficit (kPa), γ is the psychrometric constant (0.067 kPa • C −1 ), ∆ is the slope of saturation vapor pressure curve (kPa • C −1 ) and r c is canopy resistance to vapor transport (s m −1 ). The equations for these variables are shown in the Appendix B. For a full covered non-stressed canopy that transpires near to or at potential level and therefore has minimum r c = r cp (VIT point 1), the following Equation applies: For a full-covered and stressed, i.e., non-transpiring canopy with maximum r c = r cx (VIT point 2), the temperature anomaly equals: For saturated bare soil, there is no canopy resistance, i.e., r c , r cp and r cx in Equations (6) and (7), respectively, equal zero (VIT point 3). However, there is resistance to heat transfer in the boundary layer immediately above the soil surface (r s ) that should be added to r a [19] according to: Finally, dry bare soil conditions (VIT point 4) are represented by Equation (4). Due to the highly variable conditions present for dry soil mixed with f veg , Rn and G should be calculated separately, whereas ρ can be omitted: As a measure of f veg , the Normalized Difference Vegetation Index (NDVI) can be used due to its high correlation with f veg [25].
Net radiation (R n ) is the balance (difference) between incoming solar radiation R s (global radiation) and outgoing (emitted and reflected) radiation in the shortwave and longwave regions. R n can be estimated as a fraction of R s , e.g., 68% for cereal crops [36].
When measured, R n refers to a vegetated surface well supplied with water. Moran et al. (1993) [21] pointed out the importance of caution when using identical R n for bare soil. In order to improve calculations for dry bare soil conditions, R n was calculated by the Stefan-Boltzmann law, iterating the soil temperature: where T s in this equation is soil temperature ( • C); α is albedo (~0.12 for dry bare soil derived from UAV images); σ is the Stefan-Boltzmann constant σ = 5.67*10 −8 [W·m −2 ·K −4 ]; and ε s and ε s is the emissivity of soil and air, respectively (see Appendix B). The energy used for soil heating (G) is positive for warming and negative for cooling, and estimates typically range around 10% of R n , although this value may be very different for full cover and bare soil conditions. Therefore, in the present study, G was calculated with regard to a fraction of vegetation cover following Kustas and Daughtry (1990) [37]: where the NDVI was fixed to 0.9 for full-cover winter wheat and 0.2 for bare soil according to the field image estimates. Further details on the energy balance parameters calculations are provided in Appendix B.

Calculation of Seasonal and Diurnal Evapotranspiration
According to Allen et al. (1998) [27], ET p is estimated by 'correcting' the reference evapotranspiration (ET 0 ) of a well-managed grass with a coefficient (K c ) that integrates crop and climate characteristics, i.e., ET p = K c *ET 0 . In order to find ET a where the crop canopy is not full and drought stress may occur, a dual Kc model can be used by explicitly considering crop transpiration and soil evaporation: where K cb is basal crop coefficient, K s is the drought stress coefficient and K e is the evaporation coefficient of the bare soil. The units of ET a and ET 0 are mm h −1 , and the latter was calculated by the Penman-Monteith Equation [27]: where T a(hr) is the mean hourly air temperature ( • C), u 2 is the wind speed at 2 m height (m s −1 ), e 0 Ta(hr) is the saturation vapor pressure at the temperature T a (kPa) and e a is actual vapor pressure [kPa].
The crop coefficient for winter wheat in Denmark during the mid-season varies from 1.09 to 1.15 depending on the growth conditions [28]. In this study, the mid-season crop coefficient K cb full (at maximum plant cover and height) was assumed to be 1.1. For the areas without full crop cover, K cb was calculated by the following formula: where K c min is minimum K c for the bare soil assumed to be 0.15, f veg is the fraction of soil surface covered by vegetation [0.01-1] and h is plant height (m).
According to Gutman and Ignatov (1998) [38], f veg can be calculated from the NDVI: where NDVI min and NDVI max are the minimum and maximum NDVI values corresponding to bare soil and full vegetation cover, respectively. The drought stress coefficient K s was calculated using the WDI similarly to Tang, Han and Zhang (2019) [25] who used CWSI: The evaporation fraction from bare soil (K e ) is highly influenced by the fraction of soil covered by vegetation, i.e., K e depends on f veg [25] and is highly dependent on soil moisture in the top layer. In the absence of available water, K e will be close to 0. Thus, it is possible to calculate K e using the WDI: All calculations of the VIT trapezoid and ET were conducted for each flight at a halfhourly timestep at the time of the UAS flight. Further details on the parameter calculations are provided in Appendix B.

Calculation of Fraction of Transpirable Soil Water
The FTSW is the ratio between available and total transpirable soil water, with the latter being equal to the difference between FC and WP. The relation between the rate of transpiration ET a /ET p (and, correspondingly, 1-WDI) and the FTSW can be conceptualized according to Raes et al., (2009) [39] as a suite of curves differing in the level of ET p [27]. At higher ET p , a larger FTSW is needed for the crop to avoid drought stress, whereas at lower ET p , the crop can thrive at smaller FTSW. Whereas these studies considered daily values of these variables and relations, in the current study, the scheme was adopted to hourly values relevant to the diurnal WDI analysis. Analogous to Raes et al. (2009) [39], Equation (16) was iterated to fit the constant a, allowing the estimation of the FTSW from the WDI and the validation of the WDI: where c is the slope of the curve with lowest ET p (ET min ). Equation (18) can be rearranged to isolate FTSW: However, caution is needed when using Equations (16) and (17). For instance, dry soil receiving a few millimeters of rain does hold freely available water, allowing a crop to transpire at its potential rate despite FTSW being small for the entire root zone.
Using the approach to estimate FTSW makes it possible to calculate the amount of water needed to reach field capacity, i.e., the Soil Water Deficit (SWD, mm) as:

Statistical Analyses
Linear regression modelling was used to assess the relation between the WDI and soil water content and SWD. The coefficient of determination (R 2 ) was calculated to assess the strength of the relationships at 95% significance level. The daily WDI maps for the 2018 season were also correlated with the measured soil water content and the estimated SWD.
For these correlations, mean WDI per plot was used. Furthermore, WDI was validated by comparing measured WDI values to modelled ones according to Equation (16) using FTSW.
The diurnal WDI maps for 6 June 2018 and 28 June 2019 were also correlated with LWP and g s measured during the flights. For these correlations, mean WDI of 1 m 2 around the measurement locations was used. All image analysis was conducted in ArcMap (10.6.; ESRI, Redlands, CA, USA), whereas WDI, evapotranspiration maps and statistical evaluation were performed in R (ver. 3.5.1; R Core Team, 2013).

Vegetation Index/Temperature Trapezoid on Seasonal and Diurnal Scale
Parameterization of the VIT trapezoid for the 2018 season and for the diurnal flights in 2018 and 2019 is shown in Table 1. The mean "wet" reference for the crop-air temperature difference (VIT point 1) in 2018 was around −5 • C, while in 2019, it was close to 0 • C due to higher humidity (more than 60% on 28 June 2019). Lower relative humidity and higher wind speeds results in larger difference between the "wet" and "dry" baselines. The typical "dry" baseline for the crop (VIT point 2) was 4-5 • C higher than Ta (Table 1). The saturated bare soil-surface temperature difference (VIT point 3) was slightly higher due to the higher surface resistance and less transpiration compared to evapotranspiration. The largest variation in the temperature difference was seen for the "dry" soil conditions (VIT point 4) as, in the absence of available water vapor, surface temperature is defined only by accumulated incoming radiation and can be greatly reduced by wind. Therefore, under low wind speed, the aerodynamic resistance r a of soil surface to water vapor and its heat gradually increases and vice versa. The temperatures of 13 • C (on 22 May 2018, Table 1) and 25 • C (on 06 June 2018 at 12:30, Table 1) above T a correspond to conductions with higher and lower wind speed conditions, respectively. The iteration of Rn for VIT point 4 (Equation (8)) resulted in more reasonable values of T, which were otherwise underestimated, resulting in many measurement points positioned left of the dry boundary outside the trapezoid.

Water Deficit Index at Seasonal Scale and Correlation to Field Data
The VIT trapezes developed by the seasonal data in 2018 are shown in Figure 1 and the associated NDVI and WDI maps are shown in Figure 2. Due to the high evaporative demand during May and June 2018, it was not possible to follow the planned irrigation schedule. Consequently, treatment I0 was minimally irrigated, while I1 and I2 treatments suffered periodic drought stress despite irrigation (Appendix A, Table A2). The NDVI distribution indicates high variability in fc, i.e., mixed pixels and less distinct differences between treatments, which was expected with the planned irrigation schedule (Figure 2a). The flight on 22 May was conducted after irrigation. Thus, the winter wheat had lower WDI values, and most pixels were located close to the 'wet' reference in VIT (Figure 1). However, in mid-May, the winter wheat already showed signs of drought stress, as indicated by the low NDVI and high WDI (Figure 2), due to the prolonged dry spell from the end of April and throughout May, as well as limited irrigation. The grain yield in 2018 was low for all treatments, i.e., only 2.7 Mg ha −1 , 3.5 Mg ha −1 and 3.2 Mg ha −1 for I0, I1 and I2 respectively, which is almost half compared to the corresponding 2019 values of 5.9 Mg ha −1 , 5.9 Mg ha −1 and 6.0 Mg ha −1 .
On 31 May 2018, the flight was conducted when one-third of the field was already irrigated the day prior to the flight (right row) one-third of the field was irrigated the same morning as the flight (left row) and one-third of the field was not irrigated yet (middle row) (Figure 2b). The already irrigated plots transpired below potential level, thus showing early stress symptoms under the prevailing high evaporative demand. On 7 June 2018, the irrigation machine was present in the field during the flight. Taken as a whole, in 2018, insufficient irrigation under drought made treatment differences in yield rather small according to the analysis of variance (results not shown).  Table A2 in the Appendix A for irrigation time and amount).
There were significant but medium-strength correlations between the WDI and topsoil water status, with mostly stronger correlation to soil water content (SWC) at 0-20 cm compared to the soil water deficit (SWD). The

Relation between the WDI and the Fraction of Transpirable Soil Water and ET p
The relation between relative transpiration rate, i.e., ET a /ET p as predicted by (1-WDI), and FTSW in the 0-50 cm soil depth, which comprises the majority of the root zone, is shown in Figure 4a. It was possible to describe the relationship with one linear equation for the data measured on 15 and 31 May and 7 June 2018, which had similar ET p values of ca. 0.45 mm h −1 . However, the data from 22 May 2018 had a larger slope, which is likely related to the lower ET p that day. Based on Figure 4a, the relation for the lowest ET p (ET min = 0.39 mm h −1 and slope c = 2.31) is shown in Equation (17), whereby an optimal value of a = 3.31 was found by iteration. Using these parameters in Equation (17), the values (1-WDI) were derived from the measured FTSW, which fitted quite well with the average (1-WDI) per plot derived from the thermal imagery (Figure 4b).

Modelling Water Deficit Index at Diurnal Scale
The VITs derived by the diurnal data from 2018 and 2019 are shown in Figures 5 and 6, respectively, and detailed WDI maps for these diurnal flights are shown in Figures A5 and A6 in the Appendix C. In 2018, for the morning flight (10:15), certain pixel values were located outside the trapezoid, i.e., being "wetter" than the non-stressed baseline. These differences were likely due to free water present on the leaves as dew formed during the night, which was confirmed by relative humidity of almost 100% according to the meteorological data. As water evaporated, the WDI started to depict the actual plant water status, implying that the thermal data acquisition should not be conducted in the morning to ensure the absence of free water otherwise leading to skewed WDI. Similar trend was observed in 2019 with WDI increasing at noon hours.  In both years, there was no consistent response of LWP and g s to treatment (Figure 7). In 2018, the average LWP was −2.35 and −1.53 MPa for plots under low and higher irrigation, respectively, while the corresponding g s was 101 and 160 mmol m −2 s −1 . In 2019, a small plot (3 × 2 m) irrigated separately the day before the flight showed average LWP of −1.11 MPa compared to −1.43 MPa for the non-irrigated plot. The corresponding g s was more distinct, with 559 and 319 mmol m −2 s −1 , respectively. The WDI was significantly but weakly to moderately correlated with the measured leaf physiological parameters, with stronger correlation for g s compared to LWP (Figure 7).

High-Resolution Mapping of Actual Evapotranspiration at the Seasonal and Diurnal Scale
The spatiotemporal distributions of ET a for the winter wheat using the dual K c approach (Equations (10)-(15)) and the derived WDI for the 2018 season are shown in Figure 8 and for the diurnal scale in Figures 9 and 10 for 2018 and 2019, respectively. On most days, ET 0 at the time of the flight was around 0.5 mm h −1 (Table 1), with the lowest value of 0.35 mm h −1 on 22 May due to the lower incoming radiation at the time of flight. Except for the flights performed during irrigation, ET a did not exceed 0.3 mm h −1 , which pointed to low water availability and transpiration even for the full irrigation treatment. Diurnally, ET a decreased with time during the day and certain areas in the plots. Lower soil water deficit indicated a 'recovery' in evapotranspiration rate in the afternoon (Figures 9 and 10), suggesting that the ET a rate is not constant during the day and irrigation planning from a single mid-day image may be complicated.

Seasonal Dynamics of the Water Deficit Index and Correlation to Soil Water Status
Monitoring and managing drought on time is important to allow crop physiological recovery, whereas delayed response is associated with irreversible damage and yield loss [40]. The spring of 2018 was exceptionally warm and dry, and crop yields were inevitably reduced by almost half (possibly also due to heat stress) despite irrigation during the season in the I1 and I2 treatments. Modelling the relation between the WDI and FTSW in the root zone was conducted to validate the ability of the WDI to estimate the FTSW and describe the rate of ET a (Equations (16) and (17)). The model has potential to be used both ways, either deriving the FTSW using the WDI or vice versa, depending on aim and available measurements. The WDI derived from thermal imagery was significantly correlated to the WDI modelled from the FTSW (R 2 = 0.74), and knowing the FTSW along with soil FC and WP allows the assessment of the exact amount of water needed for irrigation (Equation (17)). The FTSW level needed for a crop to transpire at a potential level depends on the ET p . Therefore, future studies may explore larger variations of ET p and soil water contents by conducting UAS flights over longer periods or studying different crop systems. As noted earlier, the FTSW estimation based on WDI measurements can only be assumed to be valid during prolonged drying periods as was the case in this study during 2018. This is because a minor rainfall, which does not bring the FTSW to unity, is expected to bring the WDI to 0 due to the presence of easily available water in the root zone.
It should be mentioned that all calculations were conducted at the time of the flight in 30 min and are expressed at the hourly scale (e.g., ETa in mm h −1 ), as the typical flight duration was 20 min and mean meteorological conditions for a 30-min timeframe are the most representable. Soil moisture was the sole daily variable either measured with TDR (not exactly overlapping flight time/hour) or calculated by simple soil water balance [33] when TDR readings were done some days apart. It was assumed that WDI and ETa seasonal measurements to represent their values at the time of a flight as reliable upscale was neither possible nor objective of the study.
Under prolonged drought, winter wheat fastens its development, shortens its season, but also reduces ET a [7]. It was possible to distinguish between an irrigation that was done 1 day before the flight and the same-day irrigation (31 May 2018; Figure 2b), as winter wheat irrigated the day before already showed a higher WDI than wheat irrigated on the same day. Similarly, on 7 June 2018, when the flight was performed during irrigation, plots that were irrigated earlier on the same day already showed signs of drought stress. The WDI was higher than 0, hinting that soil at field capacity could not even fulfill the evaporative demand in this hot and dry season on the coarse sand. This has important implications for the use of WDI and it shows its dependence on soil water status (the "wet"-"dry" soil baseline). Colaizzi et al. (2003) and Köksal, (2008) [24,41] also found a relation between WDI and soil moisture. However, the WDI variation within each plot was high, and the mean WDI resulted in medium correlations to the point TDR measurements (Figure 3). Wang et al. (2018) [42] pointed out that, under dense vegetation, canopy temperature is better correlated with root-zone soil moisture than with the surface layer water content. Therefore, WDI areas with less crop cover represent the availability of the water in the topsoil and do not reflect the deeper soil water status. In this study, there was significant correlation only to the topsoil moisture in 2018 with R 2 = 0.68 on 15 May for water content and R 2 =0.6 on 2 June for water deficit. Hoffmann et al. (2016) [19] found a similar relation of surface-air temperature difference to soil water content, with R 2 values between 0.56 and 0.83. Other soil types and soil water conditions at various depths need to be explored to search for stronger correlation or a proxy for ET a under different conditions.

Water Deficit Index on the Diurnal Scale and Correlation to Leaf Water Status
On 6 June 2018, relative humidity during the night was close to 100%, which led to dew accumulation on the winter wheat leaves and lowered the temperature compared to the VIT estimation ( Figure 5). This can also be explained by the ability of the plant to replenish water during the night from the deeper levels of the soil profile. After 11:00, as the surface temperature increased, it gradually moved the pixels of the VIT range toward the "dry" WDI baseline. This led to the conclusion that UAS images should not be acquired too early in the day (e.g., 10:00) because the presence of dew on the surface of the plant may skew the WDI toward less stressed. Areas with winter wheat having more water available for transpiration recovered faster in the afternoon. In contract, plants water-limited areas remained stressed for a longer time (Figures 5 and 6).
The UAS flights are typically performed at or around noon to assess plant transpiration [43]. As shown in Figures 5 and 6, plants demonstrated stress between 12:00 and 13:00, likely due to the high air temperature triggering reduced transpiration to avoid high water loss. In addition, lower wind speed contributed to higher r s − aerodynamic resistance of heat and water vapor from the surface into the air above the canopy, leading to an increase in T c and vice versa. As the afternoon progressed, plants that are less stressed started to recover, while the others showed high levels of stress as canopy pixels moved in the VIT space from the "dry" toward the "wet" baseline ( Figures 5 and 6). The observed "partial afternoon recovery" (decrease of the WDI at 14:50, Figures A5 and A6 in the Appendix C) follows the theory that plants rehydrate during the afternoon and night. Transpiration starts in the morning, leading to lower plant water potentials and, in turn, the partial closure of the stomata. At some stage during the afternoon, soil water uptake exceeds transpiration due to high water potential gradient between plant and soil, leading to the partial recovery of transpiration and hence lower WDI as simulated by Jensen et al. (1993) [13] and shown in Figure 7 for drought-stressed barley. Therefore, afternoon UAS acquisition may be useful to find areas in the field that have higher irrigation demands.
Significant correlation to g s with R 2 up to 0.73 demonstrates the ability of the WDI to detect early crop drought stress signals (Figure 7). Stomatal closure due to increased temperature and limited water supply limits ET a and increases T c . The diurnal changes in T c are directly caused by plant transpiration [44]. In addition, g s is one of the common plant physiological indicators used to determine crop water status [13]. In 2018 and 2019, all treatments had low g s values of around 200 mmol m −2 s −1 (Figure 7a,b) compared to plot that was irrigated the day before the measurements, with g s values up to 900 mmol m −2 s −1 (Figure 7b) indicating drought stress and limited transpiration in plots with low g s . The g s and crop transpiration may vary between different cultivars of the same species [7]. Therefore, the WDI method should be further investigated as a potential phenotyping tool.
Some studies have suggested the possibility of using Tc and derived thermal indices to distinguish between water status and the predicting index from LWP values [45]. It may be difficult to derive LWP using the WDI due to the different dynamics in LWP diurnal pattern depending on the weather conditions and stress level [46]. In this study, significant but weak correlation of WDI to LWP was found (maximum R 2 = 0.56). This can be attributed to the spatial heterogeneity in the crop growth which is commonly seen in the field under same management conditions, i.e., the subtle spatial variation.

Links between Evapotranspiration and Water Deficit Index
Reference evapotranspiration calculated for full-crop conditions must be adjusted for fc analogous to the WDI (only assuming that the relationship between crop and soil is linear). This was achieved using the dual K c approach according to Allen et al. (1998) [27], which accounts for both crop transpiration and soil evaporation with regard to canopy cover and stress index. The two-source approach provides high ET a estimation accuracy [47], especially when the canopy is less dense and the soil evaporations has a larger contribution to ET a . As it can be seen from the diurnal responses (Figures 9 and 10), the ratio of ET a to ET p during the day was not constant as the WDI values changed from the morning to the afternoon. The change in the ET a /ET p ratio during the day may add uncertainties to quantifying crop water demand if using a single-shot image. It may complicate the use of the WDI for irrigation scheduling if the time of the flight for drought stress detection is not chosen correctly. On sandy soils and even under ample irrigation, water availability is limited for the noon hours when ET p reaches higher values. Previous studies have upscaled daily 'single-shot' images acquired by UAS or satellites in order to derive daily ET products, which may lead to different results depending on the time of day and instantaneous meteorological conditions [48], therefore accounting for the change in ET fraction improved by the upscaling [49]. In this study, the relation of th WDI to the FTSW clearly shows the potential use of WDI in ET a calculations since the FTSW can be used to calculate K s , as shown by Equation (10) [27]. Future studies should assess the time upscaling issue to improve estimation of crop water demand. Further refinement of the ET a calculations may be possible using nonlinear boundaries for the "wet" and "dry" baselines [50].
On the seasonal timescale, the WDI depicts the daily crop growth and supports identifying areas in the field that are prone early drought stress, prompting irrigation on a particular day. On the diurnal scale, the WDI depicts the crop behavior more precisely, which can provide more insight on the crop status and water need. Studies has revealed a 'crosstalk' between the drought stress and diurnal regulation of a range of genes differing in the vegetative and reproductive stage of crop development (e.g., [11]). At high ET p , sandy soils may not be able to fulfill crop demand to resupply water during noon heat wave. It may happen even if there is sufficient amount of irrigation present [12] and may further complicate upscaling of crop stress maps from instantaneous to daily timescales [51]. Although our study indicates solar noon and early afternoon as the best time, on the question remains of when exactly to schedule the UAS flight to capture the most representative pictures that efficiently outline the stressed areas from non-stressed. Another question remains on how the WDI can be applied in cloudy conditions since it does takes into account all atmospheric parameters [19], and whether the ability to distinguish drought-stressed plants and areas is compromised in the absence of high incoming solar radiation. This is important to know in areas where clouds occur frequently during the season. Further improvements would include the calculations of the energy balance, especially components that are changing rapidly and therefore harder to take into account (e.g., aerodynamic resistance for crop and soil surfaces). In addition, there is a need for a precise decision-making system for the irrigation need based on the WDI value thresholds.

Conclusions
This study assessed the drought stress of winter wheat grown on sandy soil in Denmark through high-resolution WDI and ET a maps and associated energy-balance and semiphysical modelling. The analysis was conducted on seasonal and diurnal variations with correlations to biophysical variables to increase the reliability and performance of drought management. Additional attention was paid to the separation of R n and G for crop and soil in order to avoid uncertainties in calculations due to mixed pixels. The conclusions can be summarized as follows.
(i) The WDI has high potential to accurately depict soil water status. It is sensitive to the topsoil water status. A semiphysical relationship was established between the FTSW in the root zone and the WDI based on ET p , which allows the direct calculation of soil water deficit from the WDI. (ii) Significant correlation of the WDI to plant stomatal conductance and leaf water potential indicates the possibility to detect early drought signals. The diurnal variation of WDI and ET a suggests that the ET a /ET p ratio varies during the day, being higher in the morning and lower at noon hours. Additionally, the largest difference between crops was seen later in the afternoon, as healthier plants started to recover faster compared to more stressed plants that demonstrated higher WDI values for a longer period. This can be used in further studies to improve the temporal upscaling of ET a with higher precision for irrigation planning. (iii) Since most of the parameters for the calculations of WDI, ET a and FTSW were derived from only a few initial variables routinely measured by the weather station, the WDI-based method allows for rapid calculation of the field energy balance and evaluation of soil and crop water status.
Future studies should assess WDI performance for different pedoclimatic and management (crops, irrigation) conditions. The WDI is based on the assumption that the relation between f veg and T a − T s from full-cover to bare soil is linear, whereas possible nonlinearity should also be studied.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author(s). The data are not publicly available due to the lack of standard, coherent and easily accessible platform for sharing data and metadata.

Appendix B. Energy Balance Components and Calculations
The following equations were used to calculate Equations (2)-(6) using the initial parameters of air temperature T a , solar radiation R s , wind speed u, air pressure p a and relative humidity RH collected by the nearby weather station and NDVI from multispectral images.
It is important to account for air density (ρ), especially with regard to the air humidity, when working in the humid climatic conditions due to the lower density of humid air compared to dry air: where ρ is the density of air (kg m −3 ), p d is the partial pressure of dry air (Pa), R d is the specific gas constant for dry air and equals 287.05 (J/kgK), p v is the partial pressure of water vapor (Pa), R v is the specific gas constant for dry air and equals 461.495 (J/kgK), and T a is air temperature (K).
Values of minimum and maximum canopy resistances r cm and r cx were obtained from the stomatal resistance (r s ) and leaf area index (LAI): where r sm and r sx are the minimum and maximum stomatal resistances (sm −1 ), respectively. For winter wheat, their values can be assumed to be r sm = 25 s m −1 and r sx = 1500 s m −1 [21]. Actual vapor pressure at hourly scale was calculated by the following equation: where e a is the actual vapor pressure (kPa), e 0 Ta(hr) is the hourly saturation vapor pressure at the temperature T a (kPa) and RH hr is the average hourly relative humidity (%). The hourly saturation vapor pressure is: where e 0 Ta(hr) is the saturation vapor pressure at the temperature T a (kPa).