**Influence of Model Grid Size on the Estimation of Surface Fluxes Using the Two Source Energy Balance Model and sUAS Imagery in Vineyards**

**Ayman Nassar 1,\*, Alfonso Torres-Rua <sup>1</sup> , William Kustas <sup>2</sup> , Hector Nieto <sup>3</sup> , Mac McKee <sup>1</sup> , Lawrence Hipps <sup>4</sup> , David Stevens <sup>1</sup> , Joseph Alfieri <sup>2</sup> , John Prueger <sup>5</sup> , Maria Mar Alsina <sup>6</sup> , Lynn McKee <sup>2</sup> , Calvin Coopmans <sup>7</sup> , Luis Sanchez <sup>6</sup> and Nick Dokoozlian <sup>6</sup>**


Received: 25 October 2019; Accepted: 10 January 2020; Published: 21 January 2020

**Abstract:** Evapotranspiration (*ET*) is a key variable for hydrology and irrigation water management, with significant importance in drought-stricken regions of the western US. This is particularly true for California, which grows much of the high-value perennial crops in the US. The advent of small Unmanned Aerial System (*sUAS*) with sensor technology similar to satellite platforms allows for the estimation of high-resolution *ET* at plant spacing scale for individual fields. However, while multiple efforts have been made to estimate *ET* from *sUAS* products, the sensitivity of *ET* models to different model grid size/resolution in complex canopies, such as vineyards, is still unknown. The variability of row spacing, canopy structure, and distance between fields makes this information necessary because additional complexity processing individual fields. Therefore, processing the entire image at a fixed resolution that is potentially larger than the plant-row separation is more efficient. From a computational perspective, there would be an advantage to running models at much coarser resolutions than the very fine native pixel size from *sUAS* imagery for operational applications. In this study, the Two-Source Energy Balance with a dual temperature (*TSEB2T*) model, which uses remotely sensed soil/substrate and canopy temperature from *sUAS* imagery, was used to estimate *ET* and identify the impact of spatial domain scale under different vine phenological conditions. The analysis relies upon high-resolution imagery collected during multiple years and times by the Utah State University *AggieAirTM sUAS* program over a commercial vineyard located near Lodi, California. This project is part of the USDA-Agricultural Research Service Grape Remote Sensing Atmospheric Profile and Evapotranspiration eXperiment (*GRAPEX*). Original spectral and thermal imagery data from *sUAS* were at 10 cm and 60 cm per pixel, respectively, and multiple spatial domain scales (3.6, 7.2, 14.4, and 30 m) were evaluated and compared against eddy covariance (*EC*) measurements. Results indicated that the *TSEB2T* model is only slightly affected in the estimation of the net radiation (*Rn*) and the soil heat flux (*G*) at different spatial resolutions, while the sensible and latent heat fluxes (*H* and *LE*, respectively) are significantly affected by coarse grid sizes. The results indicated overestimation of

*H* and underestimation of *LE* values, particularly at Landsat scale (30 m). This refers to the non-linear relationship between the land surface temperature (*LST*) and the normalized difference vegetation index (*NDVI*) at coarse model resolution. Another predominant reason for *LE* reduction in *TSEB2T* was the decrease in the aerodynamic resistance (*Ra*), which is a function of the friction velocity (*u*∗) that varies with mean canopy height and roughness length. While a small increase in grid size can be implemented, this increase should be limited to less than twice the smallest row spacing present in the *sUAS* imagery. The results also indicated that the mean *LE* at field scale is reduced by 10% to 20% at coarser resolutions, while the with-in field variability in *LE* values decreased significantly at the larger grid sizes and ranged between approximately 15% and 45%. This implies that, while the field-scale values of *LE* are fairly reliable at larger grid sizes, the with-in field variability limits its use for precision agriculture applications.

**Keywords:** evapotranspiration (*ET*); *GRAPEX*; *sUAS*; remote sensing; Two Source Energy Balance model (*TSEB*); contextual spatial domain/resolution; data aggregation; eddy covariance (*EC*)

#### **1. Introduction**

Evapotranspiration (*ET*) is a key factor in the hydrologic cycle and in irrigation demand. Conventional methods for estimating *ET*, such as lysimeters and flux towers, are limited to sampling small areas on the order of 10<sup>1</sup> to 10<sup>3</sup> m<sup>2</sup> . For that, a more efficient method is needed as *ET* varies spatially under different micro-meteorological and vegetative conditions. Accordingly, spatially distributed data are important for mapping *ET* variations over large areas, particularly in agricultural regions containing many of crop types and growth stages. In recent decades, remote sensing products from various platforms and at various spatial resolutions have been applied in modeling different environmental processes (e.g., surface energy fluxes, water and carbon balance, net primary productivity) [1]. Improved sensor systems and methods in remote sensing, and particularly the advent of small unmanned aerial systems (*sUAS*), have made these technologies a valuable source of spatial information for *ET* estimation at the canopy level. *sUAS* can offer spatial coverage with sub-meter-resolution imagery for mapping canopy and soil temperature, which are the key surface states for estimating *ET* [2]. While satellites are characterized by either coarse resolution and high temporal frequency or by high spatial resolution and low repeatability [3], *sUAS* technology, in addition to offering high-resolution data [4–6], can be described as "flexible on timing" [7]. This means that remotely sensed information can be obtained when needed or on demand using *sUAS*. For these reasons, various methods are under development to employ *sUAS* data for *ET* estimation [2].

Remote sensing is a valuable source for accessing land surface spatial information [8]. Nonetheless, spatial scaling is recognized as a challenging issue, particularly in surface-atmosphere exchange [8,9], environmental modeling, and agricultural management [10] applications and research. Previous studies by Brunsell and Gillies [11] and Giorgi [12] indicated that spatial scaling becomes more complex in cases of heterogeneous land surfaces, and homogeneity is less likely to be met in reality [13]. Various models have been developed to describe aerodynamic or energy balance fluxes, but these models assume homogeneity in terms of agricultural type, surface roughness, surface temperature, and meteorological condition [13,14]. Heat fluxes, including latent heat flux (*LE*) and sensible heat flux (*H*), are highly influenced by land surface heterogeneity [15]. Therefore, the variability in land cover within a pixel or model grid size can result in significant error in the mean pixel or grid heat flux estimation [16]. Vegetated areas with partial canopy cover will have underlying soil/substrate affecting the remotely sensed data, and hence, require models that explicitly consider the different effects of these two sources on energy exchange and sensor integration [2]. Typically, remotely sensed data at different resolutions are employed as an approximation to describe the spatial variability of the interaction between surface and atmosphere [11]. Current and future developments in remote sensing, with information spanning from sub-meters to kilometers, are making upscaling (data aggregation) a crucial issue in scientific and methodological advances. This is particularly true for understanding the physics behind climate, weather, and the surface energy balance [13,17].

In general, spatial aggregation can be performed under two different procedures: forcing inputs to a coarser resolution or aggregating the derived fluxes from initial high-resolution data (contextual spatial domain). Long et al. [18] pointed out that forcing spatial data aggregation from Landsat bands to MODIS (Moderate Resolution Imaging Spectroradiometer) resolution results in different statistical and spatial properties in *ET* estimates than at the original Landsat resolution. Study cases of *LE* resulted in inaccuracies [19,20] due to a reduction in surface variability at MODIS scale [11]. Moreover, the structure of vegetation and aerodynamic roughness influence the aggregation of turbulent fluxes and produce bias when MODIS data is used [15]. On the other hand, Bian and Butler [21] showed that low-resolution data could retain the statistical characteristics of the original data using specific aggregation techniques such as average and median. In addition, the spatial aggregation of *ET* inputs removes the effects of heterogeneity on the land surface. Still, scaling up energy fluxes from Landsat to MODIS scale is necessary in large-scale environmental models [22]. However, Landsat resolution is needed for validating modeled outputs using flux towers [23].

Several methods exist for spatial aggregation of *ET* data, but they are in the exploratory stage [24]. Ershadi et al. [14] demonstrated that *ET* results reduced by 15% when aggregating Landsat TM (Thematic Mapper) imagery by 50% using the Surface Energy Balance System (*SEBS*) model. The *ET* reduction was caused by the decrease in roughness parameterization [14]. This outcome was also supported by Brunsell and Gillies [11], who indicated that the land surface heterogeneity is highly influenced by the input forcing aggregation of Landsat TM data affecting the surface heat fluxes. In contrast, French et al. [25] found no significant difference in daily *ET* estimates when they used *METRIC* (Mapping EvapoTranspiration at high Resolution with Internalized Calibration) model and upscaled data acquired by aircraft to Landsat resolution. However, another study by Kustas and Norman [16] that used a detailed soil-vegetation atmosphere simulation model along with the thermal-based two-source energy balance model found that varying the degree of heterogeneity within a pixel, either in terms of surface roughness, moisture status, or a combination thereof, can have a significant impact on the pixel aggregated flux.

A key question related to data aggregation was raised by Su et al. [26]: "How does the level of aggregation affect surface energy fluxes as fluxes are aggregated from the resolution at which they are observed to the coarse grid cell size of the atmospheric model?". The study conducted by Guzinskia and Nieto [27] aimed to estimate *ET* using a Two Source Energy Balance (*TSEB*) model. They reported that sharpening Sentinel 3 thermal imagery at 1-km pixel resolution to higher resolution (20 m) visible/near-infrared is indicative of the main issue of the lack of fine resolution thermal-IR (InfraRed) data for input to remote sensing-based *ET* models, particularly when applied to agricultural areas. In addition, Niu et al. [28] indicated that the *TSEB* model *ET* output using *sUAS* imagery gives more reliable estimates compared to coarse-resolution data because the model can separate between canopy and soil components. Moreover, most previous studies exploring the effects of sensor resolution on modeled *ET* have used semi-empirical models (e.g., Surface Energy Balance Algorithm for Land (*SEBAL*) model) [14], while physically-based *ET* models are required to quantify changes in the water and energy exchange due to changes in fractional vegetation cover, roughness, canopy structure, phenology, etc. that are occurring at plant scale [29]. In addition, it is common knowledge that vineyards and orchard fields do not have the same row spacing. The spacing varies from 6 ft to 12 ft for vineyards [30] and from 8 ft to 18 ft for orchards [31].

In the same context as the investigations discussed above on spatial resolution and surface heterogeneity, this study investigates the impact of grid-size resolution on *LE* outputs from *TSEB* model using the component soil/substrate and canopy temperature version (*TSEB2T*) model applied to a complex agricultural canopy, namely a vineyard in California's Central Valley. The study directly quantifies the effect of sensor resolution on key *TSEB* model inputs (i.e., land surface temperature

(*LST*), leaf area index (*LAI*), canopy height (*hc*), canopy width-to-height ratio (*wc*/*hc*), and fractional cover (*fc*)) for estimating surface energy balance/*ET*. High-resolution optical and thermal data were acquired by an *sUAS* platform for vine and cover crop phenological stages at several different times during the day. In this research effort, the topics investigated include determining (a) whether the separation between canopy and soil/substrate temperature (*T<sup>c</sup>* and *T<sup>s</sup>* , respectively) using *TSEB2T* is valid for coarse spatial domains (e.g., towards Landsat scale); (b) the effect of spatial resolution of *TSEB2T* inputs on the magnitude and spatial variation of *LE*; (c) if the different spatial domain scales/pixel resolutions under study (3.6, 7.2, 14.4 and 30 m) have an impact on the magnitude of the *LE* and quantify the discrepancies as a function of resolution. *Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 4 of 35 during the day. In this research effort, the topics investigated include determining (a) whether the

#### *1.1. TSEB2T Model* separation between canopy and soil/substrate temperature (*Tc* and *Ts*, respectively) using *TSEB2T* is

and vegetation as follows:

*TSEB2T* is a physically based approach developed by Norman et al. [32] that explicitly accommodates the difference between aerodynamic and radiometric surface temperature that affect the radiative and convective exchange of energy between soil and canopy systems and the lower atmosphere. The main concept underpinning the *TSEB2T* approach is modeling of the partitioning of radiative and turbulent energy fluxes between canopy and soil systems. In this case, *H* is partitioned between soil and canopy, which is dependent mainly on *T<sup>c</sup>* and *T<sup>s</sup>* differences with the overlying atmosphere and their respective aerodynamic coupling. valid for coarse spatial domains (e.g., towards Landsat scale); (b) the effect of spatial resolution of *TSEB2T* inputs on the magnitude and spatial variation of *LE*; (c) if the different spatial domain scales/pixel resolutions under study (3.6, 7.2, 14.4 and 30 m) have an impact on the magnitude of the *LE* and quantify the discrepancies as a function of resolution. *1.1. TSEB2T Model TSEB2T* is a physically based approach developed by Norman et al. [32] that explicitly accommodates the difference between aerodynamic and radiometric surface temperature that affect

As shown in the Figure 1, the *TSEB2T* model separates the surface energy balance between soil and vegetation as follows: the radiative and convective exchange of energy between soil and canopy systems and the lower atmosphere. The main concept underpinning the *TSEB2T* approach is modeling of the partitioning

$$R\_{\rm tr} = LE + H + G\_{\rm \gamma} \tag{1}$$

$$R\_{\rm nc} = H\_{\rm c} + LE\_{\rm c} \tag{2}$$

$$R\_{\rm tts} = H\_{\rm s} + L E\_{\rm s} + G\_{\rm \prime} \tag{3}$$

where *R<sup>n</sup>* is the net radiation, *H* is the sensible heat flux, *LE* is the latent heat flux, and *G* is the soil heat flux. All units of fluxes are in W/m<sup>2</sup> . Subscripts of *c* and *s* represent the canopy and soil components, respectively. Because *T<sup>s</sup>* and *T<sup>c</sup>* can be derived from the *LST* with a high enough resolution of optical data, energy fluxes (*Rn*, *H*) can be calculated directly from the component temperatures (*T<sup>c</sup>* and *Ts*) and estimated aerodynamic resistances of canopy and soil components, while *G* is parametrized as a portion of soil net radiation (*Rns*). *LE<sup>c</sup>* and *LE<sup>s</sup>* are solved as residuals when (*T<sup>c</sup>* and *Ts*) observations are available. = + + ,(1) = + , (2) ௦ = ௦ + ௦ + , (3) where *Rn* is the net radiation, *H* is the sensible heat flux, *LE* is the latent heat flux, and *G* is the soil heat flux. All units of fluxes are in W/m2. Subscripts of *c* and *s* represent the canopy and soil components, respectively. Because *Ts* and *Tc* can be derived from the *LST* with a high enough

$$\mathbf{G} = \mathbf{c}\_{\mathbf{G}} \mathbf{R}\_{\mathbf{us}} \tag{4}$$

where *c<sup>G</sup>* is an empirical coefficient changing over the daytime [2]. while *G* is parametrized as a portion of soil net radiation (*Rns*). *LEc* and *LEs* are solved as residuals

where ீ is an empirical coefficient changing over the daytime [2].

of soil-vegetation resistance network as illustrated in Figure 1:

when (*Tc* and *Ts*) observations are available.

**Figure 1.** Schematic representation of *TSEB2T* model. **Figure 1.** Schematic representation of *TSEB2T* model.

To estimate the sensible heat flux for vegetation and canopy, Norman et al. [32] proposed a series

To estimate the sensible heat flux for vegetation and canopy, Norman et al. [32] proposed a series of soil-vegetation resistance network as illustrated in Figure 1:

$$H = H\_c + H\_s = \rho\_{air} \mathcal{C}\_p \frac{T\_{AC} - T\_A}{R\_A} = \rho\_{air} \mathcal{C}\_p \left[ \frac{T\_C - T\_{AC}}{R\_x} + \frac{T\_s - T\_{AC}}{R\_s} \right] \tag{5}$$

$$R\_A = \frac{\ln\left(\frac{z\_T - d\_0}{z\_{0H}}\right) - \Psi\_h \left(\frac{z\_T - d\_0}{L}\right) + \Psi\_h \left(\frac{z\_{0H}}{L}\right)}{\kappa' u\_\*} \tag{6}$$

where ρ*air* is the air density (kg/m<sup>3</sup> ); *C<sup>p</sup>* is the heat capacity of the air at constant pressure (J/(kg·K)); *T<sup>c</sup>* and *T<sup>s</sup>* are canopy and soil temperature (K), respectively; *TAC* is the temperature of canopy-air space (K); and *T<sup>A</sup>* is the temperature of air (K). *R<sup>A</sup>* is the aerodynamic resistance to heat transport from the soil/canopy system (s/m), *R<sup>x</sup>* is the boundary layer resistance of the canopy leaves (s/m), *R<sup>s</sup>* is the aerodynamic resistance to heat transport in the boundary layer close to the soil surface (s/m), *z<sup>T</sup>* is the measurement height for *TA*, *z*0*<sup>H</sup>* is the roughness length for heat transport, *d*<sup>0</sup> is the zero-plane displacement height (m), *L* is the Monin-Obukhov length (m), κ 0 = 0.4 is the von Karman's constant, *u*∗ is the friction velocity (m/s), and Ψ*<sup>h</sup>* is the adiabatic correction factor for the momentum.

Key factors, including *T<sup>s</sup>* and *Tc*, *LAI*, fc, *wc*/*hc*, and *hc*, are required as inputs for the *TSEB* model to parameterize the radiative and convective flux exchanges between soil/substrate and canopy. Other parameters related to micro-meteorological data are also needed to run the model. In the study conducted by Chirouze et al. [33] comparing different remote sensing *ET* models, results indicated that *TSEB* is a better model for *ET* estimation compared to others, being less sensitive to roughness parameters. This lack of sensitivity to roughness parameters was also recently verified for vineyards by Alfieri et al. [34]. The *TSEB* model has been extensively tested for years over agroecosystems [35–37], natural ecosystems [38,39], and wetlands [40,41].

*TSEB2T* was originally developed and evaluated by Kustas and Norman [42] using multiple thermal-IR radiometer viewing angles and was further refined and tested by Nieto et al. [2] applied to high resolution imagery from *sUAS* or other airborne sources. They found that *TSEB2T* gave better agreement with tower fluxes compared to other versions of *TSEB*, including *TSEB-PT* (Priestly-Taylor), *TSEB-DTD* (Dual-time-difference), and *TSEB-2T-DMS* (Data-mining sharpening of temperature). *TSEB-PT* is one version of the *TSEB* model that assumes a composite radiometric temperature (*Trad*) containing temperature contribution from the canopy and soil/substrate, which is typically provided by the radiometer. The decomposition of radiometric temperature (*Trad*) between plant canopy and soil/substrate is based on *fc*. *TSEB-DTD* is a further development of the *TSEB-PT* model described by Norman et al. [43]. The *TSEB-DTD* model is similar to the *TSEB-PT* model in that it divides the composite *Trad* into *T<sup>c</sup>* and *T<sup>s</sup>* . However, *TSEB-DTD* uses two observations of *Trad*: the first observation obtained 1.5 h after the sunrise (*Trad,*0) and the second one during the daytime (*Trad,*1). This version is less sensitive to errors in absolute radiometric surface temperature or the use of non-local air temperature observations. *TSEB-2T-DMS* partitions *T<sup>s</sup>* and *T<sup>c</sup>* using a data-mining fusion algorithm [44] to sharpen the original *LST* to be similar to the optical data, which would allow a better discrimination between *T<sup>s</sup>* and *Tc*.

The Nieto et al. [2] *TSEB2T* approach is a contextual *TSEB* that estimates *T<sup>s</sup>* and *T<sup>c</sup>* from composite *LST* imagery using the relationship between vegetation index (*VI*) and *LST* for extracting *T<sup>s</sup>* and *T<sup>c</sup>* within a spatial domain. *T<sup>s</sup>* and *T<sup>c</sup>* are calculated by averaging the temperature of pixels that are considered pure soil/substrate and pure canopy in a contextual spatial domain, namely, a two-dimensional plot of *LST* versus *VI*, such as Normalized Difference Vegetation Index (*NDVI*) (see Figure 1). That is, each pixel of the spatial domain is assigned based on *T<sup>c</sup>* and *T<sup>s</sup>* corresponding to the average temperature of the 0.6-m grids that are considered pure vegetation and bare soil, respectively. Both soil/substrate or canopy features are determined using *NDVI* threshold values (or any other vegetation index). The selection criterion for detecting the *NDVI* threshold of pure soil for bare soil interrows or, for most of the growing season, a soil senescent and cover crop stubble mixture

(substrate) (*NDVIs*) can be further supported by other sources such as *NDVI* value from a *NDVI-LAI* curve when *LAI* in the interrows is nearly zero. The pure vine canopy *NDVI* threshold (*NDVIv*) can be calculated as the mean value of pixels identified as pure vegetation in a binary (soil-vegetation) classification of a multispectral image. In cases of very dense vegetation where pure soil pixels do not exist or sparse vegetation lacking pure vegetation pixels inside the spatial domain, a linear fit between *LST* and *NDVI* can be developed where *T<sup>s</sup>* and *T<sup>c</sup>* can be estimated by previously defining the *NDVI* thresholds of canopy and bare soil (Figure 1).

#### *1.2. TSEB2T Main Inputs*

#### 1.2.1. Leaf Area Index (*LAI*)

*LAI* is one of the key inputs in *TSEB* influencing the computation of *ET* as leaves distribution is the driving factor in energy and mass exchange in this model. *LAI* is also difficult to acquire using ground-based leaf-scale measurements, due to the time-intensive effort required [45], complications using indirect methods in complex canopies, and lack of any spatial extent for mapping, even at the field scale [46]. Therefore, considerable efforts have been devoted to developing remote sensing approaches to estimate *LAI* [47].

Estimating spatial distribution of *LAI* is challenging in vineyards, with their rows of vines and interrows with little to no vegetation. A previous study conducted by Johnson [48] evaluated the *LAI-NDVI* relationship in vineyards using IKONOS satellite imagery with 1-m pixel resolution and comparing *NDVI* to ground-based *LAI* measurements. They concluded that *LAI* can be computed from *NDVI* using simple linear regression for the vineyard they studied planted with red grape in six blocks of different planting density, trellis, age, and cultivar. In addition, Johnson et al. [48] and Dobrowski et al. [49] showed that remotely sensed indices of soil and vegetation can be used to estimate *LAI*. However, a study by Fang [50] indicated that limitations exit when using vegetation indices (*VIs*) to describe the spatially distributed *LAI* due to sensitivity of the *LAI-VIs* relationship to vegetation type and substrate/soil type, and hence, will not be stable or applicable over large areas. Indeed, operational satellite retrievals of *LAI*, particularly for vineyards [51], have a level of uncertainty that could affect modeling fluxes using *TSEB*. Furthermore, canopy phenological properties (i.e., chlorophyll content and average leaf angle), along with other factors such as atmospheric scattering, soil reflectance, and the effects of mixed pixel due to a composite of soil and vegetation that changes with time and from one place to another, affect the accuracy of *LAI* estimation [47]. To improve the *LAI-VIs* relationships, numerous studies have been conducted to estimate *LAI* using statistical approaches. Artificial Neural Network (*ANN*) was very promising and is simple to use [50]; however, this method does not allow for standardization of the *LAI* estimation [52]. As described by Gonsamo and Pellikka [53], there is currently no standard or arbitrary characteristic parameters, specific vegetation types, or data sources can be employed for *LAI* estimation. Thus, researchers must develop custom models by considering the sensitivity of parameters to *LAI* within an expected range [53].

### 1.2.2. Canopy Height (*hc*)

The *h<sup>c</sup>* value is representative (mean) over the area of interest, but it can also be incorporated from spatial sources. An estimate of *h<sup>c</sup>* can be produced using high-resolution images from *sUAS* and other airborne sources processed with structure-from-motion (*SfM*) methods in Agisoft or Pix4D, among others, along with digital elevation models (*DEM*) and point clouds (*LiDAR*). The value of *h<sup>c</sup>* is required for the *TSEB2T* model to estimate surface aerodynamic roughness and radiation transmission in row crops and to calculate the foliage density, which are all required for the canopy wind attenuation model (Figure 2).

**Figure 2.** Schematic diagram for canopy *wc*/*hc* ratio. **Figure 2.** Schematic diagram for canopy *wc*/*hc* ratio.
