Tree Water Status in Apple Orchards Measured by Means of Land Surface Temperature and Vegetation Index (LST–NDVI) Trapezoidal Space Derived from Landsat 8 Satellite Images

In this study, the split window (SW) method was applied for land surface temperature (LST) retrieval using Landsat 8 in two apple orchards (Glindow, Altlandsberg). Four images were acquired during high demand of irrigation water from July to August 2018. After pre-processing images, the normalized difference vegetation index (NDVI) and LST were calculated by red, NIR, and thermal bands. The results were validated by interpolated infrared thermometer (IRT) measurements using the inverse distance weighting (IDW) method. In the next step, the temperature vegetation index (TVDI) was calculated based on the trapezoidal NDVI/LST space to determine the water status of apple trees in the case studies. Results show good agreement between interpolated LST using IRT measurements and remotely sensed LST calculation using SW in all satellite overpasses, where the absolute mean error was between 0.08 to 4.00 K and root mean square error (RMSE) values ranged between 0.71 and 4.23 K. The TVDI spatial distribution indicated that the trees suffered from water stress on 7 and 23 July and 8 August 2018 in Glindow apple orchard with the mean value of 0.69, 0.57, and 0.73, whereas in the Altlandsberg orchard on 17 August, the irrigation system compensated the water deficit as indicated by the TVDI value of 0.34. Moreover, a negative correlation between TVDI and vegetation water content (VWC) with correlation coefficient (r) of −0.81 was observed. The corresponding r for LST and VWC was equal to −0.89, which shows the inverse relation between water status and temperature-based indices. The results indicate that the LST and/or TVDI calculation using the proposed methods can be effectively applied for monitoring tree water status and support irrigation management in orchards using Landsat 8 satellite images without requiring ground measurements.


Introduction
Tree water status and the determination of appropriate irrigation measure can be supported by digital plant data as suggested in the context of precision horticulture [1,2]. Established methods for estimating water stress such as gravimetrically or electrically measured soil and plant water contents are not only time consuming, but also produce discrete local data points, which may not be representative of an entire field or orchard. To maximize water use efficiency, new irrigation management strategies are needed [3][4][5][6]. Irrigation needs in the orchard are influenced by weather, climate variables, soil, and plant response, which can be measured with thermal imaging [3,[7][8][9][10][11].
Furthermore, irrigation scheduling considering regulated deficit irrigation, is based on crop response to water stress at different growth stages [12][13][14][15]. Increased leaf temperature and decreased vegetation (leaf) water content (VWC) are signs of water stress monitored with thermal imaging [16] that can be implemented in the orchard by crane, tractor, and unmanned aerial vehicle (UAV). However, one of the key parameter in vegetation water balance studies is land surface temperature (LST). In the past, the LST was measured by infrared thermometers which was costly and time-consuming. Since the 1980s, the remote sensing (RS), particularly satellite imagery, methods have been used for LST calculations [17][18][19]. From the satellite position, "surface" means whatever is seen from the atmosphere on the surface of the earth, e.g., in an apple orchard, the tree canopy is considered as a surface, ergo the LST equals predominantly to the leaf surface temperature. An early effort for calculating LST using RS methods was made using the advanced very high resolution radiometer (AVHRR) on the NOAA 7 satellite, which acquires 1-km spatial resolution data, and applying the split window algorithm [20]. After analyzing the relevant sources of errors, the author showed that the split window method could be adopted with an accuracy of 3 • C. Among numerous applications of LST retrieval calculations, its application is still followed in estimating water status in the field [21]. In this regard, many studies have been carried out for developing and improving the LST retrieval methods [22][23][24][25][26] to estimate evapotranspiration. Moreover, the RS has been utilized for estimating vegetation indices (VI), such as the normalized difference vegetation index (NDVI) [27]. The scatter plot of these remotely sensed indices (NDVI and LST) resulted in a trapezoidal space, which has a physical concept i.e., the spatial relation (trapezoidal space) between LST and NDVI is applied to monitor the water status of crops and orchards using temperature vegetation dryness index (TVDI) [28,29]. However, the TVDI reflects a mixed signal of the soil moisture under trees and leaf water content of trees.
Due to the relatively high resolutions and the free availability of all data by the U.S. Geological Survey (USGS), the Landsat 8 launched on 11 February 2013, with 9 operational land imager (OLI) and 2 thermal infrared sensor (TIRS) bands providing the prerequisites for calculating NDVI and LST. Using the TIRS bands provided, the calculation of the LST can be performed using either two thermal bands, e.g., split window (SW) algorithm, or one band such as single channel (SC) algorithm [30]. All of the indicated algorithms are based on the radiative transfer theory (RTT) equation. The Landsat 8 data for LST retrieval calculation was approached in northwest of China [31]. The authors applied SW algorithm and validated results with moderate resolution imaging spectroradiometer (MODIS) LST products.
Although many further studies applying satellite imagery for water status estimation based on land surface temperature (LST) were carried out in field crops, the application in orchards is still in its infancy. The present study was aimed at proposing a method for TVDI estimation based on partly mixed pixels in orchard conditions by means of (i) selecting an interpolation method for gaining spatial distribution data of LST based on infrared thermometer (IRT) measurements and approaching the SW method using Landsat 8 images. Subsequently, (ii) validate the results with ground reference measurements, and (iii) estimating TVDI considering the VWC.

Study Area
In this study, two commercial apple orchards located in a fruit production area of Germany were selected: The Glindow's commercial apple orchard with longitude and latitude of the field center being 52 • 21 15.8" N 12 • 51 53.8" E (path/row: 193/23 for Landsat 8 satellite images). The cultivar of the apple (Malus × domestica Borkh.) is 'Pinova', which was harvested 24 September 2018 in the present orchard. The Altlandberg's commercial apple orchard located with longitude and latitude of the field center is 52 • 36 27.5" N 13 • 49 3.07" E (path/row: 192/23 for Landsat 8 satellite images) with the cultivar (Malus × domestica Borkh.) being 'Gala', was planted in 3.5 × 1.0 m distance. The common irrigation Sustainability 2020, 12, 70 3 of 19 system in this area is drip irrigation. Both sites have semi-humid climate and are geographically surrounded by fruit orchards (Figure 1). Sustainability 2020, 11, x FOR PEER REVIEW 3 of 22 52°36'27.5"N 13°49'3.07"E (path/row: 192/23 for Landsat 8 satellite images) with the cultivar (Malus × domestica Borkh.) being 'Gala', was planted in 3.5 × 1.0 m distance. The common irrigation system in this area is drip irrigation. Both sites have semi-humid climate and are geographically surrounded by fruit orchards (Figure 1). Daily meteorological data were recorded with weather stations (IMT280, Pessl, Austria) including air temperature, soil temperature in 5 depth, and relative humidity (RH).

Field Measurements
The LST and VWC were measured in 10 points during the Landsat 8 satellite overpasses on 7 and 23 July and 8 August 2018 between 11:00 to 12:00 in Glindow and 17 August 2018 in Altlandsberg. At each point, 20 representative leaves were sampled (n = 200) and canopy temperature was measured 20 times for each canopy using digital IRT (OS91, Omega, Norwalk, CT, USA) from sun and shade position of trees. The mean of these 20 values is considered as a representative of leaf temperature at each point.
The leaf samples were sealed in individual plastic bags, placed in a cool and dark container to avoid water losses on the return to the laboratory. Fresh weight (FW) of gathered samples was recorded immediately, and subsequently samples were oven-dried at 105 ℃ for 24 h until constant dry weight (DW) was reached. Finally, the VWC value was calculated for each sample (Equation (1)):

Landsat 8 Satellite Image Acquisition and Data Preprocessing
Landsat 8 has 9 operational land imager (OLI) spectral bands in the range 0.45 μm (blue/visible) to 2.29 μm (infrared/nonvisible) and 2 thermal infrared sensor (TIRS) bands in the high infrared range (10.0-12.5 μm). Due to the relatively high horizontal resolution of the OLI bands (30 m), the frequent sweep of the satellite over the same corridor at the earth's surface appeared in intervals of 16 days, which were free availability by the USGS (http://earthexplorer.usgs.gov). Four Landsat images were acquired during the high water demand months (July and August) of apple trees (Table 1). Daily meteorological data were recorded with weather stations (IMT280, Pessl, Austria) including air temperature, soil temperature in 5 depth, and relative humidity (RH).

Field Measurements
The LST and VWC were measured in 10 points during the Landsat 8 satellite overpasses on 7 and 23 July and 8 August 2018 between 11:00 to 12:00 in Glindow and 17 August 2018 in Altlandsberg. At each point, 20 representative leaves were sampled (n = 200) and canopy temperature was measured 20 times for each canopy using digital IRT (OS91, Omega, Norwalk, CT, USA) from sun and shade position of trees. The mean of these 20 values is considered as a representative of leaf temperature at each point.
The leaf samples were sealed in individual plastic bags, placed in a cool and dark container to avoid water losses on the return to the laboratory. Fresh weight (FW) of gathered samples was recorded immediately, and subsequently samples were oven-dried at 105 • C for 24 h until constant dry weight (DW) was reached. Finally, the VWC value was calculated for each sample (Equation (1)): (1)

Landsat 8 Satellite Image Acquisition and Data Preprocessing
Landsat 8 has 9 operational land imager (OLI) spectral bands in the range 0.45 µm (blue/visible) to 2.29 µm (infrared/nonvisible) and 2 thermal infrared sensor (TIRS) bands in the high infrared range (10.0-12.5 µm). Due to the relatively high horizontal resolution of the OLI bands (30 m), the frequent sweep of the satellite over the same corridor at the earth's surface appeared in intervals of 16 days, which were free availability by the USGS (http://earthexplorer.usgs.gov). Four Landsat images were acquired during the high water demand months (July and August) of apple trees (Table 1).
In order to preprocess the satellite images, the digital numbers (DN) of TIRS and OLI bands were converted to spectral radiance and top of atmosphere (TOA) planetary reflectance. To convert DNs to radiance the formula below was used [32]: where L λ is spectral radiance (W/(m 2 sr µm)), M L is radiance multiplicative scaling factor, Q cal is DN value of each pixel, and A L is radiance additive scaling factor. TIRS data were converted from spectral radiance to brightness temperature [32]: where BT i is top of atmosphere (TOA) brightness temperature for TIRS band i (= 10, 11) in Kelvin. The K coefficients are thermal conversion constants for Landsat 8 TIRS bands (Table 2).  The NDVI is computed by TOA planetary spectral reflectance, with solar angle correction. The DNs are converted to TOA planetary reflectance by the below equation [32]: where ρ λ andρ λ are TOA planetary spectral reflectance with and without solar angle correction (unitless), M p is reflectance multiplicative scaling factor, A p is reflectance additive scaling factor, and θ is solar elevation angle. M L , A L , M p , A p , and θ parameters in Equations (2) and (4) are extracted from the metadata file, which is downloaded with satellite images from the USGS EarthExplorer.

NDVI Calculation Using Landsat 8
After conversion of the DN to reflectance, the normalized difference vegetation index (NDVI) for the red and near infrared (NIR) bands was calculated by the following equation [33]: The value range of NDVI is −1 to 1, negative values of NDVI may point to water bodies and snow cover, while values greater than 0.5 indicate dense vegetation cover [34].

Radiative Transfer Theory (RTT) Equation and Split-Window (SW) Algorithm
A simplified radiative transfer equation expresses the apparent radiance received by a sensor [35,36]: where B i (BT i ) radiance received by channel i(=10, 11) of the sensor with brightness temperature BT i ; B i (T s ) is the ground radiance; τ i (θ) is atmospheric transmittance for channel i when viewing zenith angle is θ. The Landsat-8 TIRS is at an altitude of 705 km with a swath of 185 km. Therefore, TIRS is treated as nadir viewing, since the maximum zenith viewing angle is not more than 7.5 • , so the θ could be ignored i.e., on that angle, the effect on the atmospheric transmittance in TIRS band is negligible [37][38][39]. ε i is land surface emissivity (LSE); I ↓ i and I ↑ i are down-welling and up-welling path radiance, respectively. Regarding the Plank's law, B i (T) is calculated by: where c is light speed (2.9979 × 10 8 m/s), h is Plank constant (6.6261 × 10 −34 J·s); k is the Boltzmann constant (1.3806 × 10 −23 J/K); and λ i is the effective center wavelength of band i which is equal to 10.9 and 12 µm for band 10 and 11 TIRS; C 1 is 1.19104 × 10 8 Wµm 4 m −2 sr −1 and C 2 is 14387.7 µmK. By solving Equations (6) and (7), the LST could be estimated in Kelvin (K) for Landsat thermal bands (Equation (8)): In order to retrieve LST from Landsat 8 TIRS bands by RTT equation, τ i , I ↓ i , and I ↑ i could be analyzed by simultaneous measurements of radio sounding data during the time period of Landsat overpass. In automated applications, these types of data are not available, thus [40] simplified the RTT equation and proposed an algorithm for spilt window method by solving Equation (6) for the advanced very high resolution radiometer (AVHRR) data. Their proposed SW algorithm ( Figure 2) was applied in the present study (Equation (9)): where A 0 , A 1 and A 2 are calculated by: C and D are middle term variables that could be estimated as following: a and b coefficients are determined by linear regression between intermediate parameter (Li) and BT i .
Li is the derivative of the Plank's law function for band i(=10, 11) at brightness temperature [31]: As shown in the above equations, the A coefficients depend on LSE (ε i ) and atmospheric transmittance (τ i ).
The LSE is the ratio of radiated thermal energy by surface and a blackbody in the same temperature [41]. It can be estimated by NDVI-based emissivity method (NBEM), which was Sustainability 2020, 12, 70 6 of 19 already applied for various analyses of satellite images. NBEM uses empirical NDVI thresholds to estimate LSE from Landsat 8 [36,42]: where m and n are constant coefficients. For Landsat TIRS band 10, m, and n are 0.973 and 0.047, respectively; NDVI s = 0.2 and NDVI v = 0.5 are the thresholds of soil and full vegetation pixels [43]; ε s = 0.966 and ε v = 0.973 are soil and vegetation emissivity values [44]; P v is proportion of vegetation in mixed pixels (Equation (17)): C i is cavity effect due to the surface roughness. It is considered zero for flat surface, otherwise the value needs to be estimated [22]: where F is geometrical factor which spans the range 0 to 1, frequently set as 0.55 [45].
Since surface spectral reflections are absorbed by the atmosphere and the predominant reflector absorption is determined by the water vapor content (w) which is estimated by: where T 0 is the temperature near the surface of the earth. [46] showed that the atmospheric transmittance (τ i ) mainly depends on w for mid-latitude summer atmospheric profile. The relation between w and τ i is shown in Table 3.

Validation of LST Retrieval
For calculating LST several parameters, namely LSE, NDVI, brightness temperature (BT) and atmospheric transmittance (τi) need to be estimated. The validation was carried out with ground measurements during the time period of satellite overpass. In this regard, leaf surface temperature was measured in 10 points during the overpass to validate the LST product of Landsat 8 which is a continuous raster file, whereas the measured points are discrete value. For selecting the best interpolation method for  Table 3. Relation between atmospheric transmittance factor (τ i , i = 10, 11) and water vapor content (w) for mid-latitude summer atmospheric profile.

Validation of LST Retrieval
For calculating LST several parameters, namely LSE, NDVI, brightness temperature (BT) and atmospheric transmittance (τ i ) need to be estimated. The validation was carried out with ground measurements during the time period of satellite overpass. In this regard, leaf surface temperature was measured in 10 points during the overpass to validate the LST product of Landsat 8 which is a continuous raster file, whereas the measured points are discrete value. For selecting the best interpolation method for interpolating the raster surface, nearest neighbor (nn), inverse distance weighting (IDW), spline, and kriging models were tested. Ten observed data were divided into a calibrating set of 7 points and validation data set of 3 points. In this regard, interpolation methods were evaluated by the coefficient of determination (R 2 ) and the root mean square error (RMSE) Equations (20) and (21): where y i , y, andŷ i are observed, average of observed and interpolated temperature at point i, respectively. Moreover, the algorithm is developed for different locations (2 orchards) and time (4 individual dates).

LST-NDVI Space and TVDI Calculation
The scatter plot of LST and NDVI results in a trapezoidal space [29,47,48]. As Figure 3 shows, each part of this space shows information about canopy water status and soil moisture i.e., the dry and wet edges can be applied to obtain information about water status of the orchard. In the dry edge, a strong negative correlation can be seen between NDVI and LST which shows the water stress condition has been taken place, whereas in the wet edge, the tree is not under stress condition. The TVDI [0;1] is defined by: where a 1 and b 1 are the intercept and slope of dry edge line (Figure 3), i is the pixel number. As Figure 3 shows, for TVDI calculation, the NDVI spans from 0 to 1, whereas the range in an apple orchard is limited between 0.2 and 0.7 ( Figure 4). Moreover, each point in the trapezoidal space is related to one pixel capturing an area of 900 m 2 , thus an orchard with area less than 10 ha is not big enough to make a trapezoidal space. In this regard, the calculations were done for the whole area of a Landsat 8 satellite image and subsequently, the orchard was extracted.
To calculate the dry and wet edges, the NDVI values were divided into 10 equidistant classes. For each class, the mean value of NDVI is considered as an independent variable of the linear regression, whereas the maximum and minimum values of LST for each class are considered as depended variable for dry and wet regression lines, respectively. The scatter plot of LST and NDVI results in a trapezoidal space [29,47,48]. As Figure 3 shows, each part of this space shows information about canopy water status and soil moisture i.e., the dry and wet edges can be applied to obtain information about water status of the orchard. In the dry edge, a strong negative correlation can be seen between NDVI and LST which shows the water stress condition has been taken place, whereas in the wet edge, the tree is not under stress condition. The TVDI [0;1] is defined by: where a1 and b1 are the intercept and slope of dry edge line (Figure 3), i is the pixel number. As Figure 3 shows, for TVDI calculation, the NDVI spans from 0 to 1, whereas the range in an apple orchard is limited between 0.2 and 0.7 ( Figure 4). Moreover, each point in the trapezoidal space is related to one pixel capturing an area of 900 m 2 , thus an orchard with area less than 10 ha is not big enough to make a trapezoidal space. In this regard, the calculations were done for the whole area of a Landsat 8 satellite image and subsequently, the orchard was extracted.
To calculate the dry and wet edges, the NDVI values were divided into 10 equidistant classes. For each class, the mean value of NDVI is considered as an independent variable of the linear regression, whereas the maximum and minimum values of LST for each class are considered as depended variable for dry and wet regression lines, respectively.

Interpolation Methods to Generate LST Spatial Distribution Using IRT Ground Measurements
In order to create spatially continuous surfaces data for the LST layer using IRTs measurements, spatial interpolations were conducted. The interpolation exercise according to (a) ordinary kriging, (b) spline (Theissen and regularized method), (c) inverse distance weighting (IDW), and (d) nearest neighbor (nn) pointed out the highest accuracy of the IDW method compared to remaining interpolation methods (Table 4) with RMSE 0.21 K and 0.35 K with R 2 of 0.98 and 0.91 for calibration and cross validation on 7 July 2018 Landsat overpass, respectively. The methods were tested for Glindow and the best method-the IDW-was applied to the Altlandsberg data, as well ( Figure 5).

Split Window (SW) Algorithm to Calculate the LST Spatial Distribution Using Landsat 8 Satellite Images
In order to calculate the satellite-based LST using split window method, the atmospheric transmittance (τ i ) for bands 10 and 11of Landsat 8 was calculated (Table 5). It was found that the water vapor content (w) appeared in the range of 0.2-3.0 (g/cm 2 ) in all overpasses. As expected by increasing w, the τ i values are decreased which shows the impact of water particles on the absorption of infrared radiances measured at the thermal bands of Landsat 8 images. Table 5. Atmospheric transmittance factor (τ i ) for bands 10, 11 of Landsat 8 using water vapor content, w (g/cm 2 ) data.

Date
Location w (g/cm 2 ) τ 10 τ 11 The NDVI was used for estimation of proportion of vegetation (P v ) and LSE. The mean value of P v for all overpasses was less than 1 (P v 0.5) which marks that the orchards were not fully covered by trees canopies, thus, the LSE data should be between a percentage of LSE for bare soil and a complementary percentage of tree canopy. The value of NDVI on 7 and 23 July in Glindow is unchanged but a reduction can be seen on 8 August, whereas, no huge difference appeared between vegetation cover in these 3 dates (Figure 4). Consequently, it was expected that the NDVI values would be enhanced or similar to the values on 23 July because the leaf area of tree remains unchanged during this short time interval (23 July to 8 August), while the leaf surface and angle is assumed to change due to the water deficit. Consequently, the high temperature (Tmax = 49.5 • C) and low RH caused reduced NDVI values.

07-July
The LSE is one of indispensable inputs for LST estimation models. In the present study, the NBEM method was applied for LSE estimation, which is based on NDVI thresholds for bare soil (NDVI = 0.2) and dense vegetation cover (NDVI = 0.5). The mean value for surface emissivity was 0.98 in all overpasses (Table 6), which means for equal LST, the radiance emitted from orchard surface is 0.02 less than black body surface which is similar to the results of [49].  (15) can be estimated by linear regression between Li and BT i . As Figure 6 and Table 7 show, the a i coefficients are between −64.37 and −61.74 for bands 10 and 11, respectively whereas, the b i coefficients are about 0.45 for both thermal bands. The data were perfectly fitted to the regression line ( Figure 6) in all overpasses (R 2 = 1).  By calculating all required parameters and coefficients, the spatial distribution of LST was estimated by SW method. As Figure 7 shows, the lower temperatures are assigned to the pixels with enhanced NDVI.  By calculating all required parameters and coefficients, the spatial distribution of LST was estimated by SW method. As Figure 7 shows, the lower temperatures are assigned to the pixels with enhanced NDVI. To validate SW Landsat 8 LST retrieval, IRT observations were applied. As Figure 8 and Table 8 show, the absolute mean error of the SW method was below 4 K in all overpasses to the extent that the maximum RMSE was found on 7 July with the value of 4.23 K and absolute mean error of 3.94 K. On the other hand, the RMSE value on the 17 August overpass in Altlandsberg was only 0.71 K with absolute mean error of 0.12 K.
the absolute mean error of the SW method was below 4 K in all overpasses to the extent that the maximum RMSE was found on 7 July with the value of 4.23 K and absolute mean error of 3.94 K. On the other hand, the RMSE value on the 17 August overpass in Altlandsberg was only 0.71 K with absolute mean error of 0.12 K.
The reason for errors can be related to the very dynamic nature of LST e.g., the variation of 10 K has been reported in few centimeters or 1 K change were taken place in 1 minute [50]. As indicated in pervious parts, the ground measurements were done during the satellite transit time over the orchards to minimize these variations. All above shows that estimating the LST with 0.08-4 K error using open access satellite images is accurate enough for water status studies. It is important to note that the IDW interpolation method as well as the SW model were developed for Glindow orchard on 7 July 2018, and then applied for different dates in Glindow and Altlandsberg, as well. i.e., the application of models was validated temporally and spatially. The results proved that both models (IDW and SW) can be used for different time and location.  The reason for errors can be related to the very dynamic nature of LST e.g., the variation of 10 K has been reported in few centimeters or 1 K change were taken place in 1 min [50]. As indicated in pervious parts, the ground measurements were done during the satellite transit time over the orchards to minimize these variations. All above shows that estimating the LST with 0.08-4 K error using open access satellite images is accurate enough for water status studies. It is important to note that the IDW interpolation method as well as the SW model were developed for Glindow orchard on 7 July 2018, and then applied for different dates in Glindow and Altlandsberg, as well. i.e., the application of models was validated temporally and spatially. The results proved that both models (IDW and SW) can be used for different time and location.

LST-NDVI Trapezoidal Space for TVDI Calculation
In the next step, the equations of dry and wet edges of trapezoidal space between LST and NDVI were determined by linear regression model. The maximum and minimum values of LST for NDVI value can be used for dry and wet lines' equation using linear regression (Figure 9). By using these equations, the TVDI can be estimated. The effects of land surface emissivity (LSE) and proportion of vegetation (P v ) on LST calculation can be clearly seen in Figure 9 where data went up at NDVI = 0.2 (NDVI value for bare soil).

LST-NDVI Trapezoidal Space for TVDI Calculation
In the next step, the equations of dry and wet edges of trapezoidal space between LST and NDVI were determined by linear regression model. The maximum and minimum values of LST for NDVI value can be used for dry and wet lines' equation using linear regression (Figure 9). By using these equations, the TVDI can be estimated. The effects of land surface emissivity (LSE) and proportion of vegetation ( ) on LST calculation can be clearly seen in Figure 9 where data went up at NDVI = 0.2 (NDVI value for bare soil). As indicated in the material and methods part, in the TVDI calculations, two visible bands (NIR and red) and TIRS bands of Landsat 8 are applied. In this regard, the dry and wet edges of LST/NDVI space As indicated in the material and methods part, in the TVDI calculations, two visible bands (NIR and red) and TIRS bands of Landsat 8 are applied. In this regard, the dry and wet edges of LST/NDVI space were defined by using Equation (22), the spatial distribution of TVDI were generated ( Figure 10). The mean values of TVDI on 7 and 23 July in Glindow are 0.69 and 0.57 which show the trees are under moderate water stress and the irrigation system should start to work soon but on 8 August-as expected-the trees are under severe water stress where the mean TVDI is 0.73; therefore, the irrigation system was not working at the appropriate time and needs to start irrigation immediately. On the other hand, the drip irrigation in the Altlandsberg orchard already started to work during the measurements as well the satellite overpass, and it can be seen the TVDI is less than 0.5 (averagely 0.34), so the tree water status is in normal condition without water stress.
values of TVDI on 7 and 23 July in Glindow are 0.69 and 0.57 which show the trees are under moderate water stress and the irrigation system should start to work soon but on 8 August-as expected-the trees are under severe water stress where the mean TVDI is 0.73; therefore, the irrigation system was not working at the appropriate time and needs to start irrigation immediately. On the other hand, the drip irrigation in the Altlandsberg orchard already started to work during the measurements as well the satellite overpass, and it can be seen the TVDI is less than 0.5 (averagely 0.34), so the tree water status is in normal condition without water stress. As spatial distribution of TVDI pointed to low values in the center of the orchard and maximum values near the sides of the field on 7, 23 July 2018 in Glindow and 17 August 2018 in Altlandsberg. Here the effects of non-vegetated areas like roads and bare soils in side pixels became visible. Although the minimum values of TVDI on 8 August took place in the central pixel as well, the spatial distribution did not follow a trend like others which can be related to the high weather temperature and low relative humidity that cause an intensive dryness in the whole area without considering the land cover. The minimum value of TVDI is already equal to 0.7 which shows severe dryness.
The results illustrated that, the LST and TVDI can be applied as indicators for vegetation water content (VWC), thus the correlation between latter parameters was investigated. In order to find the relation between VWC and TVDI (and/or LST), and more importantly to validate the satellite based index, a linear regression model was developed between two parameters. As Figure 11 shows, the empirical relationship between values of measured VWC and TVDI/LST has a relatively high negative correlation coefficient. Here the effects of non-vegetated areas like roads and bare soils in side pixels became visible. Although the minimum values of TVDI on 8 August took place in the central pixel as well, the spatial distribution did not follow a trend like others which can be related to the high weather temperature and low relative humidity that cause an intensive dryness in the whole area without considering the land cover. The minimum value of TVDI is already equal to 0.7 which shows severe dryness.
The results illustrated that, the LST and TVDI can be applied as indicators for vegetation water content (VWC), thus the correlation between latter parameters was investigated. In order to find the relation between VWC and TVDI (and/or LST), and more importantly to validate the satellite based index, a linear regression model was developed between two parameters. As Figure 11 shows, the empirical relationship between values of measured VWC and TVDI/LST has a relatively high negative correlation coefficient. Relationships between VWC and temperature-based indices for different trees at different locations would enable irrigation planners to estimate VWC. The accuracy of the SW satellite-based algorithm was evaluated in studies published in the literature (e.g. [22]) which is one of main goals of the present study as well; however, it is notable to say in the case of this study, the data evaluated by ground measured data in the same time of the Landsat 8 overpass. Moreover, many studies applying satellite imagery for water status estimation based on land surface temperature (LST) were carried out in relatively large field crops (e.g. [47]), whereas, the application in apple orchards were investigated in this study. It is important to note the code of the TVDI calculation using Landsat 8 images is available on GitHub [51], the strong point of the provided algorithm is by introducing the satellite images as input to the code, several outputs including NDVI, LST, and TVDI spatial distribution maps can be obtained by the user.

Conclusions
In this study, progress has been made in combining the advantage of available open access Landsat 8 satellite images with scientific advances to estimate tree water status in apple orchards. In this regard, the TDVI was derived from land surface temperature and NDVI outlining the spatial distribution of tree water status illustrating that the satellite-based data can provide an alternative to the manual temperature readings and vegetation water content with values of absolute mean error ranging from 0.08 to 3.94 K and RMSE of 0.71 to 4.23 K between satellite-based estimated LST and measured canopy temperature for all overpasses. Moreover, the relationship between vegetation water content and TVDI (and LST) with r = −0.81 (and −0.89). It is concluded that the satellite-based methods for determining the LST, and subsequently TVDI provides advantages over ground measurements methods, because the ground observed LST are point data whereas the spatial and temporal LST changes are very dynamic. Therefore, it can be assumed that calculating the spatial distribution of LST pixel by pixel is more reliable than using interpolation methods from limited discrete point measurement. In addition, the vegetation water content can be revealed by estimated LST and TVDI, derived from thermal images for water stress detection. The remote sensing method outperforms time-consuming sampling of a large number of IRT measurements in the field. Consequently, the TVDI index using thermal bands of Landsat 8 has the potential to replace direct vegetation water content measurements and to provide a more robust measure of the tree water status and its spatial patterns of the apple production. Moreover, it can be an efficient parameter for land use change detection in mesoscale studies as well as irrigation water management planning projects.  The accuracy of the SW satellite-based algorithm was evaluated in studies published in the literature (e.g., [22]) which is one of main goals of the present study as well; however, it is notable to say in the case of this study, the data evaluated by ground measured data in the same time of the Landsat 8 overpass. Moreover, many studies applying satellite imagery for water status estimation based on land surface temperature (LST) were carried out in relatively large field crops (e.g., [47]), whereas, the application in apple orchards were investigated in this study. It is important to note the code of the TVDI calculation using Landsat 8 images is available on GitHub [51], the strong point of the provided algorithm is by introducing the satellite images as input to the code, several outputs including NDVI, LST, and TVDI spatial distribution maps can be obtained by the user.

Conclusions
In this study, progress has been made in combining the advantage of available open access Landsat 8 satellite images with scientific advances to estimate tree water status in apple orchards. In this regard, the TDVI was derived from land surface temperature and NDVI outlining the spatial distribution of tree water status illustrating that the satellite-based data can provide an alternative to the manual temperature readings and vegetation water content with values of absolute mean error ranging from 0.08 to 3.94 K and RMSE of 0.71 to 4.23 K between satellite-based estimated LST and measured canopy temperature for all overpasses. Moreover, the relationship between vegetation water content and TVDI (and LST) with r = −0.81 (and −0.89). It is concluded that the satellite-based methods for determining the LST, and subsequently TVDI provides advantages over ground measurements methods, because the ground observed LST are point data whereas the spatial and temporal LST changes are very dynamic. Therefore, it can be assumed that calculating the spatial distribution of LST pixel by pixel is more reliable than using interpolation methods from limited discrete point measurement. In addition, the vegetation water content can be revealed by estimated LST and TVDI, derived from thermal images for water stress detection. The remote sensing method outperforms time-consuming sampling of a large number of IRT measurements in the field. Consequently, the TVDI index using thermal bands of Landsat 8 has the potential to replace direct vegetation water content measurements and to provide a more robust measure of the tree water status and its spatial patterns of the apple production. Moreover, it can be an efficient parameter for land use change detection in mesoscale studies as well as irrigation water management planning projects.