Evaluation of Penman-Monteith Model Based on Sentinel-2 Data for the Estimation of Actual Evapotranspiration in Vineyards

Water scarcity is one of the most important problems of agroecosystems in Mediterranean and semiarid areas, especially for species such as vineyards that largely depend on irrigation. Actual evapotranspiration (ET) is a variable that represents water consumption of a crop, integrating climate and biophysical variables. Actual evapotranspiration models based on remote sensing data from visible bands of Sentinel-2, including Penman-Monteith–Stewart (RS-PMS) and Penman-Monteith– Leuning (RS-PML), were evaluated at different temporal scales in a Cabernet Sauvignon vineyard (Vitis vinifera L.) located in central Chile, and their performance compared with independent ET measurements from an eddy covariance system (EC) and outputs from models based on thermal infrared data from Landsat 7 and Landsat 8, such as Mapping EvapoTranspiration with high Resolution and Internalized Calibration (METRIC) and Priestley–Taylor Two-Source Model (TSEB-PT). The RS-PMS model showed the best goodness of fit for all temporal scales evaluated, especially at instantaneous and daily ET, with root mean squared error (RMSE) of 28.9 Wm−2 and 0.52 mm day−1, respectively, and Willmott agreement index (d1) values of 0.77 at instantaneous scale and 0.7 at daily scale. Additionally, both approaches of RS-PM model were evaluated incorporating a soil evaporation estimation method, one considering the soil water content (fSWC) and the other hand, using the ratio of accumulated precipitation and equivalent evaporation (fZhang), achieving the best fit at instantaneous scale for RS-PMS fSWC method with relative root mean squared error (%RMSE) of 15.2% in comparison to 58.8% of fZhang. Finally, the relevance of the RS-PMS model was highlighted in the assessment and monitoring of vineyard drip irrigation in terms of crop coefficient (Kc) estimation, which is one of the methods commonly used in irrigation planning, yielding a comparable Kc to the one obtained by the EC tower with a bias around 9%.


Introduction
Water is becoming increasingly scarce, especially for agricultural activities, and its demand for agricultural purposes is projected to increase 60% by 2050 [1], especially in irrigated crops [1][2][3][4][5]. Vineyards (Vitis vinifera L.) located in arid, semi-arid and Mediterranean climates are especially sensitive to changes in water availability which affects not only yield but also the quality of wines [6][7][8][9][10]. A key element for water management in vineyards is the accurate measurement of water use via evapotranspiration (ET) both in time and Several approaches have been proposed to improve RS-PM model, focusing on the estimation of surface conductance. In this regard, the vegetation conductance has been modeled from the LAI estimated with remote sensing data, using initially Moderate Resolution Imaging Spectroradiometer (MODIS) data that have a pixel resolution of 1000 m [40]. Despite the good agreement on magnitudes and seasonal variation in both ecosystems, the method is limited to areas where the LAI is greater than 2, since it is assumed that the canopy cover a large proportion of the soil and thus the surface conductance (g s ) is approximate to the canopy conductance (g c ).
Likewise, Leuning and collaborators [41] introduced improvements to the previous models regarding the modeling of g c and soil evaporation (RS-PML). Thus, g c is modeled as a function of maximum leaf stomatal conductance at the top of the canopy (g sx ), LAI, D a and the hyperbolic response to absorbed short-wave radiation [42], which is easily parameterized and coupled with a photosynthesis function, but does not respond to water stress, an important factor in irrigated vineyards under Mediterranean climates. Additionally, soil evaporation is assumed as a constant fraction (f ) of the equilibrium evaporation rate of soil surface, although that should be adjusted according to soil water content. However, RS-PML had good agreement when assessed in 15 sites at a spatial resolution of 1000 m, including vegetation cover of savannah, annual crops, and diverse forest classes. Authors modeled daily f as a function of soil water content [43,44], which improved the estimation of soil evaporation and thus the estimation of total ET, especially in dryland sparse vegetation.
Based on the aforementioned, modeling ET based on the Penman-Monteith equation and Sentinel-2 data in drip irrigated vineyards, the RS-PM approach should be improved incorporating soil water content data as one of the environmental variables that controls both stomatal conductance and the soil evaporation fraction. Thus, the model of stomal conductance proposed by Stewart [45] is an empirical approach based on a multiplicative environmental function from solar radiation, air temperature, vapor pressure deficit and soil water deficit, which is suitable to be coupled with the RS-PM model, hereinafter referred as Remote Sensing Penman-Monteith-Stewart model (RS-PMS). Therefore, it is hypothesized that RS-Penman-Monteith model integrated with Sentinel-2 data and coupled to a function that adjusts the stomatal conductance and soil evaporative fraction according to soil water content (RS-PMS), improves the goodness of fit of modeled actual ET. Hence, the general objective of this research was to evaluate the Penman-Monteith model for remote sensing data, based on the Leuning approach (RS-PML) and a novel approach that incorporates the Stewart stomatal conductance model (RS-PMS), on a drip irrigated vineyard located in central Chile by comparing independent data from an eddy covariance system and the METRIC and TSEB-PT models.

Study Area
The study was carried on in a drip irrigated vineyard (Vitis vinifera L. cv. Cabernet Sauvignon) during the 2017-2018 (S 1 ) and 2018-2019 (S 2 ) growing seasons. The Vineyard is in central Chile, 40 km south of the city of Santiago. The area is characterized by Mediterranean climate with annual average temperatures of 18 • C; precipitations are concentrated in winter (June-September) with an average annual total of 280 mm. The soil is characterized by Mollisol order with predominantly loamy texture within the first 0.6 m and sandy loam between 0.6 and 1 m of depth; bulk density varies between 1.38 and 1.48 g cm −3 within one meter of depth.
The vineyard was planted in 2010, rows oriented north-south with 2.5 m between rows and 1.0 m between vines. The aisles are maintained vegetation-free using mechanical and chemical weed controls; fertilization management is used to avoid nutritional deficiencies of both macro and micro elements. The amount of water applied by drip irrigation during the season was calculated as a function of reference crop evapotranspiration (ET 0 ), applying 50% of ET 0 (ET 0 × 0.5) with frequencies of 7 days throughout the season and an irrigation Remote Sens. 2021, 13, 478 4 of 25 time that varied between 8 and 10 h. The applied irrigation corresponds to the gross water demand and the total amount irrigated was measured with a flow meter, applying 481 mm in S 1 and 334 in S 2 . Additionally, the irrigation system has a single irrigation line, with two emitters per plant spaced at 0.4 m and a flow rate of 4.1 L s −1 per emitter. In each irrigation event, wetted area uniformity was observed with some cases where superficial water runoff was evidenced, which corresponded to less than 5% of the emitters. Canopy management was different in the two growing seasons. In S 1 the trellis system was vertical-shoot-positioned with three wire lines, while in S 2 the training system was without wire lines, increasing the frequency of pruning during the growing season.

Field Measurements
Ground and remote sensing data are the input of both VIS-ET and TIR-ET models (Figure 1), which are described in the following sections. Field measurements included micrometeorological, phenological and fraction of radiation transmitted by canopy. The micrometeorological data were obtained through an eddy covariance system (EC) measuring exchanges of energy and mass between the vineyard and the atmosphere [46]. Due to the prevailing wind direction during daylight hours, a west-facing EC tower was installed in the east border of the study area (  The EC was equipped with an open-path infrared gas analyzer integrated with a 3D sonic anemometer model IRGASON (Campbell Scientific Inc., Logan, UT, USA), placed 4 m above the ground surface and sampling frequency of 10 Hz. Other measuring devices in place were temperature and relative humidity HMP45C sensors (Campbell Scientific Inc., Logan, UT, USA); photosynthetically active radiation (PAR) model LI190SB (Campbell Scientific Inc., Logan, UT, USA) device; net radiometer model NR Lite2 (Kipp and Zonen B.V., The Netherlands), placed 0.5 m above the top of the vineyard canopy, which is representative of both the bare soil surface and the vegetative surface covered by the vine; two soil heat plates model HFP01-L (Campbell Scientific Inc., Logan, UT, USA), with 0.9 m of separation and 0.08 m below the soil surface. Additionally, there were two thermocouples model TCAV (Campbell Scientific Inc., Logan, UT, USA) at 0.05 m of depth and one Time Domain Reflectometry model CS616 (Campbell Scientific Inc., Logan, UT, USA) at 10 cm of depth. The data from these sensors were processed by a CR3000 datalogger (Campbell Scientific Inc., Logan, UT, USA). Traditional EC raw data quality controls were performed [47,48], including spike removal, double coordinate rotation [49], correction for frequency loss [50], air density [51] and turbulent and stationarity conditions [52]. Additionally, an analysis of energy balance closure was made, which means a linear regression between the sum of latent and sensitive heat flux (λE + H) and the subtraction between net radiation and soil heat flux (R N − G) or available energy (A) [53]. The EC data analysis was synthesized in 30-min averages and the climatological flux footprint was parameterized based on the crosswind distribution [54] in order to overlap with remote sensing data and evaluate model outputs, considering the footprint contour encompass 90% of the total flux [55][56][57].
The EC footprint was divided according to Landsat 8 pixels (green dots in Figure 2) and vines located within the pixels were randomly selected, in order to observe phenology and measure fraction of radiation transmitted by canopy (f TPAR ). Phenological observations were made every 7 days during S 1 and S 2 according to BBCH-scale [58], which assigns a decimal code to each phenological stage ranging from 00 to 99 and in vineyards were identified the beginning of Leaf development (BBCH-11), Inflorescence emergence (BBCH-53), Flowering (BBCH-61), Development of fruits (BBCH-71) and Ripening of berries (BBCH-81). The criterion used to identify a phenological stage was when 50% of the vines had reached the next phenological stage (Table 1).
Finally, f TPAR was measured at 20 points within the EC footprint as was made with phenological observation. Those points are obtained from samples made across the row in three consecutive vines each spaced by 0.3 m with a LAIPen LP 100 sensor (Photon Systems Instruments, Czech Republic), which has been used in forest canopy characterization [59][60][61]. The LAIPen LP 100 measures the gap fraction of blue band of the photosynthetically active radiation (0.4-0.5 µm), which is the most efficiently absorbed spectral band by leaf pigments and similar to LAI-2000 (LI-COR, Lincoln, NE, USA), removes radiation with wavelengths higher than 0.5 µm at which scattering is more pronounced. All measurements were made before 11:00 am in order to reduce the influence of high zenith angles. Moreover, assuming that vineyard canopy is a homogeneous turbid medium, where leaves are randomly distributed within the canopy, and individual leaf size is small when is compared with the canopy, one LAIPen LP 100 equipment measured below the canopy (R below ) and simultaneously a second equipment measured on top of the canopy (R top ).
With these measurements, f TPAR is calculated as the ratio R below /R top and applying a light extinction model [62] is estimated the Leaf Area Index (LAI E ) as an indirect measurement of LAI: where k is canopy extinction coefficient, set in a value of 0.5 in previous measurements on a Cabernet Sauvignon vineyard [63]. Then, LAI E was adjusted (LAI EA ) assuming that vineyard canopy is composed of 82% of wood (trunk, cordons and canes) by volume of surface area [64], and by inter-row space factor of 0.40, which was obtained from the ratio of vines space (1.0 m) and inter-row space (2.5 m) [65]:

Remote Sensing Data
Remote sensing data was obtained from Landsat 7, Landsat 8 and Sentinel-2 missions, with each mission having different temporal, spatial, spectral and radiometric resolution ( Table 2). Landsat 7 and Landsat 8 (Landsat 7-8) have a 16-day repeat cycle, both within visible (σ blue , σ green , σ red ) and near-infrared (σ nir ) bands have pixel resolution of 30 m, while within thermal infrared band (σ tir ) Landsat 7 has a pixel resolution of 60 m and Landasat 8 of 100 m, with a resampling to 30 m [66]. Sentinel-2 visible and near-infrared bands have a spatial resolution of 10 m with temporal frequency of 5 days at latitudes near to equator and 2-3 days at mid latitudes [67]. The remote sensing data were filtered by date during the growing season and by percentage of cloudiness lower than 10% for Landsat 7-8 and lower than 30% for Sentinel-2. In Landsat 7-8 a lower percentage of cloudiness was used, since this data source uses the TIR band for Land Surface Temperature (LST) estimation, which generates a large amount of noise when the percentage of cloudiness is higher than 10%. The Normalized Difference Vegetation Index (NDVI) was calculated from Landsat 7-8 and Sentinel-2 data (Table 3), in order to model the Leaf Area Index (LAI NDVI ). Previously the noise associated with NDVI was removed through Savitzky-Golay's algorithm [68] and therefore, a regression analysis was performed between NDVI and LAI EA to model LAI NDVI for each set of remote sensing data. In addition, the fraction of vegetation cover (f C ) was calculated as a function of LAI NDVI and the viewing angle of the sensor (θ) [69]: The LST was calculated from thermal bands of Landsat 7 (10.4-12.5 µm) and Landsat 8 (10.6-11.2 µm) according to the Single Channel Algorithm method (SC) [70,71], which is appropriate when one single thermal band of the sensor is available. Furthermore, SC method is based on surface emissivity (ε S ) and radiance (L sen ), brightness temperature (T sen ), transmissivity (τ) and downward (L d ) and upward (L u ) radiance as atmospheric variables. The atmospheric variables were obtained from the National Center for Environmental Prediction (NCEP), that is a global tropospheric analyses product providing atmospheric profiles every 6h at resolution of 1 • × 1 • [72], and it is available in the data catalog of Google Earth Engine (https://developers.google.com/earth-engine/datasets/ catalog/NCEP_RE_surface_temp#description).
The calculation of NDVI (https://code.earthengine.google.com/2f17d6fd8809f2ca1 a21448c1c61b653) and LST (https://code.earthengine.google.com/4c0f3d0af6df02bf6a8 8ba96e046f711) was supported by Google Earth Engine (GEE), a cloud-based platform for spatial data analysis [73]. GEE is based on JavaScript code for geospatial analysis and in addition provides access to several remote sensing data bases including Landsat 7, Landsat-8, Sentinel-2 and NCEP, which allow the data assimilation process with minimum computational capabilities.
The overlapping of remote sensing grid and EC climatological flux footprint were grouped in 22 pixels for Landsat 7-8 and 108 pixels for Sentinel-2, which are averaged in order to compare with 90% source area of EC climatological flux footprint [74].

Remote Sensing-Based ET Models
The Remote Sensing Penman-Monteith-Stewart model (RS-PMS) is evaluated and compared with the Remote Sensing Penman-Monteith-Leuning (RS-PML), the Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC) and the Priestley-Taylor Two Source Energy Balance (TSEB-PT). Table 4 summarizes the main features of evaluated models. Energy balance residual from H S as a function of T S , T AC , r S .

RS-PML
Leuning's approach [42] as a function of g C , g SX , LAI NDVI , D a .

•
One parameter g C model.
Four parameters g C model.

Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC)
The instantaneous λE (Wm −2 ) is calculated in METRIC [18] as a residual of the energy balance equation: where R N is net radiation (Wm −2 ); α the surface albedo; ε s the surface emissivity; R S the solar radiation (Wm −2 ); R Ld the downward longwave radiation (Wm −2 ) determined empirically by air transmissivity; R Lu is upward longwave radiation (Wm −2 ) obtained by thermal band. Therefore, shortwave and longwave net radiation are calculated following the approach suggested by [75]. The soil heat flux (G) is calculated as a fraction of R N based on field measurements, and H (W m −2 ) is calculated following the equations: where ρ a is air density (kg m −3 ); C p air specific heat at constant pressure (J kg −1 K −1 ); ∆T vertical temperature difference near the surface (K); r ah (s m −1 ) aerodynamic resistance corresponding to ∆T; and c and m are parameters empirically determined in the image. ∆T and r ah are derived from the information of two selected pixels in the evaluated scene or image, called the "cold pixel" and "hot pixel", which refer to a wet and dry surface, respectively. The λE for each of these pixels is assigned by operators, in order to derive H using Equation (4). Thus, the combined values of ∆T and r ah are iteratively obtained assuming an initial condition of neutral stability and subsequently is corrected by atmospheric stability conditions, which depends on the length of Monin-Obukhov for both the wind speed at satellite overpass and the surface roughness length [18,76,77]. As in METRIC, which estimates λE as residual of surface energy balance (Equation (4)), TSEB-PT estimates resistances of both canopy and soil surfaces, splitting the energy balance components between canopy and soil [78]: where instantaneous fluxes are expressed in Wm −2 , and subscripts refer to canopy (C) and soil (S). Here G is a fraction of R N obtained from EC measurements. Therefore, total (H), canopy (H C ) and soil (H S ) heat fluxes are given by: where T A is air temperature (K); T AC is the temperature in the air-canopy space (K); T S is soil surface temperature; r A aerodynamic resistance (s m −1 ) to heat transfer through the canopy-surface interface; r X total resistance of the canopy boundary layer (s m −1 ); and r S the soil surface boundary layer resistance (s m −1 ). Since T A is an input variable that comes from micrometeorological data; r A , r X , r S depend on canopy properties and their distribution is assumed as a serial configuration, because it is more robust for several environmental conditions and consistent with the turbulence model that is generated in vineyards inter-row space [79]; finally T C , T S and T AC are iteratively calculated considering LST from thermal band of Landsat 7-8, and assuming that fg transpires a rate based on Priestley-Taylor equation, hence instantaneous λE C is given by [80]: where α PT is the Priestley-Taylor coefficient [81] set to 1.26, γ is the psychrometric constant and s is the slope of the curve relating saturation water vapor pressure and temperature.

Remote Sensing Penman-Monteith (RS-PM)
An evaluation was made of two Penman-Monteith model approaches for remote sensing data (RS-PM). Therefore, both models are based on Penman-Monteith (PM) formulation, which takes LAI as a function of NDVI (LAI NDVI ) to estimate the available energy for latent heat flux. Thus, applying the PM equation separately to canopy (λE C ) and soil (λE S ), and assuming that λE S occurs as a fraction (f ) of the soil surface evaporation equilibrium rate [41]: where ξ = s/γ, γ is the psychrometric constant and s the slope of the curve relating saturation water vapor pressure to temperature; g a is aerodynamic conductance (m s −1 ); g c is canopy conductance (m s −1 ); D a is water vapor pressure deficit of the air (kPa) and f is a factor that modulates potential evaporation rate at the soil surface, expressed by the Priestley-Taylor equation where f = 0 when the soil is dry, and f = 1 when the soil is completely wet. In this study two approaches for the calculation of f are evaluated, on the one hand considering the ratio between the accumulated daily precipitation (P i ) and daily soil equilibrium evaporation rate for day i (E eq,s,i ) of 30 previous days [44] (f Zhang ) and on the other hand the relationship among measured soil water content (θ m ), soil water content at field capacity (θ fc ) and permanent wilting point (θ wp ) (f SWC ) [43]: The available energy absorbed by the canopy A C (Wm −2 ) and soil surface A S (Wm −2 ) depends on both available surface energy (A = R N -G) and τ as a function of LAI NDVI and the extinction coefficient (k) set in 0.5: R N (Wm −2 ) was calculated from the shortwave and longwave energy balance, which albedo (α) was estimated from Sentinel-2 data [82]; atmospheric emissivity (ε a ) was estimated from air temperature [83] and surface emissivity was considered 0.94 from ground measurements made in vineyards [84]: As in METRIC and TSEB-PT models, G is estimated as a fraction of R N based on EC data.
Moreover, RS-PML approach follows the canopy conductance model [85]: (22) where D 50 is the vapor pressure deficit in which stomatal conductance is reduced to half its maximum value; g sx is the maximum stomatal conductance at the top of the canopy; Q h is the photosynthetically active radiation (PAR) flux density at the top of the canopy; k Q , is the extinction coefficient of shortwave radiation and Q 50 is the PAR flow when stomatal conductance is half its maximum value. It should be noted that after a model sensitivity analysis [41], the parameters were set: Q 50 = 30 Wm −2 , D 50 = 0.7 kPa and k Q = 0.6. Moreover, Equation (20) is very sensitive to parameter g sx , which is specific to species and cultivar, and in this research is set as 0.0133 m s −1 [86][87][88].
The RS-PMS model is based on Equations (14)- (21) and canopy conductance (g c ) is modeled as follows [45]: where f 1 , f 2 and f 3 are functions of measured solar radiation (Rs), water vapor pressure deficit (D a ) and soil water content, which require the parameterization of k 1 and k 3 , as well as g sx . These parameters were calibrated with the even days of S 1 and S 2 and were evaluated with the odd days of those two seasons.

Daily, 7-Day and Monthly ET Extrapolation
All remote sensing-based ET models estimate instantaneous values of energy balance components in Wm −2 units. Therefore, to estimate daily ET, instantaneous λE is extrapolated using observed λE and solar radiation (R S ) ratio at time satellite overpass [30]: where ET d is total daily ET (mm d −1 ); λE so /Rs so is the ratio of latent heat to solar radiation at time satellite overpass; Rs 24 (MJ m −2 d −1 ) is average daily solar radiation; and λ is the latent heat of vaporization. Additionally, the Rs 24 ratio approach was used to extrapolate ET to a 7-day and monthly scale.

Model Assessment
The metrics to compare and evaluate the performance of remote sensing-based ET models were: • Bias represents the average difference between the observed and modeled values, so the ideal model gives the lowest value: where P i is the predicted value, corresponding to the average of pixels within the EC footprint, and O i to EC measurement. • Root mean squared error (RMSE): • Relative root mean squared error (%RMSE) which is dimensionless and expresses the error as a fraction of the measured average (O avg ), thus the model with the best performance is the one with the lowest %RMSE: • Willmott agreement index, which is a skill score that involves the variability of the observed and modeled values. Therefore, the model is perfect if P i = O i and consequently the index = 1. Conversely, if P i = O avg the index = 0 [89]:

Vineyards Seasonal Observations and Measurements
The linear regression analysis of the energy balance closure for S 1 and S 2 (Figure 3) shows that both seasons are similar in terms of slope, intercept, and coefficient of determination. The slope varies between 0.76 and 0.82, intercept around 16 Wm −2 and the coefficient of determination between 0.93 and 0.95. The pattern of LAI EA per unit of ground area, estimated from measurements of radiation transmitted by the canopy is similar during S 1 and S 2 ( Figure 5) with an initial phase of rapid growth between November and December, covering leaf development (LD) and inflorescence emergence (IE) stages, followed by a second phase where maximum values are reached during January and early February, including flowering (FW) and development of fruits (DF) stages, and finally a third phase where the LAI EA is slowly reduced from ripening of berries (RB) stage (Veraison) to the end of season (late February to April). Additionally, LAI EA values in S 2 were lower than S 1 , mainly from flowering to ripening of berries stages, possibly due to canopy management that included higher frequency of pruning in S 2 . The pattern of LAIEA per unit of ground area, estimated from measurements of radiation transmitted by the canopy is similar during S1 and S2 ( Figure 5) with an initial phase of rapid growth between November and December, covering leaf development (LD) and inflorescence emergence (IE) stages, followed by a second phase where maximum values are reached during January and early February, including flowering (FW) and development of fruits (DF) stages, and finally a third phase where the LAIEA is slowly reduced from ripening of berries (RB) stage (Veraison) to the end of season (late February to April). Additionally, LAIEA values in S2 were lower than S1, mainly from flowering to ripening of berries stages, possibly due to canopy management that included higher frequency of pruning in S2.

Biophysical Variables Based on Remote Sensing Data
The NDVI-LAI EA regression analysis using Landsat 7-8 and Sentinel-2 data are shown in Figure 6a,b respectively. In both data sets the best fit were obtained using a power function with R 2 values of 0.63 from Landsat 7-8 and 0.82 from Sentinel-2 data, suggesting that Sentinel-2 data provide better estimations of LAI NDVI and canopy properties than Landsat 7-8 data.

Biophysical Variables Based on Remote Sensing Data.
The NDVI-LAIEA regression analysis using Landsat 7-8 and Sentinel-2 data are shown in Figure 6a and 6b, respectively. In both data sets the best fit were obtained using a power function with R 2 values of 0.63 from Landsat 7-8 and 0.82 from Sentinel-2 data, suggesting that Sentinel-2 data provide better estimations of LAINDVI and canopy properties than Landsat 7-8 data.    was around 2.1 m 2 m -2 in S1 and 1.5 m 2 m -2 in S2. However, despite the differences in magnitude and %Bias, in both data sets the pattern is similar, reproducing a three-phase trend.

Assessment of Remote Sensing-Based ET Models.
Linear regression analysis of the components of the energy balance (Figure 9), shows that modeled fluxes from RS-PMS and RS-PML are the most closely grouped around the 1:1 line. However, the RS-PMS model has the best fit of instantaneous sensible (H) and latent (λE) heat fluxes with R 2 of 0.53 and 0.79, respectively. On the other hand, when comparing TIR-ET models, TSEB-PT gives the best fit of modeled λE with R 2 of 0.6, while METRIC exhibits a low R 2 of 0.3. However, comparing these models in terms of RMSE, %RMSE and d1 the difference between them is not so pronounced, having the TSEB-PT model a slightly better performance than METRIC with values of 55.

Assessment of Remote Sensing-Based ET Models
Linear regression analysis of the components of the energy balance (Figure 9), shows that modeled fluxes from RS-PMS and RS-PML are the most closely grouped around the 1:1 line. However, the RS-PMS model has the best fit of instantaneous sensible (H) and latent (λE) heat fluxes with R 2 of 0.53 and 0.79, respectively. On the other hand, when comparing TIR-ET models, TSEB-PT gives the best fit of modeled λE with R 2 of 0.6, while METRIC exhibits a low R 2 of 0.3. However, comparing these models in terms of RMSE, %RMSE and d 1 the difference between them is not so pronounced, having the TSEB-PT model a slightly better performance than METRIC with values of 55. Models based on the Penman-Monteith equation (RS-PML and RS-PMS) were evaluated with two approaches to estimate the soil evaporative fraction (fZhang and fSWC) as presented in equations 16 and 17 of section 2.4.3. In this regard, for both RS-PML and RS-PMS the approach that includes soil water content measurements (fSWC) yielded the best goodness of fit in the estimation of instantaneous λE (Table 5), where a %RMSE of 21.3% (RS-PML) and 15.2% (RS-PMS), while fZhang approach results a %RMSE of 32.2% for RS-PML and 58.8% for RS-PMS (data not shown). Likewise, the Wilmott agreement index (d1) was higher under the fSWC approach (0.65-0.77) than the fZhang approach (0.52-0.33).
The seasonal pattern of λE modeled by the TIR-ET models and compared with λE measured by the eddy covariance (λE EC) (Figure 7), show that both METRIC and TSEB-PT replicate very well the intra-seasonal variability of λE EC, however, METRIC underestimates λE in S1 (−10.2 Wm −2 ) and overestimates in S2 (16.2 Wm −2 ), while TSEB-PT underestimates both seasons (−17.6 and 1-2.7 Wm −2 ). Therefore, Figures 7g and 7h show that METRIC tends to underestimate λE in early season (November-December), while the estimation is closer to measured λE between January and March. On the other hand, Figures 7g-h show that the TSEB-PT model underestimates λE when H is overestimated and, vice versa, λE is overestimated when H is overestimated, suggesting a strong dependence of λE with the remaining components of the energy balance, particularly H.
Assessing the intra-seasonal variation of λE modeled by RS-PML and RS-PMS (Figures 8c-d), both models successfully replicate the seasonal variability of λE, especially RS-PMS with a bias of −10.8 Wm −2 in S1 and −0.1 Wm −2 in S2, while RS-PML systematically underestimates in both seasons (−19.9 and −13.3 Wm −2 ).   The seasonal pattern of λE modeled by the TIR-ET models and compared with λE measured by the eddy covariance (λE EC) (Figure 7), show that both METRIC and TSEB-PT replicate very well the intra-seasonal variability of λE EC, however, METRIC underestimates λE in S 1 (−10.2 Wm −2 ) and overestimates in S 2 (16.2 Wm −2 ), while TSEB-PT underestimates both seasons (−17.6 and 1-2.7 Wm −2 ). Therefore, Figure 7g,h show that METRIC tends to underestimate λE in early season (November-December), while the estimation is closer to measured λE between January and March. On the other hand, Figure  7g-h show that the TSEB-PT model underestimates λE when H is overestimated and, vice versa, λE is overestimated when H is overestimated, suggesting a strong dependence of λE with the remaining components of the energy balance, particularly H.
Assessing the intra-seasonal variation of λE modeled by RS-PML and RS-PMS (Figure 8c-d), both models successfully replicate the seasonal variability of λE, especially RS-PMS with a bias of −10.8 Wm −2 in S 1 and −0.1 Wm −2 in S 2 , while RS-PML systematically underestimates in both seasons (−19.9 and −13.3 Wm −2 ).

Discussion
The results obtained at instantaneous time scale both METRIC and TSEB-PT are within the range of error reported in vineyards, with a RMSE ranging from 40 to 95 Wm −2 [21,90,91] in METRIC and from 47 to 89 Wm −2 in TSEB-PT model [28,30,92]. Additionally, both METRIC and TSEB-PT showed similar performance in terms of RMSE, %RMSE and d 1 . However, TSEB-PT had slightly better performance than METRIC (Table 5). This behavior may be due to both models sharing the same input data, as reported in an evaluation of both models in herbaceous (corn, soybean, wheat) and woody (olive and vineyards) crops [23].
The slight underestimation of TSEB-PT is possibly due to an overestimation of the short and long wave radiation components and hence a systematic underestimation of R N (Figure 7c,d), this implies less energy availability for the estimation of λE and H, especially if it is taken into account that λE is a residual function of the components of the energy balance [93,94]. Conversely, an underestimation of H implies greater energy availability for λE modeling and therefore its overestimation [93,95,96] as occurs in S 2 with λE modeling by METRIC (Figure 7h).
Furthermore, LST estimation from Landsat-TIR bands require atmospheric and emissivity correction, which is a source of uncertainty in modeled H and consequently modeled λE from METRIC and TSEB-PT. However, METRIC reduces this uncertainty by including an internal calibration process, which depends on the selection of a hot (dry surface) and a cold (wet surface) pixel [15], selection done automatically throughout an exhaustive search algorithm [77]. On the other hand, TSEB-PT has greater sensitivity to LST estimation, since it is associated to estimation of T C , T AC and T S (Equations (10)- (12)) and consequently to the propagation of the error in H partition between soil (H S ) and canopy (H C ) [97].
When comparing approaches to model the soil evaporative fraction required both RS-PML and RS-PMS models, the best fit of the f SWC approach is because the energy for soil evaporation depends on the soil water content of top soil layers [41], which is captured by measuring the soil water content. However, satisfactory results have been obtained with the f Zhang approach in sparsely vegetated areas characterized by a high evaporation rate and higher frequency of precipitation [44], which evidenced that the aforementioned approach reproduces well the dynamics of soil evaporation when it depends on the precipitation of previous days. Therefore, due to the characteristics of drip irrigation, where soil evaporation is located only in the portion of soil wetted by the dripper, the f SWC approach is able to simulate more accurately the soil evaporation dynamics under such conditions. The f SWC approach consistently had the best fit in comparison to the f Zhang approach, hence hereinafter the results will be shown only for the f SWC approach for both the RS-PML and RS-PMS models, omitting the subscript f SWC .
Additionally, the results from the f SWC approach are comparable to those obtained in the evaluation of two remote sensing Penman-Monteith models (RS-PM) using visible, near-infrared and thermal infrared data from Landsat 5 and Landsat 7 platforms, yielding a higher RMSE around 39.7 Wm −2 and establishing the strong dependency between ET and soil water content, especially in irrigated vineyards and orchards, since a larger fraction of net radiation (R N ) is allocated to instantaneous latent heat flux (λE) and not to instantaneous sensible heat flux (H), as would occur on dryland vegetation covers [39].
The behavior of λE modeled by RS-PML responds to phenology, to the error associated to LAI estimation (Figure 7a,b) from remote sensing data [41], to the inherent errors in λE measurement by the eddy covariance system and to the parameterization of the maximum stomatal conductance (g sx ). The first two aspects were already discussed previously, while g sx is a species-specific parameter that varies according to phenology [41] and in this work is taken the maximum value reported in the bibliography [87]. In spite of the latter, it is possible that the selected g sx is lower than actually evidenced in the assessed seasons through measured λE, however, the g sx value satisfactorily represents the last phenological stage corresponding to Ripening of Berries (February-April) as is observed in Figure 7c,d where the bias is lower compared to the overall season.
On the other hand, the model of stomatal conductance developed by Stewart [45] provides a better description of transpiration dynamics in the vineyard, which involves the adjustment of four parameters (g SX , f 1 , f 2 and f 3 ), while the Leuning stomatal conductance approach [42] requires the adjustment of one parameter (g SX ), becoming a model that can be easily applied in a wide range of environmental conditions. Additionally, the Stewart stomatal conductance model requires soil water content measurements, which are not commonly carried out in commercial vineyards, while the Leuning stomatal conductance model requires vapor pressure deficit measurements, which are available from nearby weather stations. However, SWC measurements in commercial fields should become more frequent since water restrictions for irrigation purposes are expected to increase in the coming years. Therefore, through the widespread adoption of SWC measurements, the monitoring of available water at the root zone would be improved and the application of the RS-PMS model would be extended to different environmental and agronomic management conditions.
Analyzing different timescales where instantaneous λE was extrapolated to daily ET, 7-day and monthly ET (Table 5), the best agreement index (d 1 ) corresponds to the RS-PMS model (0.7), followed by RS-PML (0.58), TSEB-PT (0.56) and METRIC (0.5), which is consistent with respective RMSE and %RMSE values. Therefore, daily ET values modeled by METRIC have a slightly higher RMSE than those reported in vineyards [22,91,98], ranging from 0.4 to 0.8 mm day −1 . These differences possibly are due to the extrapolation method chosen, since in addition to the approach used in this research (λE/Rs ratio), there are other common methods such as ET/ET 0 [22,91] and ET/R N -G ratio [98].
In addition, daily ET error reported in vineyards by other TSEB models such as ALEXI/DisALEXI yield RMSE around 0.35-0.79 mm day −1 and %RMSE between 8 and 22% [28,30,99], which are lower than those obtained by TSEB-PT model, possibly due to the fusion of remote sensing data methodology used in ALEXI/DisALEXI models, allows a comprehensive daily characterization of vineyard water consumption.
Moreover, it is reported lower daily ET RMSE of 0.34 mm day −1 with RS-PM model that include brightness temperature data from satellite sensor [39], which improve modeled daily ET. However, both RS-PMS and RS-PML models are consistent with reported by the VIS-ET models based on multiple linear regression with RMSE around 0.62 mm day −1 [100].
The irrigation management in the evaluated vineyard takes as a criterion the ET 0 multiplied by a crop coefficient (Kc) of 0.5 (ET 0 × 0.5) during the whole season (Kc Target) (Figure 10c,d). In this context, the proposed RS-PMS model provides a method to evaluate and monitor whether the irrigation applied corresponds to the criteria selected. Therefore, Figure 9 shows the water irrigation applied in seasons S 1 and S 2 , which is around to the estimated by the criterion ET 0 × 0.5, despite irregular management in S 1 (Figure 10a) in comparison with S 2 (Figure 10b). Nonetheless, Figure 10c,d show that estimating the Kc provided by the eddy covariance system (Kc EC), which is the ratio between ET EC and ET 0 , has standard deviation of 0.04 and does not reach the Kc Target of 0.5, thus implying that the irrigation applied is lower than the estimated by the vineyard irrigation management. Hence, the Kc estimated by the RS-PMS model (Kc PMS), this means the ratio between ET RS-PMS and ET 0 , is reasonably similar to the Kc EC with a bias around 9% and standard deviation of 0.06, which would allow accurately the evaluation and monitoring of irrigation [101][102][103][104] at a spatial scale of 10 m × 10 m and a daily and weekly temporal scale, which are suitable for irrigation management. crops such as vineyards), impacting H estimation and generating uncertainties in λE, meanwhile VIS-ET models are capable to capture the heterogeneity of land surface, because of higher pixel resolution [13]. However, the data acquired by aircraft and Unmanned Airborne Vehicles (UAV) provide a high spatial resolution (< 10 m) of TIR bands, reporting better fit of the actual evapotranspiration modeled by METRIC and TSEB-PT [20,29,31,105,106], which clearly suggests that the difference in performance of the TIR-ET and VIS-ET models is given by the resolution of the input data. Additionally, under the environmental conditions and vineyard management addressed in this study (Mediterranean climate and drip irrigation), RS-PMS outperforms RS-PML, METRIC and TSEB-PT. However, several researches have highlighted that TIR-ET methods are more appropriate than VIS-ET methods under conditions of water shortage, because LST is more sensitive to vegetation water condition, while VIS-ET methods do not have a fast response to early plant water stress [15,38]. In contrast, there are other studies that have concluded that VIS-ET models, with the current specifications of the available remote sensing platforms, have better performance than TIR-ET methods due to uncertainties in LST estimation for low pixel resolution (and widely separated row crops such as vineyards), impacting H estimation and generating uncertainties in λE, meanwhile VIS-ET models are capable to capture the heterogeneity of land surface, because of higher pixel resolution [13]. However, the data acquired by aircraft and Unmanned Airborne Vehicles (UAV) provide a high spatial resolution (<10 m) of TIR bands, reporting better fit of the actual evapotranspiration modeled by METRIC and TSEB-PT [20,29,31,105,106], which clearly suggests that the difference in performance of the TIR-ET and VIS-ET models is given by the resolution of the input data.

Conclusions
The Penman-Monteith model was evaluated for the estimation of actual evapotranspiration (ET) with remote sensing data (RS-PM) at multiple temporal scales, on a vineyard located in central Chile and under drip irrigation. Therefore, two approaches of the RS-PM model were evaluated, one based on the Leuning approach (RS-PML) [41] and a novel approach based on remote sensing data that couples the RS-PM model to Stewart's stomatal conductance model [89] (RS-PMS). Both models employed Sentinel-2 remote sensing data and were compared with independent eddy covariance data and both the Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC) [18] and the Priestley-Taylor Two-Source Energy Balance (TSEB-PT) [78], which are widely used to assess actual ET and are based on infrared thermal band data from the Landsat 7 and Landsat 8 satellites, with lower spatial and temporal resolution than the Sentinel-2 visible and near-infrared band data.
The RS-PM models outperform the METRIC and TSEB-PT models, especially in the instantaneous and daily time scales where RS-PM models have the lowest root mean square error (RMSE) and the highest Wilmott agreement index (d 1 ). The higher temporal scale resolution of the input data in RS-PM models might be foremost responsible for those results since with higher resolution data obtained from aircraft and Unmanned Airborne Vehicles (UAV), both METRIC and TSEB-PT models performed similarly and even better than shown in this research for RS-PM models. On the other hand, the RS-PMS approach yields a better performance than RS-PML in all temporal scales evaluated, due to the four fitted parameters which improve the model ability to properly capture variations in soil water content, which in conjunction with the leaf area index determines the evapotranspiration rate of the vineyard managed with localized irrigation. However, the implementation of this approach under conditions different from those of this research will be subject to the availability of data to adjust the four model parameters, especially measurements of soil water content.
Finally, it was demonstrated that RS-PMS model is suitable to assess water consumption and irrigation management, since it is possible to estimate the actual Kc of the vineyard with the actual modeled ET (Kc PMS), which is one of the criteria used in irrigation management. Thus, on average the Kc PMS is biased by 9% compared to the Kc estimated with data from the Eddy covariance system (Kc EC). This is relevant taking into account it was found that during the two evaluated seasons, the vineyard's Kc was bellow of the 0.5 target Kc, equivalent to 0.33 in S 1 (2017-2018) and 0.34 in S 2 (2018-2019), representing a 32% lower than planned.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to because it is currently being used in various research projects.