Comparison of two temperature differencing methods to estimate daily evapotranspiration over a Mediterranean vineyard watershed from ASTER data

Abstract Daily evapo-transpiration (ET) was mapped at the regional extent over a Mediterranean vineyard watershed, by using ASTER imagery along with two temperature differencing methods: the Simplified Surface Energy Balance Index (S-SEBI) and the Water Deficit Index (WDI). Validation of remotely sensed estimates was conducted during almost two growth cycles (August 2007–October 2008) over seven sites that differed in soil properties, water status and canopy structure. S-SEBI and WDI were also intercompared at the watershed extent by considering ASTER imagery collected between 2002 and 2008. In order to alleviate the experimental efforts devoted to the validation exercise, ground truthing relied on in situ estimates from the HYDRUS-1D model that simulates water transfers within the vadose zone after calibration against measured soil moisture profiles. For two of the seven validation sites, the consistency of the HYDRUS-1D simulations was beforehand controlled against direct measurements with eddy covariance devices. Thus, it was shown the HYDRUS-1D simulations could be used as ground truthing. Despite the use of simple differencing methods over a complex row-structured landscape, the obtained accuracies (0.8 mm.d − 1 for S-SEBI and 1.1 mm.d − 1 for WDI) were similar to those reported in the literature for simpler canopies, and fulfilled requirements for further applications in agronomy and hydrology. WDI performed worse than S-SEBI, in spite of more determinism within the derivation of evaporative extremes used for temperature differencing. This raised the question of compromising between process description and measurement availability. Analyzing validation results suggested that among the possible factors that could affect model performance (spatial variability, soil type and color, row orientation), the first-order influence was row orientation, a property that can be characterized from very high spatial resolution remote sensing data. Finally, intercomparing S-SEBI and WDI at the watershed extent showed estimates from both models agreed within 1 mm.d − 1 , a difference similar to the model accuracies as estimated by the validation exercise. Then, time averaged maps suggested the existence of spatial patterns at the watershed extent, which may be ascribed to combined effects from soil type, soil depth and watertable level.


Introduction
modeling errors may be large when considering row structured vineyards. An alterna-48 tive possibility is using differencing methods that rely on spatial contrasts, where ET 49 is characterized from differences in surface temperature rather than from absolute val-50 ues (Menenti and Choudhury, 1993;Moran et al., 1994;Gillies et al., 1997;Roerink 51 et al., 2000). This allows to minimize errors on both data and parameterization, and to 52 benefit from variabilities captured with high spatial resolution remote sensors (Moran 53 et al., 1996;Holifield et al., 2003;Gomez et al., 2005;Sobrino et al., 2005Sobrino et al., , 200854 Mendez- Barroso et al., 2008), where the Advanced Spaceborne Thermal Emission and 55 Reflection Radiometer (ASTER) sensor uniquely offers high quality spaceborne data 56 of surface temperature (Jacob et al., 2004;Sobrino et al., 2007b;French et al., 2008;57 Sabol et al., 2009). 58 The selection of differencing methods for the remote sensing of vineyard daily ET 59 is driven by compromising between feasibility and accuracy. On the one hand, simple 60 models are attractive since they are easy to implement and they usually perform as well 61 as complex models (Timmermans et al., 2007;Kalma et al., 2008). Thus, contextual 62 models allow to avoid delicate calculations of biophysical variables, since the evapo-63 rative extremes used for the temperature differencing are beforehand estimated from 64 the spatial variability captured within thermal imagery (Kalma et al., 2008). Good per-65 formances were reported by Galleguillos et al. (2011) when mapping vineyard daily 66 In order to infer daily ET from HYDRUS-1D simulations, each of the seven val-159 idation sites was monitored for soil moisture, watertable level and vegetation canopy 160 structure. The data collection included the following items.

161
• Soil moisture profiles were sampled every 0.2 m down to 2.5 m with a Vectra  (Table 1) 165 varied according to site size and to soil heterogeneity derived from a pedological 166 map. Profiles were collected around local solar noon, biweekly and after each 167 significant (i.e. large or lengthy) rainfall. 168 • Soil horizons were characterized for depth and texture, by using in situ obser-169 vations and expert knowledge. This was performed once at all sites, for each 170 location of soil moisture profile.

171
• Watertable level was monitored down to 2.5 m, either continuously at Site 6 by 172 using automated piezometric devices, or intermittently (i.e. with a two week fre-173 quency and after each significant rainfall, around local solar noon) at other sites 174 by using manual piezometric devices. In accordance with continuous observa-175 tions at Site 6, data collected on the other sites were linearly interpolated at the 176 hourly timescale.

177
• Vegetation was monitored for canopy structure (height above soil for base and 178 top of canopy, canopy width), over one (Site 1 to 5 and 7) or more (Site 6) fields. 179 Following Trambouze (1996), these observations were performed in winter (dor-180 mant vegetation with no leave) and in summer (period of maximum vegetation 181 cover), and they were next linearly interpolated at the daily timescale. On av-182 erage over the seven validation sites, canopy height was 1.5 m with a 0.1 m 183 standard deviation (7% in relative) and canopy width was 1 m with a 0.1 m stan-184 dard deviation (10% in relative). Row orientation and spacing were measured 185 once. 186 The ground based data above-presented were collected from August 2007 to Oc-187 tober 2008, apart from meteorological data that were collected on Site 6 since 2002. 188 Each instrument was manufacturer calibrated, apart from NP device that was calibrated 189 8 to account for soil type and soil moisture. Calibration was performed against gravimetric soil moisture data at each measurement depth (soil density was estimated using 191 a Campbell DR 501 gamma probe). The calibration residual error was 0.04 m 3 .m −3 192 (15% relative error). 193 2.2.2. Remote sensing data 194 We used ASTER official products for surface reflectance, waveband emissivity and 195 radiometric temperature (Abrams, 2000). Their accuracies were 5%, 0.01 and 1.5 K 196 respectively (Thome et al., 1998;Sobrino et al., 2007b;Jacob et al., 2008;Sabol et al., 197 2009). These products resulted from 20 scenes sun-synchronously collected under  For the seven validation sites, ground truthing of daily ET was obtained from simu-220 lations of HYDRUS-1D (Simunek et al., 2008) which is a physically-based model that 221 simulates water flow in the vadose zone from Richards equation (Richards, 1931): where θ is volumetric water content, t is time, z is vertical ascending coordinate, K 223 is unsaturated hydraulic conductivity, β is the angle between the water flow and the 224 vertical axis (β = 0 • for ascending flow fall was measured at the meteorological station ( § 2.2.1). Following Riou et al. (1994),

231
ET 0 was estimated from the data collected at the meteorological station, by using the 232 Penman formulation under standard conditions that can be found in Valiantzas (2006).

233
ET 0 was then split into reference transpiration T 0 and reference evaporation E 0 , by us-

234
ing the Riou's model designed for vineyards (Riou et al., 1994;Trambouze and Voltz, 235 2001; Bsaibes, 2007). First, T 0 was calculated as a fraction of ET 0 , by assuming this 236 fraction is almost equal to the ratio of solar irradiance absorbed by vine leaves R v (ge-237 ometrically derived from solar position and canopy structure) to that absorbed by the 238 whole vineyard (1 − a v )R g : where R g is solar irradiance. Vineyard albedo a v was set to 0.2 (Bsaibes, 2007). This 240 was done for the seven validation sites whose trellis structures were similar apart from 241 row orientation ( § 2.1). Second, E 0 was determined as the residual of the ET 0 splitting 242 into T 0 and E 0 : E 0 = ET 0 − T 0 . Next, vine maximum transpiration T m was set 243 to reference transpiration T 0 , which was within the confidence interval reported by 244 Trambouze and Voltz (2001) (Neuman et al., 1974). The threshold values were set to -10 m (respectively + 0.01 m) 280 for the minimum (respectively maximum) soil potential pressure that corresponds to 281 conditions of complete dryness (respectively runoff initiation).

282
For each location where soil moisture profiles were collected with NP measure-283 ments ( The data to be used as inputs for S-SEBI and WDI are listed in Table 3. The 290 obtaining of these inputs and the resulting implementation of both models are presented 291 in the current section.

292
[  (1978) and Jackson et al. (1983), respectively. Air humidity, wind 298 speed and atmospheric pressure were linearly interpolated, since they did not depict 299 any specific temporal behavior at the day or half-day timescales.

300
The Soil Adjusted Vegetation Index (SAVI) proposed by Huete (1988) was consid-301 ered for characterizing vegetation cover fraction from ASTER images, where SAVI is 302 defined as (ρ nir are ρ red are surface reflectances over near infrared and red wavebands): The choice of this index was motivated by its ability to account for soil effects, which 304 was appropriated for studying row structured vineyards. The coefficient L (soil adjusted 305 constant) was set to 1/2 by following Huete et al. (1992) for a large variety of vegetation 306 canopies, since no value was proposed for vineyards.

311
• Albedo a s was computed as a linear combination of ASTER waveband reflectances.
Due to deficient shortwave infrared channels in 2008, we used the generic two 313 channel (red and near infrared) formulation proposed by Weiss et al. (1999), 314 where the latter was beforehand linearly corrected against the ASTER seven

321
This yielded an accuracy of about 0.20 in absolute (10% in relative).

322
• Broadband emissivity ε s was computed as a linear combination of ASTER wave-323 band emissivities, by following Ogawa et al. (2003).

324
• Atmospheric irradiance R a was estimated from measurements of air temperature 325 and vapor pressure collected at the meteorological station ( § 2.2.1), by using the 326 formulation proposed by Brutsaert (1975) for clear-sky conditions.

329
• Surface temperature T s was that provided by ASTER official products ( § 2.2.2).

330
Finally, soil heat flux G 0 was derived from ASTER imagery using the formulation 331 proposed by Clothier (1986), where G 0 is expressed as a fraction of net radiation R n : This formulation was chosen since it was recommended by Moran et al. (1994)  satellite overpass R ni . This is performed using the ratio Cd i =R nd /R ni where R nd is daily 346 net radiation. Overall, the twofold procedure yields (L is latent heat of vaporization): S-SEBI was implemented along with ASTER imagery as following ( Figure 1).

348
The temperature -albedo space was characterized by calculating, for each 10 −3 width 349 albedo class, mean albedo and extreme temperatures. Next, T max upper limit (respec-350 tively T min lower limit) was determined from the linear regression between mean albedo 351 and maximum (respectively minimum) temperature of each albedo class. As proposed 352 by Roerink et al. (2000), T max upper limit was computed by excluding albedoes be-353 low the threshold value that discriminates evaporative and radiative regimes, where the 354 threshold albedo corresponds to the maximum temperature of the concave tempera-355 ture -albedo relationship (Bastiaanssen et al., 1998). To reduce noise influence, T min 356 lower limit was computed by including all albedoes (Verstraeten et al., 2005), rather 357 than excluding classes above the threshold value (Roerink et al., 2000). To include all 358 (temperature, albedo) pairs within the outer limits, an offset was added to T max (re-359 spectively subtracted to T min ) regression line, where the offset equaled the standard 360 deviation of maximum (respectively minimum) temperature over the albedo classes.

362
Following Gomez et al. (2005), the Cd i coefficient was calibrated on a daily basis 363 using ground based measurements of net radiation collected within Site 6 ( § 2.2.1). deriving daily ET from WDI along with optical remote sensing imagery is twofold.

383
By assuming air temperature T a is homogeneous within the study area, WDI is firstly 384 computed from the differences between surface temperature T sB and extreme tempera-385 tures (maximum T sC and minimum T sA ) that correspond to the same vegetation cover 386 fraction (see Figure 2 for the setting of T sA , T sB , T sC ): where the extreme temperatures T sC and T sA are determined from energy balance cal- where C p is air isobaric specific heat that depends on T a ; γ is psychrometric constant 406 that depends on T a and atmospheric pressure P a ; ∆ is the slope of the saturated vapor 407 pressure -temperature relation that depends on T a ; V PD is air vapor pressure deficit 408 that depends on T a and air humidity e a ; R n is net radiation; G 0 is soil heat flux; r c is 409 canopy resistance to vapor transport and r a is aerodynamic resistance derived from the 410 reference formulation of Thom (1975): where k is Von Kármán constant, u is wind speed, z r is reference height, d is displace-412 ment height, z 0m and z 0h are roughness lengths for momentum and heat respectively, 413 ψ m and ψ h are stability correction functions for momentum and heat respectively, and 414 L MO is Monin-Obukhov length defined as : where ρ a is air density, u * is friction velocity, g is gravity acceleration and H is sensible 416 heat flux.

417
When implementing WDI along with ASTER imagery, the first step consisted of

425
• C p , γ, ∆ and V PD were derived from measurements of air temperature T a , air 426 humidity e a and atmospheric pressure P a at the meteorological station ( § 2.2.1).

427
Wind speed u was also derived from measurements at the meteorological station ( § 2.2.1).

429
• Determining net radiation R n and soil heat flux G 0 for each of the four vertices Vertex 4 corresponds to a dry bare soil with an infinite canopy resistance equiv-448 alent to complete stomatal closure: r c = ∞.

449
• Aerodynamic resistance r a was also estimated for each of the four vertices. First, 450 roughness length for momentum and displacement height were set to 3/100 and 451 2/3 of canopy height, respectively, by following Sene (1994) and Ortega-Farias 452 et al. (2007). Second, the ratio z om /z oh was set to 100, as recommended by 453 Giordani et al. (1996) for vineyards. The kB −1 factor was then close to five, 454 which was lower than the value of eight obtained by Verhoef et al. (1996)  The different steps of the assessment strategy, along with the involved ground based 476 and remote sensing data, are listed in Table 2. 477 We firstly controlled the consistency of HYDRUS-1D simulations, by (i)       from S-SEBI was observed for all validation sites on these three dates.

643
When comparing HYDRUS-1D estimates of daily ET against EC based retrievals, 644 the larger differences observed for Site 6 were ascribed to a larger amount of data (con-  (Seguin et al., 1999;Kalma et al., 2008). As compared to outcomes re-655 ported in the literature for lower spatial extents and summer periods only (Trambouze 656 et al., 1998;Bsaibes, 2007), the current study enlarged the assessment of HYDRUS-1D whereas S-SEBI determines these extremes from variabilities captured within thermal 707 infrared imagery -concept of contextual model as proposed by Kalma et al. (2008). with automatized frequency analysis (Delenne et al., 2008).

727
From a temporal viewpoint, we noted S-SEBI performed worse than WDI for spe-   Figure 2). This was also observed on archive image of 8 Feb 2003, and was ascribed to errors on ground based and remote sensing data as suggested by Moran et al. (1996). . 44 39 Figure 1: Typical example of the scatterplots we obtained for the surface temperature (y-ordinate) versus albedo (x-abscissa) diagram used to compute evaporative fraction with S-SEBI. Full diamonds correspond to the (temperature-albedo) pairs for minimum temperature values of each albedo class, to be used for computing the lower limit through linear regression. Empty circles correspond to the (temperature-albedo) pairs for maximum temperature values of each albedo class, to be used for computing the upper limit through linear regression. Such a scatterplot was obtained for each of the 20 ASTER scenes we considered. Figure 2: Typical example of the scatterplots we obtained for the SAVI (y-ordinate) versus T s − T a (x-abscissa) diagram used to compute WDI, where SAVI is used for characterizing vegetation cover fraction, and T s − T a is surface -air temperature gradient. Numbers are vertices that allow to determine possible extremes (in dotted lines). Segment 1-2 (respectively 3-4) corresponds to full canopies (respectively bare soils). Segment 1-3 (respectively 2-4) corresponds to wet edge with maximum evapotranspiration (respectively dry edge with negligible evapotranspiration). A, B and C are possible water status for a vineyard plot with a given vegetation cover fraction, where A is for fully wet condition, C is for complete dry condition and B is an intermediate situation that represent the pixel status. Such a scatterplot was obtained for each of the 20 ASTER scenes we considered.  Figure 2). This was also observed on archive image of 8 Feb 2003, and was ascribed to errors on ground based and remote sensing data as suggested by Moran et al. (1996). Listing of the model inputs for WDI and S-SEBI, when differentiating inputs derived from (i) ground based measurements only, (ii) remotely sensed measurements only, (iii) both ground based and remotely sensed measurements and (iv) other sources of information (e.g. bibliographic). 51 4 Summary of statistics when validating ASTER / S-SEBI and ASTER / WDI estimates of daily ET against ground based estimates from HYDRUS-1D simulations: absolute and Relative Root Mean Square Error (RMSE and RRMSE), Slope and Offset from linear regression between ASTER estimates and HYDRUS-1D simulations, and determination coefficient R2. Results for Sites 1-7 include data from both 2007 and 2008. The last two rows indicate aggregate results for all sites partitioned by year. 52 Table 1: Main characteristics of the validation sites within the Peyne watershed. The "Devices" column indicates (i) the number of locations within each site for collecting soil moisture profile from Neutron Probe (NP) device, (ii) the collection of Eddy Covariance (EC) measurements when applicable, and (iii) the collection of net radiation (Rn) measurements when applicable. Each location of NP data was equipped for piezometric measurements. EC footprint over Site 6 encompassed both sub-sites 6N and 6S, whereas Rn footprint on Site 6 encompassed sub-site 6N only. Row orientation was characterized through row azimuth, where north is origin and positive values are clockwise. A soil was considered as shallow when the depth to parent material was less than 2.5 m. An absent watertable (respectively permanent watertable) means a piezometric level below 2.5 m depth (respectively above 2.5 m depth) throughout the year. A seasonal watertable means a piezometric level above 2.5 m depth during winter and spring.  • ASTER net radiation, WDI vertexes and ET 0 . Intercomparing S-SEBI and WDI estimates of daily ET at the extent of the Peyne regional watershed 9 ASTER imageries 11 ASTER imageries Table 3: Listing of the model inputs for WDI and S-SEBI, when differentiating inputs derived from (i) ground based measurements only, (ii) remotely sensed measurements only, (iii) both ground based and remotely sensed measurements and (iv) other sources of information (e.g. bibliographic). (2) instantaneous values at ASTER overpass (11:00 UTC)