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

Abstract

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⁻², 19 W m⁻², and 0.41 mm day⁻¹, respectively. LAS overestimated H and underestimated both LE and ET by 24 W m⁻², 34 W m⁻², and 0.36 mm day⁻¹, 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² 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.

1. 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³/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,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.

2. Data

2.1. 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⁻¹, 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.

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 [23]. The vegetation height at the beginning of the growing season was 70 cm and 164 cm by the end of the season, with an average width of foliage of 70 cm.

2.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 1-min intervals. Air temperature (Ta) and atmospheric pressure (P) measurements were acquired using separate sensors that were integrated with the LAS system.

2.3. 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, 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⁻²) 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⁻²) 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 well-known 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.

2.4. Remote Sensing Data

2.4.1. 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.

2.4.2. 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.

2.5. 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 installed—adjacent 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³ cm⁻³ and 0.14 cm³ cm⁻³, 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.

3. Methods

3.1. 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, $C_n^2$. 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, $C_T^2$ and $C_q^2$, respectively, give rise to $C_n^2$. The relationship between $C_T^2$ and $C_n^2$ can be described by Equation (1) [13,38,39] as:

(1)$C_T^2=C_n^2{\left(\frac{T_a^2}{\gamma ·P}\right)}^2{\left(1+\frac{0.03}{\beta }\right)}^{-2}$

where $C_T^2$ and $C_n^2$ are in K² m^−2/3 and m^−2/3, respectively, P the atmospheric pressure (Pa), γ the refractive index coefficient for air (0.78 × 10⁻⁶ K Pa⁻¹), 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 (H_LAS) 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), $C_T^2$ can be related to the temperature scale, $T_{\ast }$, by a universal function of stability parameter $\xi =\left(z_{eff}-d\right)/L$ as:

(2)$\frac{C_T^2·{\left(z_{eff}-d\right)}^{2/3}}{T_{\ast }^2}=f_T\left(\xi \right)$

where $z_{eff}$ is the LAS effective height (m) and L is the Monin–Obukhov length. For unstable conditions (L < 0) [40,41]:

(3)$f_T=4.9{\left[1-6.1\left(\frac{z_{eff}-d}{L}\right)\right]}^{-2/3}$

The effective height, $z_{eff}$, takes into account the effects of stability conditions and slanted paths as $z_{eff}^{-2/3}·f_T\left(z_{eff}/L\right)={{\displaystyle \int }}_{u=0}^{u=1}\left(1/z^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\left(u\right)\right)f_T\left(z\left(u\right)/L\right)G\left(u\right)du$ with G(u) is the bell-shaped weighting function describing the contribution of $C_n^2$ to the total LAS signal at each point along the normalized path u [21], and $z\left(u\right)$ is the variable LAS beam height along the path (Figure 1). For slanted paths, $z\left(u\right)$ can be described relative to the maximum and minimum heights if the transmitter and/or the receiver $z_{max}$ and $z_{min}$, respectively, as $z\left(u\right)=z_{min}\left[1+\left(z_{max}/\left(z_{min}-1\right)\right)·u\right]$. The Monin–Obukhov length, L, can be estimated as:

(4)$L=\frac{T_a·u_{\ast }}{k·g·T_{\ast }}$

with k = 0.4 is the von Karman constant, $g$ the gravity (9.81 m·s⁻²), and the friction velocity, $u_{\ast }$, can be estimated using the standard Businger–Dyer flux profile [42] as:

(5)$u_{\ast }=\frac{k·U}{ln\left(\frac{z_{eff}-d}{z_0}\right)-{\psi }_m\left(\frac{z_{eff}-d}{L}\right)+{\psi }_m\left(\frac{z_0}{L}\right)}$

where U is the wind speed, ψ_m the stability correction function for the momentum transfer was calculated following [43,44] as ${\psi }_m=2ln\left[\left(1+x\right)/2\right]+ln\left[\left(1+x^2\right)/2\right]-2·arctan\left(x\right)+\pi /2$ with $x={\left[1-16\left(z_{eff}-d\right)/L\right]}^{1/4}$, d is the displacement height, and $z_0$ is the roughness length. Both d and $z_0$ were calculated as a function of crop height $h_c$ as $d=\left(2/3\right)·h_c$ and $z_0=\left(1/3\right)·h_c$, respectively [45].

H_LAS can be estimated as:

(6)$H_{LAS}=\rho ·c_p·u_{\ast }·T_{\ast }$

where $\rho$ is the air density and $c_p$ is the specific heat of air at constant pressure.

3.2. 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)$TVDI=\frac{LST -LS{T}_{min}}{LS{T}_{max}-LS{T}_{min}}$

where $LS{T}_{min}$ is the minimum LST and $LS{T}_{max}$ the maximum LST that can be estimated as a function of NDVI as $LS{T}_{max}=a+b·NVDI,$ 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].

3.3. Data Fusion

The STARFM is a data fusion approach that was developed to blend daily MODIS with the 16-day 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.

3.4. 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 H_LAS. 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:

(8)$T_R={\left[f_c·T_c^4+\left(1-f_c\right)·T_s^4\right]}^{1/4}$

with $T_R$ is the LST or sometimes referred to as radiometric surface temperature, and $T_c$ and $T_s$ are canopy and bare soil temperatures (°K). The fraction of cover fc viewed at nadir can be estimated as $f_c=1-e^{-\left(0.5·LAI\right)}$, 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].

3.5. 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.

4. Results

4.1. Observed Surface Energy Balance Closure

Observed H_EC + LE_EC 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² of 54 W m⁻² and 0.69, respectively. Adjusted sensible heat flux using the BR method (H_{EC β}) showed a mean value of 111 W m⁻² (standard deviation (SD) of 56 W m⁻²), with a 15% increase compared to that of raw H_EC of 95 W m⁻² (SD of 46 W m⁻²). Similarly, adjusted latent heat flux using the Re (LE_{EC Re}) and BR (LE_{EC β}) methods resulted in mean values of 176 (SD of 89 W m⁻²) and 159 W m⁻² (SD of 78 W m⁻²), with about 30% and 15% increase, respectively, compared to that of LE_EC of 138 W m⁻² (SD of 68 W m⁻²).

4.2. Comparison between H_LAS and H_EC

The mean and SD statistics of LAS and EC measurements of H, H_LAS, and H_EC and H_{EC β}, respectively, showed that the mean of H_LAS of 118 W m⁻² (SD of 52 W m⁻²) was ~26% higher than that of H_EC of 94 W m⁻² (SD of 56 W m⁻²) and 6% higher than that of H_{EC β} of 111 W m⁻² (SD of 46 W m⁻²). These statistics indicate that regardless of the method used to address lack of closure, H_LAS 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 H_LAS and H_EC showed some discrepancies, as indicated by an RMSD of 25 W m⁻² (R² 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 H_LAS and H_{EC β}, with an RMSD of 28 W m⁻² (R² 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 H_LAS and H_EC compared to H_LAS and H_{EC β}, 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.

4.3. H_LAS, Soil Moisture, and Vegetation Spatial Variability

4.3.1. H_LAS and β

The effects of soil moisture temporal variability on H_LAS 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 H_LAS and H_EC 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, H_LAS overestimated H_EC by a mean bias difference (MBD) of 34 W m⁻². At 1 ≥ β > 0.5, the MBD decreased to 27 W m⁻², and at β > 1, H_LAS underestimated H_EC by an MBD of −14 W m⁻². In other words, as β increased (i.e., drier surface conditions became dominant), the agreement between H_LAS and H_EC improved as the MBD consistently decreased. Also, R² of 0.76, 0.78, and 0.82 for β < 0.5, 1 ≥ β > 0.5, and β > 1, respectively, show consistently improved agreement between H_LAS and H_EC as β increased.

4.3.2. 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 H_LAS and LE_LAS.

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 NDVI_LAS and NDVI_EC, 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 NDVI_EC for each image data and time. A similar process regarding the calculation of NDVILAS was followed. The results show that both NDVI_LAS and NDVI_EC 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 NDVI_LAS was consistently, but slightly, lower than NDVI_EC (Figure 6).

A time series of footprint-integrated pairs of NDVI_LAS and NDVI_EC 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 NDVI_LAS and NDVI_EC were closely correlated (R² = 0.90), and the average NDVI_LAS was ~0.022 lower than NDVI_EC.

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 (${\sigma }_{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 ${\sigma }_{NDVI}$ 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 (${\mu }_{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 ${\mu }_{LST}$ 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.

4.3.3. H_LAS and TVDI

Estimates of TVDI

As an initial evaluation, the obtained TVDI, which ranged from 0–1, was plotted against H_EC (not shown) to make sure that it followed the known behavior, as described by [27,28,46,49]. A positive correlation between TVDI and H_EC 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 H_TSEB

The TSEB model was used to estimate H_TSEB and (Rn − G)_TSEB that can be used to evaluate both H_LAS–H_EC and H_LAS–H_TSEB against TVDI and to calculate LE_LAS, 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) (H_{TSEB RS}) and those based on the IRT (H_{TSEB IRT}) are shown in Figure 9. Two different comparisons were conducted—H_EC against H_{TSEB IRT} and H_LAS against H_{TSEB IRT}. The comparison between H_EC and H_{TSEB IRT} showed an RMSD and R² of 35 W m⁻² and 0.71, respectively, with a narrow scatter distribution around the 1:1 line. Slightly lower agreement was obtained between H_LAS and H_{TSEB IRT}, with an RMSD and R² of 51 W m⁻² and 0.64, respectively. One of the main reasons for the increased discrepancies between H_LAS and H_{TSEB IRT} was attributed to differences in the spatial representation of H_{TSEB IRT} (at a point) to that of H_LAS. Estimates of H_{TSEB RS} did not have enough data points that could be meaningfully evaluated. However, a visual inspection of Figure 9 showed that the footprint-integrated H_EC and H_LAS were closely correlated with H_{TSEB RS}. Overall, these TSEB model results are consistent with previous studies [14,35,50,56,61,62] and allowed the use of H_{TSEB IRT} to characterize surface conditions.

H_LAS and TVDI

Estimates of H_LAS, H_EC, and H_{TSEB 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 H_LAS − H_EC and H_LAS − H_{TSEB IRT} against TVDI are shown in Figure 10. The results show a negative correlation (R² = 0.55) between H_LAS − H_EC and TVDI, which indicates that as the discrepancies between H_LAS and H_EC decreased, the TVDI increased. A relatively higher correlation (R² = 0.66) was obtained between H_LAS − H_{TSEB IRT} and TVDI based on a smaller sample size. This comparison also provides a similar indication of decreased discrepancies between H_LAS and H_{TSEB IRT} with increased TVDI. In other words, both comparisons indicate that the agreement between H_LAS and both H_EC and H_{TSEB IRT} improved during higher TVDI or specifically drier surface conditions. This consistent behavior of H_LAS and TVDI suggests that TVDI can effectively be used to evaluate the effects of surface heterogeneity due to irrigation on LAS performance.

4.4. Estimates of LE and ET

4.4.1. Estimates of LE_LAS

LE_LAS was estimated as the residual of the energy balance equation as LE_LAS = Rn − G − H_LAS. Estimates of (Rn − G)_{TSEB IRT} based on the TSEB and IRT data (Section 3.4) were used to obtain relatively independent LE_LAS compared to LE_{EC Re}. Estimates of (Rn − G)_{TSEB IRT} resulted in an RMSD and R² of 34 W m⁻² and 0.93, respectively, with an underestimation indicated by an MBD of −30 W m⁻² when compared with (Rn − G)_EC (Figure 11). Two different estimates of LE_LAS were obtained: LE_{LAS TSEB IRT} based on (Rn − G)_{TSEB IRT} and LE_{LAS EC} based on (Rn − G)_EC. LE_{LAS TSEB IRT} and LE_{LAS EC} were evaluated against LE_{EC Re} and LE_{TSEB IRT}, respectively (Figure 11). These comparisons show that both estimates of LE_{LAS TSEB IRT} and LE_{LAS EC} consistently underestimated LE_{EC Re} and LE_{TSEB IRT}, respectively. The comparison between LE_{LAS EC} and LE _{TSEB IRT} showed a reasonable distribution around the 1:1 line, with an RMSD and R² of 39 W m⁻² and 0.71, respectively, and underestimation indicated by an MBD of −20 W m⁻². The comparison between LE_{LAS TSEB IRT} and LE_{EC Re} showed a narrower distribution around the 1:1 line, with an RMSD and R² of 30 W m⁻² and 0.82, respectively, and underestimation indicated by an MBD of −59 W m⁻² (Figure 11).

4.4.2. LE_LAS and TVDI

The differences LE_{LAS TSEB IRT} − LE_{EC Re} and LE_{LAS EC} − LE_{TSEB IRT} were compared with TVDI to evaluate the effects of soil moisture and vegetation growth variability on LE_LAS. The obtained results show that both differences were negatively correlated with TVDI with first comparison (i.e., LE_{LAS TSEB IRT} − LE_{EC Re} versus TVDI) resulted in R² = 0.72, and the second one (i.e., LE_{LAS EC} − LE_{TSEB IRT} versus TVDI) indicated an R² = 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 LE_{LAS TSEB IRT}, LE_{EC Re}, LE_{LAS EC}, and LE_{TSEB IRT}. Therefore, they were used to qualitatively evaluate LE_LAS and TVDI behavior. LE_{LAS EC} was used to provide a quantitative evaluation of TVDI, as discussed in Section 5.7.

4.4.3. Estimates of ET

Estimates of ET_LAS and ET_{EC Re} based on LE_LAS and LE_{EC Re}, respectively (Figure 11), show a good agreement, with an RMSD of 0.41 mm day⁻¹. An underestimation of ET_LAS to ET_{EC Re} was noticed with an MBD of −0.36 mm day⁻¹. This result was consistent with that obtained in Section 4.3, which indicated an overestimation of H_LAS to H_EC and, consequently, an underestimation of LE_LAS to LE_{EC Re}. The obtained mean values for ET_LAS and ET_{EC Re} were 2.1 and 2.5 mm day⁻¹, respectively. A better LAS performance was visually observed at ET_LAS ≤ 2 mm day⁻¹ with increased discrepancies above this value.

5. Discussion

5.1. 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 H_LAS and H_EC. 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.

5.2. H_LAS and H_EC

The obtained RMSD between H_LAS and H_EC was within acceptable levels of accuracy. Reported typical levels of accuracy of H measurements for LAS and EC were within 30 W m⁻² [40] and 50 W m⁻² [12,67,68], respectively. Also, [20] indicated that acceptable differences between H_LAS and H_EC were within 50 W m⁻². [20] compared H_LAS against H_EC over flood-irrigated olive orchards and obtained an RMSD of 36, 29, and 63 W m⁻² 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⁻² 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 H_LAS of H_EC.

5.3. H_LAS and β

Considering soil moisture effects, the use of β in Equation (1) can also partly explain some of the discrepancies between H_LAS and H_EC based on the humidity correction factor and the corresponding assumptions. The humidity correction factor [${\left(1+0.03/\beta \right)}^{-2}$] partly affects the accuracy of $C_T^2$ and H_LAS [38,39]. The humidity correction factor is less than 10% for of β ≥ 0.6 [21]. Thus, it has minimal effects on H_LAS accuracy at larger β and no need to use it at β > 1 [14,21,40]. On the other hand, the uncertainty in H_LAS can reach up to 15% and 48% at β = 0.3 and β < 0.3, respectively [40]. In other words, decreased uncertainty in H_LAS during drier conditions (i.e., H > LE and β > 1) can lead to improved agreement with H_EC 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 H_LAS and H_EC at β ≤ 0.5 increased. Equation (1) is a simplified form of Equation (A1) (Appendix B) that relates $C_n^2$, $C_T^2$, and $C_q^2$ [39], assuming that the correlation between T and q, Rtq = +1. This assumption removes the direct dependency of $C_n^2$ on $C_q^2$, and hence reduces its contribution to [${\left(1+0.03/\beta \right)}^{-2}$]. 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 H_LAS 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 H_LAS to H_EC, 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 H_LAS and H_EC 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]

5.4. 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 water-holding 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 H_LAS to H_EC and it is also consistent with the basic concepts of near surface flux exchange [55,72,73]. An evaluation of NDVI_LAS − NDVI_EC against H_LAS − H_EC showed that H_LAS − H_EC slightly increased as NDVI_LAS − NDVI_EC decreased (Figure 7b) with a low correlation (R² = 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. NDVI_LAS − NDVI_EC decreased as H_LAS − H_EC increased during relatively wet periods (increased moisture content). The typical acceptable error between H_LAS and H_EC of ~50 W m⁻² [20,40] corresponded to NDVI_LAS − NDVI_EC of ~0.018. The average of NDVI_LAS − NDVI_EC of 0.022 can partly explain the overestimation of H_LAS to H_EC of MBE of 24 W m⁻² (Table 1).

5.5. 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 ${\sigma }_{NDVI}$ during drier conditions suggested relatively more spatial variability (i.e., relatively heterogeneous) compared to higher ${\sigma }_{NDVI}$ during wetter conditions (i.e., relatively homogeneous). Also, higher ${\mu }_{LST}$ consistently corresponded with higher ${\sigma }_{NDVI}$ during drier conditions. The ${\mu }_{LST}$ 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 ${\mu }_{LST}$ 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 $NSM=\left[\left(SM-S{M}_{min}\right)/\left(S{M}_{max}-S{M}_{min}\right)\right]$, where SM is the half-hour volumetric soil moisture content, and SM_max and SM_min 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 H_LAS and H_EC 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 H_LAS. A quantitative evaluation of H_LAS based on LST–NDVI space and TVDI can further provide a more effective approach.

5.6. H_LAS 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 H_LAS and H_EC was obtained during higher TVDI during relatively drier conditions. These results suggested that H_LAS − H_EC increased during relatively wet periods as TVDI decreased (i.e., higher soil moisture contents and LE). To provide a more effective way to evaluate H_LAS, a TVDI threshold of 0.30 was identified, above which H_LAS − H_EC fell within an acceptable level of accuracy of 50 W m⁻².

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 H_LAS and LE_LAS 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.

5.7. Estimates of LE_LAS 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 point-based and a few hundered meters, respectively. Consequently, these descrepancies affected LE_{LAS TSEB IRT}, LE_{LAS EC}, and the agreement between LAS and EC estimates. An improved agreement between LE_LAS and LE_EC 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 LE_LAS over flood-irrigated winter wheat, olive yard, and drip-irrigated orange orchards that resulted in an RMSD of 49, 56, and 62 (with R² of 0.74, 0.71, and 0.68), respectively, when compared with LE_EC [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 LE_LAS over wheat, and resulted in an R² of 0.97 when compared with LE_EC. However, using (Rn − G)_EC in this manner does not provide fully independent LE_LAS and LE_EC. 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² > 0.95), with discrepancies mainly attributed to instrumentation error [80,81]. By following this approach in this study, an improved agreement between LE_{LAS EC} and LE_{EC Re} was obtained with RMSD, R², and MBD of 19 W m⁻², 0.93, and −34 W m⁻², respectively.

On the other hand, it was also noticed that the discrepancies in both comparisons (i.e., LE_{LAS TSEB IRT} vs. LE_{EC Re} and LE_{LAS EC} vs. LE_{TSEB IRT}) decreased (increased) as LE decreased (increased) during drier (wetter) conditions. The correlation between LE_{LAS TSEB IRT} − LE_{EC Re} and TVDI, as well as LE_{LAS EC} − LE_{TSEB IRT} and TVDI, suggested that both LE_{LAS TSEB IRT} − LE_{EC Re} and LE_{LAS EC} − LE_{TSEB IRT} increased during relatively higher soil moisture conditions, as characterized by decreased TVDI and increased LE [27,28,46,49]. This behavior of LE_LAS versus TVDI is consistent with the obtained results (Section 4.3) when comparing H_LAS discrepancies with TVDI.

The observed TVDI behavior in response to soil moisture variation suggested its ability to provide practical means to evaluate LE_LAS 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 LE_{LAS EC} − LE_EC to develop a threshold to evaluate LE_LAS performance. The datasets of LE_{LAS TSEB IRT} − LE_{EC Re} and LE_{LAS EC} − LE_{TSEB 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 LE_LAS and LE_{EC Re} was obtained, which can vary with the use of a larger dataset.

Measured ET_{EC Re} values throughout the growing season were ≤4 mm day⁻¹ with ≤2 mm day⁻¹ observed during periods with low soil moisture content or near the end of the season (Figure 2). Smaller discrepancies between ET_LAS and ET_{EC Re} appeared during these periods and increased right after irrigation or significant rain events when ET_EC was around or >3 mm day⁻¹. In other words, during relatively drier conditions, the LAS and EC systems showed better agreement in ET estimates.

The difference between of ET_LAS and ET_{EC Re} estimates was within 16%, which is within reported levels of accuracies. The study by [82] compared ET_LAS with those from a remote sensing-based model called SEBAL over a mixed vegetation cover and obtained an error of 17%. Other studies compared ET_EC with TSEB model estimates and reported an RMSD within 0.6 mm day⁻¹ 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 ET_LAS − ET_{EC Re} increased as TVDI decreased during wet conditions. A TVDI threshold of ~0.5 was identified, below which ET_LAS − ET_{EC Re} exceeded 15%. Additional analysis with larger dataset will need to be conducted in the future to further support this finding.

6. 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⁻², 19 W m⁻², 0.41 mm day⁻¹, respectively. LAS overestimated H and underestimated both LE and ET by 24 W m⁻², 34 W m⁻², and 0.36 mm day⁻¹, respectively.

The effects of soil moisture on H_LAS was evaluated using the Bowen ratio, β. Higher discrepancies between H_LAS and H_EC were obtained at β ≤ 0.5, but improved agreement was shown at 1 ≥ β > 0.5 and β > 1.0, with an R² 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 NDVI_LAS underestimated NDVI_EC by ~0.022, which corresponded to an overestimation of H_LAS to H_EC by ~24 W m⁻².

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 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 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.

Author Contributions

Conceptualization, H.M.E.G., J.G.-P., C.M.U.N., C.B., I.C., and A.C.; methodology, H.M.E.G. and J.G.-P.; software, H.M.E.G. and J.G.-P.; validation, H.M.E.G., J.G.-P., C.M.U.N., C.B., I.C., and A.C.; formal analysis, J.G.-P. and H.M.E.G.; investigation, J.G.-P., H.M.E.G., C.M.U.N., C.B., I.C., and A.C.; resources, J.G.-P. and H.M.E.G.; data curation, J.G.-P.; writing—original draft preparation, J.G.-P. and H.M.E.G.; writing—review and editing, H.M.E.G. and J.G.-P.; visualization, H.M.E.G. and J.G.-P.; supervision, J.G.-P.; project administration, J.G.-P.; funding acquisition, J.G.-P. and H.M.E.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was based upon work supported by the HERMANA (Herramientas para el Manejo Sostenible de Fertilización Nitrogenada y Agua) project, funded by the Spanish Ministry of Economy and Competitiveness (AGL2015-68700-R), New Mexico EPSCoR project funded by the National Science Foundation (NSF) award # IIA-1301346, and New Mexico State University.

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.

[Figure omitted. See PDF]

Figure A2

A description of the LST–NDVI space and TVDI calculations.

[Figure omitted. See PDF]

Figure A3

A description of relative soil moisture content, rain, and irrigation time series Field 7, along with Landsat overpass dates.

[Figure omitted. See PDF]

Appendix B

Appendix C

Appendix D

Figures and Table

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.

Figure 2. A plot of vineyards fraction of cover (Fc), actual evapotranspiration (ET) based on EC measurements, irrigation dates for Field 7, LAS measurements period, and Landsat overpass dates with respect to day of year (DOY) for year 2007.

Figure 3. EC surface energy balance lack of closure based on a comparison between turbulence fluxes (LE + H) and available energy (Rn − G). RMSD, root mean square difference.

Figure 4. Comparisons between HLAS and HEC β (right), HLAS and HEC Re (left), along with corresponding linear regression and 95% confidence intervals.

Figure 5. Comparison between HLAS and HEC at different β categories, with of β ≤ 0.5 (top), 1.0 ≥ β > 0.5 (middle), and β > 1.0 (bottom).

Figure 6. Normalized difference vegetation index (NDVI) (top row) and land surface temperature (LST) (middle row) based on Landsat Images along with plots of a spatial distribution of NDVI (bottom left) and LST (bottom right) of all fields.

Figure 7. Scatter plots of (a) NDVILAS and NDVIEC, and (b) NDVILAS − NDVIEC compared to HLAS − HEC.

Figure 8. Plot of LST–NDVI space, along with isolines for each Landsat overpass date.

Figure 9. Maps of the two-source energy balance (TSEB) model estimates of HTSEB (top row) and LETSEB (second row) based on Landsat, comparison between EC Re measurements and TSEB IRT estimates of Rn, G, H, and LE (third row left), scatter plot of HLAS and HTSEB IRT (third row right), comparison of footprint integrated LETSEB RS, LEEC Re, and LELAS (bottom left), comparison of HTSEB RS, HEC, and HLAS (bottom right).

Figure 10. Comparison between the error in H (HLAS − HEC) and temperature vegetation dryness index (TVDI).

Figure 11. Comparisons between (a) LELAS TSEB IRT and LEEC Re, (b) LELAS EC and LETSEB IRT, (c) (Rn − G)TSEB IRT versus (Rn − G)EC, and (d) ETLAS and ETEC Re.

Figure 12. The vineyard fields and the dominant LAS and EC footprints (left), and the weighted contribution of the vineyard fields on LAS and EC observations (right).

Table 1

Summary of sensible heat flux comparison and regression statistic.

Equation	Slope a	Intercept b (W m−2)	RMSD (W m−2)	R2	MBD * (W m−2)	No. Data Points
HLAS = a·HEC + b	0.99	25	25	0.77	24	1244
HLAS = a·HEC β + b	0.79	31	28	0.72	7	1244

* Mean bias difference (MBD).