Diurnal Dynamics of Wheat Evapotranspiration Derived from Ground-Based Thermal Imagery

The latent heat flux, one of the key components of the surface energy balance, can be inferred from remotely sensed thermal infrared data. However, discrepancies between modeled and observed evapotranspiration are large. Thermal cameras might provide a suitable tool for model evaluation under variable atmospheric conditions. Here, we evaluate the results from the Penman-Monteith, surface energy balance and Bowen ratio approaches, which estimate the diurnal course of latent heat fluxes at a ripe winter wheat stand using measured and modeled temperatures. Under overcast conditions, the models perform similarly, and radiometric image temperatures are linearly correlated with the inverted aerodynamic temperature. During clear sky conditions, the temperature of the wheat ear layer could be used to predict daytime turbulent fluxes (root mean squared error and mean absolute error: 20–35 W·m, r: 0.76–0.88), whereas spatially-averaged temperatures caused underestimation of pre-noon and overestimation of afternoon fluxes. OPEN ACCESS Remote Sens. 2014, 6 9776 Errors are dependent on the models’ ability to simulate diurnal hysteresis effects and are largest during intermittent clouds, due to the discrepancy between the timing of image capture and the time needed for the leaf-air-temperature gradient to adapt to changes in solar radiation. During such periods, we suggest using modeled surface temperatures for temporal upscaling and the validation of image data.


Introduction
Aircraft and satellite sensors operating in the thermal infrared (TIR) region of the electromagnetic spectrum provide spatially comprehensive information on the radiometric surface temperature (T rad ), which is derived from the thermal radiation emitted by the different surface components within the sensor's field of view.The diurnal and seasonal cycles of the surface temperature are affected by the relative efficiency of the components of the surface energy balance (SEB) in dissipating the available energy.Thus, SEB key components, such as the turbulent latent (λE) and sensible (H) heat flux, are frequently inferred from observations of the radiometric temperature over a wide range of spatial and temporal scales [1][2][3][4][5][6][7][8][9][10][11][12].
T rad is commonly applied as a proxy for the aerodynamic temperature (T aero ), a theoretical temperature at the thermal roughness height, which satisfies the bulk resistance formulation for the sensible heat transport [13,14].For real canopies, the quantitative relation between the radiometric and the aerodynamic temperature is, however, complex and not understood in detail, yet [11,12].Ambient conditions and canopy characteristics affect both variables differently, and they may have high spatial and temporal variability [11,13,15].Consequently, uncertainties introduce errors in H and λE estimates [12].Errors tend to be large for surfaces with partial vegetation cover, stressed vegetation and vegetation undergoing seasonal changes [8,[16][17][18].Furthermore, TIR-based studies often show a clear sky bias, since studies are predominantly carried out under clear sky conditions.Ground-based data for validating remotely sensed surface temperatures under variable cloud cover are needed [9,19].An understanding of the dynamic thermal adaptation of canopies, i.e., for the critical early and the late growth stages, will improve the interpretation of model behavior [20].Ground-based TIR cameras, allowing for a high observation frequency and for studying the spatial variability of temperatures, independently from ambient conditions, might provide a suitable tool for such studies.
The latent heat flux (λE) is one of the important SEB components where major attention is needed to improve the agreement between different climate models [21].In this study, we test the integration of temperatures derived from ground-based TIR imagery in three different single-source (which do not distinguish between soil evaporation and transpiration) evapotranspiration models (Bowen ratio, surface energy balance (SEB) and Penman-Monteith) for reproducing the diurnal course of λE.We use latent heat flux measurements provided by an eddy covariance (EC) station for validation.Models require data on net radiation, ground heat flux and standard meteorological parameters.In this study, we use micrometeorological observations from an eddy covariance station.
In particular, we aim to: (1) study the relation of the temperature of the different surface components (based on thermal infrared image, radiometer and thermocouple measurements) to the aerodynamic temperature derived from EC measurements; (2) test the suitability of the different model and temperature combinations to capture the diurnal cycle of the integrated latent heat flux from the soil and the canopy; and (3) study the effect of cloud cover on the temperatures and model results.Contrasting with radiation thermometer measurements, thermal imagery provides not only the range, but also the distribution and relative contribution of temperatures from different parts of the canopy.We hypothesize that such observations provide a deeper understanding of the performance of different algorithms under variable ambient conditions and help to explain model deviations from the reference flux measurements.

Research Site and Vegetation
The Merzenhausen research site was established by the Collaborative Research Center TR32 [22].The site is located at 93 m a.s.l.(50°55ʹ47ʹʹN, 06°17ʹ47ʹʹE) in the intensively agriculturally used area around Jülich (Germany) with a temperate climate [23].The mean annual air temperature is 10.2 °C, and the mean annual precipitation sum is 825 mm [24].Measurements are performed from 19 July to 24 July 2012.The observed area, located within a 200 m by 360 m-wide field, is cropped with winter wheat (Triticum aestivum L. cv.Tobak; sowing density: 360 plants•m −2 ).Non-destructive LAI measurements (cf.Section 2.5) indicate an average leaf area index (LAI) of 4.0 (ranging from 3.3 to 5.4) and a vegetation height of 0.87 m (ranging from 0.84 to 0.9 m) during the measurement period.The observed winter wheat is in the stage of ripening (kernel formation completed) with the flag leaves being partly senescent, indicating that the bottom of the crop canopy is senescent [25] and that the plants are reaching maturity [26].The underlying soil texture is classified as silt loam.The field was harvested on 2 August.

Experimental Set-Up
We mounted a FLIRA315 thermal imaging camera (FLIR Systems Inc., Wilsonville, OR, USA), sensitive at wavelengths of 7.5-13 μm, in a weatherproof housing at a height of 2.2 m over the soil surface.The accuracy reported for the camera is ±2 °C with a thermal sensitivity < 0.5 °C.Based on the results from preliminary field studies, the camera is tilted 36° below a horizontal line.The tilt angle was chosen to avoid the major influence of the soil on remotely sensed surface temperatures of a row-structured canopy with an erectophile leaf angle distribution.The camera's field of view (FOV, 25° × 18.8°) is pointed towards the SE, observing an area of 1 m × 1 m (Figure 1).Images (16 bit, 320 × 240 pixels) are captured at a five-minute time interval.The camera is connected to a fanless embedded box PC (Bressner Technologies, Gröbenzell, Germany) running the control software (IRControl, Automation Technology, Bad Oldesloe, Germany).

Reference Temperature Measurements
The set-up for obtaining reference leaf temperatures consists of four type T (temperature) copper-constantan thermocouples with Teflon insulation (Omega Engineering, Stamford, CT, USA) connected to a temperature data adapting device (RedLab TEMP, Meilhaus Electronic, Alling, Germany).The sensors are attached to the lower side of winter wheat leaves using semi-permeable tape.The typical and the maximum error for the thermocouple measurements (including the error for the cold junction compensation sensor) are ±0.26K and ±0.71 K, respectively.An infrared radiometer (Apogee SI-121, Apogee Instruments, Logan, UT, USA), sensitive at wavelengths of 8-14 μm (with a FOV of 18°), is mounted directly below the infrared camera housing and connected to a data logger (Combilog, Theodor Friedrichs & Co., Borgfelde, Germany).Using the same inclination as for the camera, it records the radiometric surface temperature of an area measuring 1.2 m in maximum length and 0.9 m in maximum width.The accuracy reported for the infrared radiometer is ±0.2K (for object temperatures between −30 °C and +65 °C).

Thermal Image Data Processing
We calculated the percentage area contribution of the different surface components to the total image area by means of an unobserved K-means classification (ENVI software, Research Systems Inc., Boulder, CO, USA) of manually selected daytime and nighttime images.Wheat ears (Class I) cover 30%, mixed stem and leaf surfaces (Class II) 44% and the sum of deeper shaded areas and underlying soil (Class III, describing the background) 26% of the image.We manually selected image pixels corresponding to the surface Classes I and III and, using the software, IRControl (Automation Technology, Bad Oldesloe, Germany), derive the corresponding temperature values (Class I, T rad,WE ; Class III, T rad,BG ) for each time step.The temperature corresponding to Class II (T rad,mean ) is assumed to be represented by the average image temperatures.This is motivated by the fact that stem and leaf surfaces form the largest area in the scenes and that their temperatures are usual identical to the average.Thermal infrared (8-14 μm) emissivity values between 0.93 for sparse canopies (including the effect of the underlying soil) and/or drier vegetation and 0.99 for dense, healthy vegetation canopies (with cavities trapping the radiation) have been reported [13,[27][28][29][30]. Generally, the thermal emissivity of vegetation is assumed to be higher than 0.95 [30,31].Considering the partly senescent leaves and the high fraction of dry wheat ears in the images, we apply an emissivity of 0.95 (the mean value of reported values, the lowest boundary reported for the thermal emissivity of vegetation).The applied linear spatial averaging of temperatures is less accurate than a retrieval of raw thermal emission for each pixel.However, since the manufacturer's software does not provide radiances, we are confined to this procedure.The software corrects for the so-called "background radiation" (emitted by the surrounding objects and further reflected by the object) and the radiation emitted by the atmosphere.Considering the relatively high emissivity of the canopy and, consequently, the small effect of the background radiation on temperature estimates, the "background" or "reflected temperature" was approximated using measured air temperatures.The atmospheric transmission along the observation path is calculated by providing measured air temperature (in 5 °C intervals ranging from 5 °C to 30 °C) and relative humidity values (35%, 60%, 80%, 90% and 95%) from standard meteorological measurements for each image.Differences between the raw and the temperature-and humidity-corrected temperature values range up to ±1.5 K.We extracted the surface temperatures corresponding to the surface components (T rad,WE , T rad,BG ) and the image mean (T rad,mean ), minimum (T rad,min ) and maximum temperatures (T rad,max ).For the purpose of this study, to model the turbulent latent heat flux measured by eddy covariance (EC) instruments, half-hourly average values are computed.

Ecosystem Gas Exchange, Canopy Architecture and Meteorological Data
A micrometeorological station with turbulence sensors for measuring CO 2 and H 2 O fluxes by the eddy covariance method is permanently located near the center of the wheat field, at a distance of approximately 12 m to the IR camera setup.Horizontal and vertical wind velocity are measured by a sonic anemometer (CSAT3, Campbell Scientific, Logan, UT, USA).CO 2 and H 2 O concentrations are measured using an open-path infrared gas analyzer (LI7500, LI-COR, Lincoln, NE, USA) at a measurement height of 2 m [32].Ground heat flux is measured by means of heat flux plates (Hukseflux Thermal Sensors B.V., Delft, The Netherlands) located at depths of 0.02-0.06m and repeated measurements of near-surface heat capacity density (KD2 pro, Decagon, Pullman, WA, USA) [23].Net and global radiation is measured with a net radiation sensor (NR-Lite, Kipp & Zonen, Delft, The Netherlands).Details on instruments, footprint analyses, data processing and quality filtering for estimating the turbulent latent and sensible heat fluxes are given in [24,32,33].All data are aggregated to 30-min intervals.Field measurements of the Leaf Area Index (LAI) and canopy height were performed on 23 July.Nondestructive LAI measurements are obtained using a SunScan Canopy Analysis System (Delta T Devices, Cambridge, U.K.).For obtaining representative field estimates, the data from eight plot measurements are averaged to mean values.The cloud fraction (CF) is provided by a Total Sky Imager (TSI 880, YES Inc., Turners Falls, MA, USA) located at the Jülich Observatory for Cloud Evolution (JOYCE) located approximately 8 km from the research site.To test for inconsistencies in cloud cover between both sites, we investigate TSI images for the difference, ∆CF, of CF estimates in the sector towards Merzenhausen compared to the fraction around Jülich.The cloud fraction is determined from the images at an elevation of 30°, thus with a cloud base height at ~2 km, giving information on CF at a distance of ~4 km.The 10th and 90th percentiles of ∆CF lie at −0.12 and 0.14, respectively, indicating that most of the time CF towards Merzenhausen did not deviate by more than one octa from that observed at Jülich.During intermittent cloud cover, 80% of the ∆CF lies within −0.12 to 0.32 compared to −0.02 to 0.13 for clear sky and −0.02 to 0.05 for overcast conditions.We further compare daytime half hour averages of global radiation measured at Jülich and at Merzenhausen during the observation days.Differences (∆G) can be as large as −140 W•m −2 and 220 W•m −2 for single half hour averages, in general; however, ∆G is much smaller, with 80% lying within the interval −31 to 63 W•m −2 and 50% in the interval −8 to 30 W•m −2 .Differences are, as expected, largest during intermittent cloud cover (0.2 < CF < 0.8), with 80% of ∆G within −69 to 79 W•m −2 , and smaller during overcast (−49 to 34 W•m −2 ) and clear sky (−17 to 35 W•m −2 ).Thus, radiation differences are, most of the time, much smaller than global radiation, and we therefore believe that the CF determined in Jülich is a sufficiently good estimate for Merzenhausen.

Bowen ratio
Bowen ratio method, TIR-based calculation of heat and water vapor gradients.

Penman Monteith
TIR-based parameterization of the (bulk) canopy resistance term of the PM equation.
Air density, net radiation, ground heat flux, humidity, aerodynamic resistance, air temperature, surface temperature [41][42][43] The surface energy budget, measured by the micrometeorological instruments, is usually not balanced.The evaluation of the energy balance closure was beyond the scope of our study.However, to avoid an effect of the residuum on model evaluations, we subtract it from the measured net radiation and derive the corrected net radiation Rn c (W•m −2 ): where Rn (W•m −2 ) is the measured available net radiation, λE (W•m −2 ) the latent heat flux, H the sensible heat flux (W•m −2 ) and G (W•m −2 ) the ground heat flux [42,43].The reason to correct Rn rather than any of the turbulent fluxes is that, following [44,45], the energy balance is mainly not closed due to the differences in the scales of the fluxes and the differences in the instrument footprints (and less due to measurement errors).Accordingly, one cannot decide which of the measured fluxes has to be corrected.To avoid a numeric change of the absolute values of the measured latent and sensible heat fluxes, the adjustment of the energy available to the system is considered the most reasonable way.Note that G makes a minor contribution to the overall surface energy balance with a mean absolute value of 18 W•m −2 .In the following section, the tested model approaches are described briefly.

Modeling Approach (a) "SEB"
The bulk aerodynamic or surface energy balance (SEB) method uses the bulk resistance formulation for the sensible heat flux (H, W•m −2 ), which is based on the near surface gradient of temperatures at the surface (T aero , °C) and the atmosphere (T air , °C) [46].The radiometric temperature (T rad , °C) can be used as a proxy for the aerodynamic temperature T aero [40]: where C p (1010 J•kg −1 •K −1 ) is the specific heat capacity of the air at constant pressure, ρ (kg•m −3 ) the air density and r a (s•m −1 ) the aerodynamic resistance for heat transfer from the surface to the reference height (z: 2 m).The aerodynamic resistance of the canopy is calculated from: where zm (2 m) is the height of wind measurements, zh (2 m) the height of humidity measurements, d (m) the displacement height, calculated from 1/3 × crop height, zom (m) the roughness length governing the transfer of momentum (m), calculated from 0.123 × crop height, zoh the roughness length governing the transfer of heat and water vapor, calculated from 0.1 × zom, k (0.41) the von Kármán constant and u (m•s −1 ) the wind speed.Most studies perform the stability correction of the aerodynamic resistance using the Monin-Obukhov length, calculated from measured sensible heat flux.However, to avoid the dependency of the modeling approaches on measured turbulent fluxes, we use the Richardson number (Ri) for stability correction [47,48]: For Ri < −0.008 and Ri > 0.008, the aerodynamic resistance is calculated from: where g is the acceleration of gravity (9.81 m•s −2 ), T air (Kelvin) the air temperature, T rad,mean (Kelvin) the mean remotely sensed surface temperature derived from TIR image data, z (m) the measurement height and T av (Kelvin) the average temperature of the canopy and the air temperature.Subsequently, the latent heat flux (λE, W•m −2 ) can be calculated as the residual from the energy balance: To study the agreement between the unknown aerodynamic (T aero ) and the remotely sensed surface temperature, we inverted the bulk formulation for the sensible heat flux using the sensible heat flux derived from eddy covariance measurements to compute the "inverted aerodynamic temperature" (T aero,inv ) [49].

Modeling Approach (b) "Bowen Ratio"
The Bowen ratio energy balance method is a widely used micrometeorological method to estimate latent heat fluxes because of its simplicity and robustness [50,51].At the ground level, it has been applied routinely to estimate latent and sensible heat fluxes for various land surface types, such as water, grassland, crops or forests [52].However, with the increasing availability of atmospheric water and temperature profiles at larger scales (e.g., using polar orbiting satellite systems, like the Atmospheric Infrared Sounder (AIRS)), its application to satellite-based estimates of λE becomes feasible, as demonstrated by [53].No information on surface characteristics is needed.The method is based on a rearrangement of the surface energy balance and the theory of the flux-gradient relationship [35].For the calculation, humidity and temperature observations from two different heights, net radiation and ground heat flux estimates are required.The approach is based on the assumption of the similarity of the eddy diffusivities of heat and water vapor and of a balanced energy budget.Assuming that the net horizontal advection of energy can be neglected, we apply the simplified version to compute the Bowen ratio β and the latent heat flux as follows: with ea surf as the actual vapor pressure at the surface temperature, ea the actual vapor pressure at the air temperature (kPa), γ the psychrometric constant and λ (J•kg −1 ) the latent heat of vaporization.

Modeling Approach (c) "T Modeled"
Approach (c) is similar to (a), but based on modeled surface temperatures (T modeled ) instead of T rad .We applied the modified version of the Todorovic approach [54], suggested by [38], to derive T modeled .The Todorovic model is based on the assumption that the additional energy needed by the water vapor flux to overcome the stomatal resistance (r s ) leads to a decrease in λE, which is balanced by an increase in the amount of H [37]. Thus, the additional H component describes the difference in energy needed for the water vapor flux at "wet" conditions (r s = 0) to the energy needed at actual conditions (r s > 0).The additional amount of energy leads to an increase in the canopy temperature (t, °C), causing a small increase in the sensible heat flux.The modified equation [38] includes the response of the stomatal conductance to the VPD [55]: where Δ (kPa•°C −1 ) is the slope of the temperature-saturation vapor pressure relationship, γ (kPa•°C −1 ) the psychrometric constant, VPD (kPa) the water vapor pressure deficit of the air and VPD 0 (1.5 kPa) a canopy-dependent parameter characterizing the curve of response of the stomatal conductance to the humidity deficit of the air [55][56][57].Several studies showed that stomatal closure is induced at leaf-to-air vapor pressure differences of ~1.5 kPa [58], and this value has been successfully applied for wheat canopies [38,57,59].Based on the equation for the potential evapotranspiration λE p , the potential sensible heat flux H p at potential evapotranspiration is calculated as: Accordingly, the canopy temperature at potential evapotranspiration T Cp (°C) can be derived from: The actual canopy temperature, T modeled , is computed by: Subsequently, using T modeled , the actual sensible heat flux is derived from the bulk formulation of the sensible heat flux (Equation ( 2)), and the latent heat flux is derived from the SEB (Equation ( 7)), as described for (a).

Modeling Approach (d) "Penman-Monteith"
By combining energy balance and mass transfer approaches, the Penman-Monteith (PM) equation [41,60] allows for computing evaporative fluxes from standard meteorological data.An important scaling parameter of the PM equation is the surface resistance term, r C (s•m −1 ), describing the (bulk) canopy resistance that controls the transfer of water vapor from the leaves within the canopy to the surrounding air.This term can be derived from the canopy surface temperature.For computing r C , the following equation is applied [38,61]: Subsequently, the latent heat flux is derived from the Penman-Monteith (PM) equation [39,41]:

Data Analyses and Model Error
Modeling and data analyses are performed using the R statistical programming language [62].The latent heat flux measured by the EC instruments is used as a reference data set.The root mean square error (RMSE) is based on the sum of squared errors and commonly applied in climatological studies as a measure of the goodness of predicted values compared to the reference data values.Because the RMSE is sensitive to outliers, its interpretation can be misleading [63].However, for comparison with the results from other studies, we computed the RMSE: where λE is the observed and λ the predicted latent heat flux (W•m −2 ) and n the number of observations.For estimating the overall precision of predicted latent heat fluxes, the mean absolute error (MAE), describing the absolute deviation of the modeled from measured data values with positive signs for both over-and under-estimations, is calculated as:

Ambient Conditions and Surface Temperatures
Micrometeorological data and the components of the surface energy balance (Figure 2) demonstrate the transition from cloudy conditions with occasional rainfall (19-21 July) to clear sky conditions with large temperature amplitudes (22)(23)(24).The similarity of the latent heat flux (λE) and the sensible heat flux (H) diurnal amplitudes indicate progressing senescence of the wheat plants [64] and hints at the progressing depletion of available soil water [65].Daytime evaporative fractions, defined as the ratio of λE to available energy (Rn-G) [66], range from ~0.3 to 0.5 during clear sky conditions (22-24 July) and from ~0.4 to 0.8 during cloudy conditions (19-21 July) and show the typical parabolic shape during the diurnal cycle [67].The horizontal wind speed was low to moderate (~2 m•s −1 ) at reference height z (2 m), with the highest (~5 m•s −1 ) and lowest (~ 1 m•s −1 ) maximum daytime values recorded on 19 July and 22 July, respectively (Figure 2).Dew point temperature calculations indicate that dew was likely during morning hours on 22 July.Image-based surface temperatures (T rad ) range from 3.4 °C to 36.4 °C with a mean value of ~17 °C (Table 2) and a standard deviation of ~6 °C.The infrared radiometer and the thermocouple records indicate a slightly smaller temperature range (3.9 °C to 28 °C and 4.7 °C to 30.1 °C, respectively).The statistics of T rad,WE and T rad,BG demonstrate differences in the thermal response, such as a higher thermal conductance of ripe wheat ears and dry leaves in the upper canopy layer, causing a faster response to changes in net radiation and a larger temperature amplitude, compared to the soil and deeper canopy layers (Table 2).
Pronounced daytime warming and nighttime cooling of the wheat ears (compared to T rad,BG ) are demonstrated in the sample images (Figure 1).The image-based maximum temperatures clearly exceed those recorded by thermocouples (maximum temperature differences of +10 °C between T rad,max and T leaf on July 22), whereas the spatially integrated temperatures from the radiometer are significantly lower than T leaf , i.e., at night (Figures 2 and 3).Image-based surface temperatures generally exceed air temperatures (T air ) by 2 °C-3 °C (cloudy days) up to 5 °C (clear days), but around sunset and at night (Figure 2b).spike temperatures of non-irrigated cultivars can exceed ambient air temperatures after having reached physiological maturity by up to 10 °C [68].Their position ensures that they receive a minimum of shading; thus, they are the major structures for absorbing light.Further, they have a low surface area-to-volume ratio, thus a limited ability to dissipate heat.We are, however, aware that uncertainties in the derived temperature values due to the surface thermal emissivity, changes in the reflective temperature (cf.Section 2.4) and the sun-sensor geometry have to be considered.The reflectivity of ripe wheat ears is high, resulting in a low absorption of net radiation.We account for senescent parts of the canopy by applying a comparable low emissivity value (ε = 0.95).However, the emissivity of dry leaves can be even lower [71], causing an underestimation of calculated surface temperatures.Furthermore, studies report decreasing radiometric temperatures at oblique viewing angles [72].Directional effects are, however, dependent on incident radiation (by modifying the contribution of observed sunlit and shaded leaves), the LAI and leaf angle distribution parameters [6].
Table 2. Statistics of the 30-minute interval aggregated surface temperatures (T, °C) derived from thermal imagery (T rad,mean = average; T rad,min = minimum; T rad,max = maximum; T rad,BG = temperature of the deeper shaded areas and underlying soil; T rad,WE = temperature of wheat ears), from Thermocouples 1-4 (T leaf ) and from the IR radiometer (T IRR ), of the inverted aerodynamic temperature (cf.Section 3.2) (T aero,inv ) and the modeled surface temperatures (T modeled , cf.Section 3.1.2);number of observations = 263.Despite these uncertainties, the correlation between image, point and spatially integrated temperature data is high.The highly linear correlation between T leaf and T rad,mean (Figure 3a) indicates that the assumed emissivity, temperature and humidity corrections applied during image processing (cf.Section 2.5) are correct.Comparison of T IRR with T rad,mean (Figure 3d) shows an offset of about −1 °C and a slope of ~0.96.As both measurement rely on the same physical principle, the differences should be negligible.They could be caused by small differences in the field of view and different correction algorithms (cf.Section 2.5).However, the most likely reasons are: (1) the default assumption of an emissivity of one for the derivation of T IRR ; and (2) the immanent spatial averaging of radiation in the radiometer.The latter may lead to different results than arithmetically averaging temperatures from a camera image, as thermal radiation is not linearly related to the temperature.The inconsistency between T IRR and T rad,mean leads to equivalent offsets and smaller slopes in the relation between T rad,WE and T rad,BG to T IRR (Figure 3e,f) compared to the relation to T leaf (Figure 3b,c).Linear regression fits show an increasing divergence at higher and lower temperature values (Figure 3) and for the nighttime data.T rad,WE is higher at temperatures ≥ 20 °C and T rad,BG lower at temperatures ≤ 15 °C compared to T leaf and T IRR .Cloud cover has a low impact on the relation between measured temperatures (Figure 3); however, as expected, under overcast sky conditions (CF > 0.8), the linearity is highest and the temperature range is smallest.

Figure 3.
Linear regression fit between averaged leaf temperature recorded from thermocouples (Tleaf, Figure 3a-c) and between surface T recorded from an IR radiometer (TIRR, Figure 3d-f) with remotely sensed surface temperatures (as for Table 2).The solid line represents the line of identity, while the dashed line represents the linear regression line.CF = cloud fraction; nighttime: 21.00 h-03.00 h UTC.

Radiometric Surface Temperatures versus the Inverted Aerodynamic Temperature
During overcast conditions (cloud fraction > 0.8) and for temperatures below ~20 °C, the image-based temperatures are linearly related to the inverted aerodynamic temperature (T aero,inv , cf.Section 3.1), indicating a high suitability of T rad to predict the corresponding turbulent sensible heat flux (Figure 4).T rad,WE , receiving a maximum of solar radiation, responds faster to changes in net radiation than T rad,mean and T rad,BG .As a consequence, T rad,WE is more strongly correlated to T aero,inv at temperatures > 20 °C and during clear sky conditions (CF < 0.2), whereas T rad,mean and T rad,BG are underestimating T aero,inv (Figure 4).As a consequence, for the current viewing geometry and the applied emissivity and background temperature settings (cf.Section 2.4), T rad,WE is more suitable to predict the turbulent heat fluxes, giving preference to the use of a single surface component over spatially-averaged temperatures during such periods.Contrastingly, under overcast conditions, the use of spatially-averaged temperatures, based on thermal camera imagery, performs similar.Note that awnless ears (as in this study) respond more rapidly to changes in incoming radiation [73], whereas the presence of awns can cause a different microclimate around the canopy [74,75], restricting a generalization of our observations.Scatter plots indicate the weakest correlation of T rad and T aero,inv during intermittent clouds (Figure 4) when T rad,WE and T rad,mean mostly exceed T aero,inv (cf.Section 4.2.3).This can be related to the occurrence of frequent outliers.The most likely cause for these outliers is the occurrence of sunny intervals, which cause rapid fluctuations of turbulent fluxes in response to changes in available energy.Simultaneously, a rapid increase of leaf temperatures is induced, which depends on the plant-specific time constant for thermal equilibration [76].This is, however, not captured using 30 min-averaged surface temperatures (raw data obtained at a resolution of 5 min).The timing of thermal image capture affects the suitability of derived T rad values to predict the measured flux data, which highlights the need for longer continuous observations to capture near-steady-state canopy temperatures and resulting temperature gradients [77][78][79].Consequently, especially the relation between T aero,inv and T rad,WE is nonlinear during intermittent clouds (Figure 4).

Suitability of Model-Temperature Combinations to Predict Sub-Daily λE Dynamics
The Taylor diagram [80] graphically summarizes how closely the predicted values match the observed values in terms of their correlation, RMSE and standard deviation (Figure 5).Nighttime data is excluded from the diagrams to account for the uncertainty of nocturnal EC flux measurements.The EC technique measures the turbulent transfer of heat and water vapor from the soil and the canopy integrated over the footprint area of the eddy covariance tower.Thus, the unknown contribution of soil evaporation to the total latent heat flux restricts the use of EC-based observations as validation data for evapotranspiration models, which are based on canopy surface temperatures (as in this study).Further, the uncertainty caused by the forced energy balance closure has to be taken into account (adjusted available energy, cf.Section 2, Equation ( 1)).The diagram clearly shows that measured latent heat fluxes are most accurately simulated when using T rad,WE, describing the temperature of wheat ears.These results agree with the finding on the linear correlation between T rad,WE with T aero,inv (cf.Section 3.2).The diagram compares correlation values (dotted black lines), the root-mean-square difference (solid dark red lines) and the standard variation (distance from the origin) of modeled and measured λE.Solid black line: standard deviation of measured λE.Input surface temperatures (T) as for Table 2.For model abbreviations, cf.Table 1.
On the one hand, this is related to the fast response of wheat ears to solar radiation (cf.Sections 3.1 and 3.2).On the other hand, there might be physiological reasons, related to soil water availability, for these findings, as described in the following text.During the observation period, the overall variability of soil water content is low; however, the fraction of available water (fAW) decreases to values < 50% (with fAW defined as the ratio of the difference between measured volumetric soil water content (SWCvol) and SWCvol at the wilting point to the difference between SWCvol at field capacity and SWCvol at the wilting point).For a winter wheat stand in China, [81] found that at an fAW below 50%, the plants experienced "water stress", and stomatal conductance (gc) was decreasing linearly with fAW.Reports on the overall contribution of ear activity differ [82][83][84].However, it was shown that the rate of ear photosynthesis can exceed that of the flag leaf, especially under water stress conditions [83,85].Thus, under the observed moisture and phenological conditions, ears probably remain active for a longer time period and maintain a higher water content than flag leaves [84,86,87].These studies give further reason to our findings on the suitability of ear temperatures to model turbulent fluxes during the observation period.RMSE and MAE for the complete data set and for different subsets (nighttime: 21.00 h-03.00 h UTC, daytime: 6 h-18.00 h UTC, cloud fraction ≥ 0.8, cloud fraction ≥ 0.2-0.8,cloud fraction < 0.2) are computed (Tables 3 and 4).Calculated errors using T rad,WE and T rad,mean (20-35 W•m −2 ) are well below the values published for other studies (~30-80 W•m −2 , see [3,4] for recent reviews and [88][89][90][91][92] for studies using ground-based data) and below the range of errors reported for eddy covariance (EC) measurements [42], which are here used as reference data.The highest errors (up to 90%) in latent heat flux estimates (using the surface temperature to predict the sensible heat flux) have been reported for natural semiarid areas [8].Most studies have been performed using aircraft or satellite data [4], thus restricting a comparison to our results.Further, it is worth noting that most TIR-based model evaluation studies exclude data from periods with higher fractions of cloud cover.
The results from the approaches "Bowen ratio" (b) and "T modeled" (c) show high correlation and low RMSE values and are insensitive to the input temperatures (Figure 4, Tables 3 and 4), providing robust estimates of evapotranspiration.Results suggest further investigation of the use of the Bowen ratio method for large-scale flux observations [53], despite its limitations due to the underlying assumptions [35,43].Note that we accounted for the residual of the energy balance closure, which leads to inaccuracies of the Bowen ratio method.
The approaches "Penman-Monteith" (d) and "SEB" (a) are more sensitive to the input temperature, especially during intermittent cloud cover and clear sky conditions.However, when parameterized with T rad,WE , the "Penman-Monteith" (d) and "SEB (a) algorithms show the lowest errors (RMSE and MAE ~ 20, Table 3) and the highest coefficient of determination values of the linear regression fits (r 2 = 0.83-0.85(daytime data), Figure 6).Table 3. Root mean squared error and mean absolute error (in parentheses) calculated for the difference found among eddy covariance-based measurements and model estimates of the latent heat flux using T rad,WE .;CF = cloud fraction, df = degrees of freedom.

The Impact of Cloud Cover on Model Results
RMSE and MAE are highest during intermittent clouds (cloud fraction ≥ 0.2-0.8),with the RMSE values nearly doubled.This can be related to the occurrence of frequent outliers during sunny intervals (Figure 6, cf.Section 3.2) and the limited ability of images taken at specific intervals (5 min in this study) to capture rapid fluctuations of turbulent fluxes.The differences between the MAE for different cloud conditions are lower compared to those of the RMSE (Tables 3 and 4).This indicates that the models' ability to resolve the diurnal cycle, rather than the ability to simulate the absolute flux magnitude, is diminished, the latter being more critical for temporal upscaling.Note that, resulting from the distance between the sensor and study site (cf.Section 2.5), the highest uncertainty in cloud fraction estimates occurs during intermittent cloud cover.Such inconsistencies in cloud fraction data might contribute to model errors during such periods.This case study demonstrates the complex and interrelated effects of model structure, input temperature and cloud cover on modeled fluxes.Errors in land surface models' to resolve sub-diurnal patterns of transpiration are often related to their diminished ability to simulate the natural diurnal hysteresis [93].We hypothesize that the results demonstrate the model and surface temperature-specific ability to reproduce such hysteresis effects, which are strongly affected by local cloud cover.At the observation scale, the phase angle difference between VPD and net radiation causes a phase angle difference between VPD and latent heat flux [94].Hysteresis arises from the well-known asymmetry between morning and afternoon data: for a given VPD in the morning, the stomatal conductance is higher compared to the afternoon, producing a clockwise rotation in the hysteresis curve (Figure 7).The decreasing slope at VPD values ≥ 1.5 kPa indicates the initiation of stomatal closure, giving reason to the value applied for VPD 0 in Equation (10) (cf.Section 3.3).Generally, all models capture the hysteresis effect (Figure 7).The surface temperature integrates information on stomatal closure causing a similar hysteresis effect for T rad and latent heat fluxes.The diurnal variability of this effect (and the resulting half-hourly magnitude of the hysteresis) is, however, reproduced best by those models that explicitly use VPD data (Bowen ratio, T modeled and Penman-Monteith), highlighting its importance on the partitioning of available energy during the stage of crop maturity.A sample diurnal course of the hysteresis between VPD and latent heat fluxes shows the overestimation of afternoon fluxes, i.e., for the SEB approach, followed by the Penman-Monteith equation (Figure 7).These observations explain the lower performance of the SEB and the Penman-Monteith approach during clear sky conditions (where VPD becomes relevant and stomatal closure can be initiated), as demonstrated in the RMSE and MAE statistics (Tables 3 and 4).Under prevailing moisture conditions and for the current plant physiological stage, using T rad,WE allows capturing these effects, whereas image-averaged temperatures (T rad,mean , T IRR ) underestimate pre-noon and overestimate afternoon fluxes during clear sky conditions.
The strong effect of the interrelations between model structure and cloud cover restricts the validity of our findings and highlights the need for longer time series obtained under variable cloud cover and for different plant phenological stages to study the performance of model and input temperature combinations.During the six observation days, fluxes are mainly driven by radiation.Adjusting the net radiation for energy balance closure minimizes errors for approaches that derive the latent heat flux as the residual from the energy balance (SEB approach).Therefore, error statistics shown have to be interpreted carefully.Further, note that errors that arise from the Sun-sensor geometry and the anisotropic behavior of thermal emission (i.e., during clear sky conditions) are not quantified.The fraction of the sunlit or and shaded side of the wheat ears seen by the sensor changes during the diurnal cycle [95].More generally, average climate site conditions and genotype-specific characteristics, such as the leaf water content, stomatal behavior, morphological traits, the presence of awns, canopy structure and phenological stage, strongly affect leaf and canopy temperatures.Thus, a generalization of observed temperatures and model results is strongly limited.On the other hand, such genotype-specific characteristics allow for using thermal imagery as a tool for observing plant response to stress as part of non-invasive phenotyping approaches [96].Results from our study contribute to the overall understanding of evapotranspiration model behavior for variable cloud cover conditions, when parameterized with remotely sensed surface temperatures.With increasing variability of precipitation and availability of water in the future, more comparative studies on model-specific abilities to accurately quantify evapotranspiration are needed.The integration of surface temperatures in large-scale estimates of the water balance, especially for irrigation purposes, has been used for several years.However, the availability of affordable ground-based thermal imagers and processing software has increased the possibilities of using such images to study plant response to stress and its impact on turbulent fluxes from a canopy [92].In the framework of currently developing plant phenotyping networks (e.g., the European Plant Phenotyping Network (EPPN) [97], the significance of studying canopy temperatures has been recognized.We suggest using ground-based thermal image time series of a variety of plant types and under field conditions to further evaluate the suitability of surface temperature-evapotranspiration model combinations to quantify sub-diurnal dynamics of the water balance.

Modeled Surface Temperatures
Results from the modified Todorovic approach ("T modeled" (c)) [49] outperform most approaches, i.e., with respect to modeled standard deviations (note that the position in the Taylor diagram slightly changes due to the influence of surface temperatures on the aerodynamic resistance.This indicates the superiority of modeled surface temperatures (T ,modeled ) over observed temperatures (T rad ) during certain time periods.T modeled is closely related to the inverted aerodynamic temperature; however, it overestimates this during the day (Figure 8), causing an overestimation of the daytime sensible heat flux.The linear regression fit between modeled and measured daytime (6 h-18.00 h UTC) λE shows the resulting consequent underestimation of measured latent fluxes (Figure 6, MAE ~ 25 W•m −2 ).Since T modeled are calculated from micrometeorological data with the raw values obtained at a very high temporal resolution, they are more suitable to capture nonlinear processes and short-term fluctuations.This is the most likely reason for the high agreement between modeled and measured standard deviations found for the "T modeled" (c) approach.We suggest using modeled temperatures for validating thermal infrared data obtained at a lower temporal resolution (given the availability of meteorological data in an adequate temporal resolution).

Conclusions
In this study, we demonstrate the potential and the limitations of ground-based thermal images to study the behavior of single-source evapotranspiration models under variable cloud cover, when parameterized with remotely sensed surface temperatures.At a ripe winter wheat stand, the models performed best at resolving the evapotranspiration diurnal cycle when parameterized with the temperature of the wheat ears (root mean squared error and mean absolute error between modeled and measured fluxes: 20-35 W•m −2 , r 2 : 0.76−0.88)and if integrating vapor pressure deficit data to account for natural hysteresis effects.Under clear sky conditions and for the prevailing site characteristics, which cause strong temperature differences between the different canopy layers and the soil, the simulation and the upscaling of turbulent fluxes can benefit from single-surface component temperatures.In contrast, under overcast conditions, spatially-averaged data provide similar model results.Large errors during intermittent clouds are related to the timing of thermal image capture, which affects the ability of derived surface temperatures to capture steady-state conditions of temperature gradients and causes a nonlinear relation between image-based radiometric and the aerodynamic temperature.Under such conditions, the use of discontinuous small-scale measurements to predict the complex, nonlinear processes of the surface energy balance is strongly limited, and consequently, the use of temperature-robust models is required.This further implicates that, for the temporal upscaling and the validation of remotely sensed surface temperatures, modeled temperatures based on standard meteorological data are more suitable.For future studies, we suggest using image time series obtained under a large range of environmental conditions and for a variety of plant types to further evaluate the performance of different surface temperature-model combinations to quantify sub-diurnal dynamics of the water balance and to study plant responses to abiotic stress.

Figure 2 .
Figure 2. (a) Remotely sensed surface temperature by thermal imaging (T, °C): T rad,mean = average surface T; T rad,max = maximum surface T; T leaf = averaged leaf T recorded from thermocouples; T air = air T (°C); T br(IR) = surface T recorded from an IR radiometer; T soil = soil T at −0.02 m; SWC = soil water content (Vol fraction).The dashed vertical lines represent the timing of sunrise and sunset, respectively.(b) Is as in (a), but for the temperature difference (ΔTemperature) with air temperature (°C); T rad,WE = T of wheat ears; T rad,BG = T of deeper shaded areas and underlying soil.(c) Net radiation (Rn, W•m −2 ), ground heat flux (G, W•m −2 ), eddy covariance (EC)-based sensible heat flux (H, W•m −2 ) and EC-based latent heat flux (λE, W•m −2 ).(d) Horizontal wind speed (m•s −1 ) and relative humidity (RH, %)

Figure 4 .
Figure 4. Scatter plot between image-based daytime surface temperature (T rad,mean , T rad,BG , T rad,WE ) and the inverted aerodynamic temperature T aero,inv , based on eddy covariance measurements of the sensible heat flux.CF = cloud fraction.

Figure 5 .
Figure 5.Taylor diagram displaying the statistical comparison of modeled daytime λE (for model-specific symbols, see the figure legend) to the observed daytime latent heat flux derived from eddy covariance measurements (the open circle located on the X-axis).The diagram compares correlation values (dotted black lines), the root-mean-square difference (solid dark red lines) and the standard variation (distance from the origin) of modeled and measured λE.Solid black line: standard deviation of measured λE.Input surface temperatures (T) as for Table2.For model abbreviations, cf.Table1.

Figure 6 .
Figure 6.Linear model regression fit between the daytime latent heat flux derived from eddy covariance measurements (λE EC) and modeled daytime λE.Gray dots: cloud fraction < 0.2, red triangles: cloud fraction ≥ 0.2-0.8,blue dots: cloud fraction ≥ 0.8.The solid line represents the line of identity, while the dashed line represents the linear regression line.Input temperatures (T) as for Table 2. Modeling approaches used: (a) SEB, (b) Bowen ratio, (c) T modeled, (d) Penman-Monteith.

Figure 7 .
Figure 7. Scatter plot between the vapor pressure deficit of the air (VPD) with measured λE (dotted black lines) and modeled λE (red dots) for 23 July 2012, indicating the clockwise hysteresis effect caused by the phase angle difference between net radiation and VPD.For model abbreviations, cf.Table1.
Figure 7. Scatter plot between the vapor pressure deficit of the air (VPD) with measured λE (dotted black lines) and modeled λE (red dots) for 23 July 2012, indicating the clockwise hysteresis effect caused by the phase angle difference between net radiation and VPD.For model abbreviations, cf.Table1.

1
Figure 7. Scatter plot between the vapor pressure deficit of the air (VPD) with measured λE (dotted black lines) and modeled λE (red dots) for 23 July 2012, indicating the clockwise hysteresis effect caused by the phase angle difference between net radiation and VPD.For model abbreviations, cf.Table1.
Figure 7. Scatter plot between the vapor pressure deficit of the air (VPD) with measured λE (dotted black lines) and modeled λE (red dots) for 23 July 2012, indicating the clockwise hysteresis effect caused by the phase angle difference between net radiation and VPD.For model abbreviations, cf.Table1.

Figure 8 .
Figure 8. Scatter plot between the modeled daytime canopy temperature (T modeled , cf.Section 3.3) and the inverted aerodynamic temperature (T aero,inv ), as well as the image-based temperature of wheat ears (T rad,WE ).CF = cloud fraction.

Table 1 .
Summary of the four modeling approaches: basic principle, input data and references.SEB, surface energy balance; TIR, thermal infrared; PM, Penman-Monteith.

Basic Principle Basic Input Data References SEB
[34]based approach, TIR-based calculation of H, λE derived as the residual of the SEB.Air density, net radiation, ground heat flux, aerodynamic resistance, air temperature, surface temperature[34]

Table 4 .
As in Table3, but for model estimates using T rad,mean .