2.5.4. *wc*/*h<sup>c</sup>* Ratio

*wc*/*h<sup>c</sup>* was calculated by simply dividing canopy width by canopy height at each contextual spatial domain.

#### *2.6. Goodness-of-Fit Statistics*

Evaluating the performance of the *TSEB2T* model with the *sUAS* imagery for the four different modeling grid resolutions involved comparing the estimated fluxes with measurements from the *EC* towers. Computed statistical metrics included the root mean square error (*RMSE*), the normalized root mean square error (*NRMSE*), mean absolute error (*MAE*), mean absolute percentage error (*MAPE*), and Nash–Sutcliffe efficiency coefficient (*NSE*). A value of *NSE* = 1 indicates perfect agreement between modeled and observed flux, while *NSE* approaching 0 means that the agreement is very poor, and *NSE* < 0 indicates unacceptable performance [68]. These statistical measurements are calculated as follows using *LE* as the flux:

$$RMSE = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} (LE\_{m,i} - LE\_{o,i})^2} \tag{9}$$

$$\text{NRMSE} = \frac{\text{RMSE}}{\sigma\_o} \tag{10}$$

$$MAE = \frac{\sum\_{i=1}^{n} \left| LE\_{m,i} - LE\_{o,i} \right|}{n} \tag{11}$$

$$MAPE = \frac{\sum\_{i=1}^{N} \left| \frac{LE\_{m,i} - LE\_{o,i}}{LE\_{o,i}} \right| \ast 100}{n} \tag{12}$$

$$NSE = 1 - \frac{\sum\_{i=1}^{n} \left( LE\_{m,i} - LE\_{o,i} \right)^2}{\sum\_{i=1}^{n} \left( LE\_{m,i} - \overline{LE}\_{o,i} \right)^2} \tag{13}$$

where *LEm*. denotes the modeled latent heat flux obtained from the *TSEB2T* aggregated up for the estimated flux footprint/source area, *LE<sup>o</sup>* denotes the observed values from the *EC* tower, and *n* represents the number of observations, σ*<sup>o</sup>* denotes the standard deviation of observed values.

*LE* was used for evaluating the impact of spatial resolution or grid size on modeled fluxes. At field scale, the evaluation is done using the spatial mean and coefficient of variation (*CV*) statistics. For *LE* statistical characteristics, frequency and cumulative distribution curves were used. Finally, to evaluate the effect of aggregating *LE* at 3.6 m, 7.2 m, 14.4 m, and Landsat scale, relative difference (relative error) was used. Relative difference (relative error) is defined as the root mean square error (*RMSE*) between the aggregated resolution and its reference grid size resolution of 3.6 × 3.6 m divided by the spatial mean (µ) value computed from the reference grid size (3.6 m × 3.6 m), i.e., *E<sup>r</sup>* = *RMSE*/µ [14]. Each grid value of aggregated data was compared to the *n* × *n* set of reference scale or resolution (3.6 m) grid using *E<sup>r</sup>* .

#### **3. Results and Discussion**

#### *3.1. TSEB2T Contextual Spatial Domains Validation*

#### 3.1.1. *EC* Footprint Estimation

The results of footprint analysis using the *2D* flux model developed by Kljun et al. [61] and described in Section 2.1 are shown in Figure 5 for the different *sUAS* flights.

The orientation and size of each flux footprint/source area depends on the micro-meteorological conditions at the site measured by the *EC* towers, which include the turbulence fluxes, friction velocity (*u*∗), and wind speed, which affect atmospheric stability, and canopy and *EC* measurement height, which affect the effective sampling height and wind direction that affects the orientation of the footprint. The total statistical weight of the footprint is taken to equal unity, although the actual area computed by the footprint model represents 90% of the contribution since the additional 10% essentially makes no measurable contribution. To compare the fluxes computed by the *TSEB2T* model at the different spatial resolutions with the *EC* measurements, the source area estimated by the footprint model was multiplied by the corresponding modeled fluxes (*Rn*, *H*, *LE*, *G*) using *ArcGIS10.6*. Then, a comparison between the weighted fluxes at the different spatial resolutions or grid sizes from the *TSEB2T* version of *TSEB* and *EC* measurements was performed to assess model performance.

obtained for *Rn* across different spatial domains/model grids.

the flux measurements at 3.6-m, 7.2-m, and 14.4-m modeling grid sizes, while at the 30-m resolution, the *MAE* value was around 85 W/m2. As demonstrated in Figure 6d, all values of *LE* are underestimated (below 1:1 line) with an NSE coefficient of 0.2. Furthermore, the highest *NRMSE* values were observed

**Figure 5.** Layout of 90% *EC* footprints for two towers at different times considered by this study.

#### 3.1.2. Statistical Performance

Table 3 lists the goodness-of-fit statistics between the energy fluxes using *TSEB2T* at different spatial resolutions and *EC* tower observations, while Figure 6 shows the relationship between the modeled and measured fluxes. The results indicate a significant deterioration in model performance at the 30-m grid size. A major factor that may be responsible for this poor performance in the *TSEB2T* model at 30-m resolution is that the size and dimension of the *EC* source area estimated by the footprint model cannot incorporate a representative range in the spatial variability in the fluxes at 30-m resolution. This conclusion agrees with a previous study conducted by Song et al. [69] that showed a major problem in comparing modeled and measured fluxes when there is a mismatch in pixel resolution or model grid size in the remotely sensed *ET* output and in the source area contributing to the *EC* tower measurements in a heterogeneous landscape.

Results in Table 3 also indicate that *R<sup>n</sup>* and *G* across multiple aggregated grids demonstrated a close agreement between the *TSEB2T* output and observed measurements, as indicated by lower *MAE* and *MAPE* with quite constant correlation (*R 2* ). The *MAE* and *MAPE* in the *R<sup>n</sup>* estimate at grid sizes of 3.6 m, 7.2 m, and 14.4 m accounted for less than 25 W/m<sup>2</sup> and 5%, respectively. However, at Landsat scale the *MAE* increased slightly to 29 W/m<sup>2</sup> . A similar result was obtained for *H*, where *MAE* at the finer resolutions yielded values less than 45 W/m<sup>2</sup> , while the coarser grid size of 30 m yielded a larger *MAE* of nearly 80 W/m<sup>2</sup> . As shown in Table 3, the correlation of *H* is higher than *G* and *LE*, except for 30-m resolution/model grid. This implies that the performance of the 30-m resolution is different compared to the 3.6-m, 7.2-m, and 14.4-m resolutions. The results for *LE* indicated good agreement with the flux measurements at 3.6-m, 7.2-m, and 14.4-m modeling grid sizes, while at the 30-m resolution, the *MAE* value was around 85 W/m<sup>2</sup> . As demonstrated in Figure 6d, all values of *LE* are underestimated (below 1:1 line) with an *NSE* coefficient of 0.2. Furthermore, the highest *NRMSE* values were observed for *LE*, compared with other surface fluxes, particularly at 30-m resolution. The lowest *NRMSE* was obtained for *R<sup>n</sup>* across different spatial domains/model grids.


**Table 3.** Goodness-of-fit statistics between the eddy covariance and the *TSEB2T* fluxes at different spatial scales (3.6 m, 7.2 m, 14 m, and 30 m).

**3.6 m** 

**7.2 m** 

**14.4 m** 

**Figure 5.** Layout of 90% *EC* footprints for two towers at different times considered by this study.

**Table 3.** Goodness-of-fit statistics between the eddy covariance and the TSEB2T fluxes at different

**Spatial Domain Fluxes RMSE (W/m2) NRMSE MAE (W/m2) MAPE (%) NSE R2**

*Rn* 28 0.3 25 5 0.9 0.94 *LE* 69 1.2 58 20 0.5 0.49 *H* 54 0.8 41 26 0.7 0.67 *G* 34 0.9 30 51 0.6 0.56

*Rn* 27 0.3 24 4 0.9 0.94 *LE* 66 1.2 56 19 0.5 0.53 *H* 51 0.7 36 24 0.7 0.67 *G* 33 0.8 30 50 0.6 0.58

*Rn* 25 0.3 20 4 0.9 0.95 *LE* 79 1.4 56 18 0.1 0.21 *H* 48 0.7 35 26 0.6 0.69 *G* 32 0.8 29 49 0.6 0.59

*Rn* 34 0.4 29 5 0.9 0.96 *LE* 101 1.8 86 30 0.2 0.53 *H* 93 1.3 78 67 −0.1 0.23 *G* 31 0.8 28 48 0.6 0.60

spatial scales (3.6 m, 7.2 m, 14 m, and 30 m).

**Figure 6.** Scatterplot of observed versus estimated surface fluxes using different model grid sizes/resolution with the TSEB2T model; (**a**) 3.6 m, (**b**) 7.2 m, (**c**) 14.4 m, and (**d**) 30 m. **Figure 6.** Scatterplot of observed versus estimated surface fluxes using different model grid sizes/resolution with the *TSEB2T* model; (**a**) 3.6 m, (**b**) 7.2 m, (**c**) 14.4 m, and (**d**) 30 m.

With the *TSEB2T* model and other remote sensing-based models using thermal-IR as the boundary condition, *LE* is solved as the residual component of the surface energy balance, *LE* = *R<sup>n</sup>* − *H* − *G*. Therefore, an error in the calculation of energy fluxes (*Rn*, *H*, and *G*) adversely affects the estimation of *LE*. Based on Figure 6, the *LE* estimation (or bias) is mainly influenced by the estimation of *H*. This conclusion was also reached by Kustas et al. [70], who showed the discrepancies between modeled and measured *LE* is due in large part, up to approximately 90%, to errors in modeled *H*.

### *3.2. Contextual Spatial Domain Aggregations E*ff*ects*

#### 3.2.1. The Effect of Model Grid Size on *TSEB2T* Inputs

### (a) Canopy and Soil Temperatures (*Tc*, *Ts*)

*T<sup>c</sup>* and *T<sup>s</sup>* were estimated based on a linear *LST-NDVI* relationship as described by Nieto et al. [2]. However, this relationship does not fulfill the homoscedasticity criterion when the spatial domain/resolution reaches a certain size (i.e., 30-m) as shown in Figure 7. For example, in the case of a 30-m grid size, a higher variability is observed in the *LST-NDVI* data compared with finer resolutions (3.6 m, 7.2 m, and 14.4 m). At micro-scale (e.g., 3.6 m), there are small number of pixels inside the spatial domain compared with others (7.2 m, 14.4 m, and 30 m), and exhibit an apparent linear relationship between *LST* and *NDVI*. However, at coarse resolution (e.g., 30 m), there are many more pixels, more rows of vineyard are included, and large vegetated and bare soil pixels exist inside the spatial domain. The result is a partially filled triangular shape. This indicates the relationship between *LST* and *NDVI* starts to resemble the "triangle method" [71] to estimate *ET* as the sampling domain increases.

domain increases.

*3.2. Contextual Spatial Domain Aggregations Effects* 

(a) Canopy and Soil Temperatures (*Tc*, *Ts*)

3.2.1. The Effect of Model Grid Size on *TSEB2T* Inputs

(16.85 °C) and 320 K (46.85 °C) for the *sUAS* flight in 2014.

With the *TSEB2T* model and other remote sensing-based models using thermal-IR as the boundary condition, *LE* is solved as the residual component of the surface energy balance, *LE* = *Rn* − *H* − *G*. Therefore, an error in the calculation of energy fluxes (*Rn*, *H*, and *G*) adversely affects the estimation of *LE*. Based on Figure 6, the *LE* estimation (or bias) is mainly influenced by the estimation of *H*. This conclusion was also reached by Kustas et al. [70], who showed the discrepancies between modeled and measured *LE* is due in large part, up to approximately 90%, to errors in modeled *H*.

*Tc* and *Ts* were estimated based on a linear *LST-NDVI* relationship as described by Nieto et al. [2]. However, this relationship does not fulfill the homoscedasticity criterion when the spatial domain/resolution reaches a certain size (i.e., 30-m) as shown in Figure 7. For example, in the case of a 30-m grid size, a higher variability is observed in the *LST-NDVI* data compared with finer resolutions (3.6 m, 7.2 m, and 14.4 m). At micro-scale (e.g., 3.6 m), there are small number of pixels inside the spatial domain compared with others (7.2 m, 14.4 m, and 30 m), and exhibit an apparent linear relationship between *LST* and *NDVI*. However, at coarse resolution (e.g., 30 m), there are many more pixels, more rows of vineyard are included, and large vegetated and bare soil pixels exist inside the spatial domain. The result is a partially filled triangular shape. This indicates the relationship between *LST* and *NDVI* starts to resemble the "triangle method" [71] to estimate *ET* as the sampling

**Figure 7.** The *LST-NDVI* relationship used for finding *Tc* and *Ts* as proposed by the *TSEB2T* model at different spatial domains (09 August 2014). (**a**) 3.6 m, (**b**) 7.2 m, (**c**) 14.4 m, (**d**) 30 m. L = 0.21௩\_ − 0.004௩\_ೞೠ + 0.34 ൬ <sup>൰</sup> ௩\_

Figure 8 illustrates the *T<sup>c</sup>* and *T<sup>s</sup>* maps at different resolutions, which provide an indication of the loss in spatial variability due to spatial aggregation. The ranges of *T<sup>c</sup>* and *T<sup>s</sup>* were between 290 K (16.85 ◦C) and 320 K (46.85 ◦C) for the *sUAS* flight in 2014. <sup>−</sup> 0.94 ቀ0.23൫௩\_൯ ଶ ቁ − 2.8\_ೄವ ൬ <sup>൰</sup> ௩ೄವ − 0.7 (14) *LAI* values from the *GP* model compared with the actual *LAI* field measurements showed good agreement with an *R*2 of 0.73.

(**a**) **Figure 8.** *Cont.*

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 18 of 35

**Figure 8.** Example of (**a**) canopy temperature (*Tc*) and (**b**) soil temperature (*Ts*) in Kelvin (K) at different spatial domains for 09 August 2014. **Figure 8.** Example of (**a**) canopy temperature (*Tc*) and (**b**) soil temperature (*Ts*) in Kelvin (K) at different spatial domains for 09 August 2014.

To evaluate the difference between multiple model grid sizes of *LAI* for each flight, *LAI* maps at different resolutions were estimated (see Figure 9) and statistics including the spatial mean, standard (b) Leaf Area Index (*LAI*)

deviation, and coefficient of variation (*CV*) were calculated as shown in Table 4. Figure 9 provides an indication of the loss in spatial variability in *LAI* images due to spatial aggregation. *LAI* at each contextual spatial domain/resolution was calculated using the *LAI* model (Equation (14)). Each parameter in that equation was calculated based on the pixel values inside the model grid. The ranges of *LAI* were between 0 and 2.5 for the *sUAS* flight in 2014. As illustrated in Table 4, the spatial mean With the *GP* model results, it was found that the main estimators for computing *LAI* are the mean of *NIR*/*R* ratio of the vine, area of the vine, sum of *NDVI* of the vine, standard deviation of *NIR* of the interrow, and standard deviation of *NIR*/*R* ratio of the vine. The *GP* model (Equation (14)) was applied to the remote-sensing imagery to map spatial *LAI* distribution across the study area.

$$LAI = 0.21 \text{NDVI}\_{\text{\$\upsilon\_{\text{s}}\$}} - 0.004 \text{NDVI}\_{\text{\$\upsilon\_{\text{s}}\$}} + 0.34 \left(\frac{\text{NLP}}{\text{\$\upsilon\_{\text{s}}\$}}\right)\_{\text{\tiny\upsilon\_{\text{s}}\$\text{max}}} - \frac{0.94}{\exp\left(0.23(\text{NDVI}\_{\text{\$\upsilon\_{\text{s}}\$}})^{2}\right)} - 2.8 \text{NIR}\_{\text{\tiny\upsilon\_{\text{s}}\text{STD}}} \left(\frac{\text{NIR}}{\text{\$\upsilon\_{\text{s}}\$}}\right)\_{\text{\tiny\upsilon\_{\text{STD}}}} - 0.7 \tag{14}$$

*LAI* under low *LAI* conditions using *VIs*. The frequency histogram in Figure 10 indicates the distribution of values is skewed such that the lower values are more pronounced for the flight of 2 *LAI* values from the *GP* model compared with the actual *LAI* field measurements showed good agreement with an *R* <sup>2</sup> of 0.73.

May 2016, with a non-significant change between curves from the different grid sizes, except the 30 m resolution spatial domain, which shows a higher variation compared with other scales. This behavior aligns with the decreasing *CV* values due to loss in internal or pixel variability of the *LAI* values. A similar trend of lower *CV* toward large scale (30 m) has been observed for other *TSEB2T* inputs including *hc*, *fc*, and *wc*/*hc*. To evaluate the difference between multiple model grid sizes of *LAI* for each flight, *LAI* maps at different resolutions were estimated (see Figure 9) and statistics including the spatial mean, standard deviation, and coefficient of variation (*CV*) were calculated as shown in Table 4. Figure 9 provides an indication of the loss in spatial variability in *LAI* images due to spatial aggregation. *LAI* at each contextual spatial domain/resolution was calculated using the *LAI* model (Equation (14)). Each parameter in that equation was calculated based on the pixel values inside the model grid. The ranges of *LAI* were between 0 and 2.5 for the *sUAS* flight in 2014. As illustrated in Table 4, the spatial mean value (µ) is the same across different scales, with a slight decrease in *CV*. The exception is the flight on 2 May 2016, which represents the early growing stage of the vine canopy with active/live interrow cover crop, showing a higher *CV*. Hardin and Jensen [72] also found greater uncertainty in estimating *LAI* under low *LAI* conditions using *VIs*. The frequency histogram in Figure 10 indicates the distribution of values is skewed such that the lower values are more pronounced for the flight of 2 May 2016, with a non-significant change between curves from the different grid sizes, except the 30-m resolution spatial domain, which shows a higher variation compared with other scales. This behavior aligns with the decreasing *CV* values due to loss in internal or pixel variability of the *LAI* values. A similar trend of lower *CV* toward large scale (30 m) has been observed for other *TSEB2T* inputs including *hc*, *fc*, and *wc*/*hc*.

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 19 of 35

**Figure 9.** Example of modeled *LAI* (unitless) across different spatial domains for 09 August 2014. **Figure 9.** Example of modeled *LAI* (unitless) across different spatial domains for 09 August 2014.


**Table 4.** Spatial domain effect on LAI estimation*.* **Table 4.** Spatial domain effect on *LAI* estimation.

7.2 m 0.06 0.10 1.75 14.4 m 0.06 0.10 1.66 30.0 m 0.06 0.09 1.59

**02 May 2016** 

**Figure 10.** Frequency curve of *LAI* at different times from 3.6 m and 7.2 m, 14.4 m and 30 m. **Figure 10.** Frequency curve of *LAI* at different times from 3.6 m and 7.2 m, 14.4 m and 30 m.

3.2.2. Contextual Spatial Domain Effect on Field-Scale *LE* Estimation

3.2.2. Contextual Spatial Domain Effect on Field-Scale *LE* Estimation An example of the maps of *LE* across different model grid sizes is shown in Figure 11. The maps of the energy balance components for 2014 flight at different resolutions are shown in Appendix A. The statistics (mean (*μ*) and coefficient of variation (*CV*)) for the *LE* maps at the different modeling resolutions are illustrated as bar graphs in Figures 12 and 13, respectively. For *LE*, the highest mean value is on 02 May 2016, at midafternoon. Although the grapevine canopy is fully developed by June, *LE* in May at both overpass times is higher than the acquisition in June, July, and August. However, on 3 May, the model yields the lowest *LE* values due to overcast conditions that day significantly reducing incoming solar radiation, and hence, the energy fluxes. The phenocam data (https://hrsl.ba.ars.usda.gov/awhite/CAM/) indicate the high rate of *LE* on 2 May is the result of a An example of the maps of *LE* across different model grid sizes is shown in Figure 11. The maps of the energy balance components for 2014 flight at different resolutions are shown in Appendix A. The statistics (mean (µ) and coefficient of variation (*CV*)) for the *LE* maps at the different modeling resolutions are illustrated as bar graphs in Figures 12 and 13, respectively. For *LE*, the highest mean value is on 02 May 2016, at midafternoon. Although the grapevine canopy is fully developed by June, *LE* in May at both overpass times is higher than the acquisition in June, July, and August. However, on 3 May, the model yields the lowest *LE* values due to overcast conditions that day significantly reducing incoming solar radiation, and hence, the energy fluxes. The phenocam data (https://hrsl.ba.ars.usda.gov/awhite/CAM/) indicate the high rate of *LE* on 2 May is the result of a rapidly developing vine canopy, together with a transpiring cover crop.

rapidly developing vine canopy, together with a transpiring cover crop. At a contextual spatial domain level, the magnitude of *LE* is degraded as shown in Figure 12 due to the data aggregation from the 3.6-m grid to Landsat scale (30 m). For example, the mean *LE* value from *TSEB2T* on 02 May 2016 at midafternoon was 315 W/m2 for the 3.6-m grid decreasing to 304 W/m2 for the 7.2-m grid, then decreases further to 293 W/m2 and 278 W/m2, respectively, for 14.4-m and 30-m grids. As shown in Figure 13, *CV* value slightly increases as the model grid scale/resolution size increases despite a decrease in variation of *LAI* and *LST* distribution as seen in Section 3.2.1. While *LE* degrades, the *CV* values do not show significant differences. This can be due to internal At a contextual spatial domain level, the magnitude of *LE* is degraded as shown in Figure 12 due to the data aggregation from the 3.6-m grid to Landsat scale (30 m). For example, the mean *LE* value from *TSEB2T* on 02 May 2016 at midafternoon was 315 W/m<sup>2</sup> for the 3.6-m grid decreasing to 304 W/m<sup>2</sup> for the 7.2-m grid, then decreases further to 293 W/m<sup>2</sup> and 278 W/m<sup>2</sup> , respectively, for 14.4-m and 30-m grids. As shown in Figure 13, *CV* value slightly increases as the model grid scale/resolution size increases despite a decrease in variation of *LAI* and *LST* distribution as seen in Section 3.2.1. While *LE* degrades, the *CV* values do not show significant differences. This can be due to internal *TSEB2T* compensation of the energy balance components at the different evaluated scales.

*TSEB2T* compensation of the energy balance components at the different evaluated scales.

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 21 of 35

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 21 of 35

**Figure 11.** *LE* (w/m2) aggregation at 3.6 m, 7.2 m, 14.4 m and 30 m for 09 August 2014. **Figure 11.** *LE* (W/m<sup>2</sup> ) aggregation at 3.6 m, 7.2 m, 14.4 m and 30 m for 09 August 2014.

**Figure 12.** Spatial domain effect on the mean of *LE* spatial data at different times. **Figure 12.** Spatial domain effect on the mean of *LE* spatial data at different times.

9-Aug-14 2-Jun-15 2-Jun-15 11-Jul-15

(a) (b) (c) (d)

Landsat time Landsat time Afternoon Landsat time

(d) (e) (f) (g)

11-Jul-15 2-May-16 2-May-16 3-May-16

Afternoon Afternoon Midafternoon Afternoon

3.6m 7.2m 14.4m 30m

**Figure 13.** Spatial domain effect on the coefficient of variation (*CV*) of *LE* spatial data at different

3.6m 7.2m 14.4m 30m

times.

0.00 0.10 0.20 0.30 0.40 0.50 0.60

0.00 0.10 0.20 0.30 0.40 0.50 0.60

CV

CV

Mean(LE), W/m2

11-Jul-15 2-May-16 2-May-16 3-May-16

(e) (f) (g) (h)

Afternoon Afternoon Midafternoon Afternoon

3.6m 7.2m 14.4m 30m

**Figure 13.** Spatial domain effect on the coefficient of variation (*CV*) of *LE* spatial data at different times.

#### **Figure 13.** Spatial domain effect on the coefficient of variation (*CV*) of *LE* spatial data at different times. 3.2.3. Contextual Spatial Domain Effect on *LE* Statistical Characteristics

To provide quantitative evaluation of the impact of spatial aggregation of inputs on *LE* estimation for the resulting pixel values, frequency and cumulative distribution plots for the *LE* maps are illustrated in Figure 14. This figure shows that *LE* varies at different grid sizes. The cumulative frequency distribution curves indicate that, especially at the 30-m grid size, *LE* distribution tends to have the highest cumulative values at lower *LE* range (below 300 W/m<sup>2</sup> ). A magnitude shift towards lower *LE* persists across different times, with one exception. In the case of a 30-m grid on 02 June 2015, the frequency moved up then decreased below the frequency curves of other grid sizes (3.6 m, 7.2 m, and 14.4 m). In general, the results in Figure 14 show a reduction in *LE* distribution as the scale becomes coarser. Hong et al. [22] indicated that an increase in the peak of the *LE* histogram curve spans as much 10% to 20% as a response to spatial data aggregation using *SEBAL*. In the *TSEB* model, the soil and vegetation components of the scene are treated separately, while the *SEBAL* model uses a single source approach using the composite soil/canopy temperature and is contextual defining wet and dry *ET* limits based on the hot and cold extremes in the *LST* field within the image [73]. Moreover, Ershadi et al. [14] pointed out three possible reasons behind the different results obtained from *ET* models: (a) the approach (e.g., contextual hot/cold surface temperature limits versus using absolute surface-atmosphere temperature differences) of each model to estimate *ET*, (b) the study area and eco-hydrological conditions of the surface, which may favor certain *ET* model parameterizations over others, or (c) the different models of aerodynamic resistance formulations and sensitivity to the roughness parameters.

**Figure 14.** Frequency curve (left) and cumulative frequency distribution (right) plots of instantaneous *LE* for all *sUAS* flights at 3.6 m, 7.2 m, 14.4 m, and 30 m. **Figure 14.** Frequency curve (**left**) and cumulative frequency distribution (**right**) plots of instantaneous *LE* for all *sUAS* flights at 3.6 m, 7.2 m, 14.4 m, and 30 m.

Increasing the spatial domain/resolution affects the estimation of *TSEB2T* parameters as the fine details of the surface disappear. To test these claims, *R<sup>a</sup>* (s/m) and *LST-NDVI* relationship were evaluated at different spatial domain/resolution; the latter is shown in Section 3.2.1 (a). As shown in Figure 15, there is a decreasing trend in the relative spatial mean (µ*r*) of *R<sup>a</sup>* for all flights, ranging approximately from 20% to 60%. The high variability in *R<sup>a</sup>* is related mainly to the variables that affect the friction velocity (*u*∗), which the mean canopy height and roughness length (*zoH*), which are derived from the imagery at different resolution/spatial domain. This finding is in agreement with Ershadi et al. [14] and Moran et al. [15], who indicated that the reduction of *R<sup>a</sup>* value at coarse spatial domain *Remote Sens.*  /resolution is a key factor behind the underestimation of **2019**, *11*, x FOR PEER REVIEW *LE*. 25 of 35

**Figure 15.** Variation of the relative spatial mean () of *Ra* for different flights. **Figure 15.** Variation of the relative spatial mean (µ*r*) of *R<sup>a</sup>* for different flights.

#### 3.2.4. Effects of Model Grid Size on *LE* 3.2.4. Effects of Model Grid Size on *LE*

To evaluate quantitatively the impact of model grid size via the resolution of key input data, the relative difference (relative error) () was computed using as the reference the *LE* at 3.6-m model grid size/resolution. For example, the *LE* value at the 7.2-m grid is compared to the *LE* at the 3.6-m grid size by resampling the 7.2-m grid to a 4 × 4 set of 3.6-m *LE* output which will have a uniform *LE*value at the finer resolution, and taking the difference. As illustrated in Figure 16, is calculated with the mean and percentiles (25th and 75th) for the coarser grid sizes used in the *TSEB2T* model for the different *sUAS* acquisitions. The plots demonstrate an increasing trend in as the model grid size/resolution increases/decreases. The largest value was computed for the imagery on 11 July 2015 at afternoon at nearly 45% for the Landsat resolution. In contrast to 11 July 2015, the lowest range of relative error was observed on 09 August 2014, where the ranged approximately To evaluate quantitatively the impact of model grid size via the resolution of key input data, the relative difference (relative error) (*Er*) was computed using as the reference the *LE* at 3.6-m model grid size/resolution. For example, the *LE* value at the 7.2-m grid is compared to the *LE* at the 3.6-m grid size by resampling the 7.2-m grid to a 4 × 4 set of 3.6-m *LE* output which will have a uniform *LE*-value at the finer resolution, and taking the difference. As illustrated in Figure 16, *E<sup>r</sup>* is calculated with the mean and percentiles (25th and 75th) for the coarser grid sizes used in the *TSEB2T* model for the different *sUAS* acquisitions. The plots demonstrate an increasing trend in *E<sup>r</sup>* as the model grid size/resolution increases/decreases. The largest *E<sup>r</sup>* value was computed for the imagery on 11 July 2015 at afternoon at nearly 45% for the Landsat resolution. In contrast to 11 July 2015, the lowest range of relative error was observed on 09 August 2014, where the *E<sup>r</sup>* ranged approximately between 15% for

between 20% and 40% when aggregating the Landsat data incrementally from 120 m to 960 m and using the *SEBS* model to calculate surface heat fluxes. Furthermore, Moran et al. [15] indicated that a larger error could appear (larger than 50%) in *H* estimation over a heterogeneous area due to a mix of stable and unstable conditions and the variation in aerodynamic roughness, especially for highly unstable conditions. As previously mentioned in Section 3.1, the underestimated *LE* could be influenced by overestimation in *H*, which implies that a large error is expected in the residual flux (*LE*) estimate at coarse spatial domains [70]. Furthermore, the effect of model grid size on *LE* is also visible at the 25th and 75th percentiles, which immediately increases at the 7.2-m grid size and

between 15% for the 7.2-m grid and 25% for the 30-m grid. On an average, value ranged from

the 7.2-m grid and 25% for the 30-m grid. On an average, *E<sup>r</sup>* value ranged from approximately 25% using the 7.2-m model grid size to 40% with the 30-m model resolution. continues increasing towards the 30-m resolution, providing a clear indication of increasing discrepancy with the reference grid (3.6 m) *LE* estimates.

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 26 of 35

**Figure 16.** Relative error (*Er*) at different spatial resolutions for *LE* with the triangle symbols indicating mean and light lines indicating the 25th and 75thpercentiles for the coarse grid sizes. **Figure 16.** Relative error (*Er*) at different spatial resolutions for *LE* with the triangle symbols indicating mean and light lines indicating the 25th and 75th percentiles for the coarse grid sizes.

**4. Conclusions**  The objective of this study was to assess high-resolution *LE* estimation in vineyards at different model grid sizes or resolutions, specifically 3.6 m, 7.2 m, 14.4 m, and 30 m (Landsat scale), using a physically-based *ET* model known as *TSEB2T*. The reference grid size of 3.6 m represents the finest pixel resolution that includes both vine canopy and interrow conditions, which is the resolution where the *TSEB* model algorithms of soil/substrate and canopy temperature partitioning radiation and convective energy exchange are applicable [2]. Multiple statistical measures were used to assess the effect of decreasing the spatial resolution or increasing the model grid size 2, 4, and nearly 10 times the original 3.6 m resolution. These included validation of *TSEB2T* fluxes at the different model grid sizes with the *EC* measurements, comparing *LE* spatial statistics (mean and coefficient of These results are supported by an Ershadi et al. [14] study that found the *E<sup>r</sup>* of *LE* varied between 20% and 40% when aggregating the Landsat data incrementally from 120 m to 960 m and using the *SEBS* model to calculate surface heat fluxes. Furthermore, Moran et al. [15] indicated that a larger error could appear (larger than 50%) in *H* estimation over a heterogeneous area due to a mix of stable and unstable conditions and the variation in aerodynamic roughness, especially for highly unstable conditions. As previously mentioned in Section 3.1, the underestimated *LE* could be influenced by overestimation in *H*, which implies that a large error is expected in the residual flux (*LE*) estimate at coarse spatial domains [70]. Furthermore, the effect of model grid size on *LE* is also visible at the 25th and 75th percentiles, which immediately increases at the 7.2-m grid size and continues increasing towards the 30-m resolution, providing a clear indication of increasing discrepancy with the reference grid (3.6 m) *LE* estimates.

#### resolutions using *LE* at 3.6 m grid size as the reference. The results showed that separation of *Tc* and **4. Conclusions**

*Ts*, required in *TSEB2T*, affects the *LST-NDVI* linear trend as a function of resolution of the pixels. The validation results with the flux tower measurements indicate that *Rn* and *G* discrepancies do not change across different model grid sizes, while for *H* and *LE* there is an increase in modelmeasurement differences, particularly at the 30-m resolution. This is largely caused by an overestimation in *H*, causing an underestimation in *LE* (bias), particularly at the coarsest resolution The objective of this study was to assess high-resolution *LE* estimation in vineyards at different model grid sizes or resolutions, specifically 3.6 m, 7.2 m, 14.4 m, and 30 m (Landsat scale), using a physically-based *ET* model known as *TSEB2T*. The reference grid size of 3.6 m represents the finest pixel resolution that includes both vine canopy and interrow conditions, which is the resolution where the

variation, frequency distributions) and *LE* differences over the imaged domain at the different

*TSEB* model algorithms of soil/substrate and canopy temperature partitioning radiation and convective energy exchange are applicable [2]. Multiple statistical measures were used to assess the effect of decreasing the spatial resolution or increasing the model grid size 2, 4, and nearly 10 times the original 3.6 m resolution. These included validation of *TSEB2T* fluxes at the different model grid sizes with the *EC* measurements, comparing *LE* spatial statistics (mean and coefficient of variation, frequency distributions) and *LE* differences over the imaged domain at the different resolutions using *LE* at 3.6 m grid size as the reference. The results showed that separation of *T<sup>c</sup>* and *T<sup>s</sup>* , required in *TSEB2T*, affects the *LST-NDVI* linear trend as a function of resolution of the pixels. The validation results with the flux tower measurements indicate that *R<sup>n</sup>* and *G* discrepancies do not change across different model grid sizes, while for *H* and *LE* there is an increase in model-measurement differences, particularly at the 30-m resolution. This is largely caused by an overestimation in *H*, causing an underestimation in *LE* (bias), particularly at the coarsest resolution (30-m grid size). This refers mainly to the non-linear relationship of *LST-NDVI* and the variability of *R<sup>a</sup>* due to the variables that affect the *u*<sup>∗</sup> which are the mean canopy height and roughness length, which are derived from remote sensing imagery at different spatial domain/resolution.

The effects of model grid size were evaluated at field and at grid scale using the spatial mean and coefficient of variation and relative difference, respectively. At field scale, the results show small decreases in the spatial mean over the image, ranging approximately from 10% to 20%, as the data aggregated for model grid size increased from 3.6 m to 30 m. However, the relative differences with resolution indicate a significant decrease in *LE*, ranging approximately from 25% to 45%, when aggregating the data from 3.6 m to Landsat scale (30 m). This means that, while field values of *LE* may be adequate to use, the field variability reduction limits its use for precision agriculture applications, such as identifying areas within the field under actual stress conditions or being over irrigated. These results suggest that *TSEB2T* is only applicable using imagery with high enough resolution that can readily distinguish plant canopy and soil/substrate temperatures and the modeling grid size is at a resolution where it is appropriate to apply *TSEB2T* algorithms for modeling the radiative and convective energy exchange from both the vegetation and soil substrate systems. Aggregating inputs to *TSEB2T* to multiple grid sizes of the interrow/row spacings for vineyards is not advisable, since it is likely the accuracy of surface fluxes, particularly *LE*, will deteriorate. While this study was limited to evaluating different modeling grid sizes, a future comparison with Landsat and ECOstress *ET* products is also planned, which would provide a more comprehensive scaling assessment of *ET* estimates for *sUAS*-Satellite *ET* integration. Furthermore, the effect of remote sensing resolution on the output of other *TSEB* versions such as *TSEB-PT* may be less affected and will be evaluated in a future study.

**Author Contributions:** A.N. designed, analyzed, and wrote the paper. A.T.-R. supervised the research, contributed with ideas during the interpretation of results, and reviewed the paper. W.K. supervised the research and reviewed the paper. M.M. contributed with ideas during the interpretation of results and reviewed the paper. D.S. contributed in the statistical analysis. L.H. contributed with ideas and reviewed the paper. H.N., J.A., J.P., L.M., M.M.A., C.C., L.S., and N.D. helped in collecting micrometeorological and ground measurements. All authors have read and agreed to the published version of the manuscript.

**Funding:** Funding provided by E.&J. Gallo Winery. Utah Water Research Laboratory contributed towards the acquisition and processing of the ground truth and *UAV* imagery data collected during *GRAPEX IOPs*. We would like to acknowledge the financial support for this research from NASA Applied Sciences-Water Resources Program [NNX17AF51G] and the USDA Non Assistance Cooperative Agreement 58-8042-5-092 funding. USDA is an equal opportunity provider and employer.

**Acknowledgments:** We would like to thank Aggieair Service Center team (Ian Gowing, Mark Winkelaar, and Shannon Syrstad) for their extraordinary support in this research, whose cooperation greatly improved the data collection and data processing, and the staff of Viticulture, Chemistry and Enology Division of E.&J. Gallo Winery for the assistance in the collection and processing of field data during *GRAPEX IOPs*. This project would not have been possible without the cooperation of Ernie Dosio of Pacific Agri Lands Management, along with the Sierra Loma vineyard staff, for logistical support of *GRAPEX* field and research activities. The authors would like to thank Carri Richards for editing this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

*Remote Sens.* **2020**, *12*, 342

#### **Appendix A** *Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 28 of 35 **Appendix A Appendix A**

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 28 of 35

**Figure A1.** Example of modeled *hc* (m) across different spatial domains for 09 August 2014. **Figure A1.** Example of modeled *hc* (m) across different spatial domains for 09 August 2014. **Figure A1.** Example of modeled *hc* (m) across different spatial domains for 09 August 2014.

**Figure A2.** Example of modeled *fc* across different spatial domains for 09 August 2014. **Figure A2. Figure A2.** Example of modeled Example of modeled *f fc* across different spatial domains for 09 August 2014. *<sup>c</sup>* across different spatial domains for 09 August 2014.

**LE (W/m2)**

(**b**)

**Figure A3.** Example of Modeled *wc*/*hc* different spatial domains for 09 August 2014. **Figure A3.** Example of Modeled *wc*/*hc* different spatial domains for 09 August 2014. *Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 30 of 35

**Figure A4.** *Cont.*

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 30 of 35

**Figure A4.** *Cont.*

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 31 of 35

**Figure A4.** Aggregation of surface energy fluxes across different spatial domain (3.6m, 7.2m, 14.4 m, and 30 m) on August 09, 2014: (**a**) *Rn*, (**b**) *LE*, (**c**) *H*, and (**d**) *G*. **Figure A4.** Aggregation of surface energy fluxes across different spatial domain (3.6 m, 7.2 m, 14.4 m, and 30 m) on 09 August 2014: (**a**) *Rn*, (**b**) *LE*, (**c**) *H*, and (**d**) *G*.

#### **References**

MD, USA, 2019; pp.15–17.


*Autonomous Air and Ground Sensing Systems for Agricultural Optimization and Phenotyping IV.*; SPIE: Baltimore,


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).
