Effects of Surface Heterogeneity Due to Drip Irrigation on Scintillometer Estimates of Sensible, Latent Heat Fluxes and Evapotranspiration over Vineyards

: Accurate estimates of sensible (H) and latent (LE) heat fluxes and actual evapotranspiration (ET) are required for monitoring vegetation growth and improved agricultural water management. A large aperture scintillometer (LAS) was used to provide these estimates with the objective of quantifying the effects of surface heterogeneity due to soil moisture and vegetation growth variability. The study was conducted over drip ‐ irrigated vineyards located in a semi ‐ arid region in Albacete, Spain during summer 2007. Surface heterogeneity was characterized by integrating eddy covariance (EC) observations of H, LE and ET; land surface temperature (LST) and normalized difference vegetation index (NDVI) data from Landsat and MODIS sensors; LST from an infrared thermometer (IRT); a data fusion model; and a two ‐ source surface energy balance model. The EC observations showed 16% lack of closure during unstable atmospheric conditions and was corrected using the residual method. The comparison between the LAS and EC measurements of H, LE, and ET showed root mean square difference (RMSD) of 25 W m − 2 , 19 W m − 2 , and 0.41 mm day − 1 , respectively. LAS overestimated H and underestimated both LE and ET by 24 W m − 2 , 34 W m − 2 , and 0.36 mm day − 1 , respectively. The effects of soil moisture on LAS measurement of H was evaluated using the Bowen ratio, β . Discrepancies between H LAS and H EC were higher at β ≤ 0.5 but improved at 1 ≥ β > 0.5 and β > 1.0 with R 2 of 0.76, 0.78, and 0.82, respectively. Variable vineyard growth affected LAS performance as its footprints saw lower NDVI LAS compared to that of the EC (NDVI EC ) by ~0.022. Surface heterogeneity increased during wetter periods, as characterized by the LST– NDVI space and temperature vegetation dryness index (TVDI). As TVDI increased (decreased) during drier (wetter) conditions, the discrepancies between H LAS and H EC , as well as LE LAS and LE EC Re decreased (increased). Thresholds of TVDI of 0.3, 0.25, and 0.5 were identified, above which better agreements between LAS and EC estimates of H, LE, and ET, respectively, were obtained. These findings highlight the effectiveness and ability of LAS in monitoring vegetation growth over heterogonous areas with variable soil moisture, its potential use in supporting irrigation scheduling and agricultural water management over large regions.


Introduction
The increased pressure on water resources availability and use in arid and semi-arid regions has urged many agricultural water users to adopt water conservation practices, including drip irrigation [1][2][3][4]. In Spain and generally in the European Union (EU), drip irrigation is supported by many entities, strategies, and national plans, including the EU Water Framework Directive and National Irrigation Plan [3,5,6]. One of Spain's objectives of achieving water saving of 3000 Mm 3 /year has led to an exponential increase in drip-irrigated areas that exceeded 450% between 1989-2007, with another ~19% between 2007-2015 [3,[7][8][9]. Currently, drip-irrigated areas account for ~49% of the total irrigated areas (1.79 Mha) that are mostly (~37%) covered with vineyards [3].
Vineyards have the ability to adapt to semi-arid environments and are of economic and social significance. Adoption of drip irrigation in vineyards can help in improving water management and in addressing imposed (controlled) water stress conditions. However, drip irrigation can introduce spatially variable soil moisture and vegetation growing conditions that add some challenges in obtaining accurate measurements of sensible (H), latent (LE) heat fluxes and actual evapotranspiration (ET).
Routine and accurate measurements of H, LE, and ET are required for an improved understanding of vegetation growth and water use. Eddy covariance (EC) and scintillometer systems are among the best methods that can provide such measurements, with both having pros and cons [10,11]. Some challenges can arise with the use of EC systems because they are relatively expensive, labor-intensive, requiring well-trained personnel, and provide local-scale observations [10,12]. Such local-scale observations do not adequately represent the spatial variability of the underlying surface. In many cases, large-scale observations over spatially heterogeneous surfaces are necessary. In these cases, a network of EC systems has been successfully used, which, however, adds some feasibility and management challenges.
Alternatively, advances in the scintillation method made it possible to address some of these challenges. Scintillometers can provide large-scale observations, require minimal data processing, and, technically, are relatively easy to maintain compared with a network of EC systems [13][14][15]. Scintillometers provide path-averaged measurements of H over distances of about 500 m up to 10 km. They can also provide area-averaged estimates of LE when combined with measurements of net radiation (Rn) and soil heat fluxes (G) [15][16][17]. Such large-scale observations are needed to provide more compatible observations for ground-truthing of remote sensing-based model estimates of H, LE, and ET [15,18]. Hence, scintillometers have advantages over EC systems in providing routine, efficient, and spatially representative large-scale observations that can enhance water managers' abilities to accurately monitor vegetation growth and water use [19].
Scintillometer observations are typically evaluated against those from EC systems. These evaluations often show a wide range of discrepancies depending on the underlying surface heterogeneity [14,15,17,20]. Some efforts have been made to identify the sources of these discrepancies, but there is still a need for additional evaluations. One of the main reasons for these discrepancies is attributed to differences in their represented spatial scales and surface heterogeneity.
The effects of some surface variables, including vegetation height, topography, and land cover types, on scintillometer applications were partially considered [14][15][16]21]. However, the effects of the most important variables-soil moisture and vegetation growth-have not been appropriately addressed [14,20,22]. These two variables play key roles in surface energy balance processes. Soil moisture affects partitioning of available energy (Rn − G) to H and LE. Vegetation growth affects surface roughness and ET. Some studies have accounted for surface variability due only to soil moisture, while the effects of vegetation growth on surface variability have been minimally considered or ignored. These two variables have unique combined effects in characterizing surface heterogeneity, and hence, they need to be accounted for accordingly. For such, remote sensing can effectively and appropriately characterize surface heterogeneity.
A review of recent scintillometer applications suggested that remote sensing has been minimally utilized in addressing surface heterogeneity and its effects on scintillometry. For example, [20] and [22] used only land surface temperature (LST) over surface-irrigated olive trees and natural vegetation cover, respectively, to infer the effects of soil moisture variability on scintillometer-based estimates of H. On the other hand, vegetation growth has been extensively monitored using normalized difference vegetation index (NDVI) based on surface reflectance (i.e., optical remote sensing), as in [23][24][25][26]. However, the integration of surface variability, as characterized by using remote sensing-based indices such as NDVI to explain the performance of scintillometers, has been insufficiently studied.
The combined use of LST and NDVI can provide complimentary information about surface conditions [27]. Specifically, LST and NDVI can be graphically combined to form what is called LST-NDVI space that mostly takes the shape of a triangle diagram [27,28]. This empirical triangle concept was effectively used to describe and characterize the combined effects of soil moisture and vegetation growth over most surfaces [27,28]. In an LST-NDVI space, drier (wetter) surface conditions can be represented by lower (higher) NDVI, higher (lower) LST, and vice-versa that correspond to sparsely or stressed (dense or unstressed) vegetated surfaces. An LST-NDVI space can also be numerically represented in the form of a temperature vegetation dryness index (TVDI) that allows to quantitatively characterize soil water stress, ET, and surface and root zone soil moisture [27][28][29].
Generally, scintillometer behavior relative to irrigation events has been previously qualitatively evaluated, for example, in [16,20]. However, there is still a need to quantitatively characterize the combined effects of soil moisture and vegetation growth variability to appropriately evaluate surface heterogeneity due to irrigation practices (e.g., drip irrigation) and its effects on scintillometer measurements. To address this need, this study uniquely introduced the combined use of TVDI and scintillometry.
The goal of this study was to evaluate the effects of irrigation-induced surface heterogeneity in terms of soil moisture and vegetation growth variability on large aperture scintillometer (LAS) measurements of H, LE, and ET for efficient monitoring of crop water consumption and improved irrigation water management. To achieve this goal, this study uniquely integrated an EC system, thermal and optical remote sensing data, a data fusion model, and a surface energy balance model. The objectives of the study were to: Compare EC measurements with LAS estimates; use the spatial temporal adaptive reflectance fusion model (STARFM) to develop (i.e., fuse) a daily time series of NDVI based on Landsat and MODIS images; estimates of Rn, G, H, and LE using the two-source energy balance (TSEB) model; and characterize surface heterogeneity using LST-NDVI space and TVDI. The study highlights the performance of LAS during different soil moisture conditions by providing some explanations about when and why LAS and EC measurements can closely agree.

Study Area
The experiment was carried out during the summer and fall of 2007 (1 August-31 October) in Central Spain (39°17′ N, 1°59′ W, 700 m above mean sea level-MSL) in the Albacete province, over drip-irrigated, row-oriented vineyards grown on vertical shoot positioned trellis [23]. H was measured simultaneously using an EC system and a LAS BLS900 (Scintec AG, Rottenburg, Germany). Other surface energy balance fluxes, including LE, Rn, and G, were also monitored throughout the growing season.
The study area was about 18.5 ha planted with four different 7-year-old varieties of grapes. The site was divided into 17 irrigation fields: The Cencibel variety was predominant, but Cabernet Sauvignon, Shiraz, and Merlot were also cultivated in smaller proportions. During the field campaign, other crops, including winter cereals and irrigated fruit trees, and a pine forest surrounded the site (Figure 1). The vineyards were pruned during the dormant and growing seasons to adjust the vegetation development. No-tillage or cultivation was conducted during the growing season and the application of herbicides prevented the growth of weeds. The vine spacing was 1.5 m and the row spacing 3.0 m, which were oriented approximately in the north-south direction. The drip lines were placed on the trellis below the plants, with drippers spaced 1.0 m apart. Vineyard production in 2007 was 14,000 kg ha −1 , around 7 kg/tree, which was [23] similar to previous year yields, and no infection or diseases were detected in the plants during the growing season. Figure 1. Aerial photograph of the study site with the 17 drip-irrigated vineyard fields, along with the location of the large aperture scintillometer (LAS) transmitter and receiver and eddy covariance (EC) tower. The numbers inside each polygon represent the field numbers. The LAS transmitter was located in Field 9 and the receiver was located in Field 16. The EC tower was located within Field 7 along the LAS path. The spaces between Fields 2 and 3, as well as other similar ones (~5-10 m), were for field accessibility purposes.
Irrigation was performed according to local practices, with a general rule of applying 22 mm of water to each field approximately every 12 days based on farmers' reports. However, observed irrigations over Field 7 showed that there is some variability in irrigation dates, with a total of 6 irrigations with 22 mm and a seventh one with about 11 mm, with an average of total applied irrigation depth of 143 mm throughout the growing season. The irrigation schedule over the other fields also showed similar variability (Appendix A- Figure A1). Unnecessary vegetation development increases crop water consumption, reduces soil water storage, and increases irrigation requirements. Therefore, the shoots were pruned in the experimental field around day of year (DOY) 187 and no new leaves or stems were observed after that date. So, the fractional vegetation cover (Fc) remained stable at 30%, as shown in Figure 2

LAS Measurements
The Boundary Layer Scintillometer BLS900 model from Scintec AG Rotternburg, Germany was used in this study. The LAS has an aperture diameter D = 0.15 and operates at a wavelength of 880 nm. Based on wind direction analysis, the LAS was deployed so that the wind passes dominantly perpendicular through its beam. The distance between the LAS transmitter and receiver was 640 m ( Figure 1). The LAS path was relatively horizontal, as the transmitter and receiver were at 2.8 m and 3.5 m, respectively, above the ground. The measurement was sampled at 1 Hz and averaged over 1min intervals. Air temperature (Ta) and atmospheric pressure (P) measurements were acquired using separate sensors that were integrated with the LAS system.

EC Measurements
The EC system was installed in the central part of the plot (Figure 1). The EC system consisted of a sonic anemometer CSAT-3 (Campbell Sci. Inst., Shepshed, Leicestershire LE12 9GX, UK), an infrared gas analyzer LI-7500 open system (LI-COR Inc., Lincoln, NE, USA) that records at 10 Hz (0.1 s), a humidity and temperature sensor HMP45C (Vaisala, Vantaa, Finland), and a datalogger CR5000 (Campbell Sci. Inst., UK), in which air density correction and some statistics were computed. The data were incorporated into high-frequency tables of fluxes every 30 min. The post-process was performed using the software TK2 [30], which mainly corrects the effects of separation between sensors, frequency attenuations by time averages, sensor length, and spurious values [23]. The sensors were located at 3 m above the ground-oriented northwest in the direction of prevailing winds ( Figure 1). Rn (W m −2 ) was measured using a net radiometer CNR1 (Kipp & Zonen, Delft, The Netherlands), located 4.5 m above ground. The row orientation of the vineyards' trellis introduced sunlit and shading patterns that can affect G. To account for this effect, G (W m −2 ) was measured at three positions spatially representative of the vineyard soil (i.e., the area covered by plants, at 0.75 m and 1.5 m from the plant row). G was calculated by correcting the heat fluxes measured by the flux plates by accounting for the effects of soil moisture content [31]. In each of the three G measurement positions, a Heat Fux Plate (HFP1) (Hukseflux, Delft, The Netherlands) was placed and buried at 8 cm deep along with two chrome-constantan thermocouples (TCAV, Type E) at 2 and 4 cm depths, together with a volumetric moisture reflectometer CS616 (Campbell Sci. Inst., Logan, UT, USA). A thermal infrared thermometer (IRT) was used to measure canopy and soil surface temperatures. A set of three Everest IRT (Everest Interscience Inc., Chino Hills, CA, USA) sensors accuracy of 0.3 °C. One IRT sensor was pointed at the vineyards to measure canopy temperature and the two other ones were pointed at sunlit and shaded soils to provide an average soil surface temperature. This set of canopy and soil temperatures were used to estimate LST that was used to model surface energy balance fluxes, as described in Section 3.4.

Surface Energy Balance Lack of Closure
Lack of surface energy balance fluxes closure measured by EC systems is a known issue that yet remains unresolved [32], and it appears in terms of discrepancies between turbulent heat fluxes (H + LE) and available energy (Rn − G). The lack of closure in this study was evaluated using two wellknown methods-the Bowen ratio (BR) and the Residual (Re) [32,33]. Generally, the BR method distributes the lack of closure between H and LE based on β (H/LE), and the Re method assumes this error is due to an underestimation in LE (LE = Rn − G − H) [32]. The two methods (i.e., BR and Re) were considered here because there is no consensus on a preferred one over the other. Lack of surface energy balance closure was addressed using daytime 8:30 a.m.-4:30 p.m. data during unstable atmospheric conditions.

Landsat Images
Four pairs of Landsat 5 images that consisted of LST and surface reflectance were used to characterize surface heterogeneity over the study site during the measurement period. Top of the atmosphere temperature images were atmospherically corrected to obtain at-surface LST using the radiative transfer model MODTRAN 4 [34]. These images were acquired during clear sky conditions on 4 and 27 August, and 5 and 28 September (Figure 2). Both datasets have a 30-m spatial resolution that allowed to account for field scale variability. Surface reflectance in near-infrared (NIR) and red bands were used to calculate NDVI as (NIR − Red)/(NIR + Red).
Generally, over sparsely vegetated areas with dry soil surfaces, NDVI ranges between 0-1, with low and high values correspond to bare soils and full cover dense canopies, respectively. It can reach an asymptotic value of about 0.9 over most dense canopy covers. NDVI is not strongly sensitive and indirectly related to soil moisture in the root zone. It is directly related to and can be used to describe vegetation growth conditions and some plant biophysical properties over time. It provides a delayed response in terms of vegetation growth due to water stress [28]. A well-watered transpiring vegetation has higher NDVI values compared to water-stressed less transpiring one [27,35]. So, vineyard growth can be detected with NDVI spatially and temporally. On the other hand, LST is more sensitive to soil moisture and water-stressed conditions. As a biophysical response, a surface with well transpiring vegetation has lower LST compared to one with water-stressed conditions [27,29,35,36]. The Landsat data were used to develop the LST-NDVI space and TVDI.

MODIS Images
The study used a set of 56 images of MODIS NDVI for the period between 4 August and 28 September, 2007. MODIS data have a 250-m spatial resolution that does not allow to identify field scale surface heterogeneity. Therefore, the data were used in combination with those from Landsat 5 based on a data fusion approach to develop statistically representative Landsat-like NDVI images. This process allowed filling-in the gap between Landsat overpass dates using the spatial temporal adaptive reflectance fusion model (STARFM)-a data fusion approach developed by [37]. A brief description of STARFM is provided in Section 3.3. This NDVI time series helped in characterizing vineyards growth spatial variability over time.

Root Zone Soil Moisture Content
Soil moisture content measurements were acquired in Field 7 between two rows adjacent to the EC system at multiple layers below the surface at 10, 20, 30 40, 50, and 80 cm using an EnviroScan soil water capacitance probe. Within the 3.0 m row spacing, two sets of probes were installedadjacent to the vine and near the emitters at each row. Another two sets of probes were installed at 75 cm away from each row and one set in the middle of the row (i.e., at 1.5 m from the row). The probes adjacent to, and those at 75 cm from, the vine showed a consistent soil moisture response with the irrigation events. The sensors that where located in the middle of the row, as well as those located at and below a depth of 50 cm, did not show any relevant response. The probes at shallow depths showed some response relative to significant rainfall events.
It should be noted that in this study, the probes were not calibrated to site-specific soil-water physical properties, including the range of available soil moisture between field capacity and wilting point. Measured soil hydrologic properties, including field capacity and permanent wilting point, were 0.29 cm 3 cm −3 and 0.14 cm 3 cm −3 , respectively. Detailed soil properties were provided by [23], in which they conducted root zone water balance analysis to estimate ET. While the absolute accuracy of these measurements was not necessarily needed in this study, the evolution of soil moisture response to the irrigation events was certainly important and useful to explain the behavior of LAS and EC measurements of H and LE.

Estimates of H Using LAS
A scintillometer is a device that operates at a near-infrared electromagnetic wave to measure turbulent intensity fluctuations and present them in the form of a path-averaged structure parameter of the refraction index of air, . To this purpose, an electromagnetic radiation is transmitted over a (usually) horizontal path, and the corresponding intensity fluctuations are analyzed at the receiver. Both temperature and humidity fluctuations, and , respectively, give rise to . The relationship between and can be described by Equation (1) [13,38,39] as: where and are in K 2 m −2/3 and m −2/3 , respectively, P the atmospheric pressure (Pa), γ the refractive index coefficient for air (0.78 × 10 −6 K Pa −1 ), and Ta the air temperature (°K). Equation (1) is a simplified version of Equation (A1) (Appendix B). In some cases, β can be obtained from the EC flux measurements. In this study, and to obtain LAS estimates of H (HLAS) independent of β from EC measurements, β was estimated iteratively using available (Rn − G) following the approach described by [16,40]. With the application of the Monin-Obukhov similarity theory (MOST), can be related to the temperature scale, * T , by a universal function of stability parameter ⁄ as: where is the LAS effective height (m) and L is the Monin-Obukhov length. For unstable conditions (L < 0) [40,41]: The effective height, , takes into account the effects of stability conditions and slanted paths is the bell-shaped weighting function describing the contribution of to the total LAS signal at each point along the normalized path u [21], and is the variable LAS beam height along the path ( Figure 1). For slanted paths, can be described relative to the maximum and minimum heights if the transmitter and/or the receiver and , respectively, as 1 1 ⁄ • . The Monin-Obukhov length, L, can be estimated as: with k = 0.4 is the von Karman constant, the gravity (9.81 m•s −2 ), and the friction velocity, * , can be estimated using the standard Businger-Dyer flux profile [42] as: where U is the wind speed, ψm the stability correction function for the momentum transfer was calculated following [43] and [44] as and is the roughness length. Both d and were calculated as a function of crop height ℎ as 2 3 ⁄ • ℎ and respectively [45]. HLAS can be estimated as: where is the air density and is the specific heat of air at constant pressure.

LST-NDVI Space, TVDI, and Spatial Heterogeneity
The LST-NDVI space can uniquely describe soil moisture variability and vegetation conditions over most surfaces [46]. As shown in Figure A2 (Appendix A), there is a negative correlation between LST and NDVI. In other words, as LST increases (decreases), NDVI decreases (increases). A plot of all pixel values of LST and NDVI for a given scene can reveal a triangle-shaped formulation. with upper and lower bounds representing dry and wet edges, respectively. Partially vegetated surfaces fall in the space between these two edges. A set of isolines can be developed with varying slopes (Appendix A- Figure A2) that describe soil moisture variability, H, and LE [27,28,46]. In some cases, the LST-NDVI space is represented by a trapezoid to avoid some uncertainties that appear near the lower right angle of the triangle as the isolines are set very close to each other [29]. The LST-NDVI space has been used to qualitatively describe vegetation water stress and soil moisture variability [27,29,47,48].
Quantitatively, LST and NDVI can empirically be combined to TVDI that ranges between 0-1 for wettest to driest surfaces. The TVDI concept has been utilized to directly estimate surface and root zone soil moisture content [28,46,49]. In this study, the TVDI was used as an indicator of surface heterogeneity due to soil moisture, and to partially explain the corresponding behavior of LAS measurements. The TVDI can be estimated as: (7) where is the minimum LST and the maximum LST that can be estimated as a function of NDVI as • , where a and b are linear regression coefficients of the dry edge limit line. Based on the isolines, TVDI is simply the ratio between A to B (Appendix A- Figure A2). Two sets of TVDI were calculated, one was based on LST and NDVI from the four Landsat images. The second set was based on LST from the IRT and the NDVI from STARFM over Field 7. A consistent set of LST from the IRT and NDVI from STARFM was identified based on the time of the day that would correspond to Landsat typical overpass of 12:40 p.m. The TVDI was developed using LST-Ta to account for the effects of air temperature, as recommended by [28], and the corresponding dry edge line was developed following ( Figure A1) [27,28,46].

Data Fusion
The STARFM is a data fusion approach that was developed to blend daily MODIS with the 16day Landsat reflectance data. Using statistical methods, STARFM combines spectral similarities in Landsat and MODIS image pairs with temporal changes from the 250-m MODIS data to fill-in the gap between Landsat overpass dates to provide a consistent time series of Landsat-like data at 30-m spatial resolution. The STARFM has been used in many agricultural water management applications.
It has been successfully applied with the TSEB model to provide field scale gap-filled time series of ET [50][51][52][53]. In this analysis, four pairs of Landsat and MODIS were used to develop statistically Landsat-like NDVI images from MODIS. This NDVI time series allowed to characterize vineyard growth variability within each field and within the LAS and EC footprints.

The TSEB Model
The thermal based two-source energy balance (TSEB) model [54,55] was used to provide independent estimates of H, LE, and Rn − G to evaluate HLAS. The TSEB model was applied using two different datasets-Landsat images and IRT point-based observations of LST. The TSEB model was chosen in this analysis due to its ability to provide accurate estimates of H and LE over a wide range of heterogeneous surfaces (e.g., row cropped fields) compared to other similar type of models. The TSEB model was successfully applied over heterogeneity surfaces similar to this site, including vineyard fields in California [35,50,[56][57][58]. Based on a two-source concept, the TSEB model applies the surface energy balance equation over surface components (i.e., canopy and bare soil) separately and combines the fluxes at an air-canopy interface. The main assumption of the model is that LST can be decomposed based on the fraction of cover (fc) [55] as: with is the LST or sometimes referred to as radiometric surface temperature, and and are canopy and bare soil temperatures (°K). The fraction of cover fc viewed at nadir can be estimated as 1 .

•
, where LAI is the leaf area index. A brief description of the TSEB model and its application over canopy and soil components is provided in Appendix C [54].

Footprint Model
The source areas (i.e., flux footprints) of the EC and LAS systems captured associated surface heterogeneity that contributed to measured fluxes. Flux footprints of each system partially integrated a number of fields or sub-fields ( Figure 1). A footprint analysis was carried out to estimate these source areas for each system. The size and extent of a footprint depend on many factors, including measurement height, atmospheric stability, wind speed and direction, and aerodynamic roughness length. The EC system footprints were estimated using the footprint model by [59,60]. For the LAS footprints, [59,60] model was combined with a superposition approach, as described in [14,15,20]. More details about the footprint analysis are provided in Appendix D.

Observed Surface Energy Balance Closure
Observed HEC + LEEC showed a closure ratio (CR) of 84% of Rn − G (Figure 3). In other words, an underestimation of 16% in available energy was either unmeasured or unexplained, with a root mean square difference (RMSD) and R 2 of 54 W m −2 and 0.69, respectively. Adjusted sensible heat flux using the BR method (HEC β) showed a mean value of 111 W m −2 (standard deviation (SD) of 56 W m −2 ), with a 15% increase compared to that of raw HEC of 95 W m −2 (SD of 46 W m −2 ). Similarly, adjusted latent heat flux using the Re (LEEC Re) and BR (LEEC β) methods resulted in mean values of 176 (SD of 89 W m −2 ) and 159 W m −2 (SD of 78 W m −2 ), with about 30% and 15% increase, respectively, compared to that of LEEC of 138 W m −2 (SD of 68 W m −2 ).

Comparison Between HLAS and HEC
The mean and SD statistics of LAS and EC measurements of H, HLAS, and HEC and HEC β, respectively, showed that the mean of HLAS of 118 W m −2 (SD of 52 W m −2 ) was ~26% higher than that of HEC of 94 W m −2 (SD of 56 W m −2 ) and 6% higher than that of HEC β of 111 W m −2 (SD of 46 W m −2 ). These statistics indicate that regardless of the method used to address lack of closure, HLAS overestimated those obtained by the EC system. In other words, this result suggests that the flux footprint "seen" by the LAS experienced higher H with relatively drier surface conditions compared to that "seen" by the EC, as explained in more details in Section 5. The comparison between HLAS and HEC showed some discrepancies, as indicated by an RMSD of 25 W m −2 (R 2 of 0.77) with a narrow dispersion around the 1:1 line within a 95% confidence interval (Table 1 and Figure 4). A slightly lower agreement was observed between HLAS and HEC β, with an RMSD of 28 W m −2 (R 2 of 0.72) and a wider dispersion around the 1:1 line. Also, the slope of the regression lines of 0.99 and 0.79 indicate a better agreement between HLAS and HEC compared to HLAS and HEC β, respectively. These comparisons indicate that the use of the Re method in this study to address lack of closure provided a better agreement compared to the BR method. Therefore, the analysis was completed using the Re method.

HLAS and β
The effects of soil moisture temporal variability on HLAS were evaluated using β. Generally, higher soil moisture content can lead to an increase in LE (and ET), decrease in H, and, consequently, lower β values [35]. In this study, β values were divided into three categories, including β ≤ 0.5, 1 ≥ β > 0.5, and β > 1. The first two categories generally indicate that LE (and ET) > H. The third β category indicates that H > LE. A comparison between HLAS and HEC for each β category is shown in Figure 5. The data in the third β category suggest that the site experienced relatively drier conditions compared to the other two categories. At β ≤ 0.5, HLAS overestimated HEC by a mean bias difference (MBD) of 34 W m −2 . At 1 ≥ β > 0.5, the MBD decreased to 27 W m −2 , and at β > 1, HLAS underestimated HEC by an MBD of −14 W m −2 . In other words, as β increased (i.e., drier surface conditions became dominant), the agreement between HLAS and HEC improved as the MBD consistently decreased. Also, R 2 of 0.76, 0.78, and 0.82 for β < 0.5, 1 ≥ β > 0.5, and β > 1, respectively, show consistently improved agreement between HLAS and HEC as β increased.

LST-NDVI Space and Spatial Heterogeneity
This section first shows an evaluation of LST, NDVI, and an LST-NDVI space based on the four Landsat images. Second, the LST-NDVI space was then used to develop a TVDI relationship based on the Landsat data. Third, a TVDI time series was developed based on a relatively longer dataset of LST from the IRT and NDVI from STARFM. Fourth, the TVDI time series was used as a surrogate for soil moisture to evaluate the HLAS and LELAS.

NDVI Variability
The average NDVI over each field for all four Landsat images was calculated and is shown in Figure 6. The NDVI indicates that each field experienced relatively variable growing conditions. For example, the average NDVI of all dates over Fields 2, 6, and 9 consistently show lower NDVI of 0.32, 0.30, and 0.31, compared to those over Fields 7 and 8 of 0.36 and 0.39, respectively. The corresponding average NDVI that represented, or "as seen" by, the LAS and EC footprints is shown on the x-axis in Figure 6 as points 18 and 19, and referred to herein as NDVILAS and NDVIEC, respectively. The NDVI that represented (or as seen by) the EC system was estimated by overlying (superimposing) the EC footprints over the NDVI images. Using the footprint values (which are in dimensionless units), a weighted sum of NDVI values within the footprint boundaries was calculated to estimate NDVIEC for each image data and time. A similar process regarding the calculation of NDVILAS was followed. The results show that both NDVILAS and NDVIEC of 0.35 and 0.37, respectively, were relatively higher than those of over Fields 2, 6, and 9. This pair of four NDVI values shows that NDVILAS was consistently, but slightly, lower than NDVIEC ( Figure 6).  A time series of footprint-integrated pairs of NDVILAS and NDVIEC was developed using the STARFM-based NDVI maps (Figure 7). Images with cloudy sky conditions were excluded (11 images out of 54). The results show that NDVILAS and NDVIEC were closely correlated (R 2 = 0.90), and the average NDVILAS was ~0.022 lower than NDVIEC.

LST-NDVI Space
Using the triangle concept ( Figure A2), regression lines (or isolines) were developed to represent the data for each of the four Landsat image dates ( Figure 8). The negative slope shown by all regression lines indicates that during each day, as LST increased, NDVI decreased. The standard deviation of NDVI ( ) for each day (or isoline) varied from low to high, from 0.019, 0.020, 0.021, and 0.033 from the upper-to the lower-most isoline, respectively. The suggests that there was a relatively larger spatial vineyard growth variability during wetter conditions (i.e., lower-most isoline) compared to drier ones. The mean LST ( ) varied from high to low from the upper-to the lower-most isoline from 315.9, 313.2, 307.9, and 302.2 °K, respectively. The also indicates that the surface exhibited higher LST (i.e., warmer) during drier conditions compared to wetter ones (i.e., cooler). This qualitative evaluation of LST-NDVI space was necessary to develop a TVDI time series based on the IRT data.

Estimates of TVDI
As an initial evaluation, the obtained TVDI, which ranged from 0-1, was plotted against HEC (not shown) to make sure that it followed the known behavior, as described by [27,28,46,49]. A positive correlation between TVDI and HEC was obtained. The lowest average TVDI (about 0.20) was on 28 September and corresponded to the lowest average LST (~302.2 °K). The observed behavior of the TVDI in this study agreed with the fact that higher TVDI corresponds to drier conditions [28,46,49], higher LST − Ta [29], and higher H [28,29]. Estimates of HTSEB The TSEB model was used to estimate HTSEB and (Rn − G)TSEB that can be used to evaluate both HLAS-HEC and HLAS-HTSEB against TVDI and to calculate LELAS, respectively. The TSEB model was applied using two LST datasets from the Landsat images and the IRT. Estimates of H based on the Landsat data (or remote sensing) (HTSEB RS) and those based on the IRT (HTSEB IRT) are shown in Figure  9. Two different comparisons were conducted-HEC against HTSEB IRT and HLAS against HTSEB IRT. The comparison between HEC and HTSEB IRT showed an RMSD and R 2 of 35 W m −2 and 0.71, respectively, with a narrow scatter distribution around the 1:1 line. Slightly lower agreement was obtained between HLAS and HTSEB IRT, with an RMSD and R 2 of 51 W m −2 and 0.64, respectively. One of the main reasons for the increased discrepancies between HLAS and HTSEB IRT was attributed to differences in the spatial representation of HTSEB IRT (at a point) to that of HLAS. Estimates of HTSEB RS did not have enough data points that could be meaningfully evaluated. However, a visual inspection of Figure 9 showed that the footprint-integrated HEC and HLAS were closely correlated with HTSEB RS. Overall, these TSEB model results are consistent with previous studies [14,35,50,56,61,62] and allowed the use of HTSEB IRT to characterize surface conditions.

HLAS and TVDI
Estimates of HLAS, HEC, and HTSEB IRT can be considered independent of each other-a condition that allowed to adequately evaluate soil moisture and vegetation growth variability, as characterized by the TVDI. The comparisons between HLAS − HEC and HLAS − HTSEB IRT against TVDI are shown in Figure 10. The results show a negative correlation (R 2 = 0.55) between HLAS − HEC and TVDI, which indicates that as the discrepancies between HLAS and HEC decreased, the TVDI increased. A relatively higher correlation (R 2 = 0.66) was obtained between HLAS − HTSEB IRT and TVDI based on a smaller sample size. This comparison also provides a similar indication of decreased discrepancies between HLAS and HTSEB IRT with increased TVDI. In other words, both comparisons indicate that the agreement between HLAS and both HEC and HTSEB IRT improved during higher TVDI or specifically drier surface conditions. This consistent behavior of HLAS and TVDI suggests that TVDI can effectively be used to evaluate the effects of surface heterogeneity due to irrigation on LAS performance.   (Figure 11).

LELAS and TVDI
The differences LELAS TSEB IRT − LEEC Re and LELAS EC − LETSEB IRT were compared with TVDI to evaluate the effects of soil moisture and vegetation growth variability on LELAS. The obtained results show that both differences were negatively correlated with TVDI with first comparison (i.e., LELAS TSEB IRT − LEEC Re versus TVDI) resulted in R 2 = 0.72, and the second one (i.e., LELAS EC − LETSEB IRT versus TVDI) indicated an R 2 = 0.57. These negative correlations suggest that as the two differences increased, the TVDI decreased (as soil moisture increased). These comparisons were based on a much smaller sample sizes due to the need of using consistent datasets of LELAS TSEB IRT, LEEC Re, LELAS EC, and LETSEB IRT. Therefore, they were used to qualitatively evaluate LELAS and TVDI behavior. LELAS EC was used to provide a quantitative evaluation of TVDI, as discussed in Section 5.7.

Estimates of ET
Estimates of ETLAS and ETEC Re based on LELAS and LEEC Re, respectively (Figure 11), show a good agreement, with an RMSD of 0.41 mm day −1 . An underestimation of ETLAS to ETEC Re was noticed with an MBD of −0.36 mm day −1 . This result was consistent with that obtained in Section 4.3, which indicated an overestimation of HLAS to HEC and, consequently, an underestimation of LELAS to LEEC Re. The obtained mean values for ETLAS and ETEC Re were 2.1 and 2.5 mm day −1 , respectively. A better LAS performance was visually observed at ETLAS ≤ 2 mm day −1 with increased discrepancies above this value.

Lack of Closure
Observed lack of closure of EC measurements was within the typical acceptable range of 15-20% over agricultural areas [12,32,40,56]. Closure ratios (CRs) of 0.84 and 0.86 were reported over heterogonous vineyards in California [56] and homogeneous rainfed soybean in Iowa [63][64][65]. Also, the lack of closure in this study had minimal effects on the discrepancies between HLAS and HEC. Generally, lack of closure is partly attributed to mismatch in spatial scale of measurements (i.e., Rn − G versus H + LE) and the underlying surface heterogeneity [12,32,65,66]. As indicated by [22], a combination of CR ≤ 0.80 and surface heterogeneity can increase these discrepancies. This combination highlighted the need to evaluate the effects of surface heterogeneity on LAS in this study.

HLAS and HEC
The obtained RMSD between HLAS and HEC was within acceptable levels of accuracy. Reported typical levels of accuracy of H measurements for LAS and EC were within 30 W m −2 [40] and 50 W m −2 [12,67,68], respectively. Also, [20] indicated that acceptable differences between HLAS and HEC were within 50 W m −2 .
[20] compared HLAS against HEC over flood-irrigated olive orchards and obtained an RMSD of 36, 29, and 63 W m −2 before, after, and during irrigation events, respectively, assuming heterogeneous soil moisture during irrigation and homogeneous otherwise. Also, [16] obtained an RMSD of 63 W m −2 over drip-irrigated orange orchards.
Previous studies evaluated LAS performance mainly relative to three factors: Vegetation cover uniformity, placement of EC and LAS considering their footprints, and soil moisture effects on surface heterogeneity [16,20,22,69]. In this study, the first two factors had relatively minimal effects compared to the latter. The site had homogenous cover dominantly vegetated with one vineyard variety with a constant hc and Fc during the measurement period [23]. A footprint analysis showed that the EC footprints lied mostly within Field 7 (and Field 3), and that of the LAS partially integrated nine fields, including Fields 2, 3, 4, 5, 6, 7, 8, 13, and 15. Field 7 accounted for 27% of the LAS footprints ( Figure 12). The EC footprints accounted for ~20% of, and most importantly lied entirely within, that of the LAS-such setup is considered the most ideal placement of LAS and EC systems to obtain comparable measurements [22,69]. Thus, this evaluation suggests that surface heterogeneity due to soil moisture and vegetation growth (i.e., vigor and greenness) variability can mainly explain the discrepancies between HLAS of HEC.

HLAS and β
Considering soil moisture effects, the use of β in Equation (1) can also partly explain some of the discrepancies between HLAS and HEC based on the humidity correction factor and the corresponding assumptions. The humidity correction factor [ 1 0.03⁄ ] partly affects the accuracy of and HLAS [38,39]. The humidity correction factor is less than 10% for of β ≥ 0.6 [21]. Thus, it has minimal effects on HLAS accuracy at larger β and no need to use it at β > 1 [14,21,40]. On the other hand, the uncertainty in HLAS can reach up to 15% and 48% at β = 0.3 and β < 0.3, respectively [40]. In other words, decreased uncertainty in HLAS during drier conditions (i.e., H > LE and β > 1) can lead to improved agreement with HEC compared to periods with increased soil moisture (H < LE and β < 1).
Moreover, regardless of the correction factor, its assumptions can still lead to increased discrepancies during wet conditions. The obtained results indicated that even when accounting for this factor, the discrepancies between HLAS and HEC at β ≤ 0.5 increased. Equation (1) is a simplified form of Equation (A1) (Appendix B) that relates , , and [39], assuming that the correlation between T and q, Rtq = +1. This assumption removes the direct dependency of on , and hence reduces its contribution to [ 1 0.03⁄ ]. Many studies argued the applicability of this assumption over a wide range of moisture conditions. In some cases, Rtq < 1, especially during periods with increased soil moisture content [20,38,70]. [20] found that Rtq ≤ 0.75, during which Rwt and Rwq did not follow each other, resulting in ~4% overestimation by HLAS over flood-irrigated olive orchards. Unfortunately, in this study, there were no data to evaluate the Rtq, but the study by [20] was used as a guidance, as they indicated that Rtq ≈ 0.6 for β ≤ 0.5.
The violation of this assumption, which led to overestimation of HLAS to HEC, is typical during transitions from dry to wet periods, and largest during low H values that occur during nearly neutral unstable conditions (i.e., low β values), as indicated by [71]. During such periods, the near surface layer (i.e., internal adapted layer) is not fully developed compared to very unstable conditions (i.e., during high β values). In this study, a better agreement between HLAS and HEC was clearly obtained during drier conditions (i.e., at β > 1) and a reasonable one was obtained at 1 ≥ β > 0.5. These results are consistent with previous findings by [20,40,70,71]

NDVI Variability
The NDVI provided a good indication of how healthy and variable the vineyards were growing spatially and temporally. All vineyard fields showed a consistent NDVI spatial variability over time. Such consistent NDVI behavior over time was indicative of soil temporal variability. NDVI spatial variability indicated that each vineyard field experienced physiologically variable growth compared to other fields. Spatially variable soil properties can partly explain this variable growth, since all fields were dominantly planted with one variety and there was no reported water or nutrient stress conditions. The soils in the study site were classified as Inceptisoil with a sand, clay, gravel texture [23]. These types of soils are known for their variable productivity, including texture and waterholding capacity. A similar behavior was observed over vineyards grown in California, US in fields with variable soil physical properties that caused variable growing conditions and lower yields [56].
Based on NDVI, the vineyards within the EC footprints showed relatively better growth compared to that within the LAS. This can partially explain the overestimation of HLAS to HEC and it is also consistent with the basic concepts of near surface flux exchange [55,72,73]. An evaluation of NDVILAS − NDVIEC against HLAS − HEC showed that HLAS − HEC slightly increased as NDVILAS − NDVIEC decreased ( Figure 7b) with a low correlation (R 2 = 0.17). This low correlation was partly because the study was conducted during vineyard full cover period with relatively constant Fc and minimal variation in NDVI ( Figure 2). A more effective approach in these cases is to monitor vineyards throughout the growing season. NDVILAS − NDVIEC decreased as HLAS − HEC increased during relatively wet periods (increased moisture content). The typical acceptable error between HLAS and HEC of ~50 W m −2 [20,40] (Table 1).

LST-NDVI Space and Soil Moisture
The LST-NDVI space offered an elaborate characterization of surface heterogeneity, as it combines soil moisture and vegetation growth variability based on LST and NDVI, respectively. The results indicated that during any day, some fields experienced higher LST (i.e., drier conditions) and lower NDVI compared to others. Lower during drier conditions suggested relatively more spatial variability (i.e., relatively heterogeneous) compared to higher during wetter conditions (i.e., relatively homogeneous). Also, higher consistently corresponded with higher during drier conditions. The for each isoline (or day) was used to explain and identify the status of the soil moisture of the root zone. The upper-and lower-most isolines can represent the driest and wettest edges of the triangle [27,28,46,49]. The was compared with a normalized soil moisture content in Field 7 to evaluate its temporal variability. A normalized soil moisture (NSM) content was calculated as ⁄ , where SM is the half-hour volumetric soil moisture content, and SMmax and SMmin are the maximum and minimum soil moisture contents. The soil moisture status was the lowest at the upper-most isoline during 4 August (16%), followed by 27 August (56%), 5 September (60%), and 28 September (70%) ( Figure A3). The 4 August was about four days before irrigation, 27 August was three days after irrigation, and 5 and 28 September were just one to two days after irrigation ( Figure A3).
This elaborate characterization can provide a more effective and practical approach to address some challenges in monitoring LAS performance over heterogeneous areas. Some studies characterized homogenous surface conditions based on before and after irrigation events, during which HLAS and HEC closely agreed compared to during irrigation [20,74]. However, over large areas, challenges can arise in identifying such periods due to lack of irrigation, reporting or subjectively identifying periods before and after irrigation which can lead to misleading indications. This qualitative evaluation highlighted the advantage of using the LST-NDVI space to account for soil moisture and vegetation growth variability and its appropriateness to evaluate HLAS. A quantitative evaluation of HLAS based on LST-NDVI space and TVDI can further provide a more effective approach.

HLAS and TVDI
The obtained results are consistent with the empirical triangle concept, as they showed increased TVDI that coincided with higher H during drier surface conditions (i.e., lower surface and root zone soil moisture content) [27][28][29]49,[75][76][77][78]. Also, improved agreement between HLAS and HEC was obtained during higher TVDI during relatively drier conditions. These results suggested that HLAS − HEC increased during relatively wet periods as TVDI decreased (i.e., higher soil moisture contents and LE). To provide a more effective way to evaluate HLAS, a TVDI threshold of 0.30 was identified, above which HLAS − HEC fell within an acceptable level of accuracy of 50 W m −2 . Figure 10 provides significantly useful information for practical applications. For example, thresholds of TVDI can be developed and combined with the large-scale estimates of HLAS and LELAS to manage irrigation more accurately. This analysis indicated that regardless of knowing the irrigation schedule, the TVDI was able to appropriately characterize surface conditions and LAS performance during varying soil moisture content compared to using LST or NDVI individually, as used in previous studies [20,22,74]. While these results are promising, it was important to note that adding more LST images and soil moisture measurements can provide more spatially representative TVDI.

Estimates of LELAS and ET
The resultes indicated that (Rn − G)TSEB IRT consistently underestimated (Rn − G)EC. This can partly be attributed to spatial scale mismatch between the IRT and EC systems' footrpints of 2-3 m pointbased and a few hundered meters, respectively. Consequently, these descrepancies affected LELAS TSEB IRT, LELAS EC, and the agreement between LAS and EC estimates. An improved agreement between LELAS and LEEC can be obtained when using (Rn − G)EC compared to using (Rn − G) independent from EC measurements, as indicated by some previous studies. For example, (Rn − G)IRT data were used to estimate LELAS over flood-irrigated winter wheat, olive yard, and drip-irrigated orange orchards that resulted in an RMSD of 49, 56, and 62 (with R 2 of 0.74, 0.71, and 0.68), respectively, when compared with LEEC [16]. In another study by Hoedjes et al. [79], (Rn − G)EC data were used simultaneously to address lack of energy balance closure and calculate LELAS over wheat, and resulted in an R 2 of 0.97 when compared with LEEC. However, using (Rn − G)EC in this manner does not provide fully independent LELAS and LEEC. Though, such use was not uncommon, especially in cases with data limitations [79]. Also, previous inter-comparison studies of independent measurements of Rn − G showed increased correlation (R 2 > 0.95), with discrepancies mainly attributed to instrumentation error [80,81]. By following this approach in this study, an improved agreement between LELAS EC and LEEC Re was obtained with RMSD, R 2 , and MBD of 19 W m −2 , 0.93, and −34 W m −2 , respectively.
On the other hand, it was also noticed that the discrepancies in both comparisons (i.e., LELAS TSEB IRT vs. LEEC Re and LELAS EC vs. LETSEB IRT) decreased (increased) as LE decreased (increased) during drier (wetter) conditions. The correlation between LELAS TSEB IRT − LEEC Re and TVDI, as well as LELAS EC − LETSEB IRT and TVDI, suggested that both LELAS TSEB IRT − LEEC Re and LELAS EC − LETSEB IRT increased during relatively higher soil moisture conditions, as characterized by decreased TVDI and increased LE [27,28,46,49]. This behavior of LELAS versus TVDI is consistent with the obtained results (Section 4.3) when comparing HLAS discrepancies with TVDI. The observed TVDI behavior in response to soil moisture variation suggested its ability to provide practical means to evaluate LELAS and to manage irrigation over large areas. In this study, periods with relatively higher LE occurred within ~5 days after significant wetting events (i.e., irrigation or rain). This number of days varied from one field to another, and hence, its use was challenging, not practical, and subjective. For such, the TVDI was compared with LELAS EC − LEEC to develop a threshold to evaluate LELAS performance. The datasets of LELAS TSEB IRT − LEEC Re and LELAS EC − LETSEB IRT were not considered in this case because of the limited data points. A threshold of TVDI = 0.25 was identified, above which a higher agreement between LELAS and LEEC Re was obtained, which can vary with the use of a larger dataset.
Measured ETEC Re values throughout the growing season were ≤4 mm day −1 with ≤2 mm day −1 observed during periods with low soil moisture content or near the end of the season (Figure 2). Smaller discrepancies between ETLAS and ETEC Re appeared during these periods and increased right after irrigation or significant rain events when ETEC was around or >3 mm day −1 . In other words, during relatively drier conditions, the LAS and EC systems showed better agreement in ET estimates.
The difference between of ETLAS and ETEC Re estimates was within 16%, which is within reported levels of accuracies. The study by [82] compared ETLAS with those from a remote sensing-based model called SEBAL over a mixed vegetation cover and obtained an error of 17%. Other studies compared ETEC with TSEB model estimates and reported an RMSD within 0.6 mm day −1 over agricultural sites with similar heterogeneity to this study [35,50,56]. On the other hand, suggested acceptable accuracies of ET measurements from LAS and EC were mostly within 10-15%, but in some cases, can reach up to 30% [10,11]. The difference ETLAS − ETEC Re increased as TVDI decreased during wet conditions. A TVDI threshold of ~0.5 was identified, below which ETLAS − ETEC Re exceeded 15%. Additional analysis with larger dataset will need to be conducted in the future to further support this finding.

Conclusions
This study evaluated the effects of surface heterogeneity on LAS performance due to soil moisture and vegetation growth variability caused by irrigation practices. LAS observations and estimates of H, LE, and ET were obtained over row drip-irrigated vineyards during 1 August-31 October 2007 in Albacete, Spain. The study used a suite of methods that included EC observations of H, LE, and ET; LST and NDVI from Landsat and MODIS sensors; additional LST from an IRT; a data fusion model to develop a time series of NDVI; and a surface energy balance model to estimate H and LE. The EC observations showed a lack of closure of 16% during unstable atmospheric conditions and was corrected using the Residual method. Good agreements were obtained between LAS and EC estimates of H, LE, and ET with an RMSD of 25 W m −2 , 19 W m −2 , 0.41 mm day −1 , respectively. LAS overestimated H and underestimated both LE and ET by 24 W m −2 , 34 W m −2 , and 0.36 mm day −1 , respectively.
The effects of soil moisture on HLAS was evaluated using the Bowen ratio, β. Higher discrepancies between HLAS and HEC were obtained at β ≤ 0.5, but improved agreement was shown at 1 ≥ β > 0.5 and β > 1.0, with an R 2 of 0.76, 0.78, and 0.82, respectively. Variable vineyard growth affected LAS and EC performance. The LAS footprints "saw" lower vineyard growth compared to that of the EC, as NDVILAS underestimated NDVIEC by ~0.022, which corresponded to an overestimation of HLAS to HEC by ~24 W m −2 .
The effects of surface heterogeneity on LAS performance due to the combined effects of soil moisture and vegetation growth variability were appropriately characterized using the LST-NDVI space and TVDI. Surface heterogeneity increased during wetter conditions, as indicated by the standard deviation of NDVI and mean LST. Based on its assumptions, LAS accuracy increased during drier conditions, and hence, resulted in improved LAS and EC agreements. The TVDI increased (decreased) during drier (wetter) conditions as the discrepancies between HLAS and HEC, as well as LELAS and LEEC Re decreased (increased). Thresholds of TVDI of 0.3, 0.25, and 0.5 were identified, above which improved agreement was obtained between LAS and EC estimates of H, LE, and ET, respectively. These thresholds can vary with the use of more LST-NDVI datasets and they are also region-dependent.
These findings can have important practical use, as they suggest the effectiveness and ability of LAS in monitoring vegetation growth over heterogonous areas with variable soil moisture. The combined use of LAS with LST-NDVI space and TVDI can support irrigation scheduling and agricultural water management over large regions.

Acknowledgments:
The authors would like to thank two anonymous reviewers who provided valuable comments that helped in strengthening the manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Appendix A Figure A1. A description of the vineyard fields drip irrigation schedule. The blue color represents the rainfall events. The light grey color differentiates between the months of August and September. The dark grey color represents the irrigation events over the individual fields.  Based on [38] and [39], is related to as follows: where AT and Aq are coefficients that can be calculated as a function of atmospheric pressure (p), air temperature (T), humidity (q), and optical wavelength as 0.78 • ⁄ 10 and 57.22 • 10 , cp is the specific heat at constant pressure. [38] found that during certain conditions, the factor in front of (i.e. equals 0.031 and Equation (A1) is simplified to Equation (1), as shown in Section 3.1 above.

Appendix C
Appendix C. 1

. TSEB Model Description
The surface energy balance equations over canopy and bare soil components can be described as: where subscripts s and c refer to soil and canopy components. The sensible heat flux can be estimated as: where Tac is the air temperature at an air-canopy interface, ρ the air density taken as 1.24 (kg m −3 ), cp the specific heat of air taken as 1005 (J kg −1 K −1 ), rx the total boundary layer resistance of the complete canopy leaves, rs the resistance to heat flow in the boundary layer immediately above the soil surface, and ra is the aerodynamic resistance to heat transfer. The latent heat flux for each component can be estimated as: where G can be estimated as 0.30 of Rns, αPT is the Priestly-Taylor constant taken as 1.26, fg the fraction of leaf area index (LAI) that is green, Δ the slope of the saturation vapor pressure versus temperature curve, and γ the psychrometric constant. The initial value of αPT = 1.26 is used for well-water unstressed vegetation, which can further be adjusted by following an iterative process to account for stressed vegetation. More details about the TSEB model and the methods followed to estimate the different variables can be found in many applications [35,50,56,62]. , where z is the mean plume height for diffusion from a surface source (m), zm is the measurement height, ( ) u z is the mean wind profile variables A, b and c are functions of the gamma function (), and r is the Gaussian plume model shape parameter. The diffusion in the lateral direction is assumed Gaussian [59].
To obtain the three-dimensional footprint function for the LAS, Equation (A12) was combined with the spatial weighting function of the scintillometer [14,15,20]. The LAS path was thought to consist of a series of single observation points, which, for each a single source area, was calculated using Equation (A11). Each of these source areas was then multiplied by the LAS weighting function, G (u). The summation of all individual points was normalized by the total of the source areas, which then yielded the LAS weighted footprint.