Calibration and validation of an aerodynamic method to estimate the spatial variability of sensible and latent heat fluxes over a drip-irrigated Merlot vineyard

ABSTRACT A study was carried out to calibrate and validate the aerodynamic temperature method for estimating the spatial variability of the sensible (H) and latent (LE) heat fluxes over a drip-irrigated merlot vineyard located in the Maule Region, in Chile. For this study, measurement of energy balance components and meteorological data were collected from the 2006 to 2010 growing seasons. The experimental plot was composed of a 4.25 ha of ‘Merlot’ vineyard, which was equipped with an Eddy-Covariance system and an automatic weather station. The k-fold cross-validation method was utilized to tune and validate a vineyard surface aerodynamic temperature (Taero) model, considering all of the days when Landsat scenes and ground measurements of meteorological data and surface energy balance (SEB) were available. Then, the satellite-based estimations of Taero were utilized to calculate the surface aerodynamic resistance (rah) and, subsequently, heat fluxes of H and LE. Results indicated that the estimated H and rah values were not significantly different to those measured in the vineyard (95% significance level) showing a root mean square (RMSE) and mean absolute error (MAE) between 34–29 W m−2 and 1.01–0.78 s m−1, respectively. Satellite-based computations of LE were somewhat higher than those measured at the time of satellite overpass (RMSE = 63 W m−2; MAE = 56 W m−2), presumably due to the biases embedded in the net radiation (Rn) and soil heat flux (G) computations. The proposed SEB method based on Taero is very simple to implement, presenting similar accuracies on ET mapping to those computed by complex satellite-based models.


Introduction
In semi-arid zones, sustainable management of water resources for agricultural practices such as irrigation scheduling requires the accurate determination of crop water consumption. This is because in any agricultural system the evaporation of water from the ground surface to the atmosphere is driven by the link between the soil water budget and surface energy balance (SEB), (Brutsaert 1982). Equations (1) and (2) describe the widely adopted simplified SEB: R n À G À H À LE ¼ 0; (1) where R n , LE, G, and H are the net radiation, latent heat flux, soil heat flux, and sensible heat flux, respectively (all in W m −2 ); λ is the latent heat of vaporization (J kg −1 ); ET a is the actual evapotranspiration rate (mm h −1 ). Net radiation and soil heat fluxes of the vineyard can be routinely measured using in situ net radiometers and soil heat flux plates, respectively. These radiative components of Equations (1) and (2) can also be estimated in vineyards by using simple relations, which are based on the Stefan-Boltzman Law or by empirical relations (Carrasco and Ortega-Farias 2008). Values of LE and H over vineyards are regularly measured using micrometeorological systems such as the eddy-covariance (EC), Bowen Ratio (BR) Systems, and Scintillometry (SC), among others (Hicks 1973;Oliver and Sene 1992;Heilman et al. 1994;De Bruin, Van Den Hurk, and Kohsiek 1995;Heilman et al. 1996;Green et al. 2000;Kordova-Biezuner, Mahrer, and Schwartz 2000;Ortega-Farias, Poblete-Echeverría, and Brisson 2010;Poblete-Echeverría and Ortega-Farias 2009;Balbontín-Nesvara et al. 2011;Poblete-Echeverría, Sepúlveda-Reyes, and Ortega-Farias 2014). The direct measurement of the H and LE fluxes in vineyards has been indicated as complex because they are significantly affected by the sparse characteristic of the rows that directly influence sunlight penetration and thus the microclimatic conditions of the canopy elements and the surrounding soil surface (Heilman et al. 1994;Kool et al. 2014;Palladino et al. 2013). In vineyards trained on a vertical shoot positioned system (VSP), the rows represent a heterogeneous rough surface that influences the turbulent fluxes. In field studies of drip-irrigated vineyards, Ortega-Farias, Poblete-Echeverría, and Brisson (2010), indicated that the canopy size directly influences the fractional cover (f c ) where the soil surface is the major contributor of H to the vineyard SEB. Thus, changes in H and LE fluxes are highly dependent on plant size, row orientation, row width and the trellis system.
The most important limitation of micrometeorological system measurements is that they only represent a small footprint area and they do not consider the spatial variability of LE and H (Gowda et al. 2007a;Kleissl, Hong, and Hendrickx 2009), reducing its applicability to large areas.
Since the 1990s, satellite images have been considered as a promissory technology for studying large areas and over a wide range of vegetation types and water availability (Gowda et al. 2007a;Allen et al. 2011). By inverting Equation (1), several satellite-based residual surface energy balance (RSEB) models have been developed to estimate the spatial variability of LE, H, and ET a with errors between 10% and 20% (Gowda et al. 2007a; Kalma, McVicar, and McCabe 2008;Cammalleri et al. 2012;González-Dugo et al. 2012;Petropoulos 2013). In this context, Bastiaanssen et al. (2008) indicated that the SEB Algorithm for Land (SEBAL) model estimated ET a in several rainfed and irrigated vineyards with a maximum error of 24%. Galleguillos et al. (2011) tested the Simplified SEB Index (S-SEBI) model (Roerink, Su, and Menenti 2000) with errors between 27% and 44% in the estimations of ET a in vineyards. González-Dugo et al. (2012) observed that the Two-Source Energy Balance (TSEB) model (Kustas and Norman 1997) underestimated LE by about 18% in a vineyard under low vegetation cover and in semiarid conditions. Carrasco-Benavides et al. (2014) overestimated H and LE over a drip-irrigated vineyard with errors of 16% and 17%, respectively, by using the Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC) model (Allen et al. 2007).
Equation (1) represents a reliable alternative to estimate the satellite-based values of LE, because the radiative components of the SEB (R n and G) can be routinely estimated from the visible (VIS), near-infrared (NIR), and thermal infrared (TIR) wavebands data contained in each pixel of the remotely sensed image. Uncertainties in these estimations have been found to be between 5-10% and 20-30% for R n and G, respectively (Petropoulos 2013). The most important challenge of the satellite-based RSEB models is to accurately estimate LE because it lies in the correct calculation of H (Norman, Kustas, and Humes 1995;Barbagallo, Consoli, and Russo 2009;Petropoulos 2013). Nevertheless, the satellite-based computations of R n and G for vineyards also represent a challenge because they are very sensitive to the row cropped structure (Carrasco-Benavides et al. 2014;Kool et al. 2016). Operational models such as SEBAL or METRIC (and others) require the selection of the two 'anchor' pixels (cold and hot) that make the estimation subjective and dependent on the ability of the operator in searching and isolating the most appropriate hot and cold pixels (Gowda et al. 2007b;Singh 2012, 2013;Petropoulos 2013), and/or on the availability of such extreme pixels in the satellite imagery. As an alternative, the aerodynamic model has been suggested to estimate H as a function of the bulk aerodynamic resistance (r ah , s m −1 ), surface aerodynamic temperature (T aero , K) and air temperature (T a ) (Chávez et al. 2005;Kalma, McVicar, and McCabe 2008;Chávez et al. 2010): where ρ a is the air density (kg m −3 ) of moist air; C p is the specific heat of air (1005 J kg −1 K −1 ) of dry air; T a is the air temperature (K). According to Kalma, McVicar, and McCabe (2008), Equation (3) is a reliable expression that presents several significant difficulties in its use because T aero and r ah are not easy to obtain. T aero determines the sensible heat fluxes from the surface, representing an extrapolation of T a profile down to an effective height within the canopy at which the vegetation elements of H and LE fluxes arise (Chehbouni et al. 1996). T aero and r ah can be obtained by inverting Equation (3) only when micrometeorological measures are available. From the foregoing, a simple solution to estimate H was proposed by Chávez et al. (2005), (2010)) who described an empirical method that computes H based on the estimation of T aero using a multiple linear regression function obtained from ancillary micrometeorological and remotely sensed data. This method estimates T aero by combining the surface radiometric temperature (T s ), T a and horizontal wind speed (u). In addition, the plant biomass characteristic is considered in the linear regression through the leaf area index (LAI, m 2 m −2 ). This methodology was able to estimate H over soybean and corn fields with a root mean square error (RMSE) of 30 W m −2 (Chávez et al. 2005).
Using a similar method, Chávez et al. (2010) indicated that H over a rainfed cotton was estimated with errors, ranging between −3.7% and 9.3%.
Until now, the methodology for estimating H (and then LE) based on T aero has not been applied on sparse vineyards (trained on VSP). Furthermore, considering that H is a key factor to estimate ET a under these conditions, the objective of this research was to calibrate and validate the aerodynamic method to estimate the spatial variability of H and LE over a commercial drip-irrigated Merlot vineyard.

Experimental setup
The trial was part of a 4-year experiment (from 2006 to 2010) that was carried out to collect weather and surface energy balance components on a drip-irrigated Merlot vineyard experimental plot of 4.25 ha. The whole vineyard (75 ha) is located in the Maule Region, Chile (35°25' S; 71°32' W; 125 m. a. s. l.). The experimental plot was surrounded by Merlot, Carménère, and Cabernet Sauvignon vineyards of similar characteristics ( Figure 1). The climate in the area is classified as Mediterranean semiarid with an average daily temperature of 17°C. The annual rainfall is close to 700 mm concentrated during the fall to winter period (April-August). During the spring to summer period (September-March), the rainfall is almost negligible and it represents only 2% of the yearly cumulative precipitation value. The annual cumulative short grass reference evapotranspiration, computed by the Penman-Monteith FAO 56 method (Allen et al. 1998) was approximately 1100 mm.
Merlot vines trained on a VSP system were planted in 1999 in north-south rows and placed at 1.5 m between rows and 2.5 m between vines (2667 plants ha −1 ). The vineyard canes and leaves were kept as parallelepiped where the maximum canopy height (h c ) from the soil surface reached approximately 2 m after the ripening of the cane phase. Also, the LAI and f c ranged between 0.8-1.2 m 2 m −2 and 28-31%, respectively. This vineyard was drip-irrigated using emitters with a flow capacity of 4.0 L h −1 located at the side of each vine. The soil surface in the inter-row space was free of weeds since by the end of November the soil surface was very dry. The remaining weeds were eliminated by mechanical practices to maintain the vineyard clear. To ensure that the vineyard was well irrigated, during the complete growing season, the volumetric soil water content (θ) at 0.6 m depth was weekly monitored before the irrigation event by using a portable TDR probe (6050X1, Trase System, Santa Barbara CA, USA). At the same time, the vineyard water status was monitored by measuring the midday stem water potential (Ψ s ) using a pressure chamber (model 600, PMS instruments, Albany OR, USA).

Micrometeorological measurements
Vineyard energy balance components and micrometeorological data were registered at half hour intervals using an automatic meteorological station (AMS) mounted on a tower located inside the experimental plot. In this tower, wind speed (u) and wind direction (w dir ) were monitored using a cup anemometer and a wind vane (03101-5, Young, MI, USA), respectively. Air temperature (T a ) and relative humidity (RH) were measured by a combined probe (HMP45C, Vaisala, MA, USA). The net radiation (R n ), incoming (R si ), and outgoing solar radiation (R so ) were measured by a four-component net radiometer (CNR1, Kipp & Zonen Inc., Delft, The Netherlands). All of the aforementioned sensors were placed at 4.7 m above the soil surface and data values were recorded using a data logger (CR 5000, Campbell Scientific Inc., Logan, UT, USA). A complete eddy-covariance (EC) system was installed in the same tower to provide continuous measurements of LE and H at 30 min time intervals. The EC system was composed of a fast response open path infrared gas analyser (IRGA) (LI7500, LI-COR Inc., Lincoln, NE, USA) and a threedimensional sonic anemometer (CSAT3, Campbell Scientific, Inc., Logan UT, USA) to compute LE and H, respectively. The raw data, collected by the EC system, were sampled at a high frequency of 10 Hz and stored in a PCMCIA flash memory card. Data collected were post processed off-line by using the EdiRe software package (Campbell Scientific, Inc 2008). The post processing included: the coordinate rotation (Wilczak, Oncley, and Stage 2001), sonic temperature (Schotanus, Nieuwstadt, and De Bruin 1983), and density (Webb, Pearman, and Leuning 1980) corrections. Average soil heat fluxes were recorded at 30 min intervals on an electronic data logger (CR1000, Campbell Scientific, Inc., Logan, UT). Data were obtained from eight soil heat flux plates (HFT3, Campbell Scientific Inc., Logan, UT, USA), buried at 0.08 m deep, and eight pairs of soil thermocouples (TCAV, Campbell Scientific Inc, Logan, UT, USA), buried above each flux plate at 0.02 and 0.06 m, which were regularly distributed between the inter-row and beneath the vine space. Regarding quality control, the energy balance closure was applied using the ratio of turbulent fluxes (H + LE) to the available energy (R n -G) (Ortega-Farias, Poblete-Echeverría, and Brisson 2010). Finally, values of LE and H were re-calculated using the Bowen ratio (ß = H/LE) method (Twine et al. 2000;Balbontín-Nesvara et al. 2011;Monteith and Unsworth 2013).
To explore the consistency and relationship among micrometeorological variables at the time of satellite overpass (R si , u, RH, and T a ) for different growing seasons all databases were analysed using the principal components analysis (PCA) method (Jolliffe 1990).

Landsat images: thermal waveband processing and corrections
21 Landsat scenes with less than 30% of cloud cover were available for this study (Table 1). They were downloaded from the USGS Glovis website (http://glovis.usgs. gov). The vineyard area was covered by the Landsat path/row 233/85 as a level-one terrain-corrected (L1T) product. Since Landsat 7 (ETM+) scenes contain gaps due to the scan line corrector failure (slc-off), all gaps identified inside the scene were isolated and not considered in the image processing. Each satellite scene was processed by converting the digital numbers (DN), for each pixel, to physical values following the steps indicated in the literature (Tasumi, Allen, and Trezza 2008;Allen et al. 2002;NASA 2010). Additional details on imagery processing adopted can be found in Carrasco-Benavides et al. (2012), Carrasco-Benavides et al. (2014)). All images were re-projected to the projection type WGS 84; Universal Transverse Mercator (UTM) coordinates system, zone19 South (S). All images were processed using the ERDAS Imagine processing system (Hexagon Geospatial). Land surface temperatures (LST or T s ) were computed from Landsat thermal band (bands 6 and 6.1 for Landsat 5 (TM) and ETM+, respectively). The method used to compute LST is based on the radiative transfer equation for the TIR band as follows: where L sen is the radiance measured by the Landsat sensor (W m −2 sr −1 µm −1 ); ε s is the surface thermal emissivity estimated as reported in Allen et al. (2002) ; B(T s ) is the Planck's radiance at surface temperature T s (W m −2 sr −1 µm −1 ); T s (K); L # and L " are the downwelling and upwelling atmospheric radiances, respectively (W m −2 sr −1 µm −1 ) and τ is the atmospheric transmissivity. Magnitudes indicated herein are spectral and angular dependent but they (and their symbols) have been omitted to simplify Equation (4). This has been done because these terms have been averaged in the Landsat thermal band. The computation of LST has been performed by a single-channel model where the values of B(T s ) were obtained by inverting Equation (4). Also, the radiation-surface interaction has been assumed as Lambertian (Qin, Karnieli, and Berliner 2001;Sobrino, Jiménez-Muñoz, and Paolini 2004;Sobrino et al. 2009), yielding: Finally, the LST (or T s ) was estimated following modified Planck equation: where λ is the centre wavelength of the thermal band (μm); c 1 and c 2 are Planck's function constants (c 1 = 1.19104 × 10 8 W μm 4 m −2 sr −1 and c 2 = 14387.7 μm K). The calculation of τ, L # and L " was carried out employing data of atmospheric profiles for the closest date and time of the satellite overpass. They were calculated using the Moderate Resolution Atmospheric Transmission (MODTRAN) model (Berk, Bernstein, and Robertson 1989). Considering that these atmospheric profiles are not always available for the date and time required, hindering the application of this method (Morales and  (Kalnay et al. 1996;Kistler et al. 2001;Trenberth, Fasullo, and Mackaro 2011) were also consulted. Data from the NNR database were processed as indicated in Barsi, Barker, and Schott (2003) and Barsi et al. (2005) exporting the values of the atmospheric profiles in MODTRAN format for the post processing.

Derivation of surface aerodynamic temperature and aerodynamic resistance
The surface aerodynamic temperature can be defined as the temperature at the height of the zero plane displacement (d, m) plus the momentum transfer height (Z om , m).
Since T aero and the aerodynamic resistance (r ah ) are normally unknown and not directly measured by the EC system, they can be calculated by inverting the so-called bulk aerodynamic resistance equation (Equations (3), (7), and (8)) using measurements of H, T a , and u. The aforementioned variables, along with data of crop height and wind speed height, allow us to explain the atmospheric stability through the Monin-Obukhov (MO) similarity theory (Chávez et al. 2005). The aerodynamic resistance and temperature can also be calculated iteratively as part of the surface energy balance (Dhungel et al. 2014a;Dhungel, Allen, and Trezza 2016). Briefly, the iterative method to calculate 'r ah ' is based on calculating an initial 'r ah ' for neutral atmospheric stability and then updating it using the MO method and calculating a new H. In each calculation cycle the previous step r ah (r ah-1 ) is compared with the current step r ah to evaluate if the difference in values is within 5%. If this is not the case, then a new iteration is needed until r ah values converge. The critical aspect of the convergence and iteration method occurs when the horizontal wind speed is too low (u < 1 m s −1 ). In this case, a solution for Equation (3) cannot be found. When numerical stability is not achieved, the wind speed variable values can be increased in order to reach the needed numerical convergence as indicated in Dhungel, Allen, and Trezza (2016). This adjustment approach has been successfully adopted in a recent study by Vashisht (2016) where the wind speed was fixed at 1 m s −1 when numerical stability, in the calculation of r ah , had not been achieved. In this study, a sensitivity analysis indicated that the lowest variability on an hourly evapotranspiration of grass pasture (ET o ) was around 4% when the wind speed was between 1.0 and 1.8 m s −1 .
The key issue of this approach is to accurately 'measure' H by the EC system; which allows to obtain the inverted values of T aero (T aero_inv ) and r ah (r ah_inv ): where k is the Von Karman constant (0.41); Z m , Z om , and Z oh are the horizontal wind speed measurement height, roughness length for momentum transfer, and sensible heat transfer, respectively (m). As stated above, d is the zero-plane displacement height (m); ψ m andψ h are the atmospheric stability correction factors for heat and momentum transfer, respectively. The variable L is the Monin-Obukhov atmospheric stability length or height (m), and u Ã is friction velocity (m s −1 ) at Z m , which is a function of horizontal wind speed (u, m s −1 ) measured by the 3D sonic anemometer at a height of 4.7 m above ground. The values of Z oh (= 0.1 Z om ) were computed according to the Brutsaert's method (Chávez et al. 2010) while those of d and Z om were estimated as follows (Poblete-Echeverría and Ortega-Farias 2012): Z om ¼ 0:13h c : 2.5. Calibration and validation of the T aero model in order to compute H and LE at the time of satellite overpass The proposed method to model T aero at the time of satellite overpass, consist in a multiple linear regression equation that was calibrated as follows (Chávez et al. 2005): where T aero_m is the modelled surface aerodynamic temperature ; β 0 , β 1 , β 2 , β 3 , and β 4 are empirically fitted parameters specific for the vineyard and environment conditions found in this study. Values of u, T a , and T s were measured, at the time of satellite overpass (near 1124 h at local time). Values of LAI were obtained using the Weibull function suggested by Carrasco-Benavides et al. (2014).
Due to the small size of the satellite imagery database (Table 1), the calibration and validation of the T aero_m model were executed by means of the k-fold cross-validation method (Refaeilzadeh, Tang, and Liu 2009). This process was carried-out using a script developed in R (https://www.r-project.org/). The model calibration includes the coefficient of determination (R 2 ) and mean absolute error (MAE) (Mayer and Butler 1993). For the model validation, the comparisons between the observed (T aero_inv , r ah_inv , H, and LE) and modelled (T aero_m , r ah_m , H _m , and LE _m ) values were carried-out considering the averaged fluxes of pixels that fell inside an area of 36 × 36 pixels (inside the footprint of the EC system) as indicated in Carrasco-Benavides et al. (2012); Carrasco-Benavides et al. (2014)). Estimated values of LE were calculated as residual of Equation (1) by combining the measured values of R n and G and the estimations of H (H _m ). In addition, pixel-bypixel values of R n and G were estimated to evaluate the effect of additive biases from estimations of R n and G into the computations of H and LE. In this regard, the satellitebased net radiation was calculated by adding the incoming shortwave and the outgoing longwave radiations for each pixel according to the model proposed by Tasumi, Allen, and Trezza (2008). In this model, the incoming shortwave radiation is calculated based on the surface albedo (α), which was obtained separately for each band using the atsurface reflectance TM/ETM+ sensor band, whose values were then integrated according to the intensities of the incoming solar radiation within the domain of the band. On the other hand, the longwave components were obtained from the surface emissivities calculated at each pixel according to the empirical equations developed by Tasumi (2003). Estimations of G values were calculated using Bastiaanssen (1995) approach modified by Carrasco-Benavides et al. (2014).
The statistical analysis considered the ratio (b) of the estimated to observed values (b), mean absolute error (MAE), mean absolute per cent error (MAE%) and root mean square error (RMSE) (Mayer and Butler 1993). In addition, the Student's t-test analysis was applied to check whether b was significantly different from unity at the 95% confidence level.
Finally, an example of maps for the mean values of H _m and LE _m computed at each pixel for the entire vineyard were generated by the SAGA software (Conrad et al. 2015). For the fruit set, veraison, and harvest stages, these maps were developed to evaluate the spatial variability of the H _m and LE _m inside the vineyard.

Results and discussions
Micrometeorological conditions in the vineyard were stable during the study period (Figures 2(a-d)). At the time of satellite overpass, the daily mean T a was 19.7 (±4.0)°C (Figure 2(a)) with minimum and maximum values observed during October and January, respectively. A similar trend was observed for the RH (Figure 1 2006-2007, 2007-2008, 2008-2009, and 2009-2010). Dotted boxes indicate the time of satellite overpass. of Equation (3). For low values of u close to zero, Dhungel, Allen, and Trezza (2016) suggest that the wind speed may need to be set above a specific low limit (of 1.2 or 1.3 m s −1 ) to achieve numerical stability, which could increase the evapotranspiration rates by mechanical mixing. This alternative has been suggested by the aforementioned authors to prevent the under or overestimation of fluxes if the surface energy balance iteration is ceased without convergence. In this regard, the effect of u in the convergence of Equation (3) was tested for the vineyard by varying instantaneous values of u (at 25% each time) in a range between −50% and +50%. The results demonstrated that the variations (±50%) did not substantially affect the estimation of T aero_inv and r ah_inv (data not shown). Table 2 presents the comparison among T a , T aero_inv , and T s . Air temperature at the time of satellite overpass for all selected days was between 13.52 and 24.05°C, showing differences on average of −7 (±2)°C and −15 (±3)°C versus T aero_inv and T s , respectively. Similarly, T s presented the highest temperature values (25.9-40.4°C), showing average differences of about 14 (±3)°C in relation to T a , and 8 (±2)°C in relation to T aero_inv . Moreover, aerodynamic temperature presented values between 18.7°C and 30.4°C. These results were typical and have been widely discussed by Chehbouni et al. (1996)  and Kalma, McVicar, and McCabe (2008) who indicated that differences between T a and T s as much as 10°C are expected for sparse canopies. The relationship among micrometeorological variables at the time of satellite overpass is displayed in the bi-plot graph obtained from the PCA (Figure 4). This graph exhibits the days where different weather conditions were present during the different studied years. The main composite axes (F1 and F2) account for 74.1% of the total variance allowing explaining the data behaviour in a two-dimensional graph representation. In this regard, the days and growing seasons that were projected to the left side (F1 axis) were influenced by R si and T a . For example, the T a arrow showed that the higher temperatures of the study were recorded during the period 2007-2008, especially on DOY 17. Days and growing seasons that were distributed to the right side in the same axis were influenced by u and RH (i.e. DOY 51 and 363 in the 2008-2009 and 2006-2007 seasons, respectively) presenting the highest instantaneous values of u. The second factor (F2 axis) on the upper side shows that R si , u, and RH were similar to 47% of the studied days. The rest of days and growing seasons were similar for the T a variable (lower side of F2 axis). Based on this figure, it is possible to identify that DOY 297 and DOY 337 had different weather behaviour in relation to the other studied days.
The energy closure of the EC system showed that the regression analysis through the origin resulted in a coefficient of determination and slope of 0.95 and 0.91, respectively ( Figure 5). This result indicated that the turbulent heat fluxes (H + LE) were 9% lower than the available energy (R n -G) indicating an adequate energy balance closure. Similar Table 2. Comparisons of values of air (T a ), inverted aerodynamic (T aero_inv ) and surface (T s ) temperatures at the time of satellite overpass, for all days of the year (DOY) evaluated.

Difference (°C)
Year-DOY T a (°C)   -2007, 2007-2008, 2008-2009, and 2009-2010, respectively; 'D_number' denotes the day of the year (DOY). T a is the air temperature, RH is the relative humidity, u is the wind speed and R si is the solar radiation. results were observed by Ortega-Farias, Poblete-Echeverría, and Brisson (2010) who indicated that R 2 and b were equal to 0.96 and 1.08, respectively, for an energy balance closure over a drip-irrigated Cabernet Sauvignon vineyard. In addition, in Cencibel and Cabernet Sauvignon vineyards, Balbontín-Nesvara et al. (2011) observed that the sum of hourly values of (H + LE) was lower than the available energy (R n -G) with a slope = 0.85 and R 2 = 0.94 showing a deficit in the energy balance closure of around 5%. Recently, Er-Raki et al. (2013) for a for drip-irrigated table grape ('Perlette' and 'Superior') vineyards indicated that the slope of regression forced trough the origin was between 1.04 and 1.21 (<20% of deviation). Due to the predominant lack of energy balance closure presented in our study, H and LE were corrected by forcing closure using the Bowen ratio approach, as indicated in the Materials and Method section.

Modelling and validation of the aerodynamic temperature model
The multiple linear regression model used to estimate aerodynamic temperature as a function of T s , T a , u, and LAI (from Equation (11)) is shown below (R 2 = 0.96 and MAE = 0.57°C): T aero m ¼ À 22:77 þ 24:46LAI þ 0:75T a À 0:95uþ0:20T s : The validation of T aero_m , r ah_m , H _m , and LE _m is presented in Figure 6. In Figure 6(a), it can be seen that the estimation of T aero presented trends of measured to modelled values around the 1:1 line, with almost no spreading of points around this line suggesting a good model performance. The statistical analysis presented in Table 3 indicates that there was a good agreement between the observed and estimated values of surface aerodynamic temperature with R 2 , RMSE, MAE, and MAE% equal to 0.95, 0.66°C, 0.56°C, and 2.22%, respectively. In addition, the t-test suggested that b was equal to unity (one) indicating that observed and modelled values of aerodynamic temperature were similar at the 95% confidence level. These results were comparable to those presented by Chávez et al. (2005) who indicated that a linear regression model predicted T aero for rainfed corn and soybeans (humid continental zone, Central Iowa, USA) with an RMSE = 0.90°C. Also, Chávez et al. (2010) in cotton observed that the performance of the T aero model with and without LAI was similar with RMSE ranging between 1.02°C and 2.95°C. The results presented here for the drip-irrigated vineyard are in agreement with the two aforementioned studies, demonstrating that the development of empirical linear regressions based on ancillary data is a robust procedure to estimate T aero using remote-sensing and ground-based data.
For the aerodynamic resistance, Figure 6(b) depicts that points were close (tight) to the 1:1 line. The statistical analysis (Table 3) shows values of R 2 , RMSE and MAE, and MAE % equal to 0.95, 1.01 and 0.78 s m −1 , and 2.36%, respectively. The t-test indicated that the slope was significantly equal to one suggesting that modelled and estimated values were similar. For a comparison between r ah_m and inverted values of r ah (in Equation (3)), Chávez et al. (2005) indicated that MAE and RMSE were equal to −0.16 and 2.8 s m −1 , respectively (R 2 = 0.84). It is important to highlight that r ah describes the resistance from the upward vegetation and involves friction from air flowing over vegetative surfaces (Allen et al. 1998). As a consequence, the levels of error obtained in the computations of r ah_m of this study were similar to the other aforementioned for full cover crops (Chávez et al. 2005) presumably due to the low f c of the studied vineyard (28-31%) which poorly influenced the effect of the Landsat pixel information. For H _m , Figure 6(c) shows that the cloud of points tends to distribute around the 1:1 line. The statistical analysis presented in Table 3 indicates that the estimated and measured values of H were significantly equal with RMSE = 33.62 W m −2 , MAE = 28.34 W m −2 , and MAE % = 11.37%. These values are similar to those observed by Chávez et al. (2005) who obtained RMSE values ranging from 28.3 to 33.5 W m −2 in corn and soybean. Also, Chávez et al. (2005) observed that errors about 10% in r ah would result in errors between 9% and 10% in the estimation of H. In our study, an MAE% of 2.36% in r ah_m was Table 3. Statistical analysis of observed versus modelled values of aerodynamic temperature (T aero ), aerodynamic resistance (r ah ), sensible (H), and latent (LE) heat fluxes, at the time of satellite overpass (years 2006-2010 (1) and (2) indicate that fluxes of LE were computed as residual of the surface energy balance equation utilizing measured and estimated values of net radiation and soil heat fluxes, respectively.
associated with a MAE% of 11.37% in the estimation of H over the drip-irrigated vineyard. For the same vineyard plot, Carrasco-Benavides et al. (2014) observed that the METRIC model overestimated H with errors, RMSE and MAE ranging between 13% and 16%, 51-67, and 43-57 W m −2 , respectively. For several drip-irrigated vineyards in Albacete-Spain, González-Piqueras et al. (2015) indicated that the METRIC model overestimated the sensible heat flux by about 21% with a RMSE = 55 W m −2 . Therefore, results of this study suggested that the parameterization of the aerodynamic temperature model seems to be an alternative to other remote-sensing-based algorithms that depend on selection of extreme pixels. Furthermore, it is important to indicate that the estimation of H (Equation (3)) depends on an iterative procedure of convergence based on the Monin-Obukhov (MO) atmospheric stability correction. This method has been largely analysed by Dhungel et al. (2014aDhungel et al. ( , 2014b who developed a simulation of the surface temperature (T s ) based on the Backward Averaged Accelerated Numerical Solution (BAANS) model. These authors indicated that one of the uncertain parameters that have to be specified for calculating both the surface energy balance fluxes iteratively and the Monin-Obhukov stability correction is, for instance, the surface roughness for momentum transfer (Z om ). This fact, in turn, affects the simulation of the surface temperature (T s ) by the BAANS model, according to Dhungel et al. (2014a). Thus, in Dhungel et al. (2014a) when Z om was modified (lower values) from its initial values (Z om = 0.3 m) based on LAI for desert and grassland surfaces, the BAANS simulated T s values compared better with the Landsat thermal-based T s values. However, in the case of the parameterization of T aero for the Merlot vineyard in this study, the effects of Z om in the aerodynamic resistance (r ah ) and the energy balance accounted for when using derived T aero from EC system and thus the coefficients in Equation (11) adjust or compensate for changing Z om values (and effects on T aero ). The EC system measured actual friction velocity u Ã ð Þ, which is used in the computation of actual r ah . Then, the ECbased (measured) sensible heat flux (H) is used in the inverted form of Equation (3) to solve for 'measured or inverted' T aero . Therefore, the effects of Z om , for situations when T s or T aero is much larger than the air temperature (T a ) at screen height, have been accounted for in this study. Figure 6(d) indicated that values of R n and G utilized as an input in Equation (1) affected the estimations of LE. Ground measurements of R n and G allowed to compute LE close to the 1:1 line (grey dots) (RMSE = 33 and MAE 28 W m −2 , Table 3) in comparison with the same inputs obtained for each pixel utilizing the recommended satellite-based models from the literature (black dots) (RMSE = 63 and MAE 56 W m −2 ; see Table 3). Despite the fact that the estimations of R n and G had a relatively stable spatial distribution in the pixels inside the vineyard (coefficient of variation <10%, data not shown), both variables presented deviations between 9% and 17% (data not shown).
The low performance in the estimation of LE can be associated with the difficulties to accurately estimate R n and G based only on each pixel data for the vineyard. These are composed pixels that blend data from canopy and soil, where there is a low influence of the canopy elements (f c <30%). This difficulty in the hourly estimations of R n at the vineyard, as well as in other sparse crops (Petropoulos 2013) is not surprising because similar results have also been reported by Carrasco and Ortega-Farias (2008) in a dripirrigated Cabernet Sauvignon vineyard. They found that R n values were underestimated (b ≈ 0.9; RMSE = 39-60 W m −2 ) when an R n model based on meteorological data was applied. The lack of compensation errors of H _m observed in the results computed only from satellite inputs represents a suitable explanation of the mismatching of the LE _m values. The self-compensation of errors in Equation (1), from models such as SEBAL/ METRIC, has been indicated as controversial because the selection of pairs of end/ extreme pixels (cold and hot) necessary to calibrate each image is subjective Singh 2012, 2013;Petropoulos 2013). Our results provide evidence that the available energy of a vineyard is a variable that has to be carefully calculated.

Maps of modelled H _m and LE _m
The spatial variability of the average values of H _m and LE _m for the entire vineyard (97 ha) for the fruit set, veraison, and harvest phenological stages are depicted in Figure 7. To better understand the spatial organization of the vineyard, the panchromatic satellite image from Figure 1 can be consulted. As aforementioned, during all of the evaluated days the weather conditions as well as the canopy development were similar presenting clear skies ( Figure 2). In general, the maps of overall average values of H _m (Figure 7, top maps) presented a frequency histogram skewed to the right (data not shown) where the average (±standard deviation) and mode were 285 (±30) and 302 W m −2 , respectively. The lower values of H _m were observed at non-vineyard pixels such as grass and trees. The highest values of H _m (near 400 W m −2 ) were in pixels of covered areas with buildings and bare soil. The pixels from inside the experimental area (segmented square; area ≈ 4.25 ha) showed a range of H _m between 281 and 329 W m −2 . In this case, the frequency histogram of the spatial distribution of H _m was skewed to the left (data not shown) with an average of 298 (±11) W m −2 . Overall average values of LE _m (Figure 7, bottom) for all of the phenological stages were 229 (± 37) W m −2 presenting a median of 224 W m −2 . The frequency histogram for this variable was skewed to the left (data not shown). In general, the maps of LE _m were similar presenting the highest values in pixels that fell in the small forest placed at east boundary of the vineyard (close to 365 W m −2 ). For pixels contained in the vineyard, the highest values of LE _m were computed at the Veraison stage (280-320 W m −2 ). Inside the experimental site, the overall average and median of latent heat flux were similar (227 ± 28 W m −2 ) evidencing that this variable was stable inside the aforementioned area.
Differences observed between the average values of H _m and LE _m for each phenological stage were attributable to the decreases of R si which directly influence the available energy of the vineyard. For example, the spatial variability of H _m for Fruit Set seems to have a larger value than Veraison and Harvest (Figures 7(a-c)), respectively). This was associated to the hottest days observed during fruit set. For these days the values of R si and T a were up to 30 MJ m −2 day −1 and 30°C which influenced H _m . In vineyards at Fruit Set stage, the canopy elements are normally not completely developed affecting the soil energy balance beneath and between vines. Also, soil between rows was completely dry and without weeds which produced important amount of heat. Between Veraison (from January to February) and Harvest (March) trends in average values of H _m were similar and coincident with the decreasing of daily R si and T a . At harvest, the cooling of weather conditions observed in the vineyard affected the computations of T aero_m , which is linearly dependent of T s and T air the values of the H _m pixels. The contrary was observed for average values of LE _m (Figure 7, bottom) where the highest values over the vineyard were observed at Veraison (Figure 7(e)). At this stage (between the end of January and February), the canopy elements are completely developed, resulting in a higher leaf area than Fruit Set. This increased the vine transpiration and LE in response to the high water atmospheric demand observed at Verasion.
Values of H (100-400 W m −2 ) and LE (100-300 W m −2 ) presented in this study were in the range of those indicated in previous studies (Balbontín-Nesvara et al. 2011;De Bruin, Van Den Hurk, and Kohsiek 1995;Green et al. 2000;Kool et al. 2016;Kordova-Biezuner, Mahrer, and Schwartz 2000;Macedo Pezzopane and Pedro Júnior 2003;Oliver and Sene 1992;Spano et al. 2000). These results are evidence that the proposed model successfully estimated H and LE using remote-sensing inputs as well as weather data. This result is a promising solution to the problem of estimating H and then LE in the SEB of the vineyard, considering the spatial resolution of this satellite where each pixel represents a mixture of surface features from the vineyard row and inter-row.

Conclusions
This study presented a calibration and validation of an empirical model to estimate the sensible and latent heat fluxes of a Merlot vineyard, based on the bulk aerodynamic resistance equation. Fluxes of H _m and LE _m were calculated utilizing the estimation of T aero obtained from a multiple linear regression specifically calibrated for the vineyard in central Chile. The inputs of the T aero model were: leaf area index (also spatially modelled), air temperature and wind speed measured in the vineyard. Moreover, the model incorporated land surface temperatures which were computed from the Landsat thermal band. The validation results of the estimations of T aero were similar to those measured by the EC system installed inside the experimental plot presenting RMSE = 0.66 and MAE = 0.56°C. Estimated surface aerodynamic resistance and sensible heat fluxes were also equal to those measured in the vineyard, showing RMSE and MAE of 1.01 and 0.78 s m −1 ; and 34 and 29 W m −2 , respectively. The same was observed for LE _m when a combination of ground data and satellite-based inputs were utilized (MAE% = 16%).
The calibration and validation of this method ensured its applicability to compute the fluxes of H and LE in the vineyard. This is an important step to simplify the computation of fluxes of H over large surfaces. Furthermore, errors in estimated sensible heat fluxes were lower than errors found when other remote sensing of ET methods were used (in other studies) over Merlot vineyards in central Chile. Thus, besides being a simpler method, there is evidence that the proposed T aero algorithm estimates spatial values H and LE with less error and may be adopted for routine applications.
Care should be exercised if the calibrated model developed in this study is used in a different vineyard setting/configuration and with different satellite spatial resolution than those used in this study. As indicated herein, the complex interactions among climatic factors, the soil water availability, canopy roughness, canopy architecture (individually or all together) inside the vineyard, opens new needs into the exploration of a new adjustment to the proposed model, in order to offer simple alternatives in order to reflect the spatiotemporal heterogeneity of the SEB.