Determining Evapotranspiration by Using CombinationEquation Models with Sentinel-2 Data and Comparison with Thermal-Based Energy Balance in a California Irrigated Vineyard

A new approach is proposed to derive evapotranspiration (E) and irrigation requirements by implementing the combination equation models of Penman–Monteith and Shuttleworth and Wallace with surface parameters and resistances derived from Sentinel-2 data. Surface parameters are derived from Sentinel-2 and used as an input in these models; namely: the hemispherical shortwave albedo, leaf area index and water status of the soil and canopy ensemble evaluated by using a shortwave infrared-based index. The proposed approach has been validated with data acquired during the GRAPEX (Grape Remote-sensing Atmospheric Profile and Evapotranspiration eXperiment) in California irrigated vineyards. The E products obtained with the combination equation models are evaluated by using eddy covariance flux tower measurements and are additionally compared with surface energy balance models with Landsat-7 and -8 thermal infrared data. The Shuttleworth and Wallace (S-W S-2) model provides an accuracy comparable to thermal-based methods when using local meteorological data, with daily E errors < 1 mm/day, which increased from 1 to 1.5 mm/day using meteorological forcing data from atmospheric models. The advantage of using the S-W S-2 modeling approach for monitoring ET is the high temporal revisit time of the Sentinel-2 satellites and the finer pixel resolution. These results suggest that, by integrating the thermal-based data fusion approach with the S-W S-2 modeling scheme, there is the potential to increase the frequency and reliability of satellite-based daily evapotranspiration products.


Introduction
Evapotranspiration is the key variable for determining crop water requirements and is needed for the optimal allocation of water resources for irrigation. Several methods and experiments have demonstrated that earth observation (EO) is an effective tool to derive evapotranspiration for supporting irrigation and water resources [1,2] from local to regional scales. These methods can be grouped into three main groups: (a) reflectance and vegetation index (VI) based-methods: crop coefficient and canopy parameters, such as hemispherical albedo and leaf area index LAI, are obtained by means of different analytics from reflectance or vegetation indices. These parameters are the basic inputs for the application of the widely used FAO-56 approach [3] for determining crop evapotranspiration, either by the direct calculation of the combination equation of Penman-Monteith or by using the crop coefficient and reference evapotranspiration. These methods are already implemented in operational applications for irrigation management [4] and are evolving toward the more explicit definition of the canopy conductance [5,6]; (b) thermal-based energy balance models: land surface temperature is the main input for estimating sensible heat flux and then latent heat flux as a residual term of the surface energy balance. Significant advancements have been made from the first contextual approaches using soil-vegetation-transfer models [7] to one-source models i.e., SEBAL [8], SEBS [9] and METRIC [10], and two-source models, such as TSEB [11] and ALEXI [12]. Thermal-based models have been intensively applied by using observation data from Landsat [13], which is, at the present moment, the only operational platform with medium resolution acquisitions in the thermal infrared (100 m), which are resampled to 30 m with a revisit time of 8 to 16 days depending on the site; (c) EO-driven soil water balance models: these are simulation models of water balance using EO-based input data related to crop development [14,15]. These models produce a continuous spatially distributed output, the quality of which strongly depends on the availability of reliable soil physical and hydraulic properties, as well as precipitation/irrigation inputs [16].
Compared to groups (b) and (c), (a) is particularly suitable for operational irrigation management for two main reasons: (i) there is a wide availability of medium and high resolution space-borne multispectral sensors, such as the Landsat series and Sentinel-2 of the Copernicus Program of European Space Agency, collecting data in the range of the electromagnetic spectrum used for this method, with frequent acquisitions over the same area; and (ii) the data required for achieving reliable results from soil water balance models with sufficient accuracy are often unavailable. However, it also should be noted that models in (b) are fusing and sharpening thermal observations using visible, near-IR and SWIR data from multiple satellites and developing operational evapotranspiration products that are applicable for water management and irrigation scheduling at field scales [17][18][19][20][21].
Recently, an approach based on Sentinel-2 and a combination equation with resistance terms derived by means of an empirical relationship with meteorological variables and vegetation parameters, such as the leaf area index, has been proposed [6]. However, for a given set of meteorological variables, the resistance terms in the calculation of evapotranspiration depend only on canopy parameters (i.e., LAI, fractional vegetation cover) without variations in the water status of the soil-canopy ensemble within the footprint of the meteorological dataset.
Diversely, in this study, an innovative approach for estimating crop evapotranspiration is proposed and compared with experimental data in an irrigated vineyard under Mediterranean climate conditions, aiming at the full exploitation of the Sentinel-2 sensors in terms of geometrical, temporal and radiometric resolutions. As explained in this paper, the proposed approach uses the full range of Sentinel-2 spectral observations to infer the data required for computing evapotranspiration in partial canopies by means of a modified combination equation. This modified method, based on the approach developed by Shuttleworth and Wallace [22], estimates soil and canopy contributions to the total evapotranspiration E (the sparse crop combination equation) by incorporating shortwave infrared data from Sentinel-2 for assessing the water status of the surface and modulating the resistance terms in the combination equation. Results will also be compared with two thermal-based model output sources available for the same study area.

Theory
The combination equation of Penman-Monteith for calculating evapotranspiration E is based on the simultaneous solution of the surface energy balance equation and the turbulent transport of heat and water vapor by means of resistance terms defined in different ways. The general formulation is given by the following widely used equation [23]: The different variables represented in the Equation (1)  The available energy flux density A is given by the difference between the total net radiation flux R n and soil heat flux G: with α (-) representing the spectrally integrated hemispherical albedo of the surface-for the sake of simplicity, surface albedo in the rest of the text-K ↓ representing the incoming global shortwave radiation and L* representing the net incoming long-wave radiation flux density. This latter term was calculated at the daily scale from air temperature and vapor pressure, incoming shortwave radiation and extraterrestrial solar radiation [3]. The aerodynamic resistance for sensible heat and latent heat, r a h , is given by: In this equation, h c is the canopy height, k is the von Karman's constant (0.41), z u and z T are the measurement heights for wind speed u and temperature T, d 0 is the zero-plane displacement height and z 0m and z 0h represent the roughness lengths for momentum and heat. The values of d 0 , z 0m and z 0h are taken as proportional to the canopy height h c [24], with factors 0.67, 0.123 and 0.0123, respectively.
When calculations are made for daily intervals, especially with high solar radiation, atmospheric stability corrections can be considered as negligible due to night-time compensating effects [25].
Finally, the bulk stomatal canopy resistance r s c in Equation (1) is given by: The expression of r s c in Equation (5) was first introduced by Szeicz and Long [26] and successively adopted in the formulation of the Penman-Monteith equation proposed by the FAO-56 paper [3]. The coefficient 0.5 in the denominator was introduced because, in most cases, only the upper half of crop foliage is actively contributing to the heat and vapor transfer [26]. Although the Penman-Monteith approach was originally developed for closed canopies, several hydrological studies have also applied it in the presence of canopy gaps, provided that appropriate values of resistances are introduced to account for the gap fractions. To this end, the value of r s c can be empirically determined from latent heat flux measurements by a numerical inversion of the Penman-Monteith equation [27]. Different models have been proposed to determine the canopy resistance in relation to LAI and other environmental variables [28][29][30].
The leaf resistance r leaf in Equation (5) is the reciprocal of the leaf conductance, which has been investigated in several physiological studies [31,32]. From measurements on leaves of different species, including tree and herbaceous crops, a minimum value (maximum conductance) of r leaf = 100 sm −1 has been suggested, which is also adopted in the FAO-56 approach in fully irrigated and disease-free conditions. This value has been confirmed on the basis of Fluxnet data from crop and grassland sites [33]. Assuming a constant fixed value for r leaf , i.e., 100 sm −1 , allows for the calculation of E on the basis of meteorological data and canopy parameters, such as surface albedo, LAI and crop height, which can be estimated by using remote sensing techniques. This approach has been implemented in satellite-based advisory services for irrigation management without the need for using the crop coefficient K c [2].
With an increasing water deficit in the root zone, potential stress conditions may occur due to limitations in the water uptake by roots. The variability of r leaf in relation to water stress conditions has been investigated for a Merlot grape variety in Israel that was subjected to different irrigation treatments [34], which reports measurements of the leaf resistance r leaf varying from 100 sm −1 in mild water stress to 400 sm −1 for severe stress. However, it is well known that the increase in r leaf , representing stomatal closure, occurs when the soil water content is lower than a certain threshold, by following approximately a bi-linear function, as demonstrated by observations on irrigated wheat and olive trees [35,36]. Based on these studies, the variability of the leaf resistance can be considered to vary within the range 100-400 sm −1 , depending on environmental conditions and the soil water content within the root zone.

The Sparse Canopy Combination Equation
One of the main limitations of the Penman-Monteith approach for estimating crop evapotranspiration E and irrigation requirements is the so-called "big-leaf" assumption [23], which is not valid for row or sparse crops i.e., crops that do not completely cover the soil surface. In this latter condition, the resulting value of evapotranspiration is underestimated, especially when the soil is wet, because the soil evaporation is not adequately represented.
To extend the Penman-Monteith schematization in sparse crops, Shuttleworth and Wallace [22] proposed the following formulation for explicitly partitioning canopy transpiration E c and soil evaporation E s (sparse crop combination equation): In this equation, the term D 0 represents the vapor pressure deficit at the canopy height, which can be expressed as follows: The available energy at the substrate level A s is derived from the net radiation flux density at the substrate R ns by means of the commonly used Beer's law, accordingly to the relationships: The vertical transport is controlled by two aerodynamic resistances, the first one from canopy level to the reference screen height, r a a , and the second one corresponding to the flux from the substrate to the canopy level, r a a . The values of these resistances vary linearly between two limits, ω and 0, corresponding to the closed canopy (LAI > 4) and bare substrate, as expressed by the following equations [22]: with: r a a (ω) = In Equations (13)- (16), z 0 and z 0 represent the roughness lengths for the fully closed canopy and for the substrate, respectively. The value of z 0 is proportional to the crop height by a factor of 0.05 [37,38]. The roughness length for the substrate z 0 is considered constant, with value z 0 = 0.01 m, while the coefficient n for the eddy diffusivity decay coefficient is fixed at 2.5 [22,39]. Thus, it is possible to eliminate the vapor pressure deficit at canopy height D 0 by rewriting Equation (6) in the following form: where: And the coefficients Cc and Cs coefficients in Equation (17) are given by: with: These algebraic equations allow for calculating λE in Equation (17), and the resulting value can be used back in Equation (7) to obtain D 0 and hence the two terms λE c and λE s by means of Equation (6).
In Equations (18) and (19), there are three more resistance terms, namely: the substrate resistance r s s , which regulates the evaporation from the soil and has been considered to vary between 500 (wet soil) and 2000 sm −1 (dry soil) [22]; -the bulk boundary layer resistance: -And the bulk stomatal canopy resistance, r s c , already defined in Equation (5).
In Equation (25), the coefficient 25 represents the mean boundary layer resistance. It has been shown that, due to its limited impact on the calculations, the constant value of 25 can be considered as valid for average conditions [22,39].

Linking Substrate and Canopy Resistance with SWIR Observations in the OPTRAM Approach
The detection of soil water content by means of remote sensing has been intensively investigated by using different sensors and platforms [40], but there are still strong limitations for applications at the plot scale in terms of both spatial and temporal resolution and in accounting for the effects of vegetation cover. The sensitivity of the shortwave infrared signal to vegetation water content has been confirmed, and vegetation water indices have been defined to complement the commonly applied normalized difference vegetation index, NDVI [41]. Recently, Sadeghi et al. [42] have proposed an approach (called "OPTRAM") for assessing the moisture content of the soil and vegetation ensemble based on the utilization of shortwave infrared bands of Landsat-8 and Sentinel-2. Similar to the triangle/trapezoidal method of Carlson et al. [7], which uses the domain of land surface temperature and NDVI, or fractional vegetation cover, for estimating the soil water content, Sadeghi et al. [43] define a relationship between NDVI and an index called "shortwave infrared transformed reflectance" (STR), which is given by: where ρ SWIR is the surface reflectance as detected by Sentinel-2 bands 11 and 12 (1610 and 2190 nm) and Landsat bands (similar to Sentinel-2). The two-dimensional scatter NDVI-STR for a given image (or set of images for the same area) is confined by two edges corresponding to wet and dry soil water conditions. Pixels corresponding to wet surfaces are located on the top edge of the trapezoid (wet edge), whereas pixels corresponding to dry surfaces are located on the bottom edge (dry edge, Figure 1). The boundaries of this domain are defined by plotting the pixel values of NDVI and STR for several dates (i.e., different climatic conditions) over the same area. The slope and intercepts for dry (s d , i d ) and wet (s w , i w ) edges can then be determined and the two limits in the (STR, NDVI) domain are defined by means of the following equations: A "water index" W is then calculated for each pixel from the expression: It has been shown that W, being a proxy of the soil water saturation degree, is strongly correlated with measured soil water content (R 2 > 0.8), especially when using band 12 of Sentinel-2 data [43]. However, in the approach proposed here, we intend to use W as an index to modulate the substrate and canopy resistances r s s and r leaf between wet (W = 1) and dry conditions (W = 0), corresponding to their lower and upper limit. The value of r leaf is kept constant at r leaf,min until the water index is less than a specified threshold W s ; below which, r leaf is increased linearly up to r leaf,max ( Figure 2). Although it would be quite complex to validate the resistance function depicted in the Figure 2, the underlying assumptions are consistent with other similar consolidated parameterizations i.e., for root water uptake [44,45]. The substrate resistance r s s is considered to be inversely related with W, over the entire range from 0 to 1, with limit values of 2000 and 500 sm −1 [22]. The values of r s s and r leaf resulting from this parameterization are inherently linked to the fractional vegetation cover (through NDVI) and the environmental conditions (through STR), and are finally used in the Shuttleworth-Wallace model for calculating E s and E c , as described in the Section 2.1.

Thermal-Based Energy Balance Models
The partitioning between canopy transpiration and soil evaporation has a robust conceptualization in two-source energy balance models i.e., TSEB originally proposed by Norman et al. [11]. TSEB has been intensively applied by using the thermal observations from Landsat [46]. The "two sources" approach, consisting of vegetation and soil layers, takes into account the surface heterogeneity and its influence on radiometric and aerodynamic temperatures. In the TSEB scheme, the sensible heat flux, H, is expressed as the sum of the contribution of the soil, H s , and of the canopy, H c . The vegetation directional fractional cover is used to estimate the soil and canopy temperatures, T S and T C , from the radiometric temperature by considering the viewing angle. By using an initial estimate of canopy transpiration, usually based on the Priestley-Taylor formulation [47], an iterative procedure allows for determining T C and consequently the closure of the surfaced energy balance and the calculation of the instantaneous E. Details about the application of TSEB are given by Kustas and Norman [48,49]. The daily integration is successively carried out by using a scaling factor based on the solar radiation or the evaporative fraction [50].
The accuracy of TSEB evapotranspiration estimates is strongly dependent on the accuracy of the land surface temperature inputs in relation to air temperature boundary conditions. To reduce this sensitivity, the Atmosphere-Land Exchange Inverse (ALEXI) algorithm [12] adopts a time differential correction based on the radiometric temperature from geostationary platforms providing hourly or sub-hourly thermal infrared observations. By means of a disaggregation procedure (DisALEXI), the ALEXI algorithm is extended to apply TSEB with less uncertainty caused by requiring absolute values of the land surface temperature [51] derived from sharpened Landsat-7 and Landsat-8 thermal imagery. More recently, the ALEXI/DisALEXI algorithm has been combined with a data fusion approach [52] for deriving daily maps of evapotranspiration at the spatial resolution of Landsat (30 m), as described by Cammalleri et al. [17,18]; this approach is indicated as Data Fusion in the following sections of this paper.
The TSEB/DisALEXI and the more sophisticated Data Fusion approaches have been applied for mapping evapotranspiration in irrigated vineyards in the context of the GRAPEX experiment (Grape Remote sensing Atmospheric Profile and Evapotranspiration eXperiment) [53,54]. For further details on this application of thermal-based methods in GRAPEX, the reader is referred to [55][56][57][58].
The resulting datasets, hereinafter referred to as "DisALEXI" when applying ALEXI/ DisALEXI for Landsat dates and "Data Fusion (DF)" when creating daily E maps, have been used in the present study for a cross-comparison with the proposed Shuttleworth and Wallace approach that was described in the Section 2.1.

Study Area Description and Datasets Developed in the Context of GRAPEX
The study area is the Ripperdan 720 irrigated vineyard (Vitis vinifera L.) located in the central part of California, U.S. The period of investigation considered here is April-September 2018. The site is included in the GRAPEX experiment conducted by USDA since 2013 with the aim of applying remote sensing based methods for mapping evapotranspiration, crop water use and crop stress for improving the management of irrigation in the Central Valley of California [53,54]. The climate is Mediterranean, with an average temperature between May and October of 22 • C and 12 mm of rainfall.
The  Meteorological and flux datasets for the Ripperdan site used for this study have been collected in the context of GRAPEX by the four identical eddy covariance towers installed at the S-E corner of each plot in order to match the prevailing wind conditions (Figure 3b). Wind direction data collected during the growing season from a nearby vineyard and later validated with the local flux tower wind observations indicated that nearly 95% of the upwind fetch wind was from the northwest quadrant and the upwind extent was typically within 50-75 m. Therefore, the flux towers were deployed in the southeast corner of each block to maximize the fetch and sampling within the block.
The eddy covariance (EC) systems were from Campbell Scientific, Inc. Logan Utah. The IRGASON system, measuring the water vapor/carbon dioxide together with a CSAT3 three-dimensional sonic anemometer to compute latent and sensible heat fluxes, net carbon exchange and momentum flux, were deployed at 4 m above local ground level (a.g.l.) facing due west. A four-way net radiometer (NR01 Net Radiometer Hukseflux, Delft, The Netherlands) was mounted at 4.25 m a.g.l. Within 5 m of the tower, a transect across the interrow of five soil heat flux plates (Radiation and Energy Balance Systems, Seattle, Washington, USA) at 8 cm depth, with soil thermocouples (manufactured by USDA-ARS) at 2 and 6 cm depth and a volumetric water content at 5 cm depth using Stevens Hydraprobe (HydraProbe, Stevens Water Monitoring System, Portland, OR, USA), were used to estimate soil heat flux. See Agam et al. [59] for details on soil heat flux measurements. For details on the eddy covariance flux tower measurements and post-processing, the reader is referred to Alfieri et al. [60].
E observations were available as 30 min averages and were aggregated to daily mean values for the application of the combination equation approaches, Penman-Monteith (P-M) and Shuttleworth-Wallace (S-W).
In addition, half-hourly near-surface (5 cm depth) soil moisture data are available in each block along a transect that comprises five Hydra Probes in the vine row and across the inter-, where we expect a larger reflective response of the surface in the shortwave infrared spectral region.
The cumulative irrigation applied during the study period is shown in Figure 4. There are irrigation events in late May/early June, followed by a lull until mid to late July, with significant irrigation amounts through August and early September. The cumulative volumes applied are approximately 3200 m 3 /ha for Blocks 1 to 3 and 3700 m 3 /ha for Block 4. A lesser amount of irrigation was applied to Blocks 1 and 2 during July, due to trying to induce stress conditions compared to Blocks 3 and 4 during an intensive observation period (IOP) involving the collection of the leaf level gas exchange and leaf water potential data in conjunction with airborne remote sensing imagery. See Kustas et al. [54] for a detailed description of measurements conducted during the IOPs. The measured values of evapotranspiration at the four different vineyard blocks comprising Ripperdan 720 are given in Table 1, from which one observes lower cumulative evapotranspiration at Blocks 1 and 2 as a consequence of the reduced irrigation applied. In the application of the thermal-based methods-DisALEXI and Data Fusionmeteorological data are not used as boundary conditions in the strict sense, as in the combination equation models; as such, they are less sensitive to the accuracy of the meteorological data input compared to the combination equation. For regional applications, the calculation of λE fluxes by means of DisALEXI and Data Fusion thermal approaches is based on global meteorological data provided by the Climate Forecast System Reanalysis (CFSR) of the National Center for Atmospheric Research (NCAR/UCAR).

Description of the Sentinel-2 and Landsat-7-8 Datasets and Derived Products
A main advantage of the P-M and S-W approaches proposed in this study relies on the temporal and pixel resolution of the twin platforms Sentinel-2A and 2B, herein referred to as S-2. During the five-month period from mid-April to mid-September 2018, it was possible to acquire 29 cloud-free images from Sentinel-2A and 2B over the Ripperdan site, included in the tile T10SGF (as identified in the Copernicus system). The Copernicus Open Access Hub provides free download access to Sentinel-2 products.
Top-of-atmosphere reflectance S-2 images (Level 1C) have been downloaded from the data repository of Copernicus (https://scihub.copernicus.eu/, accessed on 8 March 2020). An atmospheric correction for deriving surface reflectance coefficients in each S-2 band has been performed by using the Sen2cor processor developed by the European Space Agency and is freely available (https://step.esa.int/main/third-party-plugins-2/sen2cor/, accessed on 18 March 2020). S-2 images have a geometrical resolution of 10 m in the visible and near infrared bands and a geometrical resolution of 20 m in the red edge and shortwave (Appendix B). This latter resolution is of particular interest when compared with that of the thermal bands of Landsat-8 (100 m). For the processing carried out in this study, we have used the S-2 bands with a resolution of 10 and 20 m, these latter being successively resampled at 10 m resolution by means of the simple nearest neighbor method.
During the same period, the Landsat-7 and -8 (L-7 and L-8) had nine and seven cloudfree overpasses, respectively, along path 43 row 34; two overpass dates for each Landsat occurred on the same date of S-2 (Appendix A). However, the Landsat-7 SLC failure creates gaps over the Ripperdan blocks, limiting the usage of these images without gapfilling. Landsat surface reflectance products have been downloaded from the EROS Science Processing Architecture ESPA (https://espa.cr.usgs.gov/, accessed on 10 January 2020).
Surface reflectance coefficients from S-2 and L-7/L-8 were in good agreement, as shown by the temporal plots of NDVI in Figure 5 and of the hemispherical surface albedo in Figure 6. The values in these graphs represent the average for each of the four treatment blocks of the Ripperdan vineyard.  The albedo estimates for DisALEXI use nominal values of soil and canopy reflectances in the visible and near-infrared based on the model of Campbell and Norman [61]. In the case of Sentinel-2, the albedo has been calculated as the weighted sum of surface spectral reflectances ρ i [62]: where the weighting factors, w λ , are calculated from the mean exo-atmospheric solar irradiance in m different spectral bands (Appendix B) from the following equation: The leaf area index is the most important crop parameter in the calculation of evapotranspiration by means of the combination equation methods P-M and S-W. In the case of Landsat, LAI has been derived by using the algorithm proposed by Gao et al. [63]. For Sentinel-2, LAI values have been derived by using the procedure described in the Sentinel-2 ToolBox embedded in the ESA package SNAP [64]. More specifically, LAI is estimated from the inversion of the radiative transfer model PROSPECT + SAILH based on the Artificial Neural Network, thus taking full advantage of the spectral resolution of S-2. This approach has been validated on different crop types [65]. Table 2 provides a comparison for LAI values on the two dates during the study period, with coincident Landsat-8 and Sentinel-2 acquisitions. On the first date, the L-8 LAI is lower than S-2 LAI, whereas the opposite is true on the second date. The plots in Figure 7 show that the general temporal pattern in LAI from L-8 and S-2 is similar. Ground sampling of LAI was available in three dates but over smaller areas than represented by the 10-30 m pixel resolution of the remote sensing retrievals, so a comparison was not conducted. The decrease in June is a result of both the cover crop in the inter-row being mowed and senescent in late May and early June and July and the vine canopy undergoing structural changes and some management activities, i.e., pruning, as also shown by the NDVI decrease in mid-July. However, the LAI values for L-8 are significantly lower than S-2 in June, but closer in August and September. In general, the S-2 LAI values evidence a larger variability along the season, as well as among the four blocks, with special regard to Blocks 1 and 2 versus 3 and 4, due to the irrigation differences, particularly during June and July. This effect is also due to the different spatial resolution of the two observation datasets. The surface reflectance in Bands 11 and 12 of S-2 has been used for calculating the "shortwave infrared transformed reflectance" (STR) in Equation (21). From a theoretical point of view, the graphs in Figure 8 indicate that the values of STR based on Band 12, with a central wavelength of 2190 nm, has a larger interval-from 1.2 to 4.6-than Band 11; furthermore, there is a larger variability among the different blocks. This indicates a better sensitivity of Band 12 with respect to Band 11 of S-2 for the analyses of the present study. The temporal variation of STR is coherent with the irrigation events shown in Figure 4. We notice that, in correspondence with the significant amount of irrigation applied in August, the value of STR increases, and it decreases when there is a period without irrigation. In addition, it is worth noting that Blocks 1 and 2 correspond to lower irrigation amounts, which is reflected by the lower STR values in the time series plot.

Results
In order to determine the substrate and canopy resistances r s s and r s c by means of the approach described in the Sections 2.1 and 2.2, the NDVI-STR domain has been plotted by considering the pixels included in the S-2 image subset covering a square area of approximately 10 × 10 km 2 (~10 6 pixels) centered on the Ripperdan site (as shown in the left part of Figure 3). We have deliberately chosen to not perform a calibration of the resistances as a function of the water index W, but we instead used the standard values reported in the literature, i.e., Figure 2, which should cover a wide range of possible cases.
The entire set of S-2 images given in Appendix A has been used to produce the plot illustrated in Figure 9. The boundaries of the domain are characteristic for a given area, and they are identified by considering all of the different combinations of NDVI-STR observed over a certain scene portion. For the present study, dry and wet edges have been identified by visual inspection. This does not represent a major drawback of the methodology, since the shape of the pixel distribution on the NDVI-STR space is identical for a given location, at any time, regardless of the surface and meteorological conditions. Pixels excluded above the wet edge represent oversaturated or shadowed pixels. There are issues related to the definition of the dry and wet edges, in relation to the oversaturated pixels, which should be excluded from the NDVI-STR domain [42]. Recent research has proposed alternative ways for defining the boundaries of the NDVI-STR domain, and this is an aspect that may deserve attention for improving the methodology [66]. The values of the intercept and slope of the lines representing the wet and dry edge for the Ripperdan site are indicated in Figure 9, from which, the value of W in Equation (29) has been determined for each date.
The substrate, leaf and canopy resistances r s s , r leaf and r s c , resulting from the adoption of these limits, are plotted in Figures 10-12. In the case of the substrate resistance r s s , the values are linearly derived from W, with limits of 2000 sm −1 (dry edge) and 500 sm −1 (wet edge). In Figure 10, there is variation between 1200-1300 sm −1 until early June. From late July, r s s starts to decrease due to the irrigation water applied, in conjunction with shortwave infrared transformed reflectance STR, as illustrated in Figure 8. The leaf resistance r leaf ( Figure 11) shows a behavior similar to r s s , with a value of 200-240 sm −1 in June, rapidly decreasing from mid-July to its minimum value of 100 sm −1 at the beginning of September, when W is below the threshold of 0.6 ( Figure 2).   The combined effect of LAI and r leaf is depicted in Figure 12, showing the temporal variation of the bulk canopy resistance r s c , derived by means of Equation (5). The peak value of 1000 sm −1 is due to the low foliage density at the beginning of the season, with LAI < 0.5 (Figure 7), followed by a sharp decrease at the end of May, when LAI reaches its maximum.
The graph in Figure 13 shows the temporal variation of the average soil moisture along the transect of Block 1 with the corresponding values of the water index W derived by means of Equation (29). It is possible to notice the decrease in W in correspondence with the reduction in surface soil moisture during June, and the rapid increase with the intensification of irrigation applications from the second half of July. The value of W in this last period is much higher than during the earlier months (April), when a higher soil moisture content is recorded. A similar behavior has been observed in all four blocks. This may appear controversial, but it is explained instead by the fact that, starting from July, the vine canopy is fully developed, maintained continuously and well hydrated by the irrigation applications. The spectral response in the shortwave infrared is dependent not only on the soil but also on the vegetation water content. Hence, the water index W represents the water status of the soil-canopy ensemble and, as such, it provides useful indications for modulating surface resistances in the combination equation when vegetation is present. The calculation of E by means of the combination equation approaches, i.e., P-M and S-W S-2, has been compared with E obs from the flux towers measurements. The calculations are made using either the meteorological input from the locally installed flux towers, or using the meteorological data provided by the Climate Forecast System Reanalysis (CFSR) of the National Center for Atmospheric Research (NCAR/UCAR), similarly to DisALEXI and Data Fusion. The E values considered for the comparison correspond to pixels falling approximately within the footprint of the flux tower measurements in the N-W direction of the prevailing wind.
The results of this first comparison are summarized in the plot of Figure 14, with the corresponding statistics in Table 3. The statistical indicators have been calculated by considering the average values over the estimated footprint of the flux towers for the dates of S-2 acquisition given in Appendix A. In addition to Pearson and the determination coefficient R 2 , the root mean square error and the mean absolute error have been calculated as follows: Figure 14. Correlation between the daily evapotranspiration E measured by eddy covariance towers (E tower) and the estimates (E models) obtained by means of Penman-Monteith (P-M S-2, cross) and Shuttleworth-Wallace (S-W S-2, circles), with flux tower meteorological data and surface parameters and resistances derived from Sentinel-2. Table 3. Statistical results between daily E ( Figure 14) from P-M S-2 and S-W S-2 calculated with two different meteorological input data vs. tower measurements.  Table 3 also gives the F-statistics for the corresponding degree of freedom (n-1). The thermal-based approaches described in Section 2.3, namely DisALEXI and Data Fusion, have been applied by using Landsat-7 and -8 acquisitions (see Appendix A) and have also been compared with the measured fluxes at each individual block of the Ripperdan vineyard. The results are summarized in Table 4 and Figure 15, which are including also the S-W S-2 shown above, but disaggregated for the four management blocks.  The dates for the S-W S-2 and Data Fusion results are coincidental, with the exception of May 27, which is missing from the Data Fusion; the dates for DisALEXI are related to the Landsat overpasses. Hence, there is a difference in the number of samples for each one of the three data sets, as indicated by the degree of freedom given in Table 4.

Discussion
From the results shown in Figure 14 and Table 3, it is possible to conclude that both combination equation methods yield a reasonable agreement with the measured values of daily E. This finding confirms that the resistance modulation with OPTRAM reproduces the observed process with a satisfactory accuracy. However, although the correlation coefficients and the errors are similar for the two methods based on the combination equation, we notice that the slope of the regression is close to 1 for S-W S-2 and significantly lower for P-M S-2. This behavior was expected, considering that soil evaporation is not adequately represented in the "big-leaf" concept of the Penman-Monteith equation.
In the comparison between the Shuttleworth-Wallace and Penman-Monteith (P-M) approaches, it is interesting to observe the ratio of E S-W S-2 and E P-M plotted in Figure 16. In this graph, the continuous lines represent the variation of the ratio between E S-W S-2 and E P-M with LAI. The calculation is carried out with the two models by considering an average climatic day during the irrigation season at Ripperdan (T air = 23.6 • C; RH% = 45%; u = 1.6 ms −1 , K ↓ = 315 Wm −2 ), a value of leaf resistance r leaf of 200 sm −1 and four different values for the substrate resistance r s s . The circles are the actual values obtained for the Ripperdan, with the approach described in previous sections. As expected, the Shuttleworth and Wallace approach gives higher E for low LAI values, and it quickly converges toward Penman-Monteith for LAI > 1. The ratio S-W/P-M is greater than 1 for LAI < 1, depending on the substrate resistance r s s (inversely related to soil evaporation); the influence of r s s decreases significantly and the four lines reach an asymptotical value of the ratio, which is 1.06 in this specific case. However, for LAI = 2, the evapotranspiration estimates with Shuttleworth and Wallace are still larger than Penman-Monteith within a range of 5-20%. The S-W S-2 approach provides reliable estimates of E also in the presence of moderate water stress. During July and August, it is possible to notice a significant differentiation between the four treatment blocks, as a consequence of the different irrigation amounts applied. It is interesting to analyze the plots of Figure 12 (bulk canopy stomatal resistance) in parallel with the cumulative irrigation applied in Figure 4. Block 3, which received fairly constant irrigation, falls in the lower part of the plot of Figure 12, followed by Block 4, which has the highest total amount of irrigation but with only a minor increase during June and July. Blocks 1 and 2, where there was minimal irrigation during June and July, have higher values of r s c in response to mild hydric stress conditions. Blocks 2 and 3 correspond to the lower and upper lines in the plot of Figure 12, thus representing the maximum and minimum E modeled by S-W S-2. The corresponding values of E are represented in Figure 17 together with the daily measurements from the eddy covariance stations. The temporal behavior of Blocks 2 and 3 in the observed and modeled daily E are in good agreement with each other, with the exception of early season and end of season, around bud break and harvest in April and September, where the S-W S-2 approach is under and overestimating the observed E, respectively. This might be explained by an underestimation of the cover crop evapotranspiration in April and an overestimation of vine transpiration in September as vines begin senescing. In Figure 18, the values of the ratio between E from S-W S-2 and the reference evapotranspiration are plotted for the four Ripperdan blocks. In vineyard operations for this region, ratio values below the threshold of 0.7 (dashed line in Figure 17) are considered as an indication of water stress. During the months of July and August, Blocks 1 (diamonds) and 2 (squares) are below this threshold, as a consequence of less irrigation water applied. The reduction of E is consistent with the model predicted increase in leaf/canopy resistance for these two blocks, as shown in Figure 12. Based on the discussion above, only S-W S-2 has been considered in the comparison with thermal-based methods. From the results shown in Figure 14 and Table 4, we notice that-among the three models, S-W S-2, DisALEXI and Data Fusion-the best results are obtained by means of DisALEXI on Landsat overpass dates, corresponding to the lowest values of RMSE and MAE, even with the lowest number of data points. In all cases, the slope of the correlation is around 1, except for Blocks 1 and 2 for the Data Fusion approach, where E tends to be overestimated. The best agreement between the S-W S-2 approach and the measured fluxes in terms of correlation coefficients is observed in Block 4, whilst the lowest errors are found in Block 2.
As expected, the quality of input data affects the accuracy of results of the combination equation method S-W S-2 in all of the four blocks, as shown by the statistical indicators in Table 4. When S-W S-2 is applied by using the meteorological data from the global CFSR dataset, its performance is lower in comparison to thermal-based methods, especially in Blocks 1 and 3.

Conclusions
Based on the comparative analysis presented in this study, the S-W S-2 approach for calculating E, including the simple parameterization for the canopy and substrate resistances corrected with shortwave information, provides reasonable agreement with flux observations obtained in a variably stressed vineyard. The OPTRAM approach is able to follow the temporal evolution of the water status in the soil-vegetation ensemble, thus allowing us to modulate the resistances for the calculation of E within upper and lower limits values, which are typical for irrigated crops. The proposed method is consistent with the classical Penman-Monteith approach, with fractional vegetation cover and LAI ranging up to a fully closed canopy. Although DisALEXI outperforms S-W S-2, the main advantage of this latter approach is the spatial resolution of Sentinel-2 data (10 m) and the temporal resolution of 3-5 days. However, new advances in thermal sharpening techniques provide a way to attain the visible and near-infrared resolutions [65,66] and are being used in combining Sentinel-2 data with Sentinel-3 thermal observations to create a 20 m resolution LST product [21].
The S-W S-2 approach runs with daily averages of meteorological forcing and, therefore, does not need upscaling from instantaneous to daily values, as in the thermal-based approaches. However, this modeling approach is most accurate using local meteorological data. Using non-local weather station data leads to additional errors, which have been quantified, whereas such non-local meteorological forcing does not have a similar effect on DisALEXI results. Clearly, there are advantages and disadvantages of both thermal and shortwave-based approaches, and efforts to integrate these modeling schemes by leveraging the information and intermediate output provided by both would improve their utility in computing reliable daily E products at 10-20 m resolutions.
Further analyses are required for refining the definition of the STR-NDVI domain and the determination of both the wet and dry edges required for coefficients in Equations (27) and (28) and how these impact the model estimates of E using the S-W S-2 approach. The addition of thermal data and model output from DisALEXI to help to determine the dry and wet edges and/or better estimate the soil and canopy resistances will be explored. The availability of algorithms for deriving added-value products from S-2 data, such as LAI, the sharpened S-3 LST product and gridded meteorological data allow for the development of automated processing tools that can lead to operational applications for irrigation management and precision agriculture applications.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The satellite data used in this study are publicly available at the Copernicus Open Access Hub (https://scihub.copernicus.eu/). The CIMIS data is available after registering at this web site https://cimis.water.ca.gov/Default.aspxat. The tower flux data used in this study are available on request from the co-author W.P.Kustas (bill.kustas@usda.gov). These data are not publicly available due to research being currently conducted by other members of the GRAPEX research team.